votca 2024.2-dev
Loading...
Searching...
No Matches
part_dist.cc
Go to the documentation of this file.
1/*
2 * Copyright 2009-2021 The VOTCA Development Team (http://www.votca.org)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 *
16 */
17
18// Standard includes
19#include <fstream>
20#include <iostream>
21
22// Third party includes
23#include <boost/program_options.hpp>
24#include <boost/tokenizer.hpp>
25
26// Local VOTCA includes
29#include <votca/csg/version.h>
30
31using namespace std;
32namespace po = boost::program_options;
33using namespace votca::csg;
34using namespace votca::tools;
35
43void help_text(void) {
44 votca::csg::HelpTextHeader("csg_part_dist");
45 cout << "This program reads a topology and (set of) trajectory(ies). For "
46 "every\n"
47 "binned value of a chosen coordinate, it outputs the time-averaged "
48 "number of\n"
49 "particles, listed by particle types.\n\n";
50}
51
52void check_option(po::options_description &desc, po::variables_map &vm,
53 const string &option) {
54 if (!vm.count(option)) {
55 cout << "csg_part_dist \n\n";
56 cout << desc << endl << "parameter " << option << " is not specified\n";
57 exit(1);
58 }
59}
60
61int main(int argc, char **argv) {
62 string top_file, trj_file, ptypes_file, out_file, grid, comment,
63 coordinate = "z";
64 double min, max, step, coord;
65
66 votca::Index first_frame(0), last_frame(-1);
67
68 // Load topology+trajectory formats
71
72 // read program options
73 namespace po = boost::program_options;
74
75 // Declare the supported options.
76 po::options_description desc("Allowed options");
77
78 // let cg_engine add some program options
79 desc.add_options()("top", po::value<string>(&top_file), "topology file")(
80 "trj", po::value<string>(&trj_file), "trajectory file")(
81 "grid", po::value<string>(&grid), "output grid spacing (min:step:max)")(
82 "out", po::value<string>(&out_file),
83 "output particle distribution table")(
84 "ptypes", po::value<string>(&ptypes_file),
85 "particle types to include in the analysis\n"
86 " arg: file - particle types separated by space"
87 "\n default: all particle types")("first_frame",
88 po::value<votca::Index>(&first_frame),
89 "first frame considered for analysis")(
90 "last_frame", po::value<votca::Index>(&last_frame),
91 "last frame considered for analysis")(
92 "coord", po::value<string>(&coordinate),
93 "coordinate analyzed ('x', 'y', or 'z' (default))")(
94 "shift_com", "shift center of mass to zero")(
95 "comment", po::value<string>(&comment),
96 "store a comment in the output table")("help",
97 "produce this help message");
98
99 // now read in the command line
100 po::variables_map vm;
101 try {
102 po::store(po::parse_command_line(argc, argv, desc), vm);
103 po::notify(vm);
104 } catch (po::error &err) {
105 cout << "error parsing command line: " << err.what() << endl;
106 return -1;
107 }
108
109 // does the user want help?
110 if (vm.count("help")) {
111 help_text();
112 cout << desc << endl;
113 return 0;
114 }
115
116 check_option(desc, vm, "top");
117 check_option(desc, vm, "trj");
118 check_option(desc, vm, "out");
119 check_option(desc, vm, "grid");
120
121 if (coordinate.compare("x") != 0 && coordinate.compare("y") != 0 &&
122 coordinate.compare("z") != 0) {
123 cout << "Bad format for coordinate: " << coordinate << endl;
124 return 1;
125 }
126
127 // Parse grid
128 vector<string> toks = Tokenizer(grid, ":").ToVector();
129 if (toks.size() != 3) {
130 cout << "Wrong range format, use min:step:max\n";
131 return 1;
132 }
133 min = std::stod(toks[0]);
134 step = std::stod(toks[1]);
135 max = std::stod(toks[2]);
136 // Calculate number of bins
137 votca::Index n_bins = (votca::Index)((max - min) / (1. * step) + 1);
138 vector<votca::Index> ptypes;
139 Topology top;
140 Eigen::MatrixXi p_occ;
141 votca::Index analyzed_frames(0);
142 try {
143
144 // Load topology
145 std::unique_ptr<TopologyReader> reader = std::unique_ptr<TopologyReader>(
146 TopReaderFactory().Create(vm["top"].as<string>()));
147 if (reader == nullptr) {
148 throw std::runtime_error("input format not supported: " +
149 vm["top"].as<string>());
150 }
151
152 reader->ReadTopology(vm["top"].as<string>(), top);
153
154 // Read the particle types file and save to variable ptypes
155 if (vm.count("ptypes")) {
156 ifstream fl_ptypes;
157 fl_ptypes.open(vm["ptypes"].as<string>());
158 if (!fl_ptypes.is_open()) {
159 throw std::runtime_error("can't open " + vm["ptypes"].as<string>());
160 }
161 // build the list of particle types
162 if (fl_ptypes.eof()) {
163 throw std::runtime_error("file " + vm["ptypes"].as<string>() +
164 " is empty");
165 }
166 while (!fl_ptypes.eof()) {
167 // Not very elegant, but makes sure we don't count the same element
168 // twice
169 std::string string_tmp = "__";
170 fl_ptypes >> string_tmp;
171 if (string_tmp != "" && string_tmp != "__") {
172 ptypes.push_back(std::stol(string_tmp));
173 }
174 }
175 fl_ptypes.close();
176
177 } else {
178 // Include all particle types
179 for (const auto &mol : top.Molecules()) {
180 for (votca::Index i = 0; i < mol.BeadCount(); ++i) {
181 bool flag_found = false;
182 votca::Index part_type = std::stol(mol.getBead(i)->getType());
183 for (votca::Index ptype : ptypes) {
184 if (part_type == ptype) {
185 flag_found = true;
186 }
187 }
188 if (!flag_found) {
189 ptypes.push_back(part_type);
190 }
191 }
192 }
193 }
194
195 p_occ = Eigen::MatrixXi::Zero(ptypes.size(), ptypes.size());
196
197 // If we need to shift the center of mass, calculate the number of
198 // particles (only the ones that belong to the particle type index
199 // ptypes)
200 votca::Index n_part = 0;
201 if (vm.count("shift_com")) {
202 for (const auto &mol : top.Molecules()) {
203 for (votca::Index i = 0; i < mol.BeadCount(); ++i) {
204 votca::Index part_type = std::stol(mol.getBead(i)->getType());
205 for (votca::Index ptype : ptypes) {
206 if (part_type == ptype) {
207 ++n_part;
208 }
209 }
210 }
211 }
212 }
213
214 // Now load trajectory
215 std::unique_ptr<TrajectoryReader> trajreader =
216 std::unique_ptr<TrajectoryReader>(
217 TrjReaderFactory().Create(vm["trj"].as<string>()));
218 if (trajreader == nullptr) {
219 throw std::runtime_error("input format not supported: " +
220 vm["trj"].as<string>());
221 }
222 trajreader->Open(vm["trj"].as<string>());
223
224 // Read the trajectory. Analyze each frame to obtain
225 // particle occupancy as a function of coordinate z.
226 bool moreframes = true;
227 bool not_the_last = true;
228 while (moreframes) {
229 // Read frame
230 votca::Index frame_id = 0;
231 if (frame_id == 0) {
232 moreframes = trajreader->FirstFrame(top);
233 } else {
234 moreframes = trajreader->NextFrame(top);
235 }
236
237 // Was this the last frame we read?
238 if (last_frame == -1) {
239 not_the_last = 1;
240 } else {
241 if (frame_id <= last_frame) {
242 not_the_last = 1;
243 } else {
244 not_the_last = 0;
245 }
246 }
247
248 // Calculate new center of mass position in the direction of 'coordinate'
249 double com = 0.;
250 if (vm.count("shift_com")) {
251 for (const auto &mol : top.Molecules()) {
252 for (votca::Index i = 0; i < mol.BeadCount(); ++i) {
253 votca::Index part_type = std::stol(mol.getBead(i)->getType());
254 for (votca::Index ptype : ptypes) {
255 if (part_type == ptype) {
256 if (coordinate.compare("x") == 0) {
257 com += mol.getBead(i)->getPos().x();
258 } else if (coordinate.compare("y") == 0) {
259 com += mol.getBead(i)->getPos().y();
260 } else {
261 com += mol.getBead(i)->getPos().z();
262 }
263 }
264 }
265 }
266 }
267 com /= (double)n_part;
268 }
269
270 // Analyze frame
271 if (moreframes && frame_id >= first_frame && not_the_last) {
272 ++analyzed_frames;
273 // Loop over each atom property
274 for (const auto &mol : top.Molecules()) {
275 for (votca::Index i = 0; i < mol.BeadCount(); ++i) {
276 votca::Index part_type = std::stol(mol.getBead(i)->getType());
277 for (votca::Index j = 0; j < votca::Index(ptypes.size()); ++j) {
278 if (part_type == ptypes[j]) {
279 if (coordinate.compare("x") == 0) {
280 coord = mol.getBead(i)->getPos().x();
281 } else if (coordinate.compare("y") == 0) {
282 coord = mol.getBead(i)->getPos().y();
283 } else {
284 coord = mol.getBead(i)->getPos().z();
285 }
286
287 if (coord - com > min && coord - com < max) {
288 p_occ(j,
289 (votca::Index)std::floor((coord - com - min) / step))++;
290 }
291 }
292 }
293 }
294 }
295 }
296 ++frame_id;
297 }
298
299 trajreader->Close();
300
301 } catch (std::exception &error) {
302 cerr << "An error occured!" << endl << error.what() << endl;
303 }
304
305 // Output particle occupancy
306 try {
307 ofstream fl_out;
308 fl_out.open(vm["out"].as<string>());
309 if (!fl_out.is_open()) {
310 throw std::runtime_error("can't open " + vm["out"].as<string>());
311 }
312
313 fl_out << "#z\t" << flush;
314 for (votca::Index ptype : ptypes) {
315 fl_out << "type " << ptype << "\t" << flush;
316 }
317 fl_out << endl;
318 for (votca::Index k = 0; k < n_bins; ++k) {
319 fl_out << min + (double)k * step << "\t" << flush;
320 for (votca::Index j = 0; j < votca::Index(ptypes.size()); ++j) {
321 if (p_occ(j, k) == 0) {
322 fl_out << 0 << "\t" << flush;
323 } else {
324 fl_out << p_occ(j, k) / (1. * (double)analyzed_frames) << "\t"
325 << flush;
326 }
327 }
328 fl_out << endl;
329 }
330 fl_out.close();
331 } catch (std::exception &error) {
332 cerr << "An error occured!" << endl << error.what() << endl;
333 }
334 cout << "The table was written to " << vm["out"].as<string>() << endl;
335
336 return 0;
337}
static void RegisterPlugins(void)
topology of the whole system
Definition topology.h:81
MoleculeContainer & Molecules()
Definition topology.h:182
break string into words
Definition tokenizer.h:72
std::vector< T > ToVector()
store all words in a vector of type T, does type conversion.
Definition tokenizer.h:109
STL namespace.
FileFormatFactory< TopologyReader > & TopReaderFactory()
void HelpTextHeader(const std::string &tool_name)
Definition version.cc:42
FileFormatFactory< TrajectoryReader > & TrjReaderFactory()
Eigen::Index Index
Definition types.h:26
void help_text(void)
Definition part_dist.cc:43
int main(int argc, char **argv)
Definition part_dist.cc:61
void check_option(po::options_description &desc, po::variables_map &vm, const string &option)
Definition part_dist.cc:52