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 |
|
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(); |
72 |
|
constraintTime_ = simParams->getStatusTime(); |
73 |
|
} |
74 |
|
|
75 |
< |
constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + ".constraintForces"; |
75 |
> |
constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + |
76 |
> |
".constraintForces"; |
77 |
|
|
78 |
+ |
|
79 |
|
// create ConstraintWriter |
80 |
< |
constraintWriter_ = new ConstraintWriter(info_, constraintOutputFile_.c_str()); |
80 |
> |
constraintWriter_ = new ConstraintWriter(info_, |
81 |
> |
constraintOutputFile_.c_str()); |
82 |
|
|
83 |
|
if (!constraintWriter_){ |
84 |
|
sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n"); |
109 |
|
constraints.push_back(consPair); |
110 |
|
} |
111 |
|
} |
104 |
– |
|
112 |
|
constraintWriter_->writeConstraintForces(constraints); |
113 |
|
currConstraintTime_ += constraintTime_; |
114 |
|
} |
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 |
|
|