--- trunk/src/applications/staticProps/P2OrderParameter.cpp 2011/03/03 20:32:49 1542 +++ trunk/src/applications/staticProps/P2OrderParameter.cpp 2012/12/13 16:57:39 1819 @@ -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" @@ -50,70 +51,36 @@ namespace OpenMD { P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename, const string& sele1) - : StaticAnalyser(info, filename), doVect_(true), - selectionScript1_(sele1), evaluator1_(info), - evaluator2_(info), seleMan1_(info), seleMan2_(info) { + : StaticAnalyser(info, filename), doVect_(true), doOffset_(false), + selectionScript1_(sele1), evaluator1_(info), + evaluator2_(info), seleMan1_(info), seleMan2_(info) { setOutputName(getPrefix(filename) + ".p2"); evaluator1_.loadScriptString(sele1); - - if (!evaluator1_.isDynamic()) { - seleMan1_.setSelectionSet(evaluator1_.evaluate()); - } - } P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename, const string& sele1, const string& sele2) - : StaticAnalyser(info, filename), doVect_(false), - selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), - evaluator2_(info), seleMan1_(info), seleMan2_(info) { + : StaticAnalyser(info, filename), doVect_(false), doOffset_(false), + selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), + evaluator2_(info), seleMan1_(info), seleMan2_(info) { setOutputName(getPrefix(filename) + ".p2"); evaluator1_.loadScriptString(sele1); - evaluator2_.loadScriptString(sele2); + evaluator2_.loadScriptString(sele2); + } + + P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename, + const string& sele1, int seleOffset) + : StaticAnalyser(info, filename), doVect_(false), doOffset_(true), + seleOffset_(seleOffset), selectionScript1_(sele1), evaluator1_(info), + evaluator2_(info), seleMan1_(info), seleMan2_(info) { - 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(); - } + setOutputName(getPrefix(filename) + ".p2"); - 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)); - } + evaluator1_.loadScriptString(sele1); } void P2OrderParameter::process() { @@ -121,17 +88,19 @@ namespace OpenMD { RigidBody* rb; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; - StuntDouble* sd; - int i; - - + StuntDouble* sd1; + StuntDouble* sd2; + int ii; + int jj; + int vecCount; + DumpReader reader(info_, dumpFilename_); int nFrames = reader.getNFrames(); for (int i = 0; i < nFrames; i += step_) { reader.readFrame(i); currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); - + for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { //change the positions of atoms which belong to the rigidbodies @@ -142,41 +111,98 @@ 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(i); sd != NULL; - sd = seleMan1_.nextSelected(i)) { - 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 { + + if (doOffset_) { + + for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; + sd1 = seleMan1_.nextSelected(ii)) { + + // This will require careful rewriting if StaticProps is + // ever parallelized. For an example, see + // Thermo::getTaggedAtomPairDistance + + int sd2Index = sd1->getGlobalIndex() + seleOffset_; + sd2 = info_->getIOIndexToIntegrableObject(sd2Index); - for (vector >::iterator j = sdPairs_.begin(); - j != sdPairs_.end(); ++j) { - Vector3d vec = j->first->getPos() - j->second->getPos(); - if (usePeriodicBoundaryConditions_) - currentSnapshot_->wrapVector(vec); - vec.normalize(); - orderTensor +=outProduct(vec, vec); + Vector3d vec = sd1->getPos() - sd2->getPos(); + + if (usePeriodicBoundaryConditions_) + currentSnapshot_->wrapVector(vec); + + vec.normalize(); + + orderTensor +=outProduct(vec, vec); + vecCount++; + } + + orderTensor /= vecCount; + } else { + + seleMan2_.setSelectionSet(evaluator2_.evaluate()); + + if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) { + sprintf( painCave.errMsg, + "In frame %d, the number of selected StuntDoubles are\n" + "\tnot the 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 /= vecCount; } - - orderTensor /= sdPairs_.size(); } + 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(); Vector3d eigenvalues; Mat3x3d eigenvectors; + Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors); int which; @@ -196,29 +222,59 @@ namespace OpenMD { } RealType angle = 0.0; - + vecCount = 0; if (doVect_) { - for (sd = seleMan1_.beginSelected(i); sd != NULL; - sd = seleMan1_.nextSelected(i)) { - 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(); - if (usePeriodicBoundaryConditions_) - currentSnapshot_->wrapVector(vec); - vec.normalize(); - - angle += acos(dot(vec, director)) ; + if (doOffset_) { + + for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; + sd1 = seleMan1_.nextSelected(ii)) { + + // This will require careful rewriting if StaticProps is + // ever parallelized. For an example, see + // Thermo::getTaggedAtomPairDistance + + int sd2Index = sd1->getGlobalIndex() + seleOffset_; + sd2 = info_->getIOIndexToIntegrableObject(sd2Index); + + Vector3d vec = sd1->getPos() - sd2->getPos(); + if (usePeriodicBoundaryConditions_) + currentSnapshot_->wrapVector(vec); + vec.normalize(); + angle += acos(dot(vec, director)) ; + vecCount++; + } + angle = angle / (vecCount * NumericConstant::PI) * 180.0; + + } else { + + 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 / (vecCount * NumericConstant::PI) * 180.0; } - angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0; } OrderParam param;