--- trunk/src/applications/staticProps/P2OrderParameter.cpp 2010/05/10 17:28:26 1442 +++ 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" @@ -44,60 +45,31 @@ #include "io/DumpReader.hpp" #include "primitives/Molecule.hpp" #include "utils/NumericConstant.hpp" + +using namespace std; namespace OpenMD { - - P2OrderParameter::P2OrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) - : StaticAnalyser(info, filename), - selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), - seleMan1_(info), seleMan2_(info){ - + 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) { + 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(std::make_pair(sd1, sd2)); - } - + 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) { + setOutputName(getPrefix(filename) + ".p2"); + + evaluator1_.loadScriptString(sele1); + evaluator2_.loadScriptString(sele2); } void P2OrderParameter::process() { @@ -105,37 +77,94 @@ namespace OpenMD { RigidBody* rb; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; - + 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)) { + + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { //change the positions of atoms which belong to the rigidbodies - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { rb->updateAtoms(); - } - + } } Mat3x3d orderTensor(0.0); - for (std::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); + vecCount = 0; + + seleMan1_.setSelectionSet(evaluator1_.evaluate()); + + if (doVect_) { + + 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 /= vecCount; + + } else { + + seleMan2_.setSelectionSet(evaluator2_.evaluate()); + + 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 /= 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; @@ -155,15 +184,36 @@ namespace OpenMD { } RealType angle = 0.0; - for (std::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)) ; + vecCount = 0; + + if (doVect_) { + 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/(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; param.p2 = p2; @@ -173,20 +223,22 @@ namespace OpenMD { orderParams_.push_back(param); } - + writeP2(); - + } void P2OrderParameter::writeP2() { - std::ofstream os(getOutputFileName().c_str()); + ofstream os(getOutputFileName().c_str()); os << "#radial distribution function\n"; os<< "#selection1: (" << selectionScript1_ << ")\t"; - os << "selection2: (" << selectionScript2_ << ")\n"; + if (!doVect_) { + os << "selection2: (" << selectionScript2_ << ")\n"; + } os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n"; - for (std::size_t i = 0; i < orderParams_.size(); ++i) { + for (size_t i = 0; i < orderParams_.size(); ++i) { os << orderParams_[i].p2 << "\t" << orderParams_[i].director[0] << "\t" << orderParams_[i].director[1] << "\t"