| 1 | /* | 
| 2 | * 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 | * 1. Redistributions of source code must retain the above copyright | 
| 10 | *    notice, this list of conditions and the following disclaimer. | 
| 11 | * | 
| 12 | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 | *    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 | * | 
| 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, 234107 (2008). | 
| 39 | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 | */ | 
| 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 | #include "applications/staticProps/StaticPropsCmd.h" | 
| 53 | #include "applications/staticProps/StaticAnalyser.hpp" | 
| 54 | #include "applications/staticProps/GofR.hpp" | 
| 55 | #include "applications/staticProps/GofZ.hpp" | 
| 56 | #include "applications/staticProps/GofRZ.hpp" | 
| 57 | #include "applications/staticProps/GofRAngle.hpp" | 
| 58 | #include "applications/staticProps/GofAngle2.hpp" | 
| 59 | #include "applications/staticProps/GofXyz.hpp" | 
| 60 | #include "applications/staticProps/TwoDGofR.hpp" | 
| 61 | #include "applications/staticProps/P2OrderParameter.hpp" | 
| 62 | #include "applications/staticProps/BondOrderParameter.hpp" | 
| 63 | #include "applications/staticProps/BOPofR.hpp" | 
| 64 | #include "applications/staticProps/RippleOP.hpp" | 
| 65 | #include "applications/staticProps/SCDOrderParameter.hpp" | 
| 66 | #include "applications/staticProps/DensityPlot.hpp" | 
| 67 | #include "applications/staticProps/ObjectCount.hpp" | 
| 68 | #include "applications/staticProps/RhoZ.hpp" | 
| 69 | #include "applications/staticProps/pAngle.hpp" | 
| 70 | #include "applications/staticProps/BondAngleDistribution.hpp" | 
| 71 | #include "applications/staticProps/NanoVolume.hpp" | 
| 72 | #include "applications/staticProps/NanoLength.hpp" | 
| 73 | #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H) | 
| 74 | #include "applications/staticProps/Hxy.hpp" | 
| 75 | #endif | 
| 76 | #include "applications/staticProps/RhoR.hpp" | 
| 77 | #include "applications/staticProps/AngleR.hpp" | 
| 78 | #include "applications/staticProps/TetrahedralityParam.hpp" | 
| 79 | #include "applications/staticProps/TetrahedralityParamZ.hpp" | 
| 80 | #include "applications/staticProps/RNEMDStats.hpp" | 
| 81 |  | 
| 82 | using namespace OpenMD; | 
| 83 |  | 
| 84 | int main(int argc, char* argv[]){ | 
| 85 |  | 
| 86 |  | 
| 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 | // variable, or failing that, set it to "select all" | 
| 115 |  | 
| 116 | if (args_info.sele2_given) { | 
| 117 | sele2 = args_info.sele2_arg; | 
| 118 | } else { | 
| 119 | char* sele2Env = getenv("SELECTION1"); | 
| 120 | if (sele2Env) { | 
| 121 | sele2 = sele2Env; | 
| 122 | } else { | 
| 123 | //If sele2 is not specified, then the default behavior | 
| 124 | //should be what is already intended for sele1 | 
| 125 | sele2 = sele1; | 
| 126 | } | 
| 127 | } | 
| 128 |  | 
| 129 | bool batchMode; | 
| 130 | if (args_info.scd_given){ | 
| 131 | if (args_info.sele1_given && | 
| 132 | args_info.sele2_given && args_info.sele3_given) { | 
| 133 | batchMode = false; | 
| 134 | } 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 | 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 | simError(); | 
| 153 | } | 
| 154 | } | 
| 155 |  | 
| 156 | //parse md file and set up the system | 
| 157 | SimCreator creator; | 
| 158 | SimInfo* info = creator.createSim(dumpFileName); | 
| 159 |  | 
| 160 | 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 | } | 
| 167 | } 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 |  | 
| 173 | StaticAnalyser* analyser; | 
| 174 | if (args_info.gofr_given){ | 
| 175 | analyser= new GofR(info, dumpFileName, sele1, sele2, maxLen, | 
| 176 | args_info.nbins_arg); | 
| 177 | } else if (args_info.gofz_given) { | 
| 178 | analyser= new GofZ(info, dumpFileName, sele1, sele2, maxLen, | 
| 179 | args_info.nbins_arg); | 
| 180 | } else if (args_info.r_z_given) { | 
| 181 | analyser  = new GofRZ(info, dumpFileName, sele1, sele2, maxLen, zmaxLen, | 
| 182 | args_info.nbins_arg, args_info.nbins_z_arg); | 
| 183 | } else if (args_info.r_theta_given) { | 
| 184 | analyser  = new GofRTheta(info, dumpFileName, sele1, sele2, maxLen, | 
| 185 | args_info.nbins_arg, args_info.nanglebins_arg); | 
| 186 | } else if (args_info.r_omega_given) { | 
| 187 | analyser  = new GofROmega(info, dumpFileName, sele1, sele2, maxLen, | 
| 188 | args_info.nbins_arg, args_info.nanglebins_arg); | 
| 189 | } else if (args_info.theta_omega_given) { | 
| 190 | analyser  = new GofAngle2(info, dumpFileName, sele1, sele2, | 
| 191 | args_info.nanglebins_arg); | 
| 192 | } else if (args_info.gxyz_given) { | 
| 193 | if (args_info.refsele_given) { | 
| 194 | analyser= new GofXyz(info, dumpFileName, sele1, sele2, | 
| 195 | args_info.refsele_arg, maxLen, args_info.nbins_arg); | 
| 196 | } else { | 
| 197 | sprintf( painCave.errMsg, | 
| 198 | "--refsele must set when --gxyz is used"); | 
| 199 | painCave.severity = OPENMD_ERROR; | 
| 200 | painCave.isFatal = 1; | 
| 201 | simError(); | 
| 202 | } | 
| 203 | } else if (args_info.twodgofr_given){ | 
| 204 | if (args_info.dz_given) { | 
| 205 | analyser= new TwoDGofR(info, dumpFileName, sele1, sele2, maxLen, | 
| 206 | args_info.dz_arg, args_info.nbins_arg); | 
| 207 | } else { | 
| 208 | sprintf( painCave.errMsg, | 
| 209 | "A slab width (dz) must be specified when calculating TwoDGofR"); | 
| 210 | painCave.severity = OPENMD_ERROR; | 
| 211 | painCave.isFatal = 1; | 
| 212 | simError(); | 
| 213 | } | 
| 214 | } else if (args_info.p2_given) { | 
| 215 | if (args_info.sele1_given) { | 
| 216 | if (args_info.sele2_given) | 
| 217 | analyser  = new P2OrderParameter(info, dumpFileName, sele1, sele2); | 
| 218 | else | 
| 219 | 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 | } 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 | } 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 | args_info.rcut_arg, | 
| 237 | args_info.nbins_arg); | 
| 238 | } else { | 
| 239 | sprintf( painCave.errMsg, | 
| 240 | "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters"); | 
| 241 | painCave.severity = OPENMD_ERROR; | 
| 242 | painCave.isFatal = 1; | 
| 243 | simError(); | 
| 244 | } | 
| 245 |  | 
| 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 | } else if (args_info.tet_param_z_given) { | 
| 259 | if (args_info.rcut_given) { | 
| 260 | analyser = new TetrahedralityParamZ(info, dumpFileName, sele1, sele2, | 
| 261 | args_info.rcut_arg, | 
| 262 | args_info.nbins_arg); | 
| 263 | } 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 | } else if (args_info.bor_given){ | 
| 271 | if (args_info.rcut_given) { | 
| 272 | analyser = new BOPofR(info, dumpFileName, sele1, args_info.rcut_arg, | 
| 273 | args_info.nbins_arg, maxLen); | 
| 274 | } else { | 
| 275 | sprintf( painCave.errMsg, | 
| 276 | "A cutoff radius (rcut) must be specified when calculating Bond Order Parameters"); | 
| 277 | 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 | analyser = new BondAngleDistribution(info, dumpFileName, sele1, | 
| 284 | args_info.rcut_arg, | 
| 285 | args_info.nbins_arg); | 
| 286 | } else { | 
| 287 | sprintf( painCave.errMsg, | 
| 288 | "A cutoff radius (rcut) must be specified when calculating Bond Angle Distributions"); | 
| 289 | painCave.severity = OPENMD_ERROR; | 
| 290 | painCave.isFatal = 1; | 
| 291 | simError(); | 
| 292 | } | 
| 293 | } else if (args_info.scd_given) { | 
| 294 | if (batchMode) { | 
| 295 | analyser  = new SCDOrderParameter(info, dumpFileName, | 
| 296 | args_info.molname_arg, | 
| 297 | args_info.begin_arg, args_info.end_arg); | 
| 298 | } else{ | 
| 299 | std::string sele3 = args_info.sele3_arg; | 
| 300 | analyser  = new SCDOrderParameter(info, dumpFileName, | 
| 301 | sele1, sele2, sele3); | 
| 302 | } | 
| 303 | }else if (args_info.density_given) { | 
| 304 | analyser= new DensityPlot(info, dumpFileName, sele1, sele2, maxLen, | 
| 305 | args_info.nbins_arg); | 
| 306 | } 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 | } 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 | } else if (args_info.p_angle_given) { | 
| 315 | analyser = new pAngle(info, dumpFileName, sele1, args_info.nbins_arg); | 
| 316 | #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H) | 
| 317 | }else if (args_info.hxy_given) { | 
| 318 | analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg, | 
| 319 | args_info.nbins_y_arg, args_info.nbins_arg); | 
| 320 | #endif | 
| 321 | }else if (args_info.rho_r_given) { | 
| 322 | if (args_info.radius_given){ | 
| 323 | analyser = new RhoR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg,args_info.radius_arg); | 
| 324 | }else{ | 
| 325 | sprintf( painCave.errMsg, | 
| 326 | "A particle radius (radius) must be specified when calculating Rho(r)"); | 
| 327 | painCave.severity = OPENMD_ERROR; | 
| 328 | painCave.isFatal = 1; | 
| 329 | simError(); | 
| 330 | } | 
| 331 | } else if (args_info.hullvol_given) { | 
| 332 | analyser = new NanoVolume(info, dumpFileName, sele1); | 
| 333 | } else if (args_info.rodlength_given) { | 
| 334 | analyser = new NanoLength(info, dumpFileName, sele1); | 
| 335 | } else if (args_info.angle_r_given) { | 
| 336 | analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg); | 
| 337 | } | 
| 338 |  | 
| 339 | if (args_info.output_given) { | 
| 340 | analyser->setOutputName(args_info.output_arg); | 
| 341 | } | 
| 342 | if (args_info.step_given) { | 
| 343 | analyser->setStep(args_info.step_arg); | 
| 344 | } | 
| 345 |  | 
| 346 | analyser->process(); | 
| 347 |  | 
| 348 | delete analyser; | 
| 349 | delete info; | 
| 350 |  | 
| 351 | return 0; | 
| 352 | } |