61int main(
int argc,
char **argv) {
62 string top_file, trj_file, ptypes_file, out_file, grid, comment,
64 double min, max, step, coord;
73 namespace po = boost::program_options;
76 po::options_description desc(
"Allowed 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");
100 po::variables_map vm;
102 po::store(po::parse_command_line(argc, argv, desc), vm);
104 }
catch (po::error &err) {
105 cout <<
"error parsing command line: " << err.what() << endl;
110 if (vm.count(
"help")) {
112 cout << desc << endl;
121 if (coordinate.compare(
"x") != 0 && coordinate.compare(
"y") != 0 &&
122 coordinate.compare(
"z") != 0) {
123 cout <<
"Bad format for coordinate: " << coordinate << endl;
129 if (toks.size() != 3) {
130 cout <<
"Wrong range format, use min:step:max\n";
133 min = std::stod(toks[0]);
134 step = std::stod(toks[1]);
135 max = std::stod(toks[2]);
138 vector<votca::Index> ptypes;
140 Eigen::MatrixXi p_occ;
145 std::unique_ptr<TopologyReader> reader = std::unique_ptr<TopologyReader>(
147 if (reader ==
nullptr) {
148 throw std::runtime_error(
"input format not supported: " +
149 vm[
"top"].as<string>());
152 reader->ReadTopology(vm[
"top"].as<string>(), top);
155 if (vm.count(
"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>());
162 if (fl_ptypes.eof()) {
163 throw std::runtime_error(
"file " + vm[
"ptypes"].as<string>() +
166 while (!fl_ptypes.eof()) {
169 std::string string_tmp =
"__";
170 fl_ptypes >> string_tmp;
171 if (string_tmp !=
"" && string_tmp !=
"__") {
172 ptypes.push_back(std::stol(string_tmp));
179 for (
const auto &mol : top.
Molecules()) {
181 bool flag_found =
false;
182 votca::Index part_type = std::stol(mol.getBead(i)->getType());
184 if (part_type == ptype) {
189 ptypes.push_back(part_type);
195 p_occ = Eigen::MatrixXi::Zero(ptypes.size(), ptypes.size());
201 if (vm.count(
"shift_com")) {
202 for (
const auto &mol : top.
Molecules()) {
204 votca::Index part_type = std::stol(mol.getBead(i)->getType());
206 if (part_type == ptype) {
215 std::unique_ptr<TrajectoryReader> trajreader =
216 std::unique_ptr<TrajectoryReader>(
218 if (trajreader ==
nullptr) {
219 throw std::runtime_error(
"input format not supported: " +
220 vm[
"trj"].as<string>());
222 trajreader->Open(vm[
"trj"].as<string>());
226 bool moreframes =
true;
227 bool not_the_last =
true;
232 moreframes = trajreader->FirstFrame(top);
234 moreframes = trajreader->NextFrame(top);
238 if (last_frame == -1) {
241 if (frame_id <= last_frame) {
250 if (vm.count(
"shift_com")) {
251 for (
const auto &mol : top.
Molecules()) {
253 votca::Index part_type = std::stol(mol.getBead(i)->getType());
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();
261 com += mol.getBead(i)->getPos().z();
267 com /= (double)n_part;
271 if (moreframes && frame_id >= first_frame && not_the_last) {
274 for (
const auto &mol : top.
Molecules()) {
276 votca::Index part_type = std::stol(mol.getBead(i)->getType());
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();
284 coord = mol.getBead(i)->getPos().z();
287 if (coord - com > min && coord - com < max) {
289 (
votca::Index)std::floor((coord - com - min) / step))++;
301 }
catch (std::exception &error) {
302 cerr <<
"An error occured!" << endl << error.what() << endl;
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>());
313 fl_out <<
"#z\t" << flush;
315 fl_out <<
"type " << ptype <<
"\t" << flush;
319 fl_out << min + (double)k * step <<
"\t" << flush;
321 if (p_occ(j, k) == 0) {
322 fl_out << 0 <<
"\t" << flush;
324 fl_out << p_occ(j, k) / (1. * (double)analyzed_frames) <<
"\t"
331 }
catch (std::exception &error) {
332 cerr <<
"An error occured!" << endl << error.what() << endl;
334 cout <<
"The table was written to " << vm[
"out"].as<
string>() << endl;