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 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 <  Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doRattle_(false) {
48 >  Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10),
49 >                                  consTolerance_(1.0e-6), doRattle_(false),
50 >                                  currConstraintTime_(0.0) {
51      
52 <    if (info_->getNConstraints() > 0)
52 >    if (info_->getNGlobalConstraints() > 0)
53        doRattle_ = true;
54      
55 +    Globals* simParams = info_->getSimParams();
56  
57 <    if (info_->getSimParams()->haveDt()) {
58 <      dt_ = info_->getSimParams()->getDt();
57 >    if (simParams->haveDt()) {
58 >      dt_ = simParams->getDt();
59      } else {
60        sprintf(painCave.errMsg,
61 <              "Integrator Error: dt is not set\n");
61 >              "Rattle Error: dt is not set\n");
62        painCave.isFatal = 1;
63        simError();
64      }    
65 <    
65 >
66      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
67 +    if (simParams->haveConstraintTime()){
68 +      constraintTime_ = simParams->getConstraintTime();
69 +    } else {
70 +      constraintTime_ = simParams->getStatusTime();
71 +    }
72 +
73 +    constraintOutputFile_ = getPrefix(info_->getFinalConfigFileName()) +
74 +      ".constraintForces";
75 +
76 +    // create ConstraintWriter  
77 +    constraintWriter_ = new ConstraintWriter(info_,
78 +                                             constraintOutputFile_.c_str());  
79 +
80 +    if (!constraintWriter_){
81 +      sprintf(painCave.errMsg, "Failed to create ConstraintWriter\n");
82 +      painCave.isFatal = 1;
83 +      simError();
84 +    }
85    }
86  
87    void Rattle::constraintA() {
# Line 70 | Line 91 | namespace OpenMD {
91    void Rattle::constraintB() {
92      if (!doRattle_) return;    
93      doConstraint(&Rattle::constraintPairB);
94 +
95 +    if (currentSnapshot_->getTime() >= currConstraintTime_){
96 +      Molecule* mol;
97 +      SimInfo::MoleculeIterator mi;
98 +      ConstraintPair* consPair;
99 +      Molecule::ConstraintPairIterator cpi;
100 +      std::list<ConstraintPair*> constraints;
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 +          constraints.push_back(consPair);
107 +        }
108 +      }
109 +      constraintWriter_->writeConstraintForces(constraints);
110 +      currConstraintTime_ += constraintTime_;
111 +    }
112    }
113  
114    void Rattle::doConstraint(ConstraintPairFuncPtr func) {
# Line 82 | Line 121 | namespace OpenMD {
121      ConstraintPair* consPair;
122      Molecule::ConstraintPairIterator cpi;
123      
124 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
125 <      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
124 >    for (mol = info_->beginMolecule(mi); mol != NULL;
125 >         mol = info_->nextMolecule(mi)) {
126 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
127 >           consElem = mol->nextConstraintElem(cei)) {
128          consElem->setMoved(true);
129          consElem->setMoving(false);
130        }
# Line 97 | Line 138 | namespace OpenMD {
138  
139        //loop over every constraint pair
140  
141 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
142 <        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
141 >      for (mol = info_->beginMolecule(mi); mol != NULL;
142 >           mol = info_->nextMolecule(mi)) {
143 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
144 >             consPair = mol->nextConstraintPair(cpi)) {
145  
146  
147            //dispatch constraint algorithm
# Line 108 | Line 151 | namespace OpenMD {
151              switch(exeStatus){
152              case consFail:
153                sprintf(painCave.errMsg,
154 <                      "Constraint failure in Rattle::constrainA, Constraint Fail\n");
154 >                      "Constraint failure in Rattle::constrainA, "
155 >                      "Constraint Fail\n");
156                painCave.isFatal = 1;
157                simError();                            
158                              
159                break;
160              case consSuccess:
161 <              //constrain the pair by moving two elements
161 >              // constrain the pair by moving two elements
162                done = false;
163                consPair->getConsElem1()->setMoving(true);
164                consPair->getConsElem2()->setMoving(true);
165                break;
166              case consAlready:
167 <              //current pair is already constrained, do not need to move the elements
167 >              // current pair is already constrained, do not need to
168 >              // move the elements
169                break;
170              default:          
171 <              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
171 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
172 >                      "Error: unrecognized status");
173                painCave.isFatal = 1;
174                simError();                          
175                break;
# Line 133 | Line 179 | namespace OpenMD {
179        }//end for(iter->first())
180  
181  
182 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
183 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
182 >      for (mol = info_->beginMolecule(mi); mol != NULL;
183 >           mol = info_->nextMolecule(mi)) {
184 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
185 >             consElem = mol->nextConstraintElem(cei)) {
186            consElem->setMoved(consElem->getMoving());
187            consElem->setMoving(false);
188          }
# Line 145 | Line 193 | namespace OpenMD {
193  
194      if (!done){
195        sprintf(painCave.errMsg,
196 <              "Constraint failure in Rattle::constrainA, too many iterations: %d\n",
196 >              "Constraint failure in Rattle::constrainA, "
197 >              "too many iterations: %d\n",
198                iteration);
199        painCave.isFatal = 1;
200        simError();    
# Line 214 | Line 263 | namespace OpenMD {
263        Vector3d velB = consElem2->getVel();
264        velB -= rmb * delta;
265        consElem2->setVel(velB);
266 <
266 >      
267 >      // report the constraint force back to the constraint pair:
268 >      consPair->setConstraintForce(gab);
269        return consSuccess;
270      }
271      else
# Line 256 | Line 307 | namespace OpenMD {
307        velB -= rmb * delta;
308        consElem2->setVel(velB);
309  
310 +      // report the constraint force back to the constraint pair:
311 +      consPair->setConstraintForce(gab);
312        return consSuccess;
313      }
314      else

Diff Legend

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