| 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 | < | * [4]  Vardeman & Gezelter, in progress (2009). | 
| 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 |  |  | 
| 46 | – | #include "brains/Register.hpp" | 
| 47 |  | #include "brains/SimCreator.hpp" | 
| 48 |  | #include "brains/SimInfo.hpp" | 
| 49 |  | #include "io/DumpReader.hpp" | 
| 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 | < | //register force fields | 
| 78 | < | registerForceFields(); | 
| 79 | < |  | 
| 86 | > |  | 
| 87 |  | gengetopt_args_info args_info; | 
| 88 | < |  | 
| 88 | > |  | 
| 89 |  | //parse the command line option | 
| 90 |  | if (cmdline_parser (argc, argv, &args_info) != 0) { | 
| 91 |  | exit(1) ; | 
| 92 |  | } | 
| 93 | < |  | 
| 93 | > |  | 
| 94 |  | //get the dumpfile name | 
| 95 |  | std::string dumpFileName = args_info.input_arg; | 
| 96 |  | std::string sele1; | 
| 97 |  | std::string sele2; | 
| 98 | < | bool userSpecifiedSelect1; | 
| 92 | < | bool userSpecifiedSelect2; | 
| 93 | < |  | 
| 98 | > |  | 
| 99 |  | // check the first selection argument, or set it to the environment | 
| 100 |  | // variable, or failing that, set it to "select all" | 
| 101 | < |  | 
| 101 | > |  | 
| 102 |  | if (args_info.sele1_given) { | 
| 103 |  | sele1 = args_info.sele1_arg; | 
| 104 |  | } else { | 
| 109 |  | sele1 = "select all"; | 
| 110 |  | } | 
| 111 |  | } | 
| 112 | < |  | 
| 112 | > |  | 
| 113 |  | // check the second selection argument, or set it to the environment | 
| 114 | < | // variable, or failing that, set it to "select all" | 
| 114 | > | // variable, or failing that, set it to the first selection | 
| 115 |  |  | 
| 116 |  | if (args_info.sele2_given) { | 
| 117 |  | sele2 = args_info.sele2_arg; | 
| 118 |  | } else { | 
| 119 | < | char* sele2Env = getenv("SELECTION1"); | 
| 119 | > | char* sele2Env = getenv("SELECTION2"); | 
| 120 |  | if (sele2Env) { | 
| 121 |  | sele2 = sele2Env; | 
| 122 |  | } else { | 
| 123 | < | sele2 = "select all"; | 
| 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 |  |  | 
| 122 | – |  | 
| 123 | – | // Problems if sele1 wasn't specified, but | 
| 124 | – | // if (!args_info.scd_given) { | 
| 125 | – | //       sprintf( painCave.errMsg, | 
| 126 | – | //                "neither --sele1 option nor $SELECTION1 is set"); | 
| 127 | – | //       painCave.severity = OPENMD_ERROR; | 
| 128 | – | //       painCave.isFatal = 1; | 
| 129 | – | //       simError(); | 
| 130 | – | //     } | 
| 131 | – | //   } | 
| 132 | – |  | 
| 133 | – | // Problems if sele1 wasn't specified | 
| 134 | – |  | 
| 135 | – | //     if(!args_info.scd_given && !args_info.density_given && !args_info.slab_density_given)  { | 
| 136 | – | //       sprintf( painCave.errMsg, | 
| 137 | – | //                "neither --sele2 option nor $SELECTION1 is set"); | 
| 138 | – | //       painCave.severity = OPENMD_ERROR; | 
| 139 | – | //       painCave.isFatal = 1; | 
| 140 | – | //       simError(); | 
| 141 | – | //     } | 
| 142 | – | //   } | 
| 143 | – |  | 
| 129 |  | bool batchMode; | 
| 130 |  | if (args_info.scd_given){ | 
| 131 | < | if (args_info.sele1_given && args_info.sele2_given && args_info.sele3_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 && args_info.begin_given && args_info.end_given) { | 
| 135 | < | if (args_info.begin_arg < 0 || args_info.end_arg < 0 || args_info.begin_arg > args_info.end_arg-2) { | 
| 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"); | 
| 149 |  | " or --molname, --begin, --end are specified\n"); | 
| 150 |  | painCave.severity = OPENMD_ERROR; | 
| 151 |  | painCave.isFatal = 1; | 
| 152 | < | simError(); | 
| 165 | < |  | 
| 152 | > | simError(); | 
| 153 |  | } | 
| 154 |  | } | 
| 155 | < |  | 
| 155 | > |  | 
| 156 |  | //parse md file and set up the system | 
| 157 |  | SimCreator creator; | 
| 171 | – | std::cout << "dumpFile = " << dumpFileName << "\n"; | 
| 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; | 
| 169 | > | maxLen = std::min(std::min(hmat(0, 0), hmat(1, 1)), hmat(2, 2)) /2.0; | 
| 170 | > | zmaxLen = hmat(2,2); | 
| 171 |  | } | 
| 172 | < |  | 
| 172 | > |  | 
| 173 |  | StaticAnalyser* analyser; | 
| 174 |  | if (args_info.gofr_given){ | 
| 175 |  | analyser= new GofR(info, dumpFileName, sele1, sele2, maxLen, | 
| 176 | < | args_info.nbins_arg); | 
| 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); | 
| 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); | 
| 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); | 
| 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); | 
| 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,args_info.refsele_arg, | 
| 195 | < | maxLen, args_info.nbins_arg); | 
| 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"); | 
| 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 | < | analyser  = new P2OrderParameter(info, dumpFileName, sele1, sele2); | 
| 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); | 
| 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"); | 
| 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"); | 
| 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, args_info.rcut_arg, | 
| 284 | < | args_info.nbins_arg); | 
| 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"); | 
| 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, args_info.molname_arg, | 
| 296 | < | args_info.begin_arg, args_info.end_arg); | 
| 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, sele1, sele2, sele3); | 
| 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); | 
| 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.rnemdrt_given) { | 
| 315 | + | analyser = new RNEMDRTheta(info, dumpFileName, sele1, | 
| 316 | + | args_info.nbins_arg, args_info.nanglebins_arg); | 
| 317 | + | } else if (args_info.p_angle_given) { | 
| 318 | + | 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 | + | if (args_info.seleoffset2_given) { | 
| 325 | + | analyser  = new pAngle(info, dumpFileName, sele1, | 
| 326 | + | args_info.seleoffset_arg, | 
| 327 | + | args_info.seleoffset2_arg, | 
| 328 | + | args_info.nbins_arg); | 
| 329 | + | } else { | 
| 330 | + | analyser  = new pAngle(info, dumpFileName, sele1, | 
| 331 | + | args_info.seleoffset_arg, | 
| 332 | + | args_info.nbins_arg); | 
| 333 | + | } | 
| 334 | + | } else | 
| 335 | + | analyser  = new pAngle(info, dumpFileName, sele1, | 
| 336 | + | args_info.nbins_arg); | 
| 337 | + | } else { | 
| 338 | + | sprintf( painCave.errMsg, | 
| 339 | + | "At least one selection script (--sele1) must be specified when " | 
| 340 | + | "calculating P(angle) distributions"); | 
| 341 | + | painCave.severity = OPENMD_ERROR; | 
| 342 | + | painCave.isFatal = 1; | 
| 343 | + | simError(); | 
| 344 | + | } | 
| 345 |  | #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H) | 
| 346 |  | }else if (args_info.hxy_given) { | 
| 347 |  | analyser = new Hxy(info, dumpFileName, sele1, args_info.nbins_x_arg, | 
| 348 | < | args_info.nbins_y_arg, args_info.nbins_arg); | 
| 348 | > | args_info.nbins_y_arg, args_info.nbins_arg); | 
| 349 |  | #endif | 
| 350 |  | }else if (args_info.rho_r_given) { | 
| 351 |  | if (args_info.radius_given){ | 
| 352 |  | analyser = new RhoR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg,args_info.radius_arg); | 
| 353 |  | }else{ | 
| 354 |  | sprintf( painCave.errMsg, | 
| 355 | < | "A particle radius (radius) must be specified when calculating Rho(r)"); | 
| 355 | > | "A particle radius (radius) must be specified when calculating Rho(r)"); | 
| 356 |  | painCave.severity = OPENMD_ERROR; | 
| 357 |  | painCave.isFatal = 1; | 
| 358 |  | simError(); | 
| 359 |  | } | 
| 360 | < | }else if (args_info.hullvol_given) { | 
| 360 | > | } else if (args_info.hullvol_given) { | 
| 361 |  | analyser = new NanoVolume(info, dumpFileName, sele1); | 
| 362 | + | } else if (args_info.rodlength_given) { | 
| 363 | + | analyser = new NanoLength(info, dumpFileName, sele1); | 
| 364 | + | } else if (args_info.angle_r_given) { | 
| 365 | + | analyser = new AngleR(info, dumpFileName, sele1, maxLen,args_info.nbins_arg); | 
| 366 |  | } | 
| 367 |  |  | 
| 368 |  | if (args_info.output_given) { | 
| 371 |  | if (args_info.step_given) { | 
| 372 |  | analyser->setStep(args_info.step_arg); | 
| 373 |  | } | 
| 374 | < |  | 
| 374 | > |  | 
| 375 |  | analyser->process(); | 
| 376 | < |  | 
| 376 | > |  | 
| 377 |  | delete analyser; | 
| 378 |  | delete info; | 
| 379 | < |  | 
| 379 | > |  | 
| 380 |  | return 0; | 
| 381 |  | } | 
| 293 | – |  |