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

Diff Legend

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