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

Properties

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