--- trunk/src/constraints/Rattle.cpp 2014/04/05 20:56:01 1979 +++ trunk/src/constraints/Rattle.cpp 2014/04/15 20:36:19 1983 @@ -45,21 +45,43 @@ namespace OpenMD { #include "utils/simError.h" namespace OpenMD { - Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doRattle_(false) { + Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), + consTolerance_(1.0e-6), doRattle_(false), + currConstraintTime_(0.0) { - if (info_->getNConstraints() > 0) + if (info_->getNGlobalConstraints() > 0) doRattle_ = true; + Globals* simParams = info_->getSimParams(); - if (info_->getSimParams()->haveDt()) { - dt_ = info_->getSimParams()->getDt(); + 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() { @@ -69,6 +91,24 @@ namespace OpenMD { void Rattle::constraintB() { 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) {