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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1979 by gezelter, Sat Apr 5 20:56:01 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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   */
42  
43   #include "constraints/Rattle.hpp"
# Line 44 | 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) {
48 >  Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doRattle_(false) {
49      
50 +    if (info_->getNConstraints() > 0)
51 +      doRattle_ = true;
52 +    
53 +
54      if (info_->getSimParams()->haveDt()) {
55        dt_ = info_->getSimParams()->getDt();
56      } else {
# Line 54 | Line 59 | namespace OpenMD {
59        painCave.isFatal = 1;
60        simError();
61      }    
57 –    
62      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
63    }
64  
65    void Rattle::constraintA() {
66 <    if (info_->getNConstraints() > 0) {
67 <      doConstraint(&Rattle::constraintPairA);
64 <    }
66 >    if (!doRattle_) return;
67 >    doConstraint(&Rattle::constraintPairA);
68    }
69    void Rattle::constraintB() {
70 <    if (info_->getNConstraints() > 0) {
71 <      doConstraint(&Rattle::constraintPairB);
69 <    }
70 >    if (!doRattle_) return;    
71 >    doConstraint(&Rattle::constraintPairB);
72    }
73  
74    void Rattle::doConstraint(ConstraintPairFuncPtr func) {
75 +    if (!doRattle_) return;
76 +
77      Molecule* mol;
78      SimInfo::MoleculeIterator mi;
79      ConstraintElem* consElem;
# Line 77 | Line 81 | namespace OpenMD {
81      ConstraintPair* consPair;
82      Molecule::ConstraintPairIterator cpi;
83      
84 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
85 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
84 >    for (mol = info_->beginMolecule(mi); mol != NULL;
85 >         mol = info_->nextMolecule(mi)) {
86 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
87 >           consElem = mol->nextConstraintElem(cei)) {
88          consElem->setMoved(true);
89          consElem->setMoving(false);
90        }
# Line 92 | Line 98 | namespace OpenMD {
98  
99        //loop over every constraint pair
100  
101 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
102 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
101 >      for (mol = info_->beginMolecule(mi); mol != NULL;
102 >           mol = info_->nextMolecule(mi)) {
103 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
104 >             consPair = mol->nextConstraintPair(cpi)) {
105  
106  
107            //dispatch constraint algorithm
# Line 103 | Line 111 | namespace OpenMD {
111              switch(exeStatus){
112              case consFail:
113                sprintf(painCave.errMsg,
114 <                      "Constraint failure in Rattle::constrainA, Constraint Fail\n");
114 >                      "Constraint failure in Rattle::constrainA, "
115 >                      "Constraint Fail\n");
116                painCave.isFatal = 1;
117                simError();                            
118                              
119                break;
120              case consSuccess:
121 <              //constrain the pair by moving two elements
121 >              // constrain the pair by moving two elements
122                done = false;
123                consPair->getConsElem1()->setMoving(true);
124                consPair->getConsElem2()->setMoving(true);
125                break;
126              case consAlready:
127 <              //current pair is already constrained, do not need to move the elements
127 >              // current pair is already constrained, do not need to
128 >              // move the elements
129                break;
130              default:          
131 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
131 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
132 >                      "Error: unrecognized status");
133                painCave.isFatal = 1;
134                simError();                          
135                break;
# Line 128 | Line 139 | namespace OpenMD {
139        }//end for(iter->first())
140  
141  
142 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
143 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
142 >      for (mol = info_->beginMolecule(mi); mol != NULL;
143 >           mol = info_->nextMolecule(mi)) {
144 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
145 >             consElem = mol->nextConstraintElem(cei)) {
146            consElem->setMoved(consElem->getMoving());
147            consElem->setMoving(false);
148          }
# Line 140 | Line 153 | namespace OpenMD {
153  
154      if (!done){
155        sprintf(painCave.errMsg,
156 <              "Constraint failure in Rattle::constrainA, too many iterations: %d\n",
156 >              "Constraint failure in Rattle::constrainA, "
157 >              "too many iterations: %d\n",
158                iteration);
159        painCave.isFatal = 1;
160        simError();    
# Line 148 | Line 162 | namespace OpenMD {
162    }
163  
164    int Rattle::constraintPairA(ConstraintPair* consPair){
165 +
166      ConstraintElem* consElem1 = consPair->getConsElem1();
167      ConstraintElem* consElem2 = consPair->getConsElem2();
168  
# Line 208 | Line 223 | namespace OpenMD {
223        Vector3d velB = consElem2->getVel();
224        velB -= rmb * delta;
225        consElem2->setVel(velB);
226 <
226 >      
227 >      // report the constraint force back to the constraint pair:
228 >      consPair->setConstraintForce(gab);
229        return consSuccess;
230      }
231      else
# Line 250 | Line 267 | namespace OpenMD {
267        velB -= rmb * delta;
268        consElem2->setVel(velB);
269  
270 +      // report the constraint force back to the constraint pair:
271 +      consPair->setConstraintForce(gab);
272        return consSuccess;
273      }
274      else

Diff Legend

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