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

Comparing branches/new_design/OOPSE-4/src/integrators/NVT.cpp (file contents):
Revision 1821 by tim, Thu Dec 2 00:26:23 2004 UTC vs.
Revision 1901 by tim, Tue Jan 4 22:18:36 2005 UTC

# Line 1 | Line 1
1 + #include "integrators/IntegratorCreator.hpp"
2   #include "integrators/NVT.hpp"
3   #include "primitives/Molecule.hpp"
4   #include "utils/simError.h"
5   #include "utils/OOPSEConstant.hpp"
6 +
7   namespace oopse {
8  
9   NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
10  
11 <    Globals* globals = info_->getGlobals();
11 >    Globals* simParams = info_->getSimParams();
12  
13 <    if (globals->getUseInitXSstate()) {
13 >    if (simParams->getUseInitXSstate()) {
14          Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
15          currSnapshot->setChi(0.0);
16          currSnapshot->setIntegralOfChiDt(0.0);
17      }
18      
19 <    if (!globals->haveTargetTemp()) {
19 >    if (!simParams->haveTargetTemp()) {
20          sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
21          painCave.isFatal = 1;
22          painCave.severity = OOPSE_ERROR;
23          simError();
24      } else {
25 <        targetTemp_ = globals->getTargetTemp();
25 >        targetTemp_ = simParams->getTargetTemp();
26      }
27  
28      // We must set tauThermostat_.
29  
30 <    if (!globals->haveTauThermostat()) {
30 >    if (!simParams->haveTauThermostat()) {
31          sprintf(painCave.errMsg, "If you use the constant temperature\n"
32                                       "\tintegrator, you must set tauThermostat_.\n");
33  
# Line 33 | Line 35 | NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(inf
35          painCave.isFatal = 1;
36          simError();
37      } else {
38 <        tauThermostat_ = globals->getTauThermostat();
38 >        tauThermostat_ = simParams->getTauThermostat();
39      }
40  
41      update();
42   }
43  
44 < void NVT::update() {
44 > void NVT::doUpdate() {
45      oldVel_.resize(info_->getNIntegrableObjects());
46      oldJi_.resize(info_->getNIntegrableObjects());    
47   }
# Line 55 | Line 57 | void NVT::moveA() {
57      Vector3d pos;
58      Vector3d frc;
59  
60 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
61 <    double chi = currSnapshot->getChi();
60 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
60 >    double chi = currentSnapshot_->getChi();
61 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
62      
63      // We need the temperature at time = t for the chi update below:
64  
# Line 75 | Line 76 | void NVT::moveA() {
76  
77          // velocity half step  (use chi from previous step here):
78          //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
79 <        vel = dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
79 >        vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
80          
81          // position whole step
82          //pos[j] += dt * vel[j];
# Line 86 | Line 87 | void NVT::moveA() {
87  
88          if (integrableObject->isDirectional()) {
89  
90 <            // get and convert the torque to body frame
90 >            //convert the torque to body frame
91 >            Tb = integrableObject->lab2Body(integrableObject->getTrq());
92  
91 –            Tb = integrableObject->getTrq();
92 –            integrableObject->lab2Body(Tb);
93 –
93              // get the angular momentum, and propagate a half step
94  
95              ji = integrableObject->getJ();
# Line 105 | Line 104 | void NVT::moveA() {
104  
105      }
106      
107 <    //constraintAlgorithm->doConstrainA();
107 >    rattle->constraintA();
108  
109      // Finally, evolve chi a half step (just like a velocity) using
110      // temperature at time t, not time t+dt/2
# Line 114 | Line 113 | void NVT::moveA() {
113      chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
114      integralOfChidt += chi * dt2;
115  
116 <    currSnapshot->setChi(chi);
117 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
116 >    currentSnapshot_->setChi(chi);
117 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
118   }
119  
120   void NVT::moveB() {
# Line 133 | Line 132 | void NVT::moveB() {
132      int index;
133      // Set things up for the iteration:
134  
135 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
137 <    double chi = currSnapshot->getChi();
135 >    double chi = currentSnapshot_->getChi();
136      double oldChi = chi;
137      double  prevChi;
138 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
138 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
139  
140      index = 0;
141      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 145 | Line 143 | void NVT::moveB() {
143                 integrableObject = mol->nextIntegrableObject(j)) {
144                  oldVel_[index] = integrableObject->getVel();
145                  oldJi_[index] = integrableObject->getJ();                
146 +
147 +                ++index;    
148          }
149 <        ++index;              
149 >          
150      }
151  
152      // do the iteration:
# Line 180 | Line 180 | void NVT::moveB() {
180  
181                      // get and convert the torque to body frame
182  
183 <                    Tb = integrableObject->getTrq();
184 <                    integrableObject->lab2Body(Tb);
183 >                    Tb =  integrableObject->lab2Body(integrableObject->getTrq());
184  
185                      //for(j = 0; j < 3; j++)
186                      //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
187 <                    ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
187 >                    ji = oldJi_[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
188  
189                      integrableObject->setJ(ji);
190                  }
191 +
192 +
193 +                ++index;
194              }
195          }
196      
197  
198 <        //constraintAlgorithm->doConstrainB();
198 >        rattle->constraintB();
199  
200          if (fabs(prevChi - chi) <= chiTolerance_)
201              break;
202  
201 –        ++index;
203      }
204  
205      integralOfChidt += dt2 * chi;
206  
207 <    currSnapshot->setChi(chi);
208 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
207 >    currentSnapshot_->setChi(chi);
208 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
209   }
210  
211  
212   double NVT::calcConservedQuantity() {
213 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
214 <    double chi = currSnapshot->getChi();
215 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
213 >
214 >    double chi = currentSnapshot_->getChi();
215 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
216      double conservedQuantity;
217      double fkBT;
218      double Energy;

Diff Legend

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