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

Comparing trunk/src/integrators/NVT.cpp (file contents):
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

# Line 36 | Line 36
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).                        
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 "integrators/NVT.hpp"
# Line 51 | Line 52 | namespace OpenMD {
52      Globals* simParams = info_->getSimParams();
53  
54      if (!simParams->getUseIntialExtendedSystemState()) {
55 <      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
56 <      currSnapshot->setChi(0.0);
56 <      currSnapshot->setIntegralOfChiDt(0.0);
55 >      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
56 >      snap->setThermostat(make_pair(0.0, 0.0));
57      }
58      
59      if (!simParams->haveTargetTemp()) {
# Line 78 | Line 78 | namespace OpenMD {
78        tauThermostat_ = simParams->getTauThermostat();
79      }
80  
81 <    update();
81 >    updateSizes();
82    }
83  
84 <  void NVT::doUpdate() {
84 >  void NVT::doUpdateSizes() {
85      oldVel_.resize(info_->getNIntegrableObjects());
86 <    oldJi_.resize(info_->getNIntegrableObjects());    
86 >    oldJi_.resize(info_->getNIntegrableObjects());
87    }
88 +
89    void NVT::moveA() {
90      SimInfo::MoleculeIterator i;
91      Molecule::IntegrableObjectIterator  j;
92      Molecule* mol;
93 <    StuntDouble* integrableObject;
93 >    StuntDouble* sd;
94      Vector3d Tb;
95      Vector3d ji;
96      RealType mass;
# Line 97 | Line 98 | namespace OpenMD {
98      Vector3d pos;
99      Vector3d frc;
100  
101 <    RealType chi = currentSnapshot_->getChi();
102 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
102 <    
101 >    pair<RealType, RealType> thermostat = snap->getThermostat();
102 >
103      // We need the temperature at time = t for the chi update below:
104  
105      RealType instTemp = thermo.getTemperature();
106  
107 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
108 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
109 <           integrableObject = mol->nextIntegrableObject(j)) {
107 >    for (mol = info_->beginMolecule(i); mol != NULL;
108 >         mol = info_->nextMolecule(i)) {
109  
110 <        vel = integrableObject->getVel();
111 <        pos = integrableObject->getPos();
113 <        frc = integrableObject->getFrc();
110 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
111 >           sd = mol->nextIntegrableObject(j)) {
112  
113 <        mass = integrableObject->getMass();
113 >        vel = sd->getVel();
114 >        pos = sd->getPos();
115 >        frc = sd->getFrc();
116  
117 <        // velocity half step  (use chi from previous step here):
118 <        //vel[j] += dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - vel[j]*chi);
119 <        vel += dt2 *PhysicalConstants::energyConvert/mass*frc - dt2*chi*vel;
117 >        mass = sd->getMass();
118 >
119 >        // velocity half step (use chi from previous step here):
120 >        vel += dt2 *PhysicalConstants::energyConvert/mass*frc
121 >          - dt2*thermostat.first*vel;
122          
123          // position whole step
122 –        //pos[j] += dt * vel[j];
124          pos += dt * vel;
125  
126 <        integrableObject->setVel(vel);
127 <        integrableObject->setPos(pos);
126 >        sd->setVel(vel);
127 >        sd->setPos(pos);
128  
129 <        if (integrableObject->isDirectional()) {
129 >        if (sd->isDirectional()) {
130  
131            //convert the torque to body frame
132 <          Tb = integrableObject->lab2Body(integrableObject->getTrq());
132 >          Tb = sd->lab2Body(sd->getTrq());
133  
134            // get the angular momentum, and propagate a half step
135  
136 <          ji = integrableObject->getJ();
136 >          ji = sd->getJ();
137  
138 <          //ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
139 <          ji += dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *ji;
139 <          rotAlgo->rotate(integrableObject, ji, dt);
138 >          ji += dt2*PhysicalConstants::energyConvert*Tb
139 >            - dt2*thermostat.first *ji;
140  
141 <          integrableObject->setJ(ji);
141 >          rotAlgo_->rotate(sd, ji, dt);
142 >
143 >          sd->setJ(ji);
144          }
145        }
146  
147      }
148      
149 <    rattle->constraintA();
149 >    flucQ_->moveA();
150 >    rattle_->constraintA();
151  
152      // Finally, evolve chi a half step (just like a velocity) using
153      // temperature at time t, not time t+dt/2
154  
155 <    
156 <    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
157 <    integralOfChidt += chi * dt2;
155 >    thermostat.first += dt2 * (instTemp / targetTemp_ - 1.0)
156 >      / (tauThermostat_ * tauThermostat_);
157 >    thermostat.second += thermostat.first * dt2;
158  
159 <    currentSnapshot_->setChi(chi);
157 <    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
159 >    snap->setThermostat(thermostat);
160    }
161  
162    void NVT::moveB() {
163      SimInfo::MoleculeIterator i;
164      Molecule::IntegrableObjectIterator  j;
165      Molecule* mol;
166 <    StuntDouble* integrableObject;
166 >    StuntDouble* sd;
167      
168      Vector3d Tb;
169      Vector3d ji;    
# Line 172 | Line 174 | namespace OpenMD {
174      int index;
175      // Set things up for the iteration:
176  
177 <    RealType chi = currentSnapshot_->getChi();
178 <    RealType oldChi = chi;
177 >    pair<RealType, RealType> thermostat = snap->getThermostat();
178 >    RealType oldChi = thermostat.first;
179      RealType  prevChi;
178 –    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
180  
181      index = 0;
182 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
183 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
183 <           integrableObject = mol->nextIntegrableObject(j)) {
184 <        oldVel_[index] = integrableObject->getVel();
185 <        oldJi_[index] = integrableObject->getJ();                
182 >    for (mol = info_->beginMolecule(i); mol != NULL;
183 >         mol = info_->nextMolecule(i)) {
184  
185 +      for (sd = mol->beginIntegrableObject(j); sd != NULL;
186 +           sd = mol->nextIntegrableObject(j)) {
187 +
188 +        oldVel_[index] = sd->getVel();
189 +        
190 +        if (sd->isDirectional())
191 +          oldJi_[index] = sd->getJ();                
192 +        
193          ++index;    
194 <      }
189 <          
194 >      }          
195      }
196  
197      // do the iteration:
# Line 197 | Line 202 | namespace OpenMD {
202  
203        // evolve chi another half step using the temperature at t + dt/2
204  
205 <      prevChi = chi;
206 <      chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
205 >      prevChi = thermostat.first;
206 >      thermostat.first = oldChi + dt2 * (instTemp / targetTemp_ - 1.0)
207 >        / (tauThermostat_ * tauThermostat_);
208  
209 <      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
210 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
211 <             integrableObject = mol->nextIntegrableObject(j)) {
209 >      for (mol = info_->beginMolecule(i); mol != NULL;
210 >           mol = info_->nextMolecule(i)) {
211 >        
212 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
213 >             sd = mol->nextIntegrableObject(j)) {
214  
215 <          frc = integrableObject->getFrc();
216 <          vel = integrableObject->getVel();
215 >          frc = sd->getFrc();
216 >          vel = sd->getVel();
217  
218 <          mass = integrableObject->getMass();
218 >          mass = sd->getMass();
219  
220            // velocity half step
221 <          //for(j = 0; j < 3; j++)
222 <          //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - oldVel_[3*i + j]*chi);
223 <          vel = oldVel_[index] + dt2/mass*PhysicalConstants::energyConvert * frc - dt2*chi*oldVel_[index];
221 >
222 >          vel = oldVel_[index]
223 >            + dt2/mass*PhysicalConstants::energyConvert * frc
224 >            - dt2*thermostat.first*oldVel_[index];
225              
226 <          integrableObject->setVel(vel);
226 >          sd->setVel(vel);
227  
228 <          if (integrableObject->isDirectional()) {
228 >          if (sd->isDirectional()) {
229  
230              // get and convert the torque to body frame
231  
232 <            Tb =  integrableObject->lab2Body(integrableObject->getTrq());
232 >            Tb =  sd->lab2Body(sd->getTrq());
233  
234 <            //for(j = 0; j < 3; j++)
235 <            //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi_[3*i+j]*chi);
227 <            ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *oldJi_[index];
234 >            ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb
235 >              - dt2*thermostat.first *oldJi_[index];
236  
237 <            integrableObject->setJ(ji);
237 >            sd->setJ(ji);
238            }
239  
240  
# Line 234 | Line 242 | namespace OpenMD {
242          }
243        }
244      
245 +      rattle_->constraintB();
246  
247 <      rattle->constraintB();
239 <
240 <      if (fabs(prevChi - chi) <= chiTolerance_)
247 >      if (fabs(prevChi - thermostat.first) <= chiTolerance_)
248          break;
249  
250      }
251  
252 <    integralOfChidt += dt2 * chi;
252 >    flucQ_->moveB();
253  
254 <    currentSnapshot_->setChi(chi);
255 <    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
254 >    thermostat.second += dt2 * thermostat.first;
255 >    snap->setThermostat(thermostat);
256    }
257  
258    void NVT::resetIntegrator() {
259 <      currentSnapshot_->setChi(0.0);
253 <      currentSnapshot_->setIntegralOfChiDt(0.0);
259 >    snap->setThermostat(make_pair(0.0, 0.0));
260    }
261    
262    RealType NVT::calcConservedQuantity() {
263  
264 <    RealType chi = currentSnapshot_->getChi();
259 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264 >    pair<RealType, RealType> thermostat = snap->getThermostat();
265      RealType conservedQuantity;
266      RealType fkBT;
267      RealType Energy;
# Line 265 | Line 270 | namespace OpenMD {
270      
271      fkBT = info_->getNdf() *PhysicalConstants::kB *targetTemp_;
272  
273 <    Energy = thermo.getTotalE();
273 >    Energy = thermo.getTotalEnergy();
274  
275 <    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * PhysicalConstants::energyConvert);
275 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * thermostat.first * thermostat.first / (2.0 * PhysicalConstants::energyConvert);
276  
277 <    thermostat_potential = fkBT * integralOfChidt / PhysicalConstants::energyConvert;
277 >    thermostat_potential = fkBT * thermostat.second / PhysicalConstants::energyConvert;
278  
279      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
280  

Diff Legend

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