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 1871 by tim, Thu Dec 9 16:22:42 2004 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 + static IntegratorBuilder<NVT>* NVTCreator = new IntegratorBuilder<NVT>("NVT");
10 +
11   NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
12  
13 <    Globals* globals = info_->getGlobals();
13 >    Globals* simParams = info_->getSimParams();
14  
15 <    if (globals->getUseInitXSstate()) {
15 >    if (simParams->getUseInitXSstate()) {
16          Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
17          currSnapshot->setChi(0.0);
18          currSnapshot->setIntegralOfChiDt(0.0);
19      }
20      
21 <    if (!globals->haveTargetTemp()) {
21 >    if (!simParams->haveTargetTemp()) {
22          sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
23          painCave.isFatal = 1;
24          painCave.severity = OOPSE_ERROR;
25          simError();
26      } else {
27 <        targetTemp_ = globals->getTargetTemp();
27 >        targetTemp_ = simParams->getTargetTemp();
28      }
29  
30      // We must set tauThermostat_.
31  
32 <    if (!globals->haveTauThermostat()) {
32 >    if (!simParams->haveTauThermostat()) {
33          sprintf(painCave.errMsg, "If you use the constant temperature\n"
34                                       "\tintegrator, you must set tauThermostat_.\n");
35  
# Line 33 | Line 37 | NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(inf
37          painCave.isFatal = 1;
38          simError();
39      } else {
40 <        tauThermostat_ = globals->getTauThermostat();
40 >        tauThermostat_ = simParams->getTauThermostat();
41      }
42  
43      update();
44   }
45  
46 < void NVT::update() {
46 > void NVT::doUpdate() {
47      oldVel_.resize(info_->getNIntegrableObjects());
48      oldJi_.resize(info_->getNIntegrableObjects());    
49   }
# Line 55 | Line 59 | void NVT::moveA() {
59      Vector3d pos;
60      Vector3d frc;
61  
62 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
63 <    double chi = currSnapshot->getChi();
60 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
62 >    double chi = currentSnapshot_->getChi();
63 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
64      
65      // We need the temperature at time = t for the chi update below:
66  
# Line 75 | Line 78 | void NVT::moveA() {
78  
79          // velocity half step  (use chi from previous step here):
80          //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
81 <        vel = dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
81 >        vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
82          
83          // position whole step
84          //pos[j] += dt * vel[j];
# Line 86 | Line 89 | void NVT::moveA() {
89  
90          if (integrableObject->isDirectional()) {
91  
92 <            // get and convert the torque to body frame
92 >            //convert the torque to body frame
93 >            Tb = integrableObject->lab2Body(integrableObject->getTrq());
94  
91 –            Tb = integrableObject->getTrq();
92 –            integrableObject->lab2Body(Tb);
93 –
95              // get the angular momentum, and propagate a half step
96  
97              ji = integrableObject->getJ();
# Line 114 | Line 115 | void NVT::moveA() {
115      chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
116      integralOfChidt += chi * dt2;
117  
118 <    currSnapshot->setChi(chi);
119 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
118 >    currentSnapshot_->setChi(chi);
119 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
120   }
121  
122   void NVT::moveB() {
# Line 133 | Line 134 | void NVT::moveB() {
134      int index;
135      // Set things up for the iteration:
136  
137 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
137 <    double chi = currSnapshot->getChi();
137 >    double chi = currentSnapshot_->getChi();
138      double oldChi = chi;
139      double  prevChi;
140 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
140 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
141  
142      index = 0;
143      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 145 | Line 145 | void NVT::moveB() {
145                 integrableObject = mol->nextIntegrableObject(j)) {
146                  oldVel_[index] = integrableObject->getVel();
147                  oldJi_[index] = integrableObject->getJ();                
148 +
149 +                ++index;    
150          }
151 <        ++index;              
151 >          
152      }
153  
154      // do the iteration:
# Line 180 | Line 182 | void NVT::moveB() {
182  
183                      // get and convert the torque to body frame
184  
185 <                    Tb = integrableObject->getTrq();
184 <                    integrableObject->lab2Body(Tb);
185 >                    Tb =  integrableObject->lab2Body(integrableObject->getTrq());
186  
187                      //for(j = 0; j < 3; j++)
188                      //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
189 <                    ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
189 >                    ji = oldJi_[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
190  
191                      integrableObject->setJ(ji);
192                  }
193 +
194 +
195 +                ++index;
196              }
197          }
198      
# Line 198 | Line 202 | void NVT::moveB() {
202          if (fabs(prevChi - chi) <= chiTolerance_)
203              break;
204  
201 –        ++index;
205      }
206  
207      integralOfChidt += dt2 * chi;
208  
209 <    currSnapshot->setChi(chi);
210 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
209 >    currentSnapshot_->setChi(chi);
210 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
211   }
212  
213  
214   double NVT::calcConservedQuantity() {
215 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
216 <    double chi = currSnapshot->getChi();
217 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
215 >
216 >    double chi = currentSnapshot_->getChi();
217 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
218      double conservedQuantity;
219      double fkBT;
220      double Energy;

Diff Legend

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