ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/NPT.cpp
(Generate patch)

Comparing trunk/src/integrators/NPT.cpp (file contents):
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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   */
42  
43   #include <math.h>
# Line 62 | Line 63 | namespace OpenMD {
63   namespace OpenMD {
64  
65    NPT::NPT(SimInfo* info) :
66 <    VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
66 >    VelocityVerletIntegrator(info), etaTolerance(1e-6), chiTolerance(1e-6),
67 >    maxIterNum_(4) {
68  
69        Globals* simParams = info_->getSimParams();
70      
71        if (!simParams->getUseIntialExtendedSystemState()) {
72          Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
73 <        currSnapshot->setChi(0.0);
74 <        currSnapshot->setIntegralOfChiDt(0.0);
73 <        currSnapshot->setEta(Mat3x3d(0.0));
73 >        currSnapshot->setThermostat(make_pair(0.0, 0.0));
74 >        currSnapshot->setBarostat(Mat3x3d(0.0));
75        }
76      
77        if (!simParams->haveTargetTemp()) {
# Line 117 | Line 118 | namespace OpenMD {
118        tt2 = tauThermostat * tauThermostat;
119        tb2 = tauBarostat * tauBarostat;
120  
121 <      update();
121 >      updateSizes();
122      }
123  
124    NPT::~NPT() {
125    }
126  
127 <  void NPT::doUpdate() {
127 >  void NPT::doUpdateSizes() {
128  
129      oldPos.resize(info_->getNIntegrableObjects());
130      oldVel.resize(info_->getNIntegrableObjects());
# Line 135 | Line 136 | namespace OpenMD {
136      SimInfo::MoleculeIterator i;
137      Molecule::IntegrableObjectIterator  j;
138      Molecule* mol;
139 <    StuntDouble* integrableObject;
139 >    StuntDouble* sd;
140      Vector3d Tb, ji;
141      RealType mass;
142      Vector3d vel;
# Line 144 | Line 145 | namespace OpenMD {
145      Vector3d sc;
146      int index;
147  
148 <    chi= currentSnapshot_->getChi();
148 <    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
148 >    thermostat = snap->getThermostat();
149      loadEta();
150      
151      instaTemp =thermo.getTemperature();
# Line 153 | Line 153 | namespace OpenMD {
153      instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154      instaVol =thermo.getVolume();
155  
156 <    Vector3d  COM = info_->getCom();
156 >    Vector3d  COM = thermo.getCom();
157  
158      //evolve velocity half step
159  
160      calcVelScale();
161  
162 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
163 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
164 <           integrableObject = mol->nextIntegrableObject(j)) {
162 >    for (mol = info_->beginMolecule(i); mol != NULL;
163 >         mol = info_->nextMolecule(i)) {
164 >
165 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
166 >           sd = mol->nextIntegrableObject(j)) {
167                  
168 <        vel = integrableObject->getVel();
169 <        frc = integrableObject->getFrc();
168 >        vel = sd->getVel();
169 >        frc = sd->getFrc();
170  
171 <        mass = integrableObject->getMass();
171 >        mass = sd->getMass();
172  
173          getVelScaleA(sc, vel);
174  
175          // velocity half step  (use chi from previous step here):
176 <        //vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
176 >
177          vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
178 <        integrableObject->setVel(vel);
178 >        sd->setVel(vel);
179  
180 <        if (integrableObject->isDirectional()) {
180 >        if (sd->isDirectional()) {
181  
182            // get and convert the torque to body frame
183  
184 <          Tb = integrableObject->lab2Body(integrableObject->getTrq());
184 >          Tb = sd->lab2Body(sd->getTrq());
185  
186            // get the angular momentum, and propagate a half step
187  
188 <          ji = integrableObject->getJ();
188 >          ji = sd->getJ();
189  
190 <          //ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
191 <          ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
192 <                
193 <          rotAlgo->rotate(integrableObject, ji, dt);
190 >          ji += dt2*PhysicalConstants::energyConvert * Tb
191 >            - dt2*thermostat.first* ji;
192 >                
193 >          rotAlgo_->rotate(sd, ji, dt);
194  
195 <          integrableObject->setJ(ji);
195 >          sd->setJ(ji);
196          }
197              
198        }
199      }
200      // evolve chi and eta  half step
201  
202 <    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202 >    thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
203      
204      evolveEtaA();
205  
206      //calculate the integral of chidt
207 <    integralOfChidt += dt2 * chi;
207 >    thermostat.second += dt2 * thermostat.first;
208      
209 +    flucQ_->moveA();
210 +
211 +
212      index = 0;
213 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
214 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
215 <           integrableObject = mol->nextIntegrableObject(j)) {
216 <        oldPos[index++] = integrableObject->getPos();            
213 >    for (mol = info_->beginMolecule(i); mol != NULL;
214 >         mol = info_->nextMolecule(i)) {
215 >
216 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
217 >           sd = mol->nextIntegrableObject(j)) {
218 >
219 >        oldPos[index++] = sd->getPos();            
220 >
221        }
222      }
223      
# Line 216 | Line 225 | namespace OpenMD {
225  
226      for(int k = 0; k < maxIterNum_; k++) {
227        index = 0;
228 <      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
229 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221 <             integrableObject = mol->nextIntegrableObject(j)) {
228 >      for (mol = info_->beginMolecule(i); mol != NULL;
229 >           mol = info_->nextMolecule(i)) {
230  
231 <          vel = integrableObject->getVel();
232 <          pos = integrableObject->getPos();
231 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
232 >             sd = mol->nextIntegrableObject(j)) {
233  
234 +          vel = sd->getVel();
235 +          pos = sd->getPos();
236 +
237            this->getPosScale(pos, COM, index, sc);
238  
239            pos = oldPos[index] + dt * (vel + sc);
240 <          integrableObject->setPos(pos);    
240 >          sd->setPos(pos);    
241  
242            ++index;
243          }
244        }
245  
246 <      rattle->constraintA();
246 >      rattle_->constraintA();
247      }
248  
249      // Scale the box after all the positions have been moved:
250  
251      this->scaleSimBox();
252  
253 <    currentSnapshot_->setChi(chi);
243 <    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
253 >    snap->setThermostat(thermostat);
254  
255      saveEta();
256    }
# Line 249 | Line 259 | namespace OpenMD {
259      SimInfo::MoleculeIterator i;
260      Molecule::IntegrableObjectIterator  j;
261      Molecule* mol;
262 <    StuntDouble* integrableObject;
262 >    StuntDouble* sd;
263      int index;
264      Vector3d Tb;
265      Vector3d ji;
# Line 258 | Line 268 | namespace OpenMD {
268      Vector3d frc;
269      RealType mass;
270  
271 <
272 <    chi= currentSnapshot_->getChi();
263 <    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264 <    RealType oldChi  = chi;
271 >    thermostat = snap->getThermostat();
272 >    RealType oldChi  = thermostat.first;
273      RealType prevChi;
274  
275      loadEta();
276      
277      //save velocity and angular momentum
278      index = 0;
279 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
280 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
281 <           integrableObject = mol->nextIntegrableObject(j)) {
279 >    for (mol = info_->beginMolecule(i); mol != NULL;
280 >         mol = info_->nextMolecule(i)) {
281 >
282 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
283 >           sd = mol->nextIntegrableObject(j)) {
284                  
285 <        oldVel[index] = integrableObject->getVel();
286 <        oldJi[index] = integrableObject->getJ();
285 >        oldVel[index] = sd->getVel();
286 >
287 >        if (sd->isDirectional())
288 >           oldJi[index] = sd->getJ();
289 >
290          ++index;
291        }
292      }
# Line 286 | Line 299 | namespace OpenMD {
299        instaPress =thermo.getPressure();
300  
301        // evolve chi another half step using the temperature at t + dt/2
302 <      prevChi = chi;
303 <      chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
302 >      prevChi = thermostat.first;
303 >      thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
304  
305        //evolve eta
306        this->evolveEtaB();
307        this->calcVelScale();
308  
309        index = 0;
310 <      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
311 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 <             integrableObject = mol->nextIntegrableObject(j)) {            
310 >      for (mol = info_->beginMolecule(i); mol != NULL;
311 >           mol = info_->nextMolecule(i)) {
312  
313 <          frc = integrableObject->getFrc();
314 <          vel = integrableObject->getVel();
313 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
314 >             sd = mol->nextIntegrableObject(j)) {            
315  
316 <          mass = integrableObject->getMass();
316 >          frc = sd->getFrc();
317 >          mass = sd->getMass();
318  
319            getVelScaleB(sc, index);
320  
321            // velocity half step
322 <          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
323 <          vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
324 <          integrableObject->setVel(vel);
322 >          vel = oldVel[index]
323 >            + dt2*PhysicalConstants::energyConvert/mass* frc
324 >            - dt2*sc;
325  
326 <          if (integrableObject->isDirectional()) {
326 >          sd->setVel(vel);
327 >
328 >          if (sd->isDirectional()) {
329              // get and convert the torque to body frame
330 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
330 >            Tb = sd->lab2Body(sd->getTrq());
331  
332 <            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
333 <            ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
334 <            integrableObject->setJ(ji);
332 >            ji = oldJi[index]
333 >              + dt2*PhysicalConstants::energyConvert*Tb
334 >              - dt2*thermostat.first*oldJi[index];
335 >
336 >            sd->setJ(ji);
337            }
338  
339            ++index;
340          }
341        }
342          
343 <      rattle->constraintB();
343 >      rattle_->constraintB();
344  
345 <      if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
345 >      if ((fabs(prevChi - thermostat.first) <= chiTolerance) &&
346 >          this->etaConverged())
347          break;
348      }
349  
350      //calculate integral of chidt
351 <    integralOfChidt += dt2 * chi;
351 >    thermostat.second += dt2 * thermostat.first;
352  
353 <    currentSnapshot_->setChi(chi);
336 <    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);    
353 >    snap->setThermostat(thermostat);
354  
355 +    flucQ_->moveB();
356      saveEta();
357    }
358  
359    void NPT::resetIntegrator(){
360 <      currentSnapshot_->setChi(0.0);
361 <      currentSnapshot_->setIntegralOfChiDt(0.0);
344 <      resetEta();
360 >    snap->setThermostat(make_pair(0.0, 0.0));
361 >    resetEta();
362    }
363  
364 <
365 <    void NPT::resetEta() {
366 <      Mat3x3d etaMat(0.0);
367 <      currentSnapshot_->setEta(etaMat);    
351 <    }
352 <    
364 >  void NPT::resetEta() {
365 >    Mat3x3d etaMat(0.0);
366 >    snap->setBarostat(etaMat);    
367 >  }
368   }

Diff Legend

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