ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/constraints/Rattle.cpp
(Generate patch)

Comparing trunk/src/constraints/Rattle.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2022 by gezelter, Fri Sep 26 22:22:28 2014 UTC

# Line 45 | Line 45 | namespace OpenMD {
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 <    
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() {
# Line 70 | Line 94 | namespace OpenMD {
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) {
# Line 82 | Line 124 | namespace OpenMD {
124      ConstraintPair* consPair;
125      Molecule::ConstraintPairIterator cpi;
126      
127 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
128 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
127 >    for (mol = info_->beginMolecule(mi); mol != NULL;
128 >         mol = info_->nextMolecule(mi)) {
129 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
130 >           consElem = mol->nextConstraintElem(cei)) {
131          consElem->setMoved(true);
132          consElem->setMoving(false);
133        }
# Line 97 | Line 141 | namespace OpenMD {
141  
142        //loop over every constraint pair
143  
144 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
145 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
144 >      for (mol = info_->beginMolecule(mi); mol != NULL;
145 >           mol = info_->nextMolecule(mi)) {
146 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
147 >             consPair = mol->nextConstraintPair(cpi)) {
148  
149  
150            //dispatch constraint algorithm
# Line 108 | Line 154 | namespace OpenMD {
154              switch(exeStatus){
155              case consFail:
156                sprintf(painCave.errMsg,
157 <                      "Constraint failure in Rattle::constrainA, Constraint Fail\n");
157 >                      "Constraint failure in Rattle::constrainA, "
158 >                      "Constraint Fail\n");
159                painCave.isFatal = 1;
160                simError();                            
161                              
162                break;
163              case consSuccess:
164 <              //constrain the pair by moving two elements
164 >              // constrain the pair by moving two elements
165                done = false;
166                consPair->getConsElem1()->setMoving(true);
167                consPair->getConsElem2()->setMoving(true);
168                break;
169              case consAlready:
170 <              //current pair is already constrained, do not need to move the elements
170 >              // current pair is already constrained, do not need to
171 >              // move the elements
172                break;
173              default:          
174 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
174 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
175 >                      "Error: unrecognized status");
176                painCave.isFatal = 1;
177                simError();                          
178                break;
# Line 133 | Line 182 | namespace OpenMD {
182        }//end for(iter->first())
183  
184  
185 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
186 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
185 >      for (mol = info_->beginMolecule(mi); mol != NULL;
186 >           mol = info_->nextMolecule(mi)) {
187 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
188 >             consElem = mol->nextConstraintElem(cei)) {
189            consElem->setMoved(consElem->getMoving());
190            consElem->setMoving(false);
191          }
# Line 145 | Line 196 | namespace OpenMD {
196  
197      if (!done){
198        sprintf(painCave.errMsg,
199 <              "Constraint failure in Rattle::constrainA, too many iterations: %d\n",
199 >              "Constraint failure in Rattle::constrainA, "
200 >              "too many iterations: %d\n",
201                iteration);
202        painCave.isFatal = 1;
203        simError();    
# Line 171 | Line 223 | namespace OpenMD {
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  
# Line 214 | Line 267 | namespace OpenMD {
267        Vector3d velB = consElem2->getVel();
268        velB -= rmb * delta;
269        consElem2->setVel(velB);
270 <
270 >      
271 >      // report the constraint force back to the constraint pair:
272 >      consPair->setConstraintForce(gab);
273        return consSuccess;
274      }
275      else
# Line 256 | Line 311 | namespace OpenMD {
311        velB -= rmb * delta;
312        consElem2->setVel(velB);
313  
314 +      // report the constraint force back to the constraint pair:
315 +      consPair->setConstraintForce(gab);
316        return consSuccess;
317      }
318      else

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines