ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/StaticProps.cpp
Revision: 1999
Committed: Sat May 31 21:00:46 2014 UTC (10 years, 11 months ago) by gezelter
File size: 15741 byte(s)
Log Message:
Updating MultipoleSum to scan the cutoff radii

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 gezelter 1999 analyser = new MultipoleSum(info, dumpFileName, sele1,
249     maxLen, args_info.nbins_arg);
250 kstocke1 1522 } else if (args_info.tet_param_given) {
251     if (args_info.rcut_given) {
252     analyser = new TetrahedralityParam(info, dumpFileName, sele1,
253     args_info.rcut_arg,
254     args_info.nbins_arg);
255     } else {
256     sprintf( painCave.errMsg,
257     "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
258     painCave.severity = OPENMD_ERROR;
259     painCave.isFatal = 1;
260     simError();
261     }
262 gezelter 1782 } else if (args_info.tet_param_z_given) {
263     if (args_info.rcut_given) {
264 gezelter 1846 analyser = new TetrahedralityParamZ(info, dumpFileName, sele1, sele2,
265     args_info.rcut_arg,
266     args_info.nbins_arg);
267 gezelter 1782 } else {
268     sprintf( painCave.errMsg,
269     "A cutoff radius (rcut) must be specified when calculating Tetrahedrality Parameters");
270     painCave.severity = OPENMD_ERROR;
271     painCave.isFatal = 1;
272     simError();
273     }
274 gezelter 1992 } else if (args_info.ior_given){
275 gezelter 1513 if (args_info.rcut_given) {
276 gezelter 1992 analyser = new IcosahedralOfR(info, dumpFileName, sele1,
277     args_info.rcut_arg,
278     args_info.nbins_arg, maxLen);
279     } else {
280     sprintf( painCave.errMsg,
281     "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
282     painCave.severity = OPENMD_ERROR;
283     painCave.isFatal = 1;
284     simError();
285     }
286     } else if (args_info.for_given){
287     if (args_info.rcut_given) {
288     analyser = new FCCOfR(info, dumpFileName, sele1, args_info.rcut_arg,
289 kstocke1 1522 args_info.nbins_arg, maxLen);
290 chuckv 1445 } else {
291 gezelter 1513 sprintf( painCave.errMsg,
292 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters");
293 gezelter 1513 painCave.severity = OPENMD_ERROR;
294     painCave.isFatal = 1;
295     simError();
296     }
297     } else if (args_info.bad_given){
298     if (args_info.rcut_given) {
299 gezelter 1846 analyser = new BondAngleDistribution(info, dumpFileName, sele1,
300     args_info.rcut_arg,
301 kstocke1 1522 args_info.nbins_arg);
302 gezelter 1513 } else {
303     sprintf( painCave.errMsg,
304 kstocke1 1522 "A cutoff radius (rcut) must be specified when calculating Bond Angle Distributions");
305 gezelter 1513 painCave.severity = OPENMD_ERROR;
306     painCave.isFatal = 1;
307     simError();
308 gezelter 1846 }
309 gezelter 1513 } else if (args_info.scd_given) {
310     if (batchMode) {
311 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
312     args_info.molname_arg,
313 kstocke1 1522 args_info.begin_arg, args_info.end_arg);
314 gezelter 1513 } else{
315     std::string sele3 = args_info.sele3_arg;
316 gezelter 1846 analyser = new SCDOrderParameter(info, dumpFileName,
317     sele1, sele2, sele3);
318 gezelter 1513 }
319     }else if (args_info.density_given) {
320     analyser= new DensityPlot(info, dumpFileName, sele1, sele2, maxLen,
321 kstocke1 1522 args_info.nbins_arg);
322 gezelter 1513 } else if (args_info.count_given) {
323     analyser = new ObjectCount(info, dumpFileName, sele1 );
324     } else if (args_info.slab_density_given) {
325     analyser = new RhoZ(info, dumpFileName, sele1, args_info.nbins_arg);
326 gezelter 1879 } else if (args_info.rnemdz_given) {
327     analyser = new RNEMDZ(info, dumpFileName, sele1, args_info.nbins_arg);
328     } else if (args_info.rnemdr_given) {
329     analyser = new RNEMDR(info, dumpFileName, sele1, args_info.nbins_arg);
330 gezelter 1953 } else if (args_info.rnemdrt_given) {
331     analyser = new RNEMDRTheta(info, dumpFileName, sele1,
332     args_info.nbins_arg, args_info.nanglebins_arg);
333 gezelter 1994 } else if (args_info.nitrile_given) {
334     analyser = new NitrileFrequencyMap(info, dumpFileName, sele1,
335     args_info.nbins_arg);
336 gezelter 1513 } else if (args_info.p_angle_given) {
337 gezelter 1979 if (args_info.sele1_given) {
338     if (args_info.sele2_given)
339     analyser = new pAngle(info, dumpFileName, sele1, sele2,
340     args_info.nbins_arg);
341     else
342 gezelter 1991 if (args_info.seleoffset_given) {
343     if (args_info.seleoffset2_given) {
344     analyser = new pAngle(info, dumpFileName, sele1,
345     args_info.seleoffset_arg,
346     args_info.seleoffset2_arg,
347     args_info.nbins_arg);
348     } else {
349     analyser = new pAngle(info, dumpFileName, sele1,
350     args_info.seleoffset_arg,
351     args_info.nbins_arg);
352     }
353     } else
354 gezelter 1979 analyser = new pAngle(info, dumpFileName, sele1,
355     args_info.nbins_arg);
356     } else {
357     sprintf( painCave.errMsg,
358     "At least one selection script (--sele1) must be specified when "
359     "calculating P(angle) distributions");
360     painCave.severity = OPENMD_ERROR;
361     painCave.isFatal = 1;
362     simError();
363     }
364 gezelter 957 #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
365 gezelter 1513 }else if (args_info.hxy_given) {
366     analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg,
367 kstocke1 1522 args_info.nbins_y_arg, args_info.nbins_arg);
368 gezelter 956 #endif
369 gezelter 1513 }else if (args_info.rho_r_given) {
370     if (args_info.radius_given){
371     analyser = new RhoR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg,args_info.radius_arg);
372     }else{
373     sprintf( painCave.errMsg,
374 kstocke1 1522 "A particle radius (radius) must be specified when calculating Rho(r)");
375 gezelter 1513 painCave.severity = OPENMD_ERROR;
376     painCave.isFatal = 1;
377     simError();
378 chuckv 1091 }
379 kstocke1 1522 } else if (args_info.hullvol_given) {
380 gezelter 1513 analyser = new NanoVolume(info, dumpFileName, sele1);
381 gezelter 1585 } else if (args_info.rodlength_given) {
382     analyser = new NanoLength(info, dumpFileName, sele1);
383 kstocke1 1522 } else if (args_info.angle_r_given) {
384     analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg);
385 gezelter 1513 }
386 gezelter 1846
387 gezelter 1513 if (args_info.output_given) {
388     analyser->setOutputName(args_info.output_arg);
389     }
390     if (args_info.step_given) {
391     analyser->setStep(args_info.step_arg);
392     }
393 gezelter 1846
394 gezelter 1513 analyser->process();
395 kstocke1 1522
396 gezelter 1513 delete analyser;
397     delete info;
398 gezelter 1846
399 gezelter 1513 return 0;
400 tim 310 }

Properties

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