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 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 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 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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).                        
40   */
41  
42   #include <math.h>
# Line 46 | Line 46
46   #include "integrators/NPT.hpp"
47   #include "math/SquareMatrix3.hpp"
48   #include "primitives/Molecule.hpp"
49 < #include "utils/OOPSEConstant.hpp"
49 > #include "utils/PhysicalConstants.hpp"
50   #include "utils/simError.h"
51  
52   // Basic isotropic thermostating and barostating via the Melchionna
# Line 59 | Line 59
59   //
60   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61  
62 < namespace oopse {
62 > namespace OpenMD {
63  
64 < NPT::NPT(SimInfo* info) :
64 >  NPT::NPT(SimInfo* info) :
65      VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
66  
67 <    Globals* simParams = info_->getSimParams();
67 >      Globals* simParams = info_->getSimParams();
68      
69 <    if (!simParams->getUseInitXSstate()) {
69 >      if (!simParams->getUseIntialExtendedSystemState()) {
70          Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
71          currSnapshot->setChi(0.0);
72          currSnapshot->setIntegralOfChiDt(0.0);
73          currSnapshot->setEta(Mat3x3d(0.0));
74 <    }
74 >      }
75      
76 <    if (!simParams->haveTargetTemp()) {
76 >      if (!simParams->haveTargetTemp()) {
77          sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
78          painCave.isFatal = 1;
79 <        painCave.severity = OOPSE_ERROR;
79 >        painCave.severity = OPENMD_ERROR;
80          simError();
81 <    } else {
81 >      } else {
82          targetTemp = simParams->getTargetTemp();
83 <    }
83 >      }
84  
85 <    // We must set tauThermostat
86 <    if (!simParams->haveTauThermostat()) {
85 >      // We must set tauThermostat
86 >      if (!simParams->haveTauThermostat()) {
87          sprintf(painCave.errMsg, "If you use the constant temperature\n"
88 <                                     "\tintegrator, you must set tauThermostat_.\n");
88 >                "\tintegrator, you must set tauThermostat.\n");
89  
90 <        painCave.severity = OOPSE_ERROR;
90 >        painCave.severity = OPENMD_ERROR;
91          painCave.isFatal = 1;
92          simError();
93 <    } else {
93 >      } else {
94          tauThermostat = simParams->getTauThermostat();
95 <    }
95 >      }
96  
97 <    if (!simParams->haveTargetPressure()) {
97 >      if (!simParams->haveTargetPressure()) {
98          sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
99 <                                     "   without a targetPressure!\n");
99 >                "   without a targetPressure!\n");
100  
101          painCave.isFatal = 1;
102          simError();
103 <    } else {
103 >      } else {
104          targetPressure = simParams->getTargetPressure();
105 <    }
105 >      }
106      
107 <    if (!simParams->haveTauBarostat()) {
107 >      if (!simParams->haveTauBarostat()) {
108          sprintf(painCave.errMsg,
109                  "If you use the NPT integrator, you must set tauBarostat.\n");
110 <        painCave.severity = OOPSE_ERROR;
110 >        painCave.severity = OPENMD_ERROR;
111          painCave.isFatal = 1;
112          simError();
113 <    } else {
113 >      } else {
114          tauBarostat = simParams->getTauBarostat();
115 <    }
115 >      }
116      
117 <    tt2 = tauThermostat * tauThermostat;
118 <    tb2 = tauBarostat * tauBarostat;
117 >      tt2 = tauThermostat * tauThermostat;
118 >      tb2 = tauBarostat * tauBarostat;
119  
120 <    update();
121 < }
120 >      update();
121 >    }
122  
123 < NPT::~NPT() {
124 < }
123 >  NPT::~NPT() {
124 >  }
125  
126 < void NPT::doUpdate() {
126 >  void NPT::doUpdate() {
127  
128      oldPos.resize(info_->getNIntegrableObjects());
129      oldVel.resize(info_->getNIntegrableObjects());
130      oldJi.resize(info_->getNIntegrableObjects());
131  
132 < }
132 >  }
133  
134 < void NPT::moveA() {
134 >  void NPT::moveA() {
135      SimInfo::MoleculeIterator i;
136      Molecule::IntegrableObjectIterator  j;
137      Molecule* mol;
138      StuntDouble* integrableObject;
139      Vector3d Tb, ji;
140 <    double mass;
140 >    RealType mass;
141      Vector3d vel;
142      Vector3d pos;
143      Vector3d frc;
# Line 150 | Line 150 | void NPT::moveA() {
150      
151      instaTemp =thermo.getTemperature();
152      press = thermo.getPressureTensor();
153 <    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
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();
# Line 160 | Line 160 | void NPT::moveA() {
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)) {
163 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
164 >           integrableObject = mol->nextIntegrableObject(j)) {
165                  
166 <            vel = integrableObject->getVel();
167 <            frc = integrableObject->getFrc();
166 >        vel = integrableObject->getVel();
167 >        frc = integrableObject->getFrc();
168  
169 <            mass = integrableObject->getMass();
169 >        mass = integrableObject->getMass();
170  
171 <            getVelScaleA(sc, vel);
171 >        getVelScaleA(sc, vel);
172  
173 <            // velocity half step  (use chi from previous step here):
174 <            //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
175 <            vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
176 <            integrableObject->setVel(vel);
173 >        // velocity half step  (use chi from previous step here):
174 >        //vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
175 >        vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
176 >        integrableObject->setVel(vel);
177  
178 <            if (integrableObject->isDirectional()) {
178 >        if (integrableObject->isDirectional()) {
179  
180 <                // get and convert the torque to body frame
180 >          // get and convert the torque to body frame
181  
182 <                Tb = integrableObject->lab2Body(integrableObject->getTrq());
182 >          Tb = integrableObject->lab2Body(integrableObject->getTrq());
183  
184 <                // get the angular momentum, and propagate a half step
184 >          // get the angular momentum, and propagate a half step
185  
186 <                ji = integrableObject->getJ();
186 >          ji = integrableObject->getJ();
187  
188 <                //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
189 <                ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
188 >          //ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
189 >          ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
190                  
191 <                rotAlgo->rotate(integrableObject, ji, dt);
191 >          rotAlgo->rotate(integrableObject, ji, dt);
192  
193 <                integrableObject->setJ(ji);
194 <            }
193 >          integrableObject->setJ(ji);
194 >        }
195              
196 <        }
196 >      }
197      }
198      // evolve chi and eta  half step
199  
# Line 206 | Line 206 | void NPT::moveA() {
206      
207      index = 0;
208      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
209 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
210 <               integrableObject = mol->nextIntegrableObject(j)) {
211 <            oldPos[index++] = integrableObject->getPos();            
212 <        }
209 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
210 >           integrableObject = mol->nextIntegrableObject(j)) {
211 >        oldPos[index++] = integrableObject->getPos();            
212 >      }
213      }
214      
215      //the first estimation of r(t+dt) is equal to  r(t)
216  
217      for(int k = 0; k < maxIterNum_; k++) {
218 <        index = 0;
219 <        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
220 <            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221 <                   integrableObject = mol->nextIntegrableObject(j)) {
218 >      index = 0;
219 >      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
220 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221 >             integrableObject = mol->nextIntegrableObject(j)) {
222  
223 <                vel = integrableObject->getVel();
224 <                pos = integrableObject->getPos();
223 >          vel = integrableObject->getVel();
224 >          pos = integrableObject->getPos();
225  
226 <                this->getPosScale(pos, COM, index, sc);
226 >          this->getPosScale(pos, COM, index, sc);
227  
228 <                pos = oldPos[index] + dt * (vel + sc);
229 <                integrableObject->setPos(pos);    
228 >          pos = oldPos[index] + dt * (vel + sc);
229 >          integrableObject->setPos(pos);    
230  
231 <                ++index;
232 <           }
233 <        }
231 >          ++index;
232 >        }
233 >      }
234  
235 <        rattle->constraintA();
235 >      rattle->constraintA();
236      }
237  
238      // Scale the box after all the positions have been moved:
# Line 243 | Line 243 | void NPT::moveA() {
243      currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
244  
245      saveEta();
246 < }
246 >  }
247  
248 < void NPT::moveB(void) {
248 >  void NPT::moveB(void) {
249      SimInfo::MoleculeIterator i;
250      Molecule::IntegrableObjectIterator  j;
251      Molecule* mol;
# Line 256 | Line 256 | void NPT::moveB(void) {
256      Vector3d sc;
257      Vector3d vel;
258      Vector3d frc;
259 <    double mass;
259 >    RealType mass;
260  
261  
262      chi= currentSnapshot_->getChi();
263      integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264 <    double oldChi  = chi;
265 <    double prevChi;
264 >    RealType oldChi  = chi;
265 >    RealType prevChi;
266  
267      loadEta();
268      
269      //save velocity and angular momentum
270      index = 0;
271      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
272 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
273 <               integrableObject = mol->nextIntegrableObject(j)) {
272 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
273 >           integrableObject = mol->nextIntegrableObject(j)) {
274                  
275 <            oldVel[index] = integrableObject->getVel();
276 <            oldJi[index] = integrableObject->getJ();
277 <            ++index;
278 <        }
275 >        oldVel[index] = integrableObject->getVel();
276 >        oldJi[index] = integrableObject->getJ();
277 >        ++index;
278 >      }
279      }
280  
281      // do the iteration:
282      instaVol =thermo.getVolume();
283  
284      for(int k = 0; k < maxIterNum_; k++) {
285 <        instaTemp =thermo.getTemperature();
286 <        instaPress =thermo.getPressure();
285 >      instaTemp =thermo.getTemperature();
286 >      instaPress =thermo.getPressure();
287  
288 <        // evolve chi another half step using the temperature at t + dt/2
289 <        prevChi = chi;
290 <        chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
288 >      // evolve chi another half step using the temperature at t + dt/2
289 >      prevChi = chi;
290 >      chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
291  
292 <        //evolve eta
293 <        this->evolveEtaB();
294 <        this->calcVelScale();
292 >      //evolve eta
293 >      this->evolveEtaB();
294 >      this->calcVelScale();
295  
296 <        index = 0;
297 <        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
298 <            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 <                   integrableObject = mol->nextIntegrableObject(j)) {            
296 >      index = 0;
297 >      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
298 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 >             integrableObject = mol->nextIntegrableObject(j)) {            
300  
301 <                frc = integrableObject->getFrc();
302 <                vel = integrableObject->getVel();
301 >          frc = integrableObject->getFrc();
302 >          vel = integrableObject->getVel();
303  
304 <                mass = integrableObject->getMass();
304 >          mass = integrableObject->getMass();
305  
306 <                getVelScaleB(sc, index);
306 >          getVelScaleB(sc, index);
307  
308 <                // velocity half step
309 <                //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
310 <                vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
311 <                integrableObject->setVel(vel);
308 >          // velocity half step
309 >          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
310 >          vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
311 >          integrableObject->setVel(vel);
312  
313 <                if (integrableObject->isDirectional()) {
314 <                    // get and convert the torque to body frame
315 <                    Tb = integrableObject->lab2Body(integrableObject->getTrq());
313 >          if (integrableObject->isDirectional()) {
314 >            // get and convert the torque to body frame
315 >            Tb = integrableObject->lab2Body(integrableObject->getTrq());
316  
317 <                    //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
318 <                    ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
319 <                    integrableObject->setJ(ji);
320 <                }
317 >            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
318 >            ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
319 >            integrableObject->setJ(ji);
320 >          }
321  
322 <                ++index;
323 <            }
324 <        }
322 >          ++index;
323 >        }
324 >      }
325          
326 <        rattle->constraintB();
326 >      rattle->constraintB();
327  
328 <        if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
329 <            break;
328 >      if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
329 >        break;
330      }
331  
332      //calculate integral of chidt
# Line 336 | Line 336 | void NPT::moveB(void) {
336      currentSnapshot_->setIntegralOfChiDt(integralOfChidt);    
337  
338      saveEta();
339 < }
339 >  }
340  
341 +  void NPT::resetIntegrator(){
342 +      currentSnapshot_->setChi(0.0);
343 +      currentSnapshot_->setIntegralOfChiDt(0.0);
344 +      resetEta();
345 +  }
346 +
347 +
348 +    void NPT::resetEta() {
349 +      Mat3x3d etaMat(0.0);
350 +      currentSnapshot_->setEta(etaMat);    
351 +    }
352 +    
353   }

Diff Legend

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