ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/StaticProps.cpp
Revision: 2049
Committed: Tue Jan 6 21:44:10 2015 UTC (10 years, 3 months ago) by gezelter
File size: 18843 byte(s)
Log Message:
Added Geometric HBond StaticProps module.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42
43 #include <iostream>
44 #include <fstream>
45 #include <string>
46
47 #include "brains/SimCreator.hpp"
48 #include "brains/SimInfo.hpp"
49 #include "io/DumpReader.hpp"
50 #include "utils/simError.h"
51
52 #include "applications/staticProps/StaticPropsCmd.h"
53 #include "applications/staticProps/StaticAnalyser.hpp"
54 #include "applications/staticProps/GofR.hpp"
55 #include "applications/staticProps/GofZ.hpp"
56 #include "applications/staticProps/GofRZ.hpp"
57 #include "applications/staticProps/GofRAngle.hpp"
58 #include "applications/staticProps/GofAngle2.hpp"
59 #include "applications/staticProps/GofXyz.hpp"
60 #include "applications/staticProps/TwoDGofR.hpp"
61 #include "applications/staticProps/P2OrderParameter.hpp"
62 #include "applications/staticProps/BondOrderParameter.hpp"
63 #include "applications/staticProps/BOPofR.hpp"
64 #include "applications/staticProps/RippleOP.hpp"
65 #include "applications/staticProps/SCDOrderParameter.hpp"
66 #include "applications/staticProps/DensityPlot.hpp"
67 #include "applications/staticProps/ObjectCount.hpp"
68 #include "applications/staticProps/RhoZ.hpp"
69 #include "applications/staticProps/pAngle.hpp"
70 #include "applications/staticProps/BondAngleDistribution.hpp"
71 #include "applications/staticProps/NanoVolume.hpp"
72 #include "applications/staticProps/NanoLength.hpp"
73 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
74 #include "applications/staticProps/Hxy.hpp"
75 #endif
76 #include "applications/staticProps/RhoR.hpp"
77 #include "applications/staticProps/AngleR.hpp"
78 #include "applications/staticProps/TetrahedralityParam.hpp"
79 #include "applications/staticProps/TetrahedralityParamZ.hpp"
80 #include "applications/staticProps/TetrahedralityParamXYZ.hpp"
81 #include "applications/staticProps/RNEMDStats.hpp"
82 #include "applications/staticProps/NitrileFrequencyMap.hpp"
83 #include "applications/staticProps/MultipoleSum.hpp"
84 #include "applications/staticProps/SurfaceDiffusion.hpp"
85 #include "applications/staticProps/HBondGeometric.hpp"
86
87 using namespace OpenMD;
88
89 int main(int argc, char* argv[]){
90
91
92 gengetopt_args_info args_info;
93
94 //parse the command line option
95 if (cmdline_parser (argc, argv, &args_info) != 0) {
96 exit(1) ;
97 }
98
99 //get the dumpfile name
100 std::string dumpFileName = args_info.input_arg;
101 std::string sele1;
102 std::string sele2;
103 std::string sele3;
104
105 // check the first selection argument, or set it to the environment
106 // variable, or failing that, set it to "select all"
107
108 if (args_info.sele1_given) {
109 sele1 = args_info.sele1_arg;
110 } else {
111 char* sele1Env= getenv("SELECTION1");
112 if (sele1Env) {
113 sele1 = sele1Env;
114 } else {
115 sele1 = "select all";
116 }
117 }
118
119 // check the second selection argument, or set it to the environment
120 // variable, or failing that, set it to the first selection
121
122 if (args_info.sele2_given) {
123 sele2 = args_info.sele2_arg;
124 } else {
125 char* sele2Env = getenv("SELECTION2");
126 if (sele2Env) {
127 sele2 = sele2Env;
128 } else {
129 //If sele2 is not specified, then the default behavior
130 //should be what is already intended for sele1
131 sele2 = sele1;
132 }
133 }
134
135 // check the third selection argument, which is only set if
136 // requested by the user
137
138 if (args_info.sele3_given) sele3 = args_info.sele3_arg;
139
140 bool batchMode;
141 if (args_info.scd_given){
142 if (args_info.sele1_given &&
143 args_info.sele2_given && args_info.sele3_given) {
144 batchMode = false;
145 } else if (args_info.molname_given &&
146 args_info.begin_given && args_info.end_given) {
147 if (args_info.begin_arg < 0 ||
148 args_info.end_arg < 0 || args_info.begin_arg > args_info.end_arg-2) {
149 sprintf( painCave.errMsg,
150 "below conditions are not satisfied:\n"
151 "0 <= begin && 0<= end && begin <= end-2\n");
152 painCave.severity = OPENMD_ERROR;
153 painCave.isFatal = 1;
154 simError();
155 }
156 batchMode = true;
157 } else{
158 sprintf( painCave.errMsg,
159 "either --sele1, --sele2, --sele3 are specified,"
160 " or --molname, --begin, --end are specified\n");
161 painCave.severity = OPENMD_ERROR;
162 painCave.isFatal = 1;
163 simError();
164 }
165 }
166
167 //parse md file and set up the system
168 SimCreator creator;
169 SimInfo* info = creator.createSim(dumpFileName);
170
171 RealType maxLen;
172 RealType zmaxLen;
173 if (args_info.length_given) {
174 maxLen = args_info.length_arg;
175 if (args_info.zlength_given){
176 zmaxLen = args_info.zlength_arg;
177 }
178 } else {
179 Mat3x3d hmat = info->getSnapshotManager()->getCurrentSnapshot()->getHmat();
180 maxLen = std::min(std::min(hmat(0, 0), hmat(1, 1)), hmat(2, 2)) /2.0;
181 zmaxLen = hmat(2,2);
182 }
183
184 StaticAnalyser* analyser;
185 if (args_info.gofr_given){
186 analyser= new GofR(info, dumpFileName, sele1, sele2, maxLen,
187 args_info.nbins_arg);
188 } else if (args_info.gofz_given) {
189 analyser= new GofZ(info, dumpFileName, sele1, sele2, maxLen,
190 args_info.nbins_arg);
191 } else if (args_info.r_z_given) {
192 analyser = new GofRZ(info, dumpFileName, sele1, sele2, maxLen, zmaxLen,
193 args_info.nbins_arg, args_info.nbins_z_arg);
194 } else if (args_info.r_theta_given) {
195 if (args_info.sele3_given)
196 analyser = new GofRTheta(info, dumpFileName, sele1, sele2, sele3, maxLen,
197 args_info.nbins_arg, args_info.nanglebins_arg);
198 else
199 analyser = new GofRTheta(info, dumpFileName, sele1, sele2, maxLen,
200 args_info.nbins_arg, args_info.nanglebins_arg);
201 } else if (args_info.r_omega_given) {
202 if (args_info.sele3_given)
203 analyser = new GofROmega(info, dumpFileName, sele1, sele2, sele3, maxLen,
204 args_info.nbins_arg, args_info.nanglebins_arg);
205 else
206 analyser = new GofROmega(info, dumpFileName, sele1, sele2, maxLen,
207 args_info.nbins_arg, args_info.nanglebins_arg);
208
209 } else if (args_info.theta_omega_given) {
210 if (args_info.sele3_given)
211 analyser = new GofAngle2(info, dumpFileName, sele1, sele2, sele3,
212 args_info.nanglebins_arg);
213 else
214 analyser = new GofAngle2(info, dumpFileName, sele1, sele2,
215 args_info.nanglebins_arg);
216 } else if (args_info.gxyz_given) {
217 if (args_info.refsele_given) {
218 analyser= new GofXyz(info, dumpFileName, sele1, sele2,
219 args_info.refsele_arg, maxLen, args_info.nbins_arg);
220 } else {
221 sprintf( painCave.errMsg,
222 "--refsele must set when --gxyz is used");
223 painCave.severity = OPENMD_ERROR;
224 painCave.isFatal = 1;
225 simError();
226 }
227 } else if (args_info.twodgofr_given){
228 if (args_info.dz_given) {
229 analyser= new TwoDGofR(info, dumpFileName, sele1, sele2, maxLen,
230 args_info.dz_arg, args_info.nbins_arg);
231 } else {
232 sprintf( painCave.errMsg,
233 "A slab width (dz) must be specified when calculating TwoDGofR");
234 painCave.severity = OPENMD_ERROR;
235 painCave.isFatal = 1;
236 simError();
237 }
238 } else if (args_info.p2_given) {
239 if (args_info.sele1_given) {
240 if (args_info.sele2_given)
241 analyser = new P2OrderParameter(info, dumpFileName, sele1, sele2);
242 else
243 if (args_info.seleoffset_given)
244 analyser = new P2OrderParameter(info, dumpFileName, sele1,
245 args_info.seleoffset_arg);
246 else
247 analyser = new P2OrderParameter(info, dumpFileName, sele1);
248 } else {
249 sprintf( painCave.errMsg,
250 "At least one selection script (--sele1) must be specified when calculating P2 order parameters");
251 painCave.severity = OPENMD_ERROR;
252 painCave.isFatal = 1;
253 simError();
254 }
255 } else if (args_info.rp2_given){
256 analyser = new RippleOP(info, dumpFileName, sele1, sele2);
257 } else if (args_info.bo_given){
258 if (args_info.rcut_given) {
259 analyser = new BondOrderParameter(info, dumpFileName, sele1,
260 args_info.rcut_arg,
261 args_info.nbins_arg);
262 } else {
263 sprintf( painCave.errMsg,
264 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
265 painCave.severity = OPENMD_ERROR;
266 painCave.isFatal = 1;
267 simError();
268 }
269 } else if (args_info.multipole_given){
270 analyser = new MultipoleSum(info, dumpFileName, sele1,
271 maxLen, args_info.nbins_arg);
272 } else if (args_info.tet_param_given) {
273 if (args_info.rcut_given) {
274 analyser = new TetrahedralityParam(info, dumpFileName, sele1,
275 args_info.rcut_arg,
276 args_info.nbins_arg);
277 } else {
278 sprintf( painCave.errMsg,
279 "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
280 painCave.severity = OPENMD_ERROR;
281 painCave.isFatal = 1;
282 simError();
283 }
284 } else if (args_info.tet_param_z_given) {
285 if (args_info.rcut_given) {
286 analyser = new TetrahedralityParamZ(info, dumpFileName, sele1, sele2,
287 args_info.rcut_arg,
288 args_info.nbins_arg);
289 } else {
290 sprintf( painCave.errMsg,
291 "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
292 painCave.severity = OPENMD_ERROR;
293 painCave.isFatal = 1;
294 simError();
295 }
296 } else if (args_info.tet_param_xyz_given) {
297 if (!args_info.rcut_given) {
298 sprintf( painCave.errMsg,
299 "A cutoff radius (rcut) must be specified when calculating"
300 " Tetrahedrality Parameters");
301 painCave.severity = OPENMD_ERROR;
302 painCave.isFatal = 1;
303 simError();
304 }
305 if (!args_info.voxelSize_given) {
306 sprintf( painCave.errMsg,
307 "A voxel size must be specified when calculating"
308 " volume-resolved Tetrahedrality Parameters");
309 painCave.severity = OPENMD_ERROR;
310 painCave.isFatal = 1;
311 simError();
312 }
313 if (!args_info.gaussWidth_given) {
314 sprintf( painCave.errMsg,
315 "A gaussian width must be specified when calculating"
316 " volume-resolved Tetrahedrality Parameters");
317 painCave.severity = OPENMD_ERROR;
318 painCave.isFatal = 1;
319 simError();
320 }
321 analyser = new TetrahedralityParamXYZ(info, dumpFileName, sele1, sele2,
322 args_info.rcut_arg,
323 args_info.voxelSize_arg,
324 args_info.gaussWidth_arg);
325 } else if (args_info.ior_given){
326 if (args_info.rcut_given) {
327 analyser = new IcosahedralOfR(info, dumpFileName, sele1,
328 args_info.rcut_arg,
329 args_info.nbins_arg, maxLen);
330 } else {
331 sprintf( painCave.errMsg,
332 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
333 painCave.severity = OPENMD_ERROR;
334 painCave.isFatal = 1;
335 simError();
336 }
337 } else if (args_info.for_given){
338 if (args_info.rcut_given) {
339 analyser = new FCCOfR(info, dumpFileName, sele1, args_info.rcut_arg,
340 args_info.nbins_arg, maxLen);
341 } else {
342 sprintf( painCave.errMsg,
343 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
344 painCave.severity = OPENMD_ERROR;
345 painCave.isFatal = 1;
346 simError();
347 }
348 } else if (args_info.bad_given){
349 if (args_info.rcut_given) {
350 analyser = new BondAngleDistribution(info, dumpFileName, sele1,
351 args_info.rcut_arg,
352 args_info.nbins_arg);
353 } else {
354 sprintf( painCave.errMsg,
355 "A cutoff radius (rcut) must be specified when calculating Bond Angle Distributions");
356 painCave.severity = OPENMD_ERROR;
357 painCave.isFatal = 1;
358 simError();
359 }
360 } else if (args_info.scd_given) {
361 if (batchMode) {
362 analyser = new SCDOrderParameter(info, dumpFileName,
363 args_info.molname_arg,
364 args_info.begin_arg, args_info.end_arg);
365 } else{
366 analyser = new SCDOrderParameter(info, dumpFileName,
367 sele1, sele2, sele3);
368 }
369 }else if (args_info.density_given) {
370 analyser= new DensityPlot(info, dumpFileName, sele1, sele2, maxLen,
371 args_info.nbins_arg);
372 } else if (args_info.count_given) {
373 analyser = new ObjectCount(info, dumpFileName, sele1 );
374 } else if (args_info.slab_density_given) {
375 analyser = new RhoZ(info, dumpFileName, sele1, args_info.nbins_arg);
376 } else if (args_info.rnemdz_given) {
377 analyser = new RNEMDZ(info, dumpFileName, sele1, args_info.nbins_arg);
378 } else if (args_info.rnemdr_given) {
379 analyser = new RNEMDR(info, dumpFileName, sele1, args_info.nbins_arg);
380 } else if (args_info.rnemdrt_given) {
381 analyser = new RNEMDRTheta(info, dumpFileName, sele1,
382 args_info.nbins_arg, args_info.nanglebins_arg);
383 } else if (args_info.nitrile_given) {
384 analyser = new NitrileFrequencyMap(info, dumpFileName, sele1,
385 args_info.nbins_arg);
386 } else if (args_info.p_angle_given) {
387 if (args_info.sele1_given) {
388 if (args_info.sele2_given)
389 analyser = new pAngle(info, dumpFileName, sele1, sele2,
390 args_info.nbins_arg);
391 else
392 if (args_info.seleoffset_given) {
393 if (args_info.seleoffset2_given) {
394 analyser = new pAngle(info, dumpFileName, sele1,
395 args_info.seleoffset_arg,
396 args_info.seleoffset2_arg,
397 args_info.nbins_arg);
398 } else {
399 analyser = new pAngle(info, dumpFileName, sele1,
400 args_info.seleoffset_arg,
401 args_info.nbins_arg);
402 }
403 } else
404 analyser = new pAngle(info, dumpFileName, sele1,
405 args_info.nbins_arg);
406 } else {
407 sprintf( painCave.errMsg,
408 "At least one selection script (--sele1) must be specified when "
409 "calculating P(angle) distributions");
410 painCave.severity = OPENMD_ERROR;
411 painCave.isFatal = 1;
412 simError();
413 }
414 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
415 }else if (args_info.hxy_given) {
416 analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg,
417 args_info.nbins_y_arg, args_info.nbins_arg);
418 #endif
419 }else if (args_info.surfDiffusion_given){
420 analyser = new SurfaceDiffusion(info, dumpFileName, sele1, maxLen);
421 }else if (args_info.rho_r_given) {
422 if (args_info.radius_given){
423 analyser = new RhoR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg,args_info.radius_arg);
424 }else{
425 sprintf( painCave.errMsg,
426 "A particle radius (radius) must be specified when calculating Rho(r)");
427 painCave.severity = OPENMD_ERROR;
428 painCave.isFatal = 1;
429 simError();
430 }
431 } else if (args_info.hullvol_given) {
432 analyser = new NanoVolume(info, dumpFileName, sele1);
433 } else if (args_info.rodlength_given) {
434 analyser = new NanoLength(info, dumpFileName, sele1);
435 } else if (args_info.angle_r_given) {
436 analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg);
437 } else if (args_info.hbond_given){
438 if (args_info.rcut_given) {
439 if (args_info.thetacut_given) {
440
441 analyser = new HBondGeometric(info, dumpFileName, sele1, sele2,
442 args_info.rcut_arg,
443 args_info.thetacut_arg,
444 args_info.nbins_arg);
445 } else {
446 sprintf( painCave.errMsg,
447 "A cutoff angle (thetacut) must be specified when calculating Hydrogen Bonding Statistics");
448 painCave.severity = OPENMD_ERROR;
449 painCave.isFatal = 1;
450 simError();
451 }
452 } else {
453 sprintf( painCave.errMsg,
454 "A cutoff radius (rcut) must be specified when calculating Hydrogen Bonding Statistics");
455 painCave.severity = OPENMD_ERROR;
456 painCave.isFatal = 1;
457 simError();
458 }
459 }
460
461 if (args_info.output_given) {
462 analyser->setOutputName(args_info.output_arg);
463 }
464 if (args_info.step_given) {
465 analyser->setStep(args_info.step_arg);
466 }
467
468 analyser->process();
469
470 delete analyser;
471 delete info;
472
473 return 0;
474 }

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date