--- trunk/src/applications/staticProps/P2OrderParameter.cpp 2006/10/18 21:58:48 1078 +++ trunk/src/applications/staticProps/P2OrderParameter.cpp 2011/10/18 13:44:44 1657 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,15 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [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). */ #include "applications/staticProps/P2OrderParameter.hpp" @@ -44,48 +44,65 @@ #include "io/DumpReader.hpp" #include "primitives/Molecule.hpp" #include "utils/NumericConstant.hpp" -namespace oopse { +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); + + 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) { + 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 = OOPSE_ERROR; + 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 = OOPSE_ERROR; + 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 = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); - + } - + int i; int j; StuntDouble* sd1; @@ -94,10 +111,8 @@ namespace oopse { sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) { - sdPairs_.push_back(std::make_pair(sd1, sd2)); - } - - + sdPairs_.push_back(make_pair(sd1, sd2)); + } } void P2OrderParameter::process() { @@ -105,6 +120,8 @@ namespace oopse { RigidBody* rb; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; + StuntDouble* sd; + int i, ii; DumpReader reader(info_, dumpFilename_); int nFrames = reader.getNFrames(); @@ -112,30 +129,54 @@ namespace oopse { 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); + + 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); + vec.normalize(); + orderTensor += outProduct(vec, vec); + } + } + + orderTensor /= seleMan1_.getSelectionCount(); + + } 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(); + orderTensor +=outProduct(vec, vec); + } + + orderTensor /= sdPairs_.size(); } - orderTensor /= sdPairs_.size(); + orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity(); Vector3d eigenvalues; Mat3x3d eigenvectors; + Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors); int which; @@ -155,43 +196,66 @@ namespace oopse { } 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)) ; + RealType p4 = 0.0; + + if (doVect_) { + for (sd = seleMan1_.beginSelected(ii); sd != NULL; + sd = seleMan1_.nextSelected(ii)) { + if (sd->isDirectional()) { + Vector3d vec = sd->getA().getColumn(2); + vec.normalize(); + RealType myDot = dot(vec, director); + angle += acos(myDot); + p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0; + } + } + angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0; + p4 /= (seleMan1_.getSelectionCount()); + } 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(); + RealType myDot = dot(vec, director); + angle += acos(myDot); + p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0; + } + angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0; + p4 /= (seleMan1_.getSelectionCount()); } - angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0; OrderParam param; param.p2 = p2; param.director = director; param.angle = angle; + param.p4 = p4; 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"; - os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n"; + if (!doVect_) { + os << "selection2: (" << selectionScript2_ << ")\n"; + } + os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\t\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" << orderParams_[i].director[2] << "\t" - << orderParams_[i].angle << "\n"; + << orderParams_[i].angle << "\t" + << orderParams_[i].p4 << "\n"; }