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

Comparing trunk/src/constraints/Shake.cpp (file contents):
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1989 by gezelter, Fri Apr 18 20:07:09 2014 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 45 | Line 45 | namespace OpenMD {
45   #include "utils/simError.h"
46   namespace OpenMD {
47  
48 <  Shake::Shake(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doShake_(false) {
48 >  Shake::Shake(SimInfo* info) : info_(info), maxConsIteration_(10),
49 >                                consTolerance_(1.0e-6), doShake_(false),
50 >                                currConstraintTime_(0.0)  {
51      
52 <    if (info_->getNConstraints() > 0)
52 >    if (info_->getNGlobalConstraints() > 0)
53        doShake_ = true;
54      
55 +    if (!doShake_) return;
56 +
57 +    Globals* simParams = info_->getSimParams();
58 +
59      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
60 +    if (simParams->haveConstraintTime()){
61 +      constraintTime_ = simParams->getConstraintTime();
62 +    } else {
63 +      constraintTime_ = simParams->getStatusTime();
64 +    }
65 +
66 +    constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) +
67 +      ".constraintForces";
68 +
69 +    // create ConstraintWriter  
70 +    constraintWriter_ = new ConstraintWriter(info_,
71 +                                             constraintOutputFile_.c_str());  
72 +
73 +    if (!constraintWriter_){
74 +      sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n");
75 +      painCave.isFatal = 1;
76 +      simError();
77 +    }
78    }
79    
80    void Shake::constraintR() {
# Line 60 | Line 84 | namespace OpenMD {
84    void Shake::constraintF() {
85      if (!doShake_) return;    
86      doConstraint(&Shake::constraintPairF);
87 +
88 +    if (currentSnapshot_->getTime() >= currConstraintTime_){
89 +      Molecule* mol;
90 +      SimInfo::MoleculeIterator mi;
91 +      ConstraintPair* consPair;
92 +      Molecule::ConstraintPairIterator cpi;
93 +      std::list<ConstraintPair*> constraints;
94 +      for (mol = info_->beginMolecule(mi); mol != NULL;
95 +           mol = info_->nextMolecule(mi)) {
96 +        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
97 +             consPair = mol->nextConstraintPair(cpi)) {
98 +          
99 +          constraints.push_back(consPair);
100 +        }
101 +      }
102 +      
103 +      constraintWriter_->writeConstraintForces(constraints);
104 +      currConstraintTime_ += constraintTime_;
105 +    }
106    }
107  
108    void Shake::doConstraint(ConstraintPairFuncPtr func) {
# Line 72 | Line 115 | namespace OpenMD {
115      ConstraintPair* consPair;
116      Molecule::ConstraintPairIterator cpi;
117      
118 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
119 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
118 >    for (mol = info_->beginMolecule(mi); mol != NULL;
119 >         mol = info_->nextMolecule(mi)) {
120 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
121 >           consElem = mol->nextConstraintElem(cei)) {
122          consElem->setMoved(true);
123          consElem->setMoving(false);
124        }
# Line 87 | Line 132 | namespace OpenMD {
132  
133        //loop over every constraint pair
134  
135 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
136 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
135 >      for (mol = info_->beginMolecule(mi); mol != NULL;
136 >           mol = info_->nextMolecule(mi)) {
137 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
138 >             consPair = mol->nextConstraintPair(cpi)) {
139 >          
140  
93 –
141            //dispatch constraint algorithm
142            if(consPair->isMoved()) {
143              int exeStatus = (this->*func)(consPair);
# Line 98 | Line 145 | namespace OpenMD {
145              switch(exeStatus){
146              case consFail:
147                sprintf(painCave.errMsg,
148 <                      "Constraint failure in Shake::constrainA, Constraint Fail\n");
148 >                      "Constraint failure in Shake::constrainA, "
149 >                      "Constraint Fail\n");
150                painCave.isFatal = 1;
151                simError();                            
152                              
153                break;
154              case consSuccess:
155 <              //constrain the pair by moving two elements
155 >              // constrain the pair by moving two elements
156                done = false;
157                consPair->getConsElem1()->setMoving(true);
158                consPair->getConsElem2()->setMoving(true);
159                break;
160              case consAlready:
161 <              //current pair is already constrained, do not need to move the elements
161 >              // current pair is already constrained, do not need to
162 >              // move the elements
163                break;
164              default:          
165 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
165 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
166 >                      "Error: unrecognized status");
167                painCave.isFatal = 1;
168                simError();                          
169                break;
# Line 123 | Line 173 | namespace OpenMD {
173        }//end for(iter->first())
174  
175  
176 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
177 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
176 >      for (mol = info_->beginMolecule(mi); mol != NULL;
177 >           mol = info_->nextMolecule(mi)) {
178 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
179 >             consElem = mol->nextConstraintElem(cei)) {
180            consElem->setMoved(consElem->getMoving());
181            consElem->setMoving(false);
182          }
# Line 135 | Line 187 | namespace OpenMD {
187  
188      if (!done){
189        sprintf(painCave.errMsg,
190 <              "Constraint failure in Shake::constrainA, too many iterations: %d\n",
190 >              "Constraint failure in Shake::constrainA, "
191 >              "too many iterations: %d\n",
192                iteration);
193        painCave.isFatal = 1;
194        simError();    
# Line 195 | Line 248 | namespace OpenMD {
248        //set atom2's position
249        posB -= rmb * delta;
250        consElem2->setPos(posB);
251 +
252 +      // report the constraint force back to the constraint pair:
253 +      consPair->setConstraintForce(gab);
254        return consSuccess;
255      }
256      else
# Line 238 | Line 294 | namespace OpenMD {
294        consElem1->setFrc(frcA);
295        consElem2->setFrc(frcB);
296  
297 +      // report the constraint force back to the constraint pair:
298 +      consPair->setConstraintForce(gab);
299        return consSuccess;
300      }
301      else

Diff Legend

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