--- trunk/src/constraints/Rattle.cpp 2006/05/17 21:51:42 963 +++ trunk/src/constraints/Rattle.cpp 2014/09/26 22:22:28 2022 @@ -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,39 +28,95 @@ * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "constraints/Rattle.hpp" #include "primitives/Molecule.hpp" #include "utils/simError.h" -namespace oopse { +namespace OpenMD { - Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6) { + Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), + consTolerance_(1.0e-6), doRattle_(false), + currConstraintTime_(0.0) { - if (info_->getSimParams()->haveDt()) { - dt_ = info_->getSimParams()->getDt(); + if (info_->getNGlobalConstraints() > 0) + doRattle_ = true; + + if (!doRattle_) return; + + Globals* simParams = info_->getSimParams(); + + if (simParams->haveDt()) { + dt_ = simParams->getDt(); } else { sprintf(painCave.errMsg, - "Integrator Error: dt is not set\n"); + "Rattle Error: dt is not set\n"); painCave.isFatal = 1; simError(); } - + currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); + if (simParams->haveConstraintTime()){ + constraintTime_ = simParams->getConstraintTime(); + } else { + constraintTime_ = simParams->getStatusTime(); + } + + constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + + ".constraintForces"; + + + // create ConstraintWriter + constraintWriter_ = new ConstraintWriter(info_, + constraintOutputFile_.c_str()); + + if (!constraintWriter_){ + sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n"); + painCave.isFatal = 1; + simError(); + } } void Rattle::constraintA() { - if (info_->getNConstraints() > 0) { - doConstraint(&Rattle::constraintPairA); - } + if (!doRattle_) return; + doConstraint(&Rattle::constraintPairA); } void Rattle::constraintB() { - if (info_->getNConstraints() > 0) { - doConstraint(&Rattle::constraintPairB); + if (!doRattle_) return; + doConstraint(&Rattle::constraintPairB); + + if (currentSnapshot_->getTime() >= currConstraintTime_){ + Molecule* mol; + SimInfo::MoleculeIterator mi; + ConstraintPair* consPair; + Molecule::ConstraintPairIterator cpi; + std::list constraints; + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; + consPair = mol->nextConstraintPair(cpi)) { + + constraints.push_back(consPair); + } + } + constraintWriter_->writeConstraintForces(constraints); + currConstraintTime_ += constraintTime_; } } void Rattle::doConstraint(ConstraintPairFuncPtr func) { + if (!doRattle_) return; + Molecule* mol; SimInfo::MoleculeIterator mi; ConstraintElem* consElem; @@ -77,8 +124,10 @@ namespace oopse { ConstraintPair* consPair; Molecule::ConstraintPairIterator cpi; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (consElem = mol->beginConstraintElem(cei); consElem != NULL; + consElem = mol->nextConstraintElem(cei)) { consElem->setMoved(true); consElem->setMoving(false); } @@ -92,8 +141,10 @@ namespace oopse { //loop over every constraint pair - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; + consPair = mol->nextConstraintPair(cpi)) { //dispatch constraint algorithm @@ -103,22 +154,25 @@ namespace oopse { switch(exeStatus){ case consFail: sprintf(painCave.errMsg, - "Constraint failure in Rattle::constrainA, Constraint Fail\n"); + "Constraint failure in Rattle::constrainA, " + "Constraint Fail\n"); painCave.isFatal = 1; simError(); break; case consSuccess: - //constrain the pair by moving two elements + // constrain the pair by moving two elements done = false; consPair->getConsElem1()->setMoving(true); consPair->getConsElem2()->setMoving(true); break; case consAlready: - //current pair is already constrained, do not need to move the elements + // current pair is already constrained, do not need to + // move the elements break; default: - sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status"); + sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() " + "Error: unrecognized status"); painCave.isFatal = 1; simError(); break; @@ -128,8 +182,10 @@ namespace oopse { }//end for(iter->first()) - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (consElem = mol->beginConstraintElem(cei); consElem != NULL; + consElem = mol->nextConstraintElem(cei)) { consElem->setMoved(consElem->getMoving()); consElem->setMoving(false); } @@ -140,7 +196,8 @@ namespace oopse { if (!done){ sprintf(painCave.errMsg, - "Constraint failure in Rattle::constrainA, too many iterations: %d\n", + "Constraint failure in Rattle::constrainA, " + "too many iterations: %d\n", iteration); painCave.isFatal = 1; simError(); @@ -148,6 +205,7 @@ namespace oopse { } int Rattle::constraintPairA(ConstraintPair* consPair){ + ConstraintElem* consElem1 = consPair->getConsElem1(); ConstraintElem* consElem2 = consPair->getConsElem2(); @@ -165,13 +223,14 @@ namespace oopse { RealType rabsq = consPair->getConsDistSquare(); RealType diffsq = rabsq - pabsq; + // the original rattle code from alan tidesley if (fabs(diffsq) > (consTolerance_ * rabsq * 2)){ Vector3d oldPosA = consElem1->getPrevPos(); Vector3d oldPosB = consElem2->getPrevPos(); - Vector3d rab = oldPosA - oldPosB; + Vector3d rab = oldPosA - oldPosB; currentSnapshot_->wrapVector(rab); @@ -208,7 +267,9 @@ namespace oopse { Vector3d velB = consElem2->getVel(); velB -= rmb * delta; consElem2->setVel(velB); - + + // report the constraint force back to the constraint pair: + consPair->setConstraintForce(gab); return consSuccess; } else @@ -250,6 +311,8 @@ namespace oopse { velB -= rmb * delta; consElem2->setVel(velB); + // report the constraint force back to the constraint pair: + consPair->setConstraintForce(gab); return consSuccess; } else