ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/NPT.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/integrators/NPT.cpp (file contents):
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC vs.
Revision 1822 by tim, Thu Dec 2 02:08:29 2004 UTC

# Line 3 | Line 3
3   #include "brains/SimInfo.hpp"
4   #include "brains/Thermo.hpp"
5   #include "integrators/NPT.hpp"
6 + #include "math/SquareMatrix3.hpp"
7 + #include "primitives/Molecule.hpp"
8 + #include "utils/OOPSEConstant.hpp"
9   #include "utils/simError.h"
10  
8 –
11   // Basic isotropic thermostating and barostating via the Melchionna
12   // modification of the Hoover algorithm:
13   //
# Line 27 | Line 29 | NPT::NPT(SimInfo* info) :
29          Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
30          currSnapshot->setChi(0.0);
31          currSnapshot->setIntegralOfChiDt(0.0);
32 +        currSnapshot->setEta(Mat3x3d(0.0));
33      }
34      
35      if (!globals->haveTargetTemp()) {
# Line 79 | Line 82 | void NPT::update() {
82   NPT::~NPT() {
83   }
84  
85 < void NPT::update() {
85 > void NPT::doUpdate() {
86 >    VelocityVerletIntegrator::update();
87      oldPos.resize(info_->getNIntegrableObjects());
88      oldVel.resize(info_->getNIntegrableObjects());
89      oldJi.resize(info_->getNIntegrableObjects());
90      // We need NkBT a lot, so just set it here: This is the RAW number
91      // of integrableObjects, so no subtraction or addition of constraints or
92      // orientational degrees of freedom:
93 <    NkBT = info_->getNGlobalIntegrableObjects()*kB *targetTemp;
93 >    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
94  
95      // fkBT is used because the thermostat operates on more degrees of freedom
96      // than the barostat (when there are particles with orientational degrees
97      // of freedom).  
98 <    fkBT = info_->getNdf()*kB *targetTemp;    
98 >    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
99   }
100  
101   void NPT::moveA() {
102 <    typename SimInfo::MoleculeIterator i;
103 <    typename Molecule::IntegrableObjectIterator  j;
102 >    SimInfo::MoleculeIterator i;
103 >    Molecule::IntegrableObjectIterator  j;
104      Molecule* mol;
105      StuntDouble* integrableObject;
106      Vector3d Tb, ji;
# Line 105 | Line 109 | void NPT::moveA() {
109      Vector3d pos;
110      Vector3d frc;
111      Vector3d sc;
108 –    Vector3d COM;
112      int index;
113 +
114 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115 +    chi= currSnapshot->getChi();
116 +    integralOfChidt = currSnapshot->getIntegralOfChiDt();
117 +    loadEta();
118      
119 <    instaTemp = tStats->getTemperature();
120 <    tStats->getPressureTensor(press);
121 <    instaPress = p_convert * (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
122 <    instaVol = tStats->getVolume();
119 >    instaTemp =thermo.getTemperature();
120 >    press = thermo.getPressureTensor();
121 >    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
122 >    instaVol =thermo.getVolume();
123  
124 <    tStats->getCOM(COM);
124 >    Vector3d  COM = info_->getCom();
125  
126      //evolve velocity half step
127  
# Line 131 | Line 139 | void NPT::moveA() {
139              getVelScaleA(sc, vel);
140  
141              // velocity half step  (use chi from previous step here):
142 <            //vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
143 <            vel += dt2*eConvert/mass* frc - dt2*sc;
142 >            //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
143 >            vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
144              integrableObject->setVel(vel);
145  
146              if (integrableObject->isDirectional()) {
# Line 146 | Line 154 | void NPT::moveA() {
154  
155                  ji = integrableObject->getJ();
156  
157 <                //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
158 <                ji += dt2*eConvert * Tb - dt2*chi* ji;
157 >                //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
158 >                ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
159                  
160 <                this->rotationPropagation(integrableObject, ji);
160 >                rotAlgo->rotate(integrableObject, ji, dt);
161  
162                  integrableObject->setJ(ji);
163              }
# Line 158 | Line 166 | void NPT::moveA() {
166      }
167      // evolve chi and eta  half step
168  
169 <    evolveChiA();
169 >    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
170 >    
171      evolveEtaA();
172  
173      //calculate the integral of chidt
# Line 198 | Line 207 | void NPT::moveA() {
207      // Scale the box after all the positions have been moved:
208  
209      this->scaleSimBox();
210 +
211 +    currSnapshot->setChi(chi);
212 +    currSnapshot->setIntegralOfChiDt(integralOfChidt);
213 +
214 +    saveEta();
215   }
216  
217   void NPT::moveB(void) {
218 <    typename SimInfo::MoleculeIterator i;
219 <    typename Molecule::IntegrableObjectIterator  j;
218 >    SimInfo::MoleculeIterator i;
219 >    Molecule::IntegrableObjectIterator  j;
220      Molecule* mol;
221      StuntDouble* integrableObject;
222      int index;
# Line 212 | Line 226 | void NPT::moveB(void) {
226      Vector3d vel;
227      Vector3d frc;
228      double mass;
229 +
230 +
231 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
232 +    chi= currSnapshot->getChi();
233 +    integralOfChidt = currSnapshot->getIntegralOfChiDt();
234 +    double oldChi  = chi;
235 +    double prevChi;
236  
237 +    loadEta();
238 +    
239      //save velocity and angular momentum
240      index = 0;
241      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 220 | Line 243 | void NPT::moveB(void) {
243                 integrableObject = mol->nextIntegrableObject(j)) {
244                  
245              oldVel[index] = integrableObject->getVel();
246 <            oldJi[3 * i + j] = integrableObject->getJ();
246 >            oldJi[index] = integrableObject->getJ();
247              ++index;
248          }
249      }
250  
251      // do the iteration:
252 <    instaVol = tStats->getVolume();
252 >    instaVol =thermo.getVolume();
253  
254      for(int k = 0; k < maxIterNum_; k++) {
255 <        instaTemp = tStats->getTemperature();
256 <        instaPress = tStats->getPressure();
255 >        instaTemp =thermo.getTemperature();
256 >        instaPress =thermo.getPressure();
257  
258          // evolve chi another half step using the temperature at t + dt/2
259 <        this->evolveChiB();
259 >        prevChi = chi;
260 >        chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
261 >
262 >        //evolve eta
263          this->evolveEtaB();
264          this->calcVelScale();
265  
# Line 241 | Line 267 | void NPT::moveB(void) {
267          for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
268              for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
269                     integrableObject = mol->nextIntegrableObject(j)) {            
270 <                integrableObject->getFrc(frc);
270 >
271 >                frc = integrableObject->getFrc();
272                  vel = integrableObject->getVel();
273  
274                  mass = integrableObject->getMass();
275  
276 <                getVelScaleB(sc, i);
276 >                getVelScaleB(sc, index);
277  
278                  // velocity half step
279 <                //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * eConvert - sc[j]);
280 <                vel = oldVel[index] + dt2*eConvert/mass* frc - dt2*sc;
279 >                //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
280 >                vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
281                  integrableObject->setVel(vel);
282  
283                  if (integrableObject->isDirectional()) {
# Line 258 | Line 285 | void NPT::moveB(void) {
285                      Tb = integrableObject->getTrq();
286                      integrableObject->lab2Body(Tb);
287  
288 <                    //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
289 <                    ji = oldJi[index] + dt2*eConvert*Tb - dt2*chi*oldJi[index];
288 >                    //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
289 >                    ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
290                      integrableObject->setJ(ji);
291                  }
292  
# Line 269 | Line 296 | void NPT::moveB(void) {
296          
297          //constraintAlgorithm->doConstrainB();
298  
299 <        if (this->chiConverged() && this->etaConverged())
299 >        if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
300              break;
301      }
302  
303      //calculate integral of chidt
304      integralOfChidt += dt2 * chi;
278 – }
305  
306 < void NPT::evolveChiA() {
307 <    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
282 <    oldChi = chi;
283 < }
306 >    currSnapshot->setChi(chi);
307 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);    
308  
309 < void NPT::evolveChiB() {
286 <    prevChi = chi;
287 <    chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
309 >    saveEta();
310   }
311  
290 – bool NPT::chiConverged() {
291 –    return (fabs(prevChi - chi) <= chiTolerance);
312   }
293 –
294 – }

Diff Legend

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