ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/StaticProps.cpp
Revision: 1998
Committed: Fri May 30 19:48:35 2014 UTC (10 years, 11 months ago) by gezelter
File size: 15966 byte(s)
Log Message:
Added the MultipoleSum

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 1879 #include "applications/staticProps/RNEMDStats.hpp"
81 gezelter 1994 #include "applications/staticProps/NitrileFrequencyMap.hpp"
82 gezelter 1998 #include "applications/staticProps/MultipoleSum.hpp"
83 gezelter 1879
84 gezelter 1390 using namespace OpenMD;
85 tim 310
86     int main(int argc, char* argv[]){
87    
88 gezelter 1513
89     gengetopt_args_info args_info;
90    
91     //parse the command line option
92     if (cmdline_parser (argc, argv, &args_info) != 0) {
93     exit(1) ;
94     }
95    
96     //get the dumpfile name
97     std::string dumpFileName = args_info.input_arg;
98     std::string sele1;
99     std::string sele2;
100    
101     // check the first selection argument, or set it to the environment
102     // variable, or failing that, set it to "select all"
103    
104     if (args_info.sele1_given) {
105     sele1 = args_info.sele1_arg;
106     } else {
107     char* sele1Env= getenv("SELECTION1");
108     if (sele1Env) {
109     sele1 = sele1Env;
110     } else {
111     sele1 = "select all";
112     }
113     }
114    
115     // check the second selection argument, or set it to the environment
116 gezelter 1937 // variable, or failing that, set it to the first selection
117 gezelter 1513
118     if (args_info.sele2_given) {
119     sele2 = args_info.sele2_arg;
120     } else {
121 gezelter 1937 char* sele2Env = getenv("SELECTION2");
122 gezelter 1513 if (sele2Env) {
123     sele2 = sele2Env;
124     } else {
125 gezelter 1819 //If sele2 is not specified, then the default behavior
126     //should be what is already intended for sele1
127 jmichalk 1785 sele2 = sele1;
128 gezelter 1513 }
129     }
130 gezelter 1846
131 gezelter 1513 bool batchMode;
132     if (args_info.scd_given){
133 gezelter 1846 if (args_info.sele1_given &&
134     args_info.sele2_given && args_info.sele3_given) {
135 gezelter 1513 batchMode = false;
136 gezelter 1846 } else if (args_info.molname_given &&
137     args_info.begin_given && args_info.end_given) {
138     if (args_info.begin_arg < 0 ||
139     args_info.end_arg < 0 || args_info.begin_arg > args_info.end_arg-2) {
140 gezelter 1513 sprintf( painCave.errMsg,
141     "below conditions are not satisfied:\n"
142     "0 <= begin && 0<= end && begin <= end-2\n");
143     painCave.severity = OPENMD_ERROR;
144     painCave.isFatal = 1;
145     simError();
146     }
147     batchMode = true;
148     } else{
149     sprintf( painCave.errMsg,
150     "either --sele1, --sele2, --sele3 are specified,"
151     " or --molname, --begin, --end are specified\n");
152     painCave.severity = OPENMD_ERROR;
153     painCave.isFatal = 1;
154 gezelter 1846 simError();
155 gezelter 1513 }
156     }
157    
158     //parse md file and set up the system
159     SimCreator creator;
160     SimInfo* info = creator.createSim(dumpFileName);
161 tim 310
162 gezelter 1513 RealType maxLen;
163     RealType zmaxLen;
164     if (args_info.length_given) {
165     maxLen = args_info.length_arg;
166     if (args_info.zlength_given){
167     zmaxLen = args_info.zlength_arg;
168 chuckv 1445 }
169 gezelter 1513 } else {
170     Mat3x3d hmat = info->getSnapshotManager()->getCurrentSnapshot()->getHmat();
171     maxLen = std::min(std::min(hmat(0, 0), hmat(1, 1)), hmat(2, 2)) /2.0;
172     zmaxLen = hmat(2,2);
173     }
174 kstocke1 1522
175 gezelter 1513 StaticAnalyser* analyser;
176     if (args_info.gofr_given){
177     analyser= new GofR(info, dumpFileName, sele1, sele2, maxLen,
178 kstocke1 1522 args_info.nbins_arg);
179 gezelter 1513 } else if (args_info.gofz_given) {
180     analyser= new GofZ(info, dumpFileName, sele1, sele2, maxLen,
181 kstocke1 1522 args_info.nbins_arg);
182 gezelter 1513 } else if (args_info.r_z_given) {
183     analyser = new GofRZ(info, dumpFileName, sele1, sele2, maxLen, zmaxLen,
184 kstocke1 1522 args_info.nbins_arg, args_info.nbins_z_arg);
185 gezelter 1513 } else if (args_info.r_theta_given) {
186     analyser = new GofRTheta(info, dumpFileName, sele1, sele2, maxLen,
187 kstocke1 1522 args_info.nbins_arg, args_info.nanglebins_arg);
188 gezelter 1513 } else if (args_info.r_omega_given) {
189     analyser = new GofROmega(info, dumpFileName, sele1, sele2, maxLen,
190 kstocke1 1522 args_info.nbins_arg, args_info.nanglebins_arg);
191 gezelter 1513 } else if (args_info.theta_omega_given) {
192     analyser = new GofAngle2(info, dumpFileName, sele1, sele2,
193 kstocke1 1522 args_info.nanglebins_arg);
194 gezelter 1513 } else if (args_info.gxyz_given) {
195     if (args_info.refsele_given) {
196 gezelter 1846 analyser= new GofXyz(info, dumpFileName, sele1, sele2,
197     args_info.refsele_arg, maxLen, args_info.nbins_arg);
198 gezelter 1039 } else {
199 gezelter 1513 sprintf( painCave.errMsg,
200 kstocke1 1522 "--refsele must set when --gxyz is used");
201 gezelter 1513 painCave.severity = OPENMD_ERROR;
202     painCave.isFatal = 1;
203     simError();
204 tim 310 }
205 gezelter 1513 } else if (args_info.twodgofr_given){
206     if (args_info.dz_given) {
207     analyser= new TwoDGofR(info, dumpFileName, sele1, sele2, maxLen,
208 kstocke1 1522 args_info.dz_arg, args_info.nbins_arg);
209 chuckv 1445 } else {
210 gezelter 1513 sprintf( painCave.errMsg,
211 kstocke1 1522 "A slab width (dz) must be specified when calculating TwoDGofR");
212 gezelter 1513 painCave.severity = OPENMD_ERROR;
213     painCave.isFatal = 1;
214     simError();
215 kstocke1 1522 }
216 gezelter 1513 } else if (args_info.p2_given) {
217 gezelter 1542 if (args_info.sele1_given) {
218     if (args_info.sele2_given)
219     analyser = new P2OrderParameter(info, dumpFileName, sele1, sele2);
220     else
221 gezelter 1819 if (args_info.seleoffset_given)
222     analyser = new P2OrderParameter(info, dumpFileName, sele1,
223     args_info.seleoffset_arg);
224     else
225     analyser = new P2OrderParameter(info, dumpFileName, sele1);
226 gezelter 1542 } else {
227     sprintf( painCave.errMsg,
228     "At least one selection script (--sele1) must be specified when calculating P2 order parameters");
229     painCave.severity = OPENMD_ERROR;
230     painCave.isFatal = 1;
231     simError();
232     }
233 gezelter 1513 } else if (args_info.rp2_given){
234     analyser = new RippleOP(info, dumpFileName, sele1, sele2);
235     } else if (args_info.bo_given){
236     if (args_info.rcut_given) {
237     analyser = new BondOrderParameter(info, dumpFileName, sele1,
238 kstocke1 1522 args_info.rcut_arg,
239     args_info.nbins_arg);
240 gezelter 1513 } else {
241     sprintf( painCave.errMsg,
242 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
243 gezelter 1513 painCave.severity = OPENMD_ERROR;
244     painCave.isFatal = 1;
245     simError();
246 tim 544 }
247 gezelter 1998 } else if (args_info.multipole_given){
248     if (args_info.rcut_given) {
249     analyser = new MultipoleSum(info, dumpFileName, sele1,
250     args_info.rcut_arg);
251     } else {
252     sprintf( painCave.errMsg,
253     "A cutoff radius (rcut) must be specified when calculating Multipole Sums");
254     painCave.severity = OPENMD_ERROR;
255     painCave.isFatal = 1;
256     simError();
257     }
258 kstocke1 1522
259     } else if (args_info.tet_param_given) {
260     if (args_info.rcut_given) {
261     analyser = new TetrahedralityParam(info, dumpFileName, sele1,
262     args_info.rcut_arg,
263     args_info.nbins_arg);
264     } else {
265     sprintf( painCave.errMsg,
266     "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
267     painCave.severity = OPENMD_ERROR;
268     painCave.isFatal = 1;
269     simError();
270     }
271 gezelter 1782 } else if (args_info.tet_param_z_given) {
272     if (args_info.rcut_given) {
273 gezelter 1846 analyser = new TetrahedralityParamZ(info, dumpFileName, sele1, sele2,
274     args_info.rcut_arg,
275     args_info.nbins_arg);
276 gezelter 1782 } else {
277     sprintf( painCave.errMsg,
278     "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
279     painCave.severity = OPENMD_ERROR;
280     painCave.isFatal = 1;
281     simError();
282     }
283 gezelter 1992 } else if (args_info.ior_given){
284 gezelter 1513 if (args_info.rcut_given) {
285 gezelter 1992 analyser = new IcosahedralOfR(info, dumpFileName, sele1,
286     args_info.rcut_arg,
287     args_info.nbins_arg, maxLen);
288     } else {
289     sprintf( painCave.errMsg,
290     "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
291     painCave.severity = OPENMD_ERROR;
292     painCave.isFatal = 1;
293     simError();
294     }
295     } else if (args_info.for_given){
296     if (args_info.rcut_given) {
297     analyser = new FCCOfR(info, dumpFileName, sele1, args_info.rcut_arg,
298 kstocke1 1522 args_info.nbins_arg, maxLen);
299 chuckv 1445 } else {
300 gezelter 1513 sprintf( painCave.errMsg,
301 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
302 gezelter 1513 painCave.severity = OPENMD_ERROR;
303     painCave.isFatal = 1;
304     simError();
305     }
306     } else if (args_info.bad_given){
307     if (args_info.rcut_given) {
308 gezelter 1846 analyser = new BondAngleDistribution(info, dumpFileName, sele1,
309     args_info.rcut_arg,
310 kstocke1 1522 args_info.nbins_arg);
311 gezelter 1513 } else {
312     sprintf( painCave.errMsg,
313 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Angle Distributions");
314 gezelter 1513 painCave.severity = OPENMD_ERROR;
315     painCave.isFatal = 1;
316     simError();
317 gezelter 1846 }
318 gezelter 1513 } else if (args_info.scd_given) {
319     if (batchMode) {
320 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
321     args_info.molname_arg,
322 kstocke1 1522 args_info.begin_arg, args_info.end_arg);
323 gezelter 1513 } else{
324     std::string sele3 = args_info.sele3_arg;
325 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
326     sele1, sele2, sele3);
327 gezelter 1513 }
328     }else if (args_info.density_given) {
329     analyser= new DensityPlot(info, dumpFileName, sele1, sele2, maxLen,
330 kstocke1 1522 args_info.nbins_arg);
331 gezelter 1513 } else if (args_info.count_given) {
332     analyser = new ObjectCount(info, dumpFileName, sele1 );
333     } else if (args_info.slab_density_given) {
334     analyser = new RhoZ(info, dumpFileName, sele1, args_info.nbins_arg);
335 gezelter 1879 } else if (args_info.rnemdz_given) {
336     analyser = new RNEMDZ(info, dumpFileName, sele1, args_info.nbins_arg);
337     } else if (args_info.rnemdr_given) {
338     analyser = new RNEMDR(info, dumpFileName, sele1, args_info.nbins_arg);
339 gezelter 1953 } else if (args_info.rnemdrt_given) {
340     analyser = new RNEMDRTheta(info, dumpFileName, sele1,
341     args_info.nbins_arg, args_info.nanglebins_arg);
342 gezelter 1994 } else if (args_info.nitrile_given) {
343     analyser = new NitrileFrequencyMap(info, dumpFileName, sele1,
344     args_info.nbins_arg);
345 gezelter 1513 } else if (args_info.p_angle_given) {
346 gezelter 1979 if (args_info.sele1_given) {
347     if (args_info.sele2_given)
348     analyser = new pAngle(info, dumpFileName, sele1, sele2,
349     args_info.nbins_arg);
350     else
351 gezelter 1991 if (args_info.seleoffset_given) {
352     if (args_info.seleoffset2_given) {
353     analyser = new pAngle(info, dumpFileName, sele1,
354     args_info.seleoffset_arg,
355     args_info.seleoffset2_arg,
356     args_info.nbins_arg);
357     } else {
358     analyser = new pAngle(info, dumpFileName, sele1,
359     args_info.seleoffset_arg,
360     args_info.nbins_arg);
361     }
362     } else
363 gezelter 1979 analyser = new pAngle(info, dumpFileName, sele1,
364     args_info.nbins_arg);
365     } else {
366     sprintf( painCave.errMsg,
367     "At least one selection script (--sele1) must be specified when "
368     "calculating P(angle) distributions");
369     painCave.severity = OPENMD_ERROR;
370     painCave.isFatal = 1;
371     simError();
372     }
373 gezelter 957 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
374 gezelter 1513 }else if (args_info.hxy_given) {
375     analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg,
376 kstocke1 1522 args_info.nbins_y_arg, args_info.nbins_arg);
377 gezelter 956 #endif
378 gezelter 1513 }else if (args_info.rho_r_given) {
379     if (args_info.radius_given){
380     analyser = new RhoR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg,args_info.radius_arg);
381     }else{
382     sprintf( painCave.errMsg,
383 kstocke1 1522 "A particle radius (radius) must be specified when calculating Rho(r)");
384 gezelter 1513 painCave.severity = OPENMD_ERROR;
385     painCave.isFatal = 1;
386     simError();
387 chuckv 1091 }
388 kstocke1 1522 } else if (args_info.hullvol_given) {
389 gezelter 1513 analyser = new NanoVolume(info, dumpFileName, sele1);
390 gezelter 1585 } else if (args_info.rodlength_given) {
391     analyser = new NanoLength(info, dumpFileName, sele1);
392 kstocke1 1522 } else if (args_info.angle_r_given) {
393     analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg);
394 gezelter 1513 }
395 gezelter 1846
396 gezelter 1513 if (args_info.output_given) {
397     analyser->setOutputName(args_info.output_arg);
398     }
399     if (args_info.step_given) {
400     analyser->setStep(args_info.step_arg);
401     }
402 gezelter 1846
403 gezelter 1513 analyser->process();
404 kstocke1 1522
405 gezelter 1513 delete analyser;
406     delete info;
407 gezelter 1846
408 gezelter 1513 return 0;
409 tim 310 }

Properties

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