OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
StaticProps.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include <fstream>
46#include <iostream>
47#include <memory>
48#include <string>
49
50#include "StaticPropsCmd.hpp"
51#include "applications/staticProps/BOPofR.hpp"
52#include "applications/staticProps/BondAngleDistribution.hpp"
53#include "applications/staticProps/BondOrderParameter.hpp"
54#include "applications/staticProps/DensityPlot.hpp"
55#include "applications/staticProps/GofAngle2.hpp"
56#include "applications/staticProps/GofR.hpp"
57#include "applications/staticProps/GofRAngle.hpp"
58#include "applications/staticProps/GofRAngle2.hpp"
59#include "applications/staticProps/GofRZ.hpp"
60#include "applications/staticProps/GofXyz.hpp"
61#include "applications/staticProps/GofZ.hpp"
62#include "applications/staticProps/KirkwoodBuff.hpp"
63#include "applications/staticProps/NanoLength.hpp"
64#include "applications/staticProps/NanoVolume.hpp"
65#include "applications/staticProps/ObjectCount.hpp"
66#include "applications/staticProps/P2OrderParameter.hpp"
67#include "applications/staticProps/P2R.hpp"
68#include "applications/staticProps/PipeDensity.hpp"
69#include "applications/staticProps/RhoZ.hpp"
70#include "applications/staticProps/RippleOP.hpp"
71#include "applications/staticProps/SCDOrderParameter.hpp"
72#include "applications/staticProps/StaticAnalyser.hpp"
73#include "applications/staticProps/TwoDGofR.hpp"
74#include "applications/staticProps/pAngle.hpp"
75#include "brains/SimCreator.hpp"
76#include "brains/SimInfo.hpp"
77#include "io/DumpReader.hpp"
78#include "utils/simError.h"
79#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
80#include "applications/staticProps/Hxy.hpp"
81#endif
82#include "applications/staticProps/AngleR.hpp"
83#include "applications/staticProps/ChargeDensityZ.hpp"
84#include "applications/staticProps/ChargeHistogram.hpp"
85#include "applications/staticProps/ChargeR.hpp"
86#include "applications/staticProps/ChargeZ.hpp"
87#include "applications/staticProps/CoordinationNumber.hpp"
88#include "applications/staticProps/CurrentDensity.hpp"
89#include "applications/staticProps/DensityHistogram.hpp"
90#include "applications/staticProps/DipoleOrientation.hpp"
91#include "applications/staticProps/Field.hpp"
92#include "applications/staticProps/HBondGeometric.hpp"
93#include "applications/staticProps/HBondR.hpp"
94#include "applications/staticProps/HBondRvol.hpp"
95#include "applications/staticProps/HBondZ.hpp"
96#include "applications/staticProps/HBondZvol.hpp"
97#include "applications/staticProps/Kirkwood.hpp"
98#include "applications/staticProps/MassDensityR.hpp"
99#include "applications/staticProps/MassDensityZ.hpp"
100#include "applications/staticProps/MomentumHistogram.hpp"
101#include "applications/staticProps/MultipoleSum.hpp"
102#include "applications/staticProps/NitrileFrequencyMap.hpp"
103#include "applications/staticProps/NumberR.hpp"
104#include "applications/staticProps/NumberZ.hpp"
105#include "applications/staticProps/OrderParameterProbZ.hpp"
106#include "applications/staticProps/PositionZ.hpp"
107#include "applications/staticProps/PotDiff.hpp"
108#include "applications/staticProps/RNEMDStats.hpp"
109#include "applications/staticProps/RhoR.hpp"
110#include "applications/staticProps/SurfaceDiffusion.hpp"
111#include "applications/staticProps/TetrahedralityHBMatrix.hpp"
112#include "applications/staticProps/TetrahedralityParam.hpp"
113#include "applications/staticProps/TetrahedralityParamDens.hpp"
114#include "applications/staticProps/TetrahedralityParamR.hpp"
115#include "applications/staticProps/TetrahedralityParamXYZ.hpp"
116#include "applications/staticProps/TetrahedralityParamZ.hpp"
117#include "applications/staticProps/VelocityZ.hpp"
118
119using namespace OpenMD;
120
121int main(int argc, char* argv[]) {
122 gengetopt_args_info args_info;
123
124 // parse the command line option
125 if (cmdline_parser(argc, argv, &args_info) != 0) { exit(1); }
126
127 // get the dumpfile name
128 std::string dumpFileName = args_info.input_arg;
129 std::string sele1;
130 std::string sele2;
131 std::string sele3;
132 std::string comsele;
133
134 // check the first selection argument, or set it to the environment
135 // variable, or failing that, set it to "select all"
136
137 if (args_info.sele1_given) {
138 sele1 = args_info.sele1_arg;
139 } else {
140 char* sele1Env = getenv("SELECTION1");
141 if (sele1Env) {
142 sele1 = sele1Env;
143 } else {
144 sele1 = "select all";
145 }
146 }
147
148 // check the second selection argument, or set it to the environment
149 // variable, or failing that, set it to the first selection
150
151 if (args_info.sele2_given) {
152 sele2 = args_info.sele2_arg;
153 } else {
154 char* sele2Env = getenv("SELECTION2");
155 if (sele2Env) {
156 sele2 = sele2Env;
157 } else {
158 // If sele2 is not specified, then the default behavior
159 // should be what is already intended for sele1
160 sele2 = sele1;
161 }
162 }
163
164 // check the third selection argument, which is only set if
165 // requested by the user
166
167 if (args_info.sele3_given) sele3 = args_info.sele3_arg;
168
169 // check the comsele selection argument, which is only set if
170 // requested by the user
171
172 if (args_info.comsele_given) comsele = args_info.comsele_arg;
173
174 bool batchMode(false);
175 if (args_info.scd_given) {
176 if (args_info.sele1_given && args_info.sele2_given &&
177 args_info.sele3_given) {
178 batchMode = false;
179 } else if (args_info.molname_given && args_info.begin_given &&
180 args_info.end_given) {
181 if (args_info.begin_arg < 0 || args_info.end_arg < 0 ||
182 args_info.begin_arg > args_info.end_arg - 2) {
183 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
184 "below conditions are not satisfied:\n"
185 "0 <= begin && 0<= end && begin <= end-2\n");
186 painCave.severity = OPENMD_ERROR;
187 painCave.isFatal = 1;
188 simError();
189 }
190 batchMode = true;
191 } else {
192 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
193 "either --sele1, --sele2, --sele3 are specified,"
194 " or --molname, --begin, --end are specified\n");
195 painCave.severity = OPENMD_ERROR;
196 painCave.isFatal = 1;
197 simError();
198 }
199 }
200
201 // parse md file and set up the system
202 SimCreator creator;
203 SimInfo* info = creator.createSim(dumpFileName);
204
205 // important step until we move this into SimCreator:
206 info->update();
207
208 // convert privilegedAxis to corresponding integer
209 // x axis -> 0
210 // y axis -> 1
211 // z axis -> 2 (default)
212
213 int privilegedAxis;
214 switch (args_info.privilegedAxis_arg) {
215 case privilegedAxis_arg_x:
216 privilegedAxis = 0;
217 break;
218 case privilegedAxis_arg_y:
219 privilegedAxis = 1;
220 break;
221 case privilegedAxis_arg_z:
222 default:
223 privilegedAxis = 2;
224 break;
225 }
226
227 int privilegedAxis2;
228 switch (args_info.privilegedAxis2_arg) {
229 case privilegedAxis2_arg_x:
230 privilegedAxis2 = 0;
231 break;
232 case privilegedAxis2_arg_y:
233 privilegedAxis2 = 1;
234 break;
235 case privilegedAxis2_arg_z:
236 default:
237 privilegedAxis2 = 2;
238 break;
239 }
240
241 RealType maxLen;
242 RealType zmaxLen(0.0);
243 if (args_info.length_given) {
244 maxLen = args_info.length_arg;
245 if (args_info.zlength_given) { zmaxLen = args_info.zlength_arg; }
246 } else {
247 Mat3x3d hmat = info->getSnapshotManager()->getCurrentSnapshot()->getHmat();
248 // The maximum length for radial distribution functions is actually half
249 // the smallest box length
250 maxLen = std::min(std::min(hmat(0, 0), hmat(1, 1)), hmat(2, 2)) / 2.0;
251 // This should be the extent of the privileged axis:
252 zmaxLen = hmat(privilegedAxis, privilegedAxis);
253 }
254
255 int nanglebins, nrbins;
256 // in case we override nbins with nrbins:
257 if (args_info.nrbins_given) {
258 nrbins = args_info.nrbins_arg;
259 } else {
260 nrbins = args_info.nbins_arg;
261 }
262 // in case we override nbins with nanglebins:
263 if (args_info.nanglebins_given) {
264 nanglebins = args_info.nanglebins_arg;
265 } else {
266 nanglebins = args_info.nbins_arg;
267 }
268
269 RealType binWidth = args_info.binWidth_arg;
270
271 // override default vander waals radius for fictious atoms in a model
272 RealType vRadius;
273 if (args_info.v_radius_given) {
274 vRadius = args_info.v_radius_arg;
275 } else {
276 vRadius = 1.52;
277 }
278
279 int momentum_type;
280 switch (args_info.momentum_arg) {
281 case momentum_arg_P:
282 momentum_type = 0;
283 break;
284 case momentum_arg_J:
285 default:
286 momentum_type = 1;
287 break;
288 }
289
290 int momentum_comp;
291 switch (args_info.component_arg) {
292 case component_arg_x:
293 momentum_comp = 0;
294 break;
295 case component_arg_y:
296 momentum_comp = 1;
297 break;
298 case component_arg_z:
299 default:
300 momentum_comp = 2;
301 break;
302 }
303
304 std::unique_ptr<StaticAnalyser> analyser {nullptr};
305
306 if (args_info.gofr_given) {
307 analyser = std::make_unique<GofR>(info, dumpFileName, sele1, sele2, maxLen,
308 nrbins);
309 } else if (args_info.gofz_given) {
310 analyser =
311 std::make_unique<GofZ>(info, dumpFileName, sele1, sele2, maxLen,
312 zmaxLen, args_info.nbins_arg, privilegedAxis);
313 } else if (args_info.r_z_given) {
314 analyser = std::make_unique<GofRZ>(info, dumpFileName, sele1, sele2, maxLen,
315 zmaxLen, nrbins, args_info.nbins_z_arg,
316 privilegedAxis);
317 } else if (args_info.r_theta_given) {
318 if (args_info.sele3_given)
319 analyser = std::make_unique<GofRTheta>(info, dumpFileName, sele1, sele2,
320 sele3, maxLen, nrbins, nanglebins);
321 else
322 analyser = std::make_unique<GofRTheta>(info, dumpFileName, sele1, sele2,
323 maxLen, nrbins, nanglebins);
324 } else if (args_info.r_omega_given) {
325 if (args_info.sele3_given)
326 analyser = std::make_unique<GofROmega>(info, dumpFileName, sele1, sele2,
327 sele3, maxLen, nrbins, nanglebins);
328 else
329 analyser = std::make_unique<GofROmega>(info, dumpFileName, sele1, sele2,
330 maxLen, nrbins, nanglebins);
331 } else if (args_info.theta_omega_given) {
332 if (args_info.sele3_given)
333 analyser = std::make_unique<GofAngle2>(info, dumpFileName, sele1, sele2,
334 sele3, nanglebins);
335 else
336 analyser = std::make_unique<GofAngle2>(info, dumpFileName, sele1, sele2,
337 nanglebins);
338 } else if (args_info.r_theta_omega_given) {
339 if (args_info.sele3_given)
340 analyser = std::make_unique<GofRAngle2>(
341 info, dumpFileName, sele1, sele2, sele3, maxLen, nrbins, nanglebins);
342 else
343 analyser = std::make_unique<GofRAngle2>(info, dumpFileName, sele1, sele2,
344 maxLen, nrbins, nanglebins);
345 } else if (args_info.gxyz_given) {
346 if (args_info.refsele_given) {
347 analyser = std::make_unique<GofXyz>(info, dumpFileName, sele1, sele2,
348 args_info.refsele_arg, maxLen,
349 args_info.nbins_arg);
350 } else {
351 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
352 "--refsele must set when --gxyz is used");
353 painCave.severity = OPENMD_ERROR;
354 painCave.isFatal = 1;
355 simError();
356 }
357 } else if (args_info.twodgofr_given) {
358 if (args_info.dz_given) {
359 analyser = std::make_unique<TwoDGofR>(info, dumpFileName, sele1, sele2,
360 maxLen, args_info.dz_arg, nrbins);
361 } else {
362 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
363 "A slab width (dz) must be specified when calculating TwoDGofR");
364 painCave.severity = OPENMD_ERROR;
365 painCave.isFatal = 1;
366 simError();
367 }
368 } else if (args_info.kirkwood_buff_given) {
369 if (args_info.sele1_given && args_info.sele2_given) {
370 analyser = std::make_unique<KirkwoodBuff>(info, dumpFileName, sele1,
371 sele2, maxLen, nrbins);
372 } else {
373 snprintf(
374 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
375 "Two selection scripts (--sele1 and --sele2) must be specified when "
376 "calculating Kirkwood Buff integrals");
377 painCave.severity = OPENMD_ERROR;
378 painCave.isFatal = 1;
379 simError();
380 }
381 } else if (args_info.p2_given) {
382 if (args_info.sele1_given) {
383 if (args_info.sele2_given) {
384 analyser = std::make_unique<P2OrderParameter>(info, dumpFileName, sele1,
385 sele2);
386 } else if (args_info.seleoffset_given) {
387 analyser = std::make_unique<P2OrderParameter>(info, dumpFileName, sele1,
388 args_info.seleoffset_arg);
389 } else {
390 analyser =
391 std::make_unique<P2OrderParameter>(info, dumpFileName, sele1);
392 }
393 } else {
394 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
395 "At least one selection script (--sele1) must be specified when "
396 "calculating P2 order parameters");
397 painCave.severity = OPENMD_ERROR;
398 painCave.isFatal = 1;
399 simError();
400 }
401 } else if (args_info.rp2_given) {
402 analyser = std::make_unique<RippleOP>(info, dumpFileName, sele1, sele2);
403 } else if (args_info.bo_given) {
404 if (args_info.rcut_given) {
405 analyser = std::make_unique<BondOrderParameter>(
406 info, dumpFileName, sele1, args_info.rcut_arg, args_info.nbins_arg);
407 } else {
408 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
409 "A cutoff radius (rcut) must be specified when calculating Bond "
410 "Order "
411 "Parameters");
412 painCave.severity = OPENMD_ERROR;
413 painCave.isFatal = 1;
414 simError();
415 }
416 } else if (args_info.multipole_given) {
417 analyser = std::make_unique<MultipoleSum>(info, dumpFileName, sele1, maxLen,
418 args_info.nbins_arg);
419
420 } else if (args_info.tet_param_given) {
421 if (args_info.rcut_given) {
422 analyser = std::make_unique<TetrahedralityParam>(
423 info, dumpFileName, sele1, args_info.rcut_arg, args_info.nbins_arg);
424 } else {
425 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
426 "A cutoff radius (rcut) must be specified when calculating "
427 "Tetrahedrality "
428 "Parameters");
429 painCave.severity = OPENMD_ERROR;
430 painCave.isFatal = 1;
431 simError();
432 }
433 } else if (args_info.tet_param_z_given) {
434 if (args_info.rcut_given) {
435 analyser = std::make_unique<TetrahedralityParamZ>(
436 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
437 args_info.nbins_arg, privilegedAxis);
438 } else {
439 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
440 "A cutoff radius (rcut) must be specified when calculating "
441 "Tetrahedrality "
442 "Parameters");
443 painCave.severity = OPENMD_ERROR;
444 painCave.isFatal = 1;
445 simError();
446 }
447 } else if (args_info.tet_param_r_given) {
448 if (args_info.rcut_given) {
449 if (args_info.sele3_given) {
450 analyser = std::make_unique<TetrahedralityParamR>(
451 info, dumpFileName, sele1, sele2, sele3, args_info.rcut_arg, maxLen,
452 nrbins);
453 } else {
454 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
455 "Selection3 (--sele3) must be given when calculating "
456 "Tetrahedrality Parameter Qk(r)");
457 painCave.severity = OPENMD_ERROR;
458 painCave.isFatal = 1;
459 simError();
460 }
461 } else {
462 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
463 "A cutoff radius (rcut) must be specified when calculating "
464 "Tetrahedrality "
465 "Parameters");
466 painCave.severity = OPENMD_ERROR;
467 painCave.isFatal = 1;
468 simError();
469 }
470 } else if (args_info.tet_param_dens_given) {
471 if (args_info.rcut_given) {
472 analyser = std::make_unique<TetrahedralityParamDens>(
473 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
474 args_info.nbins_arg);
475 } else {
476 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
477 "A cutoff radius (rcut) must be specified when calculating "
478 "Tetrahedrality "
479 "Parameters");
480 painCave.severity = OPENMD_ERROR;
481 painCave.isFatal = 1;
482 simError();
483 }
484 } else if (args_info.tet_hb_given) {
485 if (args_info.rcut_given) {
486 analyser = std::make_unique<TetrahedralityHBMatrix>(
487 info, dumpFileName, sele1, args_info.rcut_arg, args_info.OOcut_arg,
488 args_info.thetacut_arg, args_info.OHcut_arg, args_info.nbins_arg);
489 } else {
490 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
491 "A cutoff radius (rcut) must be specified when calculating "
492 " Tetrahedrality Hydrogen Bonding Matrix");
493 painCave.severity = OPENMD_ERROR;
494 painCave.isFatal = 1;
495 simError();
496 }
497 } else if (args_info.tet_param_xyz_given) {
498 if (!args_info.rcut_given) {
499 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
500 "A cutoff radius (rcut) must be specified when calculating"
501 " Tetrahedrality Parameters");
502 painCave.severity = OPENMD_ERROR;
503 painCave.isFatal = 1;
504 simError();
505 }
506 if (!args_info.voxelSize_given) {
507 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
508 "A voxel size must be specified when calculating"
509 " volume-resolved Tetrahedrality Parameters");
510 painCave.severity = OPENMD_ERROR;
511 painCave.isFatal = 1;
512 simError();
513 }
514 if (!args_info.gaussWidth_given) {
515 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
516 "A gaussian width must be specified when calculating"
517 " volume-resolved Tetrahedrality Parameters");
518 painCave.severity = OPENMD_ERROR;
519 painCave.isFatal = 1;
520 simError();
521 }
522 analyser = std::make_unique<TetrahedralityParamXYZ>(
523 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
524 args_info.voxelSize_arg, args_info.gaussWidth_arg);
525 } else if (args_info.ior_given) {
526 if (args_info.rcut_given) {
527 analyser = std::make_unique<IcosahedralOfR>(
528 info, dumpFileName, sele1, args_info.rcut_arg, nrbins, maxLen);
529 } else {
530 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
531 "A cutoff radius (rcut) must be specified when calculating Bond "
532 "Order "
533 "Parameters");
534 painCave.severity = OPENMD_ERROR;
535 painCave.isFatal = 1;
536 simError();
537 }
538 } else if (args_info.for_given) {
539 if (args_info.rcut_given) {
540 analyser = std::make_unique<FCCOfR>(info, dumpFileName, sele1,
541 args_info.rcut_arg, nrbins, maxLen);
542 } else {
543 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
544 "A cutoff radius (rcut) must be specified when calculating Bond "
545 "Order "
546 "Parameters");
547 painCave.severity = OPENMD_ERROR;
548 painCave.isFatal = 1;
549 simError();
550 }
551 } else if (args_info.bad_given) {
552 if (args_info.rcut_given) {
553 analyser = std::make_unique<BondAngleDistribution>(
554 info, dumpFileName, sele1, args_info.rcut_arg, args_info.nbins_arg);
555 } else {
556 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
557 "A cutoff radius (rcut) must be specified when calculating Bond "
558 "Angle "
559 "Distributions");
560 painCave.severity = OPENMD_ERROR;
561 painCave.isFatal = 1;
562 simError();
563 }
564 } else if (args_info.scd_given) {
565 if (batchMode) {
566 analyser = std::make_unique<SCDOrderParameter>(
567 info, dumpFileName, args_info.molname_arg, args_info.begin_arg,
568 args_info.end_arg);
569 } else {
570 analyser = std::make_unique<SCDOrderParameter>(info, dumpFileName, sele1,
571 sele2, sele3);
572 }
573 } else if (args_info.density_given) {
574 analyser = std::make_unique<DensityPlot>(info, dumpFileName, sele1, sele2,
575 maxLen, args_info.nbins_arg);
576 } else if (args_info.count_given) {
577 analyser = std::make_unique<ObjectCount>(info, dumpFileName, sele1);
578 } else if (args_info.slab_density_given) {
579 analyser = std::make_unique<RhoZ>(info, dumpFileName, sele1,
580 args_info.nbins_arg, privilegedAxis);
581 } else if (args_info.eam_density_given) {
582 analyser = std::make_unique<DensityHistogram>(info, dumpFileName, sele1,
583 args_info.nbins_arg);
584 } else if (args_info.momentum_distribution_given) {
585 analyser = std::make_unique<MomentumHistogram>(
586 info, dumpFileName, sele1, args_info.nbins_arg, momentum_type,
587 momentum_comp);
588 } else if (args_info.net_charge_given) {
589 analyser = std::make_unique<ChargeHistogram>(info, dumpFileName, sele1,
590 args_info.nbins_arg);
591 } else if (args_info.current_density_given) {
592 analyser = std::make_unique<CurrentDensity>(
593 info, dumpFileName, sele1, args_info.nbins_arg, privilegedAxis);
594 } else if (args_info.chargez_given) {
595 analyser = std::make_unique<ChargeZ>(info, dumpFileName, sele1,
596 args_info.nbins_arg, privilegedAxis);
597 } else if (args_info.charger_given) {
598 analyser = std::make_unique<ChargeR>(info, dumpFileName, sele1, maxLen,
599 args_info.nbins_arg);
600 } else if (args_info.numberz_given) {
601 analyser = std::make_unique<NumberZ>(info, dumpFileName, sele1,
602 args_info.nbins_arg, privilegedAxis);
603 } else if (args_info.numberr_given) {
604 analyser = std::make_unique<NumberR>(info, dumpFileName, sele1, maxLen,
605 args_info.nbins_arg);
606 } else if (args_info.massdensityz_given) {
607 analyser = std::make_unique<MassDensityZ>(
608 info, dumpFileName, sele1, args_info.nbins_arg, privilegedAxis);
609 } else if (args_info.massdensityr_given) {
610 analyser = std::make_unique<MassDensityR>(info, dumpFileName, sele1, maxLen,
611 args_info.nbins_arg);
612 } else if (args_info.charge_density_z_given) {
613 analyser = std::make_unique<ChargeDensityZ>(
614 info, dumpFileName, sele1, args_info.nbins_arg, vRadius,
615 args_info.atom_name_arg, args_info.gen_xyz_flag, privilegedAxis);
616 } else if (args_info.countz_given) {
617 analyser = std::make_unique<PositionZ>(info, dumpFileName, sele1,
618 args_info.nbins_arg, privilegedAxis);
619 } else if (args_info.pipe_density_given) {
620 switch (privilegedAxis) {
621 case 0:
622 analyser = std::make_unique<PipeDensity>(
623 info, dumpFileName, sele1, args_info.nbins_y_arg,
624 args_info.nbins_z_arg, privilegedAxis);
625 break;
626 case 1:
627 analyser = std::make_unique<PipeDensity>(
628 info, dumpFileName, sele1, args_info.nbins_z_arg,
629 args_info.nbins_x_arg, privilegedAxis);
630 break;
631 case 2:
632 default:
633 analyser = std::make_unique<PipeDensity>(
634 info, dumpFileName, sele1, args_info.nbins_x_arg,
635 args_info.nbins_y_arg, privilegedAxis);
636 break;
637 }
638 } else if (args_info.rnemdz_given) {
639 analyser = std::make_unique<RNEMDZ>(info, dumpFileName, sele1,
640 args_info.nbins_arg, privilegedAxis);
641 } else if (args_info.rnemdr_given) {
642 analyser = std::make_unique<RNEMDR>(info, dumpFileName, sele1, comsele,
643 nrbins, binWidth);
644 } else if (args_info.rnemdrt_given) {
645 analyser = std::make_unique<RNEMDRTheta>(info, dumpFileName, sele1, comsele,
646 nrbins, binWidth, nanglebins);
647 } else if (args_info.nitrile_given) {
648 analyser = std::make_unique<NitrileFrequencyMap>(info, dumpFileName, sele1,
649 args_info.nbins_arg);
650 } else if (args_info.p_angle_given) {
651 if (args_info.sele1_given) {
652 if (args_info.sele2_given)
653 analyser = std::make_unique<pAngle>(info, dumpFileName, sele1, sele2,
654 args_info.nbins_arg);
655 else if (args_info.seleoffset_given) {
656 if (args_info.seleoffset2_given) {
657 analyser = std::make_unique<pAngle>(
658 info, dumpFileName, sele1, args_info.seleoffset_arg,
659 args_info.seleoffset2_arg, args_info.nbins_arg);
660 } else {
661 analyser = std::make_unique<pAngle>(info, dumpFileName, sele1,
662 args_info.seleoffset_arg,
663 args_info.nbins_arg);
664 }
665 } else
666 analyser = std::make_unique<pAngle>(info, dumpFileName, sele1,
667 args_info.nbins_arg);
668 } else {
669 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
670 "At least one selection script (--sele1) must be specified when "
671 "calculating P(angle) distributions");
672 painCave.severity = OPENMD_ERROR;
673 painCave.isFatal = 1;
674 simError();
675 }
676#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
677 } else if (args_info.hxy_given) {
678 analyser = std::make_unique<Hxy>(
679 info, dumpFileName, sele1, args_info.nbins_x_arg, args_info.nbins_y_arg,
680 args_info.nbins_z_arg, args_info.nbins_arg);
681#endif
682 } else if (args_info.cn_given || args_info.scn_given || args_info.gcn_given) {
683 if (args_info.rcut_given) {
684 if (args_info.cn_given) {
685 analyser = std::make_unique<CoordinationNumber>(
686 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
687 args_info.nbins_arg);
688 } else if (args_info.scn_given) {
689 analyser =
690 std::make_unique<SCN>(info, dumpFileName, sele1, sele2,
691 args_info.rcut_arg, args_info.nbins_arg);
692 } else if (args_info.gcn_given) {
693 analyser =
694 std::make_unique<GCN>(info, dumpFileName, sele1, sele2,
695 args_info.rcut_arg, args_info.nbins_arg);
696 }
697 } else {
698 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
699 "A cutoff radius (rcut) must be specified when calculating\n"
700 "\t Coordination Numbers");
701 painCave.severity = OPENMD_ERROR;
702 painCave.isFatal = 1;
703 simError();
704 }
705 } else if (args_info.surfDiffusion_given) {
706 analyser =
707 std::make_unique<SurfaceDiffusion>(info, dumpFileName, sele1, maxLen);
708 } else if (args_info.rho_r_given) {
709 if (args_info.radius_given) {
710 analyser = std::make_unique<RhoR>(info, dumpFileName, sele1, maxLen,
711 nrbins, args_info.radius_arg);
712 } else {
713 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
714 "A particle radius (radius) must be specified when calculating "
715 "Rho(r)");
716 painCave.severity = OPENMD_ERROR;
717 painCave.isFatal = 1;
718 simError();
719 }
720 } else if (args_info.hullvol_given) {
721 analyser = std::make_unique<NanoVolume>(info, dumpFileName, sele1);
722 } else if (args_info.rodlength_given) {
723 analyser = std::make_unique<NanoLength>(info, dumpFileName, sele1);
724 } else if (args_info.angle_r_given) {
725 if (args_info.sele1_given) {
726 if (args_info.sele2_given)
727 analyser = std::make_unique<AngleR>(info, dumpFileName, sele1, sele2,
728 maxLen, nrbins);
729 else if (args_info.seleoffset_given) {
730 if (args_info.seleoffset2_given) {
731 analyser = std::make_unique<AngleR>(
732 info, dumpFileName, sele1, args_info.seleoffset_arg,
733 args_info.seleoffset2_arg, maxLen, nrbins);
734 } else {
735 analyser = std::make_unique<AngleR>(info, dumpFileName, sele1,
736 args_info.seleoffset_arg, maxLen,
737 nrbins);
738 }
739 } else
740 analyser =
741 std::make_unique<AngleR>(info, dumpFileName, sele1, maxLen, nrbins);
742 } else {
743 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
744 "At least one selection script (--sele1) must be specified when "
745 "calculating Angle(r) values");
746 painCave.severity = OPENMD_ERROR;
747 painCave.isFatal = 1;
748 simError();
749 }
750 } else if (args_info.p2r_given) {
751 if (args_info.sele1_given) {
752 if (args_info.sele2_given)
753 analyser = std::make_unique<P2R>(info, dumpFileName, sele1, sele2,
754 args_info.nbins_arg);
755 else if (args_info.seleoffset_given) {
756 if (args_info.seleoffset2_given) {
757 analyser = std::make_unique<P2R>(
758 info, dumpFileName, sele1, args_info.seleoffset_arg,
759 args_info.seleoffset2_arg, args_info.nbins_arg);
760 } else {
761 analyser = std::make_unique<P2R>(info, dumpFileName, sele1,
762 args_info.seleoffset_arg,
763 args_info.nbins_arg);
764 }
765 } else
766 analyser = std::make_unique<P2R>(info, dumpFileName, sele1,
767 args_info.nbins_arg);
768 } else {
769 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
770 "At least one selection script (--sele1) must be specified when "
771 "calculating P2R values");
772 painCave.severity = OPENMD_ERROR;
773 painCave.isFatal = 1;
774 simError();
775 }
776 } else if (args_info.hbond_given) {
777 if (args_info.rcut_given) {
778 if (args_info.thetacut_given) {
779 analyser = std::make_unique<HBondGeometric>(
780 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
781 args_info.thetacut_arg, args_info.nbins_arg);
782 } else {
783 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
784 "A cutoff angle (thetacut) must be specified when calculating "
785 "Hydrogen "
786 "Bonding Statistics");
787 painCave.severity = OPENMD_ERROR;
788 painCave.isFatal = 1;
789 simError();
790 }
791 } else {
792 snprintf(
793 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
794 "A cutoff radius (rcut) must be specified when calculating Hydrogen "
795 "Bonding Statistics");
796 painCave.severity = OPENMD_ERROR;
797 painCave.isFatal = 1;
798 simError();
799 }
800
801 } else if (args_info.hbondz_given) {
802 if (args_info.rcut_given) {
803 if (args_info.thetacut_given) {
804 analyser = std::make_unique<HBondZ>(
805 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
806 args_info.thetacut_arg, args_info.nbins_arg);
807 } else {
808 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
809 "A cutoff angle (thetacut) must be specified when calculating "
810 "Hydrogen "
811 "Bonding Statistics");
812 painCave.severity = OPENMD_ERROR;
813 painCave.isFatal = 1;
814 simError();
815 }
816 } else {
817 snprintf(
818 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
819 "A cutoff radius (rcut) must be specified when calculating Hydrogen "
820 "Bonding Statistics");
821 painCave.severity = OPENMD_ERROR;
822 painCave.isFatal = 1;
823 simError();
824 }
825 } else if (args_info.hbondzvol_given) {
826 if (args_info.rcut_given) {
827 if (args_info.thetacut_given) {
828 analyser = std::make_unique<HBondZvol>(
829 info, dumpFileName, sele1, sele2, args_info.rcut_arg,
830 args_info.thetacut_arg, args_info.nbins_arg);
831 } else {
832 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
833 "A cutoff angle (thetacut) must be specified when calculating "
834 "Hydrogen "
835 "Bonding Statistics");
836 painCave.severity = OPENMD_ERROR;
837 painCave.isFatal = 1;
838 simError();
839 }
840 } else {
841 snprintf(
842 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
843 "A cutoff radius (rcut) must be specified when calculating Hydrogen "
844 "Bonding Statistics");
845 painCave.severity = OPENMD_ERROR;
846 painCave.isFatal = 1;
847 simError();
848 }
849 } else if (args_info.hbondr_given) {
850 if (args_info.rcut_given) {
851 if (args_info.thetacut_given) {
852 analyser = std::make_unique<HBondR>(
853 info, dumpFileName, sele1, sele2, sele3, args_info.rcut_arg, maxLen,
854 args_info.thetacut_arg, args_info.nrbins_arg);
855 } else {
856 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
857 "A cutoff angle (thetacut) must be specified when calculating "
858 "Hydrogen "
859 "Bonding Statistics");
860 painCave.severity = OPENMD_ERROR;
861 painCave.isFatal = 1;
862 simError();
863 }
864 } else {
865 snprintf(
866 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
867 "A cutoff radius (rcut) must be specified when calculating Hydrogen "
868 "Bonding Statistics");
869 painCave.severity = OPENMD_ERROR;
870 painCave.isFatal = 1;
871 simError();
872 }
873 } else if (args_info.hbondrvol_given) {
874 if (args_info.rcut_given) {
875 if (args_info.thetacut_given) {
876 analyser = std::make_unique<HBondRvol>(
877 info, dumpFileName, sele1, sele2, sele3, args_info.rcut_arg, maxLen,
878 args_info.thetacut_arg, args_info.nrbins_arg);
879 } else {
880 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
881 "A cutoff angle (thetacut) must be specified when calculating "
882 "Hydrogen "
883 "Bonding Statistics");
884 painCave.severity = OPENMD_ERROR;
885 painCave.isFatal = 1;
886 simError();
887 }
888 } else {
889 snprintf(
890 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
891 "A cutoff radius (rcut) must be specified when calculating Hydrogen "
892 "Bonding Statistics");
893 painCave.severity = OPENMD_ERROR;
894 painCave.isFatal = 1;
895 simError();
896 }
897 } else if (args_info.potDiff_given) {
898 analyser = std::make_unique<PotDiff>(info, dumpFileName, sele1);
899 } else if (args_info.kirkwood_given) {
900 analyser = std::make_unique<Kirkwood>(info, dumpFileName, sele1, sele2,
901 maxLen, nrbins);
902 } else if (args_info.kirkwoodQ_given) {
903 analyser = std::make_unique<KirkwoodQuadrupoles>(info, dumpFileName, sele1,
904 sele2, maxLen, nrbins);
905 } else if (args_info.densityfield_given) {
906 analyser = std::make_unique<DensityField>(info, dumpFileName, sele1,
907 args_info.voxelSize_arg);
908 } else if (args_info.velocityfield_given) {
909 analyser = std::make_unique<VelocityField>(info, dumpFileName, sele1,
910 args_info.voxelSize_arg);
911 } else if (args_info.velocityZ_given) {
912 switch (privilegedAxis) {
913 case 0:
914 if (privilegedAxis2 == 1) {
915 analyser = std::make_unique<VelocityZ>(
916 info, dumpFileName, sele1, args_info.nbins_x_arg,
917 args_info.nbins_y_arg, privilegedAxis, privilegedAxis2);
918 } else if (privilegedAxis2 == 2) {
919 analyser = std::make_unique<VelocityZ>(
920 info, dumpFileName, sele1, args_info.nbins_x_arg,
921 args_info.nbins_z_arg, privilegedAxis, privilegedAxis2);
922 }
923 break;
924 case 1:
925 if (privilegedAxis2 == 0) {
926 analyser = std::make_unique<VelocityZ>(
927 info, dumpFileName, sele1, args_info.nbins_y_arg,
928 args_info.nbins_x_arg, privilegedAxis, privilegedAxis2);
929 } else if (privilegedAxis2 == 2) {
930 analyser = std::make_unique<VelocityZ>(
931 info, dumpFileName, sele1, args_info.nbins_y_arg,
932 args_info.nbins_z_arg, privilegedAxis, privilegedAxis2);
933 }
934 break;
935 case 2:
936 default:
937 if (privilegedAxis2 == 0) {
938 analyser = std::make_unique<VelocityZ>(
939 info, dumpFileName, sele1, args_info.nbins_z_arg,
940 args_info.nbins_x_arg, privilegedAxis, privilegedAxis2);
941 } else if (privilegedAxis2 == 1) {
942 analyser = std::make_unique<VelocityZ>(
943 info, dumpFileName, sele1, args_info.nbins_z_arg,
944 args_info.nbins_y_arg, privilegedAxis, privilegedAxis2);
945 }
946 break;
947 }
948 } else if (args_info.dipole_orientation_given) {
949 if (args_info.dipoleX_given && args_info.dipoleY_given &&
950 args_info.dipoleZ_given)
951 analyser = std::make_unique<DipoleOrientation>(
952 info, dumpFileName, sele1, args_info.dipoleX_arg,
953 args_info.dipoleY_arg, args_info.dipoleZ_arg, args_info.nbins_arg,
954 privilegedAxis);
955 else {
956 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
957 "Dipole components must be provided.");
958 painCave.severity = OPENMD_ERROR;
959 painCave.isFatal = 1;
960 simError();
961 }
962
963 } else if (args_info.order_prob_given) {
964 if (args_info.dipoleX_given && args_info.dipoleY_given &&
965 args_info.dipoleZ_given)
966 analyser = std::make_unique<OrderParameterProbZ>(
967 info, dumpFileName, sele1, args_info.dipoleX_arg,
968 args_info.dipoleY_arg, args_info.dipoleZ_arg, args_info.nbins_arg,
969 privilegedAxis);
970 else {
971 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
972 "Dipole components must be provided.");
973 painCave.severity = OPENMD_ERROR;
974 painCave.isFatal = 1;
975 simError();
976 }
977 }
978
979 if (analyser != NULL) {
980 if (args_info.output_given) {
981 analyser->setOutputName(args_info.output_arg);
982 }
983 if (args_info.step_given) { analyser->setStep(args_info.step_arg); }
984
985 analyser->process();
986 } else {
987 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
988 "StaticProps: No Analyser was created, nothing to do!");
989 painCave.severity = OPENMD_ERROR;
990 painCave.isFatal = 1;
991 simError();
992 }
993
994 delete info;
995
996 return 0;
997}
PotDiff.hpp
StaticAnalyser for Potential Energy changes with charges turned off.
StaticPropsCmd.hpp
The header file for the command line option parser generated by GNU Gengetopt version 2....
OpenMD::SimCreator
The only responsibility of SimCreator is to parse the meta-data file and create a SimInfo instance ba...
Definition SimCreator.hpp:71
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
OpenMD::SimInfo::update
void update()
update
Definition SimInfo.cpp:697
OpenMD::SimInfo::getSnapshotManager
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
OpenMD::Snapshot::getHmat
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
OpenMD::SnapshotManager::getCurrentSnapshot
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Definition SnapshotManager.hpp:89
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60
gengetopt_args_info
Where the command line options are stored.
Definition Dump2XYZCmd.hpp:39