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

# User Rev Content
1 gezelter 507 /*
2 tim 310 * 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 gezelter 1390 * 1. Redistributions of source code must retain the above copyright
10 tim 310 * notice, this list of conditions and the following disclaimer.
11     *
12 gezelter 1390 * 2. Redistributions in binary form must reproduce the above copyright
13 tim 310 * 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 gezelter 1390 *
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 gezelter 1879 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 gezelter 1782 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 tim 310 */
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 tim 311 #include "applications/staticProps/StaticPropsCmd.h"
53 tim 543 #include "applications/staticProps/StaticAnalyser.hpp"
54 tim 311 #include "applications/staticProps/GofR.hpp"
55 xsun 1213 #include "applications/staticProps/GofZ.hpp"
56 gezelter 1440 #include "applications/staticProps/GofRZ.hpp"
57 tim 311 #include "applications/staticProps/GofRAngle.hpp"
58     #include "applications/staticProps/GofAngle2.hpp"
59     #include "applications/staticProps/GofXyz.hpp"
60 gezelter 1454 #include "applications/staticProps/TwoDGofR.hpp"
61 tim 543 #include "applications/staticProps/P2OrderParameter.hpp"
62 gezelter 1039 #include "applications/staticProps/BondOrderParameter.hpp"
63 chuckv 1128 #include "applications/staticProps/BOPofR.hpp"
64 xsun 980 #include "applications/staticProps/RippleOP.hpp"
65 tim 544 #include "applications/staticProps/SCDOrderParameter.hpp"
66 tim 545 #include "applications/staticProps/DensityPlot.hpp"
67 gezelter 1513 #include "applications/staticProps/ObjectCount.hpp"
68 tim 840 #include "applications/staticProps/RhoZ.hpp"
69 gezelter 1413 #include "applications/staticProps/pAngle.hpp"
70 chuckv 1180 #include "applications/staticProps/BondAngleDistribution.hpp"
71     #include "applications/staticProps/NanoVolume.hpp"
72 gezelter 1585 #include "applications/staticProps/NanoLength.hpp"
73 gezelter 957 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
74 xsun 955 #include "applications/staticProps/Hxy.hpp"
75 gezelter 956 #endif
76 chuckv 1091 #include "applications/staticProps/RhoR.hpp"
77 kstocke1 1522 #include "applications/staticProps/AngleR.hpp"
78     #include "applications/staticProps/TetrahedralityParam.hpp"
79 gezelter 1782 #include "applications/staticProps/TetrahedralityParamZ.hpp"
80 gezelter 2015 #include "applications/staticProps/TetrahedralityParamXYZ.hpp"
81 gezelter 1879 #include "applications/staticProps/RNEMDStats.hpp"
82 gezelter 1994 #include "applications/staticProps/NitrileFrequencyMap.hpp"
83 gezelter 1998 #include "applications/staticProps/MultipoleSum.hpp"
84 jmichalk 2031 #include "applications/staticProps/SurfaceDiffusion.hpp"
85 gezelter 2049 #include "applications/staticProps/HBondGeometric.hpp"
86 gezelter 1879
87 gezelter 1390 using namespace OpenMD;
88 tim 310
89     int main(int argc, char* argv[]){
90    
91 gezelter 1513
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 gezelter 2023 std::string sele3;
104 gezelter 1513
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 gezelter 1937 // variable, or failing that, set it to the first selection
121 gezelter 1513
122     if (args_info.sele2_given) {
123     sele2 = args_info.sele2_arg;
124     } else {
125 gezelter 1937 char* sele2Env = getenv("SELECTION2");
126 gezelter 1513 if (sele2Env) {
127     sele2 = sele2Env;
128     } else {
129 gezelter 1819 //If sele2 is not specified, then the default behavior
130     //should be what is already intended for sele1
131 jmichalk 1785 sele2 = sele1;
132 gezelter 1513 }
133     }
134 gezelter 1846
135 gezelter 2023 // 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 gezelter 1513 bool batchMode;
141     if (args_info.scd_given){
142 gezelter 1846 if (args_info.sele1_given &&
143     args_info.sele2_given && args_info.sele3_given) {
144 gezelter 1513 batchMode = false;
145 gezelter 1846 } 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 gezelter 1513 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 gezelter 1846 simError();
164 gezelter 1513 }
165     }
166    
167     //parse md file and set up the system
168     SimCreator creator;
169     SimInfo* info = creator.createSim(dumpFileName);
170 tim 310
171 gezelter 1513 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 chuckv 1445 }
178 gezelter 1513 } 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 kstocke1 1522
184 gezelter 1513 StaticAnalyser* analyser;
185     if (args_info.gofr_given){
186     analyser= new GofR(info, dumpFileName, sele1, sele2, maxLen,
187 kstocke1 1522 args_info.nbins_arg);
188 gezelter 1513 } else if (args_info.gofz_given) {
189     analyser= new GofZ(info, dumpFileName, sele1, sele2, maxLen,
190 kstocke1 1522 args_info.nbins_arg);
191 gezelter 1513 } else if (args_info.r_z_given) {
192     analyser = new GofRZ(info, dumpFileName, sele1, sele2, maxLen, zmaxLen,
193 kstocke1 1522 args_info.nbins_arg, args_info.nbins_z_arg);
194 gezelter 1513 } else if (args_info.r_theta_given) {
195 gezelter 2023 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 gezelter 1513 } else if (args_info.r_omega_given) {
202 gezelter 2023 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 gezelter 1513 } else if (args_info.theta_omega_given) {
210 gezelter 2023 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 gezelter 1513 } else if (args_info.gxyz_given) {
217     if (args_info.refsele_given) {
218 gezelter 1846 analyser= new GofXyz(info, dumpFileName, sele1, sele2,
219     args_info.refsele_arg, maxLen, args_info.nbins_arg);
220 gezelter 1039 } else {
221 gezelter 1513 sprintf( painCave.errMsg,
222 kstocke1 1522 "--refsele must set when --gxyz is used");
223 gezelter 1513 painCave.severity = OPENMD_ERROR;
224     painCave.isFatal = 1;
225     simError();
226 tim 310 }
227 gezelter 1513 } else if (args_info.twodgofr_given){
228     if (args_info.dz_given) {
229     analyser= new TwoDGofR(info, dumpFileName, sele1, sele2, maxLen,
230 kstocke1 1522 args_info.dz_arg, args_info.nbins_arg);
231 chuckv 1445 } else {
232 gezelter 1513 sprintf( painCave.errMsg,
233 kstocke1 1522 "A slab width (dz) must be specified when calculating TwoDGofR");
234 gezelter 1513 painCave.severity = OPENMD_ERROR;
235     painCave.isFatal = 1;
236     simError();
237 kstocke1 1522 }
238 gezelter 1513 } else if (args_info.p2_given) {
239 gezelter 1542 if (args_info.sele1_given) {
240     if (args_info.sele2_given)
241     analyser = new P2OrderParameter(info, dumpFileName, sele1, sele2);
242     else
243 gezelter 1819 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 gezelter 1542 } 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 gezelter 1513 } 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 kstocke1 1522 args_info.rcut_arg,
261     args_info.nbins_arg);
262 gezelter 1513 } else {
263     sprintf( painCave.errMsg,
264 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
265 gezelter 1513 painCave.severity = OPENMD_ERROR;
266     painCave.isFatal = 1;
267     simError();
268 tim 544 }
269 gezelter 1998 } else if (args_info.multipole_given){
270 gezelter 1999 analyser = new MultipoleSum(info, dumpFileName, sele1,
271     maxLen, args_info.nbins_arg);
272 kstocke1 1522 } 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 gezelter 1782 } else if (args_info.tet_param_z_given) {
285     if (args_info.rcut_given) {
286 gezelter 1846 analyser = new TetrahedralityParamZ(info, dumpFileName, sele1, sele2,
287     args_info.rcut_arg,
288     args_info.nbins_arg);
289 gezelter 1782 } 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 gezelter 2015 } 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 gezelter 1992 } else if (args_info.ior_given){
326 gezelter 1513 if (args_info.rcut_given) {
327 gezelter 1992 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 kstocke1 1522 args_info.nbins_arg, maxLen);
341 chuckv 1445 } else {
342 gezelter 1513 sprintf( painCave.errMsg,
343 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
344 gezelter 1513 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 gezelter 1846 analyser = new BondAngleDistribution(info, dumpFileName, sele1,
351     args_info.rcut_arg,
352 kstocke1 1522 args_info.nbins_arg);
353 gezelter 1513 } else {
354     sprintf( painCave.errMsg,
355 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Angle Distributions");
356 gezelter 1513 painCave.severity = OPENMD_ERROR;
357     painCave.isFatal = 1;
358     simError();
359 gezelter 1846 }
360 gezelter 1513 } else if (args_info.scd_given) {
361     if (batchMode) {
362 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
363     args_info.molname_arg,
364 kstocke1 1522 args_info.begin_arg, args_info.end_arg);
365 gezelter 1513 } else{
366 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
367     sele1, sele2, sele3);
368 gezelter 1513 }
369     }else if (args_info.density_given) {
370     analyser= new DensityPlot(info, dumpFileName, sele1, sele2, maxLen,
371 kstocke1 1522 args_info.nbins_arg);
372 gezelter 1513 } 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 gezelter 1879 } 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 gezelter 1953 } else if (args_info.rnemdrt_given) {
381     analyser = new RNEMDRTheta(info, dumpFileName, sele1,
382     args_info.nbins_arg, args_info.nanglebins_arg);
383 gezelter 1994 } else if (args_info.nitrile_given) {
384     analyser = new NitrileFrequencyMap(info, dumpFileName, sele1,
385     args_info.nbins_arg);
386 gezelter 1513 } else if (args_info.p_angle_given) {
387 gezelter 1979 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 gezelter 1991 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 gezelter 1979 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 gezelter 957 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
415 gezelter 1513 }else if (args_info.hxy_given) {
416     analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg,
417 kstocke1 1522 args_info.nbins_y_arg, args_info.nbins_arg);
418 gezelter 956 #endif
419 jmichalk 2031 }else if (args_info.surfDiffusion_given){
420     analyser = new SurfaceDiffusion(info, dumpFileName, sele1, maxLen);
421 gezelter 1513 }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 kstocke1 1522 "A particle radius (radius) must be specified when calculating Rho(r)");
427 gezelter 1513 painCave.severity = OPENMD_ERROR;
428     painCave.isFatal = 1;
429     simError();
430 chuckv 1091 }
431 kstocke1 1522 } else if (args_info.hullvol_given) {
432 gezelter 1513 analyser = new NanoVolume(info, dumpFileName, sele1);
433 gezelter 1585 } else if (args_info.rodlength_given) {
434     analyser = new NanoLength(info, dumpFileName, sele1);
435 kstocke1 1522 } else if (args_info.angle_r_given) {
436     analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg);
437 gezelter 2049 } 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 gezelter 1513 }
460 gezelter 1846
461 gezelter 1513 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 gezelter 1846
468 gezelter 1513 analyser->process();
469 kstocke1 1522
470 gezelter 1513 delete analyser;
471     delete info;
472 gezelter 1846
473 gezelter 1513 return 0;
474 tim 310 }

Properties

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