--- trunk/src/applications/staticProps/P2OrderParameter.cpp 2011/05/11 16:32:48 1558 +++ trunk/src/applications/staticProps/P2OrderParameter.cpp 2012/12/12 20:37:29 1818 @@ -36,7 +36,8 @@ * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "applications/staticProps/P2OrderParameter.hpp" @@ -57,10 +58,6 @@ namespace OpenMD { setOutputName(getPrefix(filename) + ".p2"); evaluator1_.loadScriptString(sele1); - - if (!evaluator1_.isDynamic()) { - seleMan1_.setSelectionSet(evaluator1_.evaluate()); - } } P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename, @@ -72,47 +69,7 @@ namespace OpenMD { setOutputName(getPrefix(filename) + ".p2"); evaluator1_.loadScriptString(sele1); - evaluator2_.loadScriptString(sele2); - - if (!evaluator1_.isDynamic()) { - seleMan1_.setSelectionSet(evaluator1_.evaluate()); - }else { - sprintf( painCave.errMsg, - "--sele1 must be static selection\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - if (!evaluator2_.isDynamic()) { - seleMan2_.setSelectionSet(evaluator2_.evaluate()); - }else { - sprintf( painCave.errMsg, - "--sele2 must be static selection\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) { - sprintf( painCave.errMsg, - "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - - } - - int i; - int j; - StuntDouble* sd1; - StuntDouble* sd2; - for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j); - sd1 != NULL && sd2 != NULL; - sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) { - - sdPairs_.push_back(make_pair(sd1, sd2)); - } + evaluator2_.loadScriptString(sele2); } void P2OrderParameter::process() { @@ -120,9 +77,12 @@ namespace OpenMD { RigidBody* rb; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; - StuntDouble* sd; - int i, ii; - + StuntDouble* sd1; + StuntDouble* sd2; + int ii; + int jj; + int vecCount; + DumpReader reader(info_, dumpFilename_); int nFrames = reader.getNFrames(); @@ -140,37 +100,65 @@ namespace OpenMD { } Mat3x3d orderTensor(0.0); + vecCount = 0; + seleMan1_.setSelectionSet(evaluator1_.evaluate()); + if (doVect_) { - if (evaluator1_.isDynamic()) - seleMan1_.setSelectionSet(evaluator1_.evaluate()); - - for (sd = seleMan1_.beginSelected(ii); sd != NULL; - sd = seleMan1_.nextSelected(ii)) { - if (sd->isDirectional()) { - Vector3d vec = sd->getA().getColumn(2); + for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; + sd1 = seleMan1_.nextSelected(ii)) { + if (sd1->isDirectional()) { + Vector3d vec = sd1->getA().getColumn(2); vec.normalize(); orderTensor += outProduct(vec, vec); + vecCount++; } } - orderTensor /= seleMan1_.getSelectionCount(); + orderTensor /= vecCount; } else { + + seleMan2_.setSelectionSet(evaluator2_.evaluate()); - for (vector >::iterator j = sdPairs_.begin(); - j != sdPairs_.end(); ++j) { - Vector3d vec = j->first->getPos() - j->second->getPos(); + if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) { + sprintf( painCave.errMsg, + "In frame %d, the number of selected StuntDoubles are not\n" + "\tthe same in --sele1 and sele2\n", i); + painCave.severity = OPENMD_INFO; + painCave.isFatal = 0; + simError(); + } + + for (sd1 = seleMan1_.beginSelected(ii), + sd2 = seleMan2_.beginSelected(jj); + sd1 != NULL && sd2 != NULL; + sd1 = seleMan1_.nextSelected(ii), + sd2 = seleMan2_.nextSelected(jj)) { + + Vector3d vec = sd1->getPos() - sd2->getPos(); + if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(vec); + vec.normalize(); + orderTensor +=outProduct(vec, vec); + vecCount++; } - orderTensor /= sdPairs_.size(); + orderTensor /= vecCount; } + if (vecCount == 0) { + sprintf( painCave.errMsg, + "In frame %d, the number of selected vectors was zero.\n" + "\tThis will not give a meaningful order parameter.", i); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity(); @@ -196,27 +184,35 @@ namespace OpenMD { } RealType angle = 0.0; + vecCount = 0; if (doVect_) { - for (sd = seleMan1_.beginSelected(ii); sd != NULL; - sd = seleMan1_.nextSelected(ii)) { - if (sd->isDirectional()) { - Vector3d vec = sd->getA().getColumn(2); + for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; + sd1 = seleMan1_.nextSelected(ii)) { + if (sd1->isDirectional()) { + Vector3d vec = sd1->getA().getColumn(2); vec.normalize(); angle += acos(dot(vec, director)); + vecCount++; } } - angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0; + angle = angle/(vecCount*NumericConstant::PI)*180.0; } else { - for (vector >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) { - Vector3d vec = j->first->getPos() - j->second->getPos(); + for (sd1 = seleMan1_.beginSelected(ii), + sd2 = seleMan2_.beginSelected(jj); + sd1 != NULL && sd2 != NULL; + sd1 = seleMan1_.nextSelected(ii), + sd2 = seleMan2_.nextSelected(jj)) { + + Vector3d vec = sd1->getPos() - sd2->getPos(); if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(vec); vec.normalize(); angle += acos(dot(vec, director)) ; + vecCount++; } - angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0; + angle = angle / (vecCount * NumericConstant::PI) * 180.0; } OrderParam param;