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

Diff Legend

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