ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/constraints/Rattle.cpp
(Generate patch)

Comparing trunk/src/constraints/Rattle.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 44 | Line 44 | namespace oopse {
44   #include "utils/simError.h"
45   namespace oopse {
46  
47 < Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6) {
47 >  Rattle::Rattle(SimInfo* info) : info_(info), maxConsIteration_(10), consTolerance_(1.0e-6) {
48      
49      if (info_->getSimParams()->haveDt()) {
50 <        dt_ = info_->getSimParams()->getDt();
50 >      dt_ = info_->getSimParams()->getDt();
51      } else {
52 <            sprintf(painCave.errMsg,
53 <                    "Integrator Error: dt is not set\n");
54 <            painCave.isFatal = 1;
55 <            simError();
52 >      sprintf(painCave.errMsg,
53 >              "Integrator Error: dt is not set\n");
54 >      painCave.isFatal = 1;
55 >      simError();
56      }    
57      
58      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
59 < }
59 >  }
60  
61 < void Rattle::constraintA() {
61 >  void Rattle::constraintA() {
62      if (info_->getNConstraints() > 0) {
63 <        doConstraint(&Rattle::constraintPairA);
63 >      doConstraint(&Rattle::constraintPairA);
64      }
65 < }
66 < void Rattle::constraintB() {
65 >  }
66 >  void Rattle::constraintB() {
67      if (info_->getNConstraints() > 0) {
68 <        doConstraint(&Rattle::constraintPairB);
68 >      doConstraint(&Rattle::constraintPairB);
69      }
70 < }
70 >  }
71  
72 < void Rattle::doConstraint(ConstraintPairFuncPtr func) {
72 >  void Rattle::doConstraint(ConstraintPairFuncPtr func) {
73      Molecule* mol;
74      SimInfo::MoleculeIterator mi;
75      ConstraintElem* consElem;
# Line 78 | Line 78 | void Rattle::doConstraint(ConstraintPairFuncPtr func)
78      Molecule::ConstraintPairIterator cpi;
79      
80      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
81 <        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
82 <            consElem->setMoved(true);
83 <            consElem->setMoving(false);
84 <        }
81 >      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
82 >        consElem->setMoved(true);
83 >        consElem->setMoving(false);
84 >      }
85      }
86      
87      //main loop of constraint algorithm
88      bool done = false;
89      int iteration = 0;
90      while(!done && iteration < maxConsIteration_){
91 <        done = true;
91 >      done = true;
92  
93 <        //loop over every constraint pair
93 >      //loop over every constraint pair
94  
95 <        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
96 <            for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
95 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
96 >        for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
97  
98  
99 <                //dispatch constraint algorithm
100 <                if(consPair->isMoved()) {
101 <                    int exeStatus = (this->*func)(consPair);
99 >          //dispatch constraint algorithm
100 >          if(consPair->isMoved()) {
101 >            int exeStatus = (this->*func)(consPair);
102  
103 <                    switch(exeStatus){
104 <                        case consFail:
105 <                            sprintf(painCave.errMsg,
106 <                            "Constraint failure in Rattle::constrainA, Constraint Fail\n");
107 <                            painCave.isFatal = 1;
108 <                            simError();                            
103 >            switch(exeStatus){
104 >            case consFail:
105 >              sprintf(painCave.errMsg,
106 >                      "Constraint failure in Rattle::constrainA, Constraint Fail\n");
107 >              painCave.isFatal = 1;
108 >              simError();                            
109                              
110 <                            break;
111 <                        case consSuccess:
112 <                            //constrain the pair by moving two elements
113 <                            done = false;
114 <                            consPair->getConsElem1()->setMoving(true);
115 <                            consPair->getConsElem2()->setMoving(true);
116 <                            break;
117 <                        case consAlready:
118 <                            //current pair is already constrained, do not need to move the elements
119 <                            break;
120 <                        default:          
121 <                            sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
122 <                            painCave.isFatal = 1;
123 <                            simError();                          
124 <                            break;
125 <                    }      
126 <                }
127 <            }
128 <        }//end for(iter->first())
110 >              break;
111 >            case consSuccess:
112 >              //constrain the pair by moving two elements
113 >              done = false;
114 >              consPair->getConsElem1()->setMoving(true);
115 >              consPair->getConsElem2()->setMoving(true);
116 >              break;
117 >            case consAlready:
118 >              //current pair is already constrained, do not need to move the elements
119 >              break;
120 >            default:          
121 >              sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstrain() Error: unrecognized status");
122 >              painCave.isFatal = 1;
123 >              simError();                          
124 >              break;
125 >            }      
126 >          }
127 >        }
128 >      }//end for(iter->first())
129  
130  
131 <        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
132 <            for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
133 <                consElem->setMoved(consElem->getMoving());
134 <                consElem->setMoving(false);
135 <            }
136 <        }
131 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
132 >        for (consElem = mol->beginConstraintElem(cei); consElem != NULL; consElem = mol->nextConstraintElem(cei)) {
133 >          consElem->setMoved(consElem->getMoving());
134 >          consElem->setMoving(false);
135 >        }
136 >      }
137  
138 <        iteration++;
138 >      iteration++;
139      }//end while
140  
141      if (!done){
# Line 145 | Line 145 | void Rattle::doConstraint(ConstraintPairFuncPtr func)
145        painCave.isFatal = 1;
146        simError();    
147      }
148 < }
148 >  }
149  
150 < int Rattle::constraintPairA(ConstraintPair* consPair){
151 <  ConstraintElem* consElem1 = consPair->getConsElem1();
152 <  ConstraintElem* consElem2 = consPair->getConsElem2();
150 >  int Rattle::constraintPairA(ConstraintPair* consPair){
151 >    ConstraintElem* consElem1 = consPair->getConsElem1();
152 >    ConstraintElem* consElem2 = consPair->getConsElem2();
153  
154 <  Vector3d posA = consElem1->getPos();
155 <  Vector3d posB = consElem2->getPos();
154 >    Vector3d posA = consElem1->getPos();
155 >    Vector3d posB = consElem2->getPos();
156  
157 <  Vector3d pab = posA -posB;  
158 <
159 <  //periodic boundary condition
157 >    Vector3d pab = posA -posB;  
158  
159 <  currentSnapshot_->wrapVector(pab);
159 >    //periodic boundary condition
160  
161 <  double pabsq = pab.lengthSquare();
161 >    currentSnapshot_->wrapVector(pab);
162  
163 <  double rabsq = consPair->getConsDistSquare();
166 <  double diffsq = rabsq - pabsq;
163 >    RealType pabsq = pab.lengthSquare();
164  
165 <  // the original rattle code from alan tidesley
166 <  if (fabs(diffsq) > (consTolerance_ * rabsq * 2)){
165 >    RealType rabsq = consPair->getConsDistSquare();
166 >    RealType diffsq = rabsq - pabsq;
167 >
168 >    // the original rattle code from alan tidesley
169 >    if (fabs(diffsq) > (consTolerance_ * rabsq * 2)){
170      
171 <    Vector3d oldPosA = consElem1->getPrevPos();
172 <    Vector3d oldPosB = consElem2->getPrevPos();      
171 >      Vector3d oldPosA = consElem1->getPrevPos();
172 >      Vector3d oldPosB = consElem2->getPrevPos();      
173  
174 <    Vector3d rab = oldPosA - oldPosB;    
174 >      Vector3d rab = oldPosA - oldPosB;    
175  
176 <    currentSnapshot_->wrapVector(rab);
176 >      currentSnapshot_->wrapVector(rab);
177  
178 <    double rpab = dot(rab, pab);
179 <    double rpabsq = rpab * rpab;
178 >      RealType rpab = dot(rab, pab);
179 >      RealType rpabsq = rpab * rpab;
180  
181 <    if (rpabsq < (rabsq * -diffsq)){
182 <      return consFail;
183 <    }
181 >      if (rpabsq < (rabsq * -diffsq)){
182 >        return consFail;
183 >      }
184  
185 <    double rma = 1.0 / consElem1->getMass();
186 <    double rmb = 1.0 / consElem2->getMass();
185 >      RealType rma = 1.0 / consElem1->getMass();
186 >      RealType rmb = 1.0 / consElem2->getMass();
187  
188 <    double gab = diffsq / (2.0 * (rma + rmb) * rpab);
188 >      RealType gab = diffsq / (2.0 * (rma + rmb) * rpab);
189  
190 <    Vector3d delta = rab * gab;
190 >      Vector3d delta = rab * gab;
191  
192 <    //set atom1's position
193 <    posA += rma * delta;    
194 <    consElem1->setPos(posA);
192 >      //set atom1's position
193 >      posA += rma * delta;    
194 >      consElem1->setPos(posA);
195  
196 <    //set atom2's position
197 <    posB -= rmb * delta;
198 <    consElem2->setPos(posB);
196 >      //set atom2's position
197 >      posB -= rmb * delta;
198 >      consElem2->setPos(posB);
199  
200 <    delta /= dt_;
200 >      delta /= dt_;
201      
202 <    //set atom1's velocity
203 <    Vector3d velA = consElem1->getVel();
204 <    velA += rma * delta;
205 <    consElem1->setVel(velA);
202 >      //set atom1's velocity
203 >      Vector3d velA = consElem1->getVel();
204 >      velA += rma * delta;
205 >      consElem1->setVel(velA);
206  
207 <    //set atom2's velocity
208 <    Vector3d velB = consElem2->getVel();
209 <    velB -= rmb * delta;
210 <    consElem2->setVel(velB);
207 >      //set atom2's velocity
208 >      Vector3d velB = consElem2->getVel();
209 >      velB -= rmb * delta;
210 >      consElem2->setVel(velB);
211  
212 <    return consSuccess;
213 <  }
214 <  else
215 <    return consAlready;
212 >      return consSuccess;
213 >    }
214 >    else
215 >      return consAlready;
216    
217 < }
217 >  }
218  
219  
220 < int Rattle::constraintPairB(ConstraintPair* consPair){
220 >  int Rattle::constraintPairB(ConstraintPair* consPair){
221      ConstraintElem* consElem1 = consPair->getConsElem1();
222      ConstraintElem* consElem2 = consPair->getConsElem2();
223  
# Line 234 | Line 234 | int Rattle::constraintPairB(ConstraintPair* consPair){
234  
235      currentSnapshot_->wrapVector(rab);
236  
237 <    double rma = 1.0 / consElem1->getMass();
238 <    double rmb = 1.0 / consElem2->getMass();
237 >    RealType rma = 1.0 / consElem1->getMass();
238 >    RealType rmb = 1.0 / consElem2->getMass();
239  
240 <    double rvab = dot(rab, dv);
240 >    RealType rvab = dot(rab, dv);
241  
242 <    double gab = -rvab / ((rma + rmb) * consPair->getConsDistSquare());
242 >    RealType gab = -rvab / ((rma + rmb) * consPair->getConsDistSquare());
243  
244      if (fabs(gab) > consTolerance_){
245        Vector3d delta = rab * gab;
# Line 255 | Line 255 | int Rattle::constraintPairB(ConstraintPair* consPair){
255      else
256        return consAlready;
257  
258 < }
258 >  }
259  
260   }

Diff Legend

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