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

Comparing:
branches/development/src/constraints/Shake.cpp (file contents), Revision 1748 by gezelter, Wed Jun 6 19:54:09 2012 UTC vs.
trunk/src/constraints/Shake.cpp (file contents), Revision 1982 by gezelter, Tue Apr 15 12:12:23 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 47 | Line 47 | namespace OpenMD {
47  
48    Shake::Shake(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doShake_(false) {
49      
50 <    if (info_->getNConstraints() > 0)
50 >    if (info_->getNGlobalConstraints() > 0)
51        doShake_ = true;
52      
53 +    Globals* simParams = info_->getSimParams();
54 +
55      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
56 +    if (simParams->haveConstraintTime()){
57 +      constraintTime_ = simParams->getConstraintTime();
58 +    } else {
59 +      constraintTime_ = simParams->getStatusTime();
60 +    }
61 +
62 +    constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + ".constraintForces";
63 +
64 +    // create ConstraintWriter  
65 +    constraintWriter_ = new ConstraintWriter(info_, constraintOutputFile_.c_str());  
66 +
67 +    if (!constraintWriter_){
68 +      sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n");
69 +      painCave.isFatal = 1;
70 +      simError();
71 +    }
72    }
73    
74    void Shake::constraintR() {
# Line 60 | Line 78 | namespace OpenMD {
78    void Shake::constraintF() {
79      if (!doShake_) return;    
80      doConstraint(&Shake::constraintPairF);
81 +
82 +    if (currentSnapshot_->getTime() >= currConstraintTime_){
83 +      Molecule* mol;
84 +      SimInfo::MoleculeIterator mi;
85 +      ConstraintPair* consPair;
86 +      Molecule::ConstraintPairIterator cpi;
87 +      std::list<ConstraintPair*> constraints;
88 +      for (mol = info_->beginMolecule(mi); mol != NULL;
89 +           mol = info_->nextMolecule(mi)) {
90 +        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
91 +             consPair = mol->nextConstraintPair(cpi)) {
92 +          
93 +          constraints.push_back(consPair);
94 +        }
95 +      }
96 +      
97 +      constraintWriter_->writeConstraintForces(constraints);
98 +      currConstraintTime_ += constraintTime_;
99 +    }
100    }
101  
102    void Shake::doConstraint(ConstraintPairFuncPtr func) {
# Line 72 | Line 109 | namespace OpenMD {
109      ConstraintPair* consPair;
110      Molecule::ConstraintPairIterator cpi;
111      
112 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
113 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
112 >    for (mol = info_->beginMolecule(mi); mol != NULL;
113 >         mol = info_->nextMolecule(mi)) {
114 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
115 >           consElem = mol->nextConstraintElem(cei)) {
116          consElem->setMoved(true);
117          consElem->setMoving(false);
118        }
# Line 87 | Line 126 | namespace OpenMD {
126  
127        //loop over every constraint pair
128  
129 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
130 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
129 >      for (mol = info_->beginMolecule(mi); mol != NULL;
130 >           mol = info_->nextMolecule(mi)) {
131 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
132 >             consPair = mol->nextConstraintPair(cpi)) {
133 >          
134  
93 –
135            //dispatch constraint algorithm
136            if(consPair->isMoved()) {
137              int exeStatus = (this->*func)(consPair);
# Line 98 | Line 139 | namespace OpenMD {
139              switch(exeStatus){
140              case consFail:
141                sprintf(painCave.errMsg,
142 <                      "Constraint failure in Shake::constrainA, Constraint Fail\n");
142 >                      "Constraint failure in Shake::constrainA, "
143 >                      "Constraint Fail\n");
144                painCave.isFatal = 1;
145                simError();                            
146                              
147                break;
148              case consSuccess:
149 <              //constrain the pair by moving two elements
149 >              // constrain the pair by moving two elements
150                done = false;
151                consPair->getConsElem1()->setMoving(true);
152                consPair->getConsElem2()->setMoving(true);
153                break;
154              case consAlready:
155 <              //current pair is already constrained, do not need to move the elements
155 >              // current pair is already constrained, do not need to
156 >              // move the elements
157                break;
158              default:          
159 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
159 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
160 >                      "Error: unrecognized status");
161                painCave.isFatal = 1;
162                simError();                          
163                break;
# Line 123 | Line 167 | namespace OpenMD {
167        }//end for(iter->first())
168  
169  
170 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
171 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
170 >      for (mol = info_->beginMolecule(mi); mol != NULL;
171 >           mol = info_->nextMolecule(mi)) {
172 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
173 >             consElem = mol->nextConstraintElem(cei)) {
174            consElem->setMoved(consElem->getMoving());
175            consElem->setMoving(false);
176          }
# Line 135 | Line 181 | namespace OpenMD {
181  
182      if (!done){
183        sprintf(painCave.errMsg,
184 <              "Constraint failure in Shake::constrainA, too many iterations: %d\n",
184 >              "Constraint failure in Shake::constrainA, "
185 >              "too many iterations: %d\n",
186                iteration);
187        painCave.isFatal = 1;
188        simError();    
# Line 195 | Line 242 | namespace OpenMD {
242        //set atom2's position
243        posB -= rmb * delta;
244        consElem2->setPos(posB);
245 +
246 +      // report the constraint force back to the constraint pair:
247 +      consPair->setConstraintForce(gab);
248        return consSuccess;
249      }
250      else
# Line 238 | Line 288 | namespace OpenMD {
288        consElem1->setFrc(frcA);
289        consElem2->setFrc(frcB);
290  
291 +      // report the constraint force back to the constraint pair:
292 +      consPair->setConstraintForce(gab);
293        return consSuccess;
294      }
295      else

Diff Legend

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