| 45 |
|
#include "utils/simError.h" |
| 46 |
|
namespace OpenMD { |
| 47 |
|
|
| 48 |
< |
Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doRattle_(false) { |
| 48 |
> |
Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), |
| 49 |
> |
consTolerance_(1.0e-6), doRattle_(false), |
| 50 |
> |
currConstraintTime_(0.0) { |
| 51 |
|
|
| 52 |
< |
if (info_->getNConstraints() > 0) |
| 52 |
> |
if (info_->getNGlobalConstraints() > 0) |
| 53 |
|
doRattle_ = true; |
| 54 |
|
|
| 55 |
+ |
if (!doRattle_) return; |
| 56 |
|
|
| 57 |
< |
if (info_->getSimParams()->haveDt()) { |
| 58 |
< |
dt_ = info_->getSimParams()->getDt(); |
| 57 |
> |
Globals* simParams = info_->getSimParams(); |
| 58 |
> |
|
| 59 |
> |
if (simParams->haveDt()) { |
| 60 |
> |
dt_ = simParams->getDt(); |
| 61 |
|
} else { |
| 62 |
|
sprintf(painCave.errMsg, |
| 63 |
< |
"Integrator Error: dt is not set\n"); |
| 63 |
> |
"Rattle Error: dt is not set\n"); |
| 64 |
|
painCave.isFatal = 1; |
| 65 |
|
simError(); |
| 66 |
|
} |
| 67 |
+ |
|
| 68 |
|
currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 69 |
+ |
if (simParams->haveConstraintTime()){ |
| 70 |
+ |
constraintTime_ = simParams->getConstraintTime(); |
| 71 |
+ |
} else { |
| 72 |
+ |
constraintTime_ = simParams->getStatusTime(); |
| 73 |
+ |
} |
| 74 |
+ |
|
| 75 |
+ |
constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + |
| 76 |
+ |
".constraintForces"; |
| 77 |
+ |
|
| 78 |
+ |
|
| 79 |
+ |
// create ConstraintWriter |
| 80 |
+ |
constraintWriter_ = new ConstraintWriter(info_, |
| 81 |
+ |
constraintOutputFile_.c_str()); |
| 82 |
+ |
|
| 83 |
+ |
if (!constraintWriter_){ |
| 84 |
+ |
sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n"); |
| 85 |
+ |
painCave.isFatal = 1; |
| 86 |
+ |
simError(); |
| 87 |
+ |
} |
| 88 |
|
} |
| 89 |
|
|
| 90 |
|
void Rattle::constraintA() { |
| 94 |
|
void Rattle::constraintB() { |
| 95 |
|
if (!doRattle_) return; |
| 96 |
|
doConstraint(&Rattle::constraintPairB); |
| 97 |
+ |
|
| 98 |
+ |
if (currentSnapshot_->getTime() >= currConstraintTime_){ |
| 99 |
+ |
Molecule* mol; |
| 100 |
+ |
SimInfo::MoleculeIterator mi; |
| 101 |
+ |
ConstraintPair* consPair; |
| 102 |
+ |
Molecule::ConstraintPairIterator cpi; |
| 103 |
+ |
std::list<ConstraintPair*> constraints; |
| 104 |
+ |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 105 |
+ |
mol = info_->nextMolecule(mi)) { |
| 106 |
+ |
for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; |
| 107 |
+ |
consPair = mol->nextConstraintPair(cpi)) { |
| 108 |
+ |
|
| 109 |
+ |
constraints.push_back(consPair); |
| 110 |
+ |
} |
| 111 |
+ |
} |
| 112 |
+ |
constraintWriter_->writeConstraintForces(constraints); |
| 113 |
+ |
currConstraintTime_ += constraintTime_; |
| 114 |
+ |
} |
| 115 |
|
} |
| 116 |
|
|
| 117 |
|
void Rattle::doConstraint(ConstraintPairFuncPtr func) { |
| 223 |
|
RealType rabsq = consPair->getConsDistSquare(); |
| 224 |
|
RealType diffsq = rabsq - pabsq; |
| 225 |
|
|
| 226 |
+ |
|
| 227 |
|
// the original rattle code from alan tidesley |
| 228 |
|
if (fabs(diffsq) > (consTolerance_ * rabsq * 2)){ |
| 229 |
|
|
| 230 |
|
Vector3d oldPosA = consElem1->getPrevPos(); |
| 231 |
|
Vector3d oldPosB = consElem2->getPrevPos(); |
| 232 |
|
|
| 233 |
< |
Vector3d rab = oldPosA - oldPosB; |
| 233 |
> |
Vector3d rab = oldPosA - oldPosB; |
| 234 |
|
|
| 235 |
|
currentSnapshot_->wrapVector(rab); |
| 236 |
|
|
