ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/StaticProps.cpp
Revision: 1846
Committed: Thu Jan 31 15:55:47 2013 UTC (12 years, 3 months ago) by gezelter
File size: 13287 byte(s)
Log Message:
Bug fix in Tetrahedrality parameter Z exposed another bug in SimInfo.

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

Properties

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