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 1983 by gezelter, Tue Apr 15 20:36:19 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 +    Globals* simParams = info_->getSimParams();
56 +
57      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
58 +    if (simParams->haveConstraintTime()){
59 +      constraintTime_ = simParams->getConstraintTime();
60 +    } else {
61 +      constraintTime_ = simParams->getStatusTime();
62 +    }
63 +
64 +    constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) +
65 +      ".constraintForces";
66 +
67 +    // create ConstraintWriter  
68 +    constraintWriter_ = new ConstraintWriter(info_,
69 +                                             constraintOutputFile_.c_str());  
70 +
71 +    if (!constraintWriter_){
72 +      sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n");
73 +      painCave.isFatal = 1;
74 +      simError();
75 +    }
76    }
77    
78    void Shake::constraintR() {
# Line 60 | Line 82 | namespace OpenMD {
82    void Shake::constraintF() {
83      if (!doShake_) return;    
84      doConstraint(&Shake::constraintPairF);
85 +
86 +    if (currentSnapshot_->getTime() >= currConstraintTime_){
87 +      Molecule* mol;
88 +      SimInfo::MoleculeIterator mi;
89 +      ConstraintPair* consPair;
90 +      Molecule::ConstraintPairIterator cpi;
91 +      std::list<ConstraintPair*> constraints;
92 +      for (mol = info_->beginMolecule(mi); mol != NULL;
93 +           mol = info_->nextMolecule(mi)) {
94 +        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
95 +             consPair = mol->nextConstraintPair(cpi)) {
96 +          
97 +          constraints.push_back(consPair);
98 +        }
99 +      }
100 +      
101 +      constraintWriter_->writeConstraintForces(constraints);
102 +      currConstraintTime_ += constraintTime_;
103 +    }
104    }
105  
106    void Shake::doConstraint(ConstraintPairFuncPtr func) {
# Line 72 | Line 113 | namespace OpenMD {
113      ConstraintPair* consPair;
114      Molecule::ConstraintPairIterator cpi;
115      
116 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
117 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
116 >    for (mol = info_->beginMolecule(mi); mol != NULL;
117 >         mol = info_->nextMolecule(mi)) {
118 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
119 >           consElem = mol->nextConstraintElem(cei)) {
120          consElem->setMoved(true);
121          consElem->setMoving(false);
122        }
# Line 87 | Line 130 | namespace OpenMD {
130  
131        //loop over every constraint pair
132  
133 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
134 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
133 >      for (mol = info_->beginMolecule(mi); mol != NULL;
134 >           mol = info_->nextMolecule(mi)) {
135 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
136 >             consPair = mol->nextConstraintPair(cpi)) {
137 >          
138  
93 –
139            //dispatch constraint algorithm
140            if(consPair->isMoved()) {
141              int exeStatus = (this->*func)(consPair);
# Line 98 | Line 143 | namespace OpenMD {
143              switch(exeStatus){
144              case consFail:
145                sprintf(painCave.errMsg,
146 <                      "Constraint failure in Shake::constrainA, Constraint Fail\n");
146 >                      "Constraint failure in Shake::constrainA, "
147 >                      "Constraint Fail\n");
148                painCave.isFatal = 1;
149                simError();                            
150                              
151                break;
152              case consSuccess:
153 <              //constrain the pair by moving two elements
153 >              // constrain the pair by moving two elements
154                done = false;
155                consPair->getConsElem1()->setMoving(true);
156                consPair->getConsElem2()->setMoving(true);
157                break;
158              case consAlready:
159 <              //current pair is already constrained, do not need to move the elements
159 >              // current pair is already constrained, do not need to
160 >              // move the elements
161                break;
162              default:          
163 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
163 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
164 >                      "Error: unrecognized status");
165                painCave.isFatal = 1;
166                simError();                          
167                break;
# Line 123 | Line 171 | namespace OpenMD {
171        }//end for(iter->first())
172  
173  
174 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
175 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
174 >      for (mol = info_->beginMolecule(mi); mol != NULL;
175 >           mol = info_->nextMolecule(mi)) {
176 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
177 >             consElem = mol->nextConstraintElem(cei)) {
178            consElem->setMoved(consElem->getMoving());
179            consElem->setMoving(false);
180          }
# Line 135 | Line 185 | namespace OpenMD {
185  
186      if (!done){
187        sprintf(painCave.errMsg,
188 <              "Constraint failure in Shake::constrainA, too many iterations: %d\n",
188 >              "Constraint failure in Shake::constrainA, "
189 >              "too many iterations: %d\n",
190                iteration);
191        painCave.isFatal = 1;
192        simError();    
# Line 195 | Line 246 | namespace OpenMD {
246        //set atom2's position
247        posB -= rmb * delta;
248        consElem2->setPos(posB);
249 +
250 +      // report the constraint force back to the constraint pair:
251 +      consPair->setConstraintForce(gab);
252        return consSuccess;
253      }
254      else
# Line 238 | Line 292 | namespace OpenMD {
292        consElem1->setFrc(frcA);
293        consElem2->setFrc(frcB);
294  
295 +      // report the constraint force back to the constraint pair:
296 +      consPair->setConstraintForce(gab);
297        return consSuccess;
298      }
299      else

Diff Legend

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