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 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).          
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_->getSimParams()->haveDt()) {
51 <      dt_ = info_->getSimParams()->getDt();
50 >    if (info_->getNGlobalConstraints() > 0)
51 >      doRattle_ = true;
52 >    
53 >    Globals* simParams = info_->getSimParams();
54 >
55 >    if (simParams->haveDt()) {
56 >      dt_ = simParams->getDt();
57      } else {
58        sprintf(painCave.errMsg,
59 <              "Integrator Error: dt is not set\n");
59 >              "Rattle Error: dt is not set\n");
60        painCave.isFatal = 1;
61        simError();
62      }    
63 <    
63 >
64      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
65 +    if (simParams->haveConstraintTime()){
66 +      constraintTime_ = simParams->getConstraintTime();
67 +    } else {
68 +      constraintTime_ = simParams->getStatusTime();
69 +    }
70 +
71 +    constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) + ".constraintForces";
72 +
73 +    // create ConstraintWriter  
74 +    constraintWriter_ = new ConstraintWriter(info_, constraintOutputFile_.c_str());  
75 +
76 +    if (!constraintWriter_){
77 +      sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n");
78 +      painCave.isFatal = 1;
79 +      simError();
80 +    }
81    }
82  
83    void Rattle::constraintA() {
84 <    if (info_->getNConstraints() > 0) {
85 <      doConstraint(&Rattle::constraintPairA);
64 <    }
84 >    if (!doRattle_) return;
85 >    doConstraint(&Rattle::constraintPairA);
86    }
87    void Rattle::constraintB() {
88 <    if (info_->getNConstraints() > 0) {
89 <      doConstraint(&Rattle::constraintPairB);
88 >    if (!doRattle_) return;    
89 >    doConstraint(&Rattle::constraintPairB);
90 >
91 >    if (currentSnapshot_->getTime() >= currConstraintTime_){
92 >      Molecule* mol;
93 >      SimInfo::MoleculeIterator mi;
94 >      ConstraintPair* consPair;
95 >      Molecule::ConstraintPairIterator cpi;
96 >      std::list<ConstraintPair*> constraints;
97 >      for (mol = info_->beginMolecule(mi); mol != NULL;
98 >           mol = info_->nextMolecule(mi)) {
99 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
100 >             consPair = mol->nextConstraintPair(cpi)) {
101 >          
102 >          constraints.push_back(consPair);
103 >        }
104 >      }
105 >      
106 >      constraintWriter_->writeConstraintForces(constraints);
107 >      currConstraintTime_ += constraintTime_;
108      }
109    }
110  
111    void Rattle::doConstraint(ConstraintPairFuncPtr func) {
112 +    if (!doRattle_) return;
113 +
114      Molecule* mol;
115      SimInfo::MoleculeIterator mi;
116      ConstraintElem* consElem;
# Line 77 | Line 118 | namespace OpenMD {
118      ConstraintPair* consPair;
119      Molecule::ConstraintPairIterator cpi;
120      
121 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
122 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
121 >    for (mol = info_->beginMolecule(mi); mol != NULL;
122 >         mol = info_->nextMolecule(mi)) {
123 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
124 >           consElem = mol->nextConstraintElem(cei)) {
125          consElem->setMoved(true);
126          consElem->setMoving(false);
127        }
# Line 92 | Line 135 | namespace OpenMD {
135  
136        //loop over every constraint pair
137  
138 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
139 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
138 >      for (mol = info_->beginMolecule(mi); mol != NULL;
139 >           mol = info_->nextMolecule(mi)) {
140 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
141 >             consPair = mol->nextConstraintPair(cpi)) {
142  
143  
144            //dispatch constraint algorithm
# Line 103 | Line 148 | namespace OpenMD {
148              switch(exeStatus){
149              case consFail:
150                sprintf(painCave.errMsg,
151 <                      "Constraint failure in Rattle::constrainA, Constraint Fail\n");
151 >                      "Constraint failure in Rattle::constrainA, "
152 >                      "Constraint Fail\n");
153                painCave.isFatal = 1;
154                simError();                            
155                              
156                break;
157              case consSuccess:
158 <              //constrain the pair by moving two elements
158 >              // constrain the pair by moving two elements
159                done = false;
160                consPair->getConsElem1()->setMoving(true);
161                consPair->getConsElem2()->setMoving(true);
162                break;
163              case consAlready:
164 <              //current pair is already constrained, do not need to move the elements
164 >              // current pair is already constrained, do not need to
165 >              // move the elements
166                break;
167              default:          
168 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
168 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
169 >                      "Error: unrecognized status");
170                painCave.isFatal = 1;
171                simError();                          
172                break;
# Line 128 | Line 176 | namespace OpenMD {
176        }//end for(iter->first())
177  
178  
179 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
180 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
179 >      for (mol = info_->beginMolecule(mi); mol != NULL;
180 >           mol = info_->nextMolecule(mi)) {
181 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
182 >             consElem = mol->nextConstraintElem(cei)) {
183            consElem->setMoved(consElem->getMoving());
184            consElem->setMoving(false);
185          }
# Line 140 | Line 190 | namespace OpenMD {
190  
191      if (!done){
192        sprintf(painCave.errMsg,
193 <              "Constraint failure in Rattle::constrainA, too many iterations: %d\n",
193 >              "Constraint failure in Rattle::constrainA, "
194 >              "too many iterations: %d\n",
195                iteration);
196        painCave.isFatal = 1;
197        simError();    
# Line 148 | Line 199 | namespace OpenMD {
199    }
200  
201    int Rattle::constraintPairA(ConstraintPair* consPair){
202 +
203      ConstraintElem* consElem1 = consPair->getConsElem1();
204      ConstraintElem* consElem2 = consPair->getConsElem2();
205  
# Line 208 | Line 260 | namespace OpenMD {
260        Vector3d velB = consElem2->getVel();
261        velB -= rmb * delta;
262        consElem2->setVel(velB);
263 <
263 >      
264 >      // report the constraint force back to the constraint pair:
265 >      consPair->setConstraintForce(gab);
266        return consSuccess;
267      }
268      else
# Line 250 | Line 304 | namespace OpenMD {
304        velB -= rmb * delta;
305        consElem2->setVel(velB);
306  
307 +      // report the constraint force back to the constraint pair:
308 +      consPair->setConstraintForce(gab);
309        return consSuccess;
310      }
311      else

Diff Legend

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