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 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    Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6), doRattle_(false) {
49      
50 <    if (info_->getNConstraints() > 0)
50 >    if (info_->getNGlobalConstraints() > 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");
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() {
# Line 70 | Line 87 | namespace OpenMD {
87    void Rattle::constraintB() {
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) {
# Line 82 | 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 97 | 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 108 | 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 133 | 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 145 | 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 214 | 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 256 | 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