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 1773 by tim, Mon Nov 22 20:55:52 2004 UTC vs.
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC

# Line 2 | Line 2 | NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(inf
2  
3   namespace oopse {
4  
5 < NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){
6 <    GenericData *data;
7 <    DoubleGenericData *chiValue;
8 <    DoubleGenericData *integralOfChidtValue;
5 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
6  
7 <    chiValue = NULL;
11 <    integralOfChidtValue = NULL;
7 >    Globals* globals = info_->getGlobals();
8  
9 <    chi = 0.0;
10 <    have_tau_thermostat = 0;
11 <    have_target_temp = 0;
12 <    have_chi_tolerance = 0;
13 <    integralOfChidt = 0.0;
9 >    if (globals->getUseInitXSstate()) {
10 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
11 >        currSnapshot->setChi(0.0);
12 >        currSnapshot->setIntegralOfChiDt(0.0);
13 >    }
14 >    
15 >    if (!globals->haveTargetTemp()) {
16 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
17 >        painCave.isFatal = 1;
18 >        painCave.severity = OOPSE_ERROR;
19 >        simError();
20 >    } else {
21 >        targetTemp_ = globals->getTargetTemp();
22 >    }
23  
24 <    if (theInfo->useInitXSstate) {
24 >    // We must set tauThermostat_.
25  
26 <        // retrieve chi and integralOfChidt from simInfo
27 <        data = info->getPropertyByName(CHIVALUE_ID);
26 >    if (!globals->haveTauThermostat()) {
27 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
28 >                                     "\tintegrator, you must set tauThermostat_.\n");
29  
30 <        if (data) {
31 <            chiValue = dynamic_cast<DoubleGenericData *>(data);
32 <        }
33 <
34 <        data = info->getPropertyByName(INTEGRALOFCHIDT_ID);
29 <
30 <        if (data) {
31 <            integralOfChidtValue = dynamic_cast<DoubleGenericData *>(data);
32 <        }
33 <
34 <        // chi and integralOfChidt should appear by pair
35 <        if (chiValue && integralOfChidtValue) {
36 <            chi = chiValue->getData();
37 <            integralOfChidt = integralOfChidtValue->getData();
38 <        }
30 >        painCave.severity = OOPSE_ERROR;
31 >        painCave.isFatal = 1;
32 >        simError();
33 >    } else {
34 >        tauThermostat_ = globals->getTauThermostat();
35      }
36  
37 <    oldVel_.resize(info_->getNIntegrableObjects());
42 <    oldJi_.resize(info_->getNIntegrableObjects());
37 >    update();
38   }
39  
40   NVT::~NVT() {
41   }
42  
43 + void NVT::update() {
44 +    oldVel_.resize(info_->getNIntegrableObjects());
45 +    oldJi_.resize(info_->getNIntegrableObjects());    
46 + }
47   void NVT::moveA() {
48      typename SimInfo::MoleculeIterator i;
49      typename Molecule::IntegrableObjectIterator  j;
# Line 74 | Line 73 | void NVT::moveA() {
73          mass = integrableObject->getMass();
74  
75          // velocity half step  (use chi from previous step here):
76 <        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
77 <        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
76 >        //vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
77 >        vel = dt2 *eConvert/mass*frc - dt2*chi*vel;
78          
79          // position whole step
80          //pos[j] += dt * vel[j];
81 <        pos += fullStep * vel;
81 >        pos += dt * vel;
82  
83          integrableObject->setVel(vel);
84          integrableObject->setPos(pos);
# Line 95 | Line 94 | void NVT::moveA() {
94  
95              ji = integrableObject->getJ();
96  
97 <            //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
98 <            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
97 >            //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
98 >            ji += dt2*eConvert*Tb - dt2*ch *ji;
99              this->rotationPropagation(integrableObject, ji);
100  
101              integrableObject->setJ(ji);
# Line 105 | Line 104 | void NVT::moveA() {
104  
105      }
106      
107 <    if (nConstrained)
109 <        constrainA();
107 >    //constraintAlgorithm->doConstrainA();
108  
109      // Finally, evolve chi a half step (just like a velocity) using
110      // temperature at time t, not time t+dt/2
# Line 115 | Line 113 | void NVT::moveA() {
113      double chi = currSnapshot->getChi();
114      double integralOfChidt = currSnapshot->getIntegralOfChiDt();
115      
116 <    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
117 <    integralOfChidt += chi * halfStep;
116 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
117 >    integralOfChidt += chi * dt2;
118  
119      currSnapshot->setChi(chi);
120      currSnapshot->setIntegralOfChiDt(integralOfChidt);
# Line 162 | Line 160 | void NVT::moveB() {
160          // evolve chi another half step using the temperature at t + dt/2
161  
162          prevChi = chi;
163 <        chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
163 >        chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
164  
165          for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
166              for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 175 | Line 173 | void NVT::moveB() {
173  
174                  // velocity half step
175                  //for(j = 0; j < 3; j++)
176 <                //    vel[j] = oldVel_[3*i+j] + halfStep * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
177 <                vel = oldVel_ + halfStep/mass*eConvert * frc - halfStep*chi*oldVel_[index];
176 >                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
177 >                vel = oldVel_ + dt2/mass*eConvert * frc - dt2*chi*oldVel_[index];
178              
179                  integrableObject->setVel(vel);
180  
# Line 188 | Line 186 | void NVT::moveB() {
186                      integrableObject->lab2Body(Tb);
187  
188                      //for(j = 0; j < 3; j++)
189 <                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
190 <                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];
189 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
190 >                    ji += dt2*eConvert*Tb - dt2*ch *oldJi_[index];
191  
192                      integrableObject->setJ(ji);
193                  }
# Line 197 | Line 195 | void NVT::moveB() {
195          }
196      
197  
198 <        if (nConstrained)
201 <            constrainB();
198 >        //constraintAlgorithm->doConstrainB();
199  
200          if (fabs(prevChi - chi) <= chiTolerance)
201              break;
# Line 206 | Line 203 | void NVT::moveB() {
203          ++index;
204      }
205  
206 <    integralOfChidt += halfStep * chi;
206 >    integralOfChidt += dt2 * chi;
207  
208      currSnapshot->setChi(chi);
209      currSnapshot->setIntegralOfChiDt(integralOfChidt);
210   }
211  
215 – template <typename T>
216 – int NVT<T>::readyCheck() {
212  
213 <    //check parent's readyCheck() first
219 <    if (T::readyCheck() == -1)
220 <        return -1;
221 <
222 <    // First check to see if we have a target temperature.
223 <    // Not having one is fatal.
224 <
225 <    if (!have_target_temp) {
226 <        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
227 <        painCave.isFatal = 1;
228 <        painCave.severity = OOPSE_ERROR;
229 <        simError();
230 <        return -1;
231 <    }
232 <
233 <    // We must set tauThermostat_.
234 <
235 <    if (!have_tau_thermostat) {
236 <        sprintf(painCave.errMsg, "If you use the constant temperature\n"
237 <                                     "\tintegrator, you must set tauThermostat_.\n");
238 <
239 <        painCave.severity = OOPSE_ERROR;
240 <        painCave.isFatal = 1;
241 <        simError();
242 <        return -1;
243 <    }
244 <
245 <    if (!have_chi_tolerance) {
246 <        sprintf(painCave.errMsg, "In NVT integrator: setting chi tolerance to 1e-6\n");
247 <        chiTolerance = 1e - 6;
248 <        have_chi_tolerance = 1;
249 <        painCave.severity = OOPSE_INFO;
250 <        painCave.isFatal = 0;
251 <        simError();
252 <    }
253 <
254 <    return 1;
255 < }
256 <
257 < template <typename T>
258 < double NVT<T>::getConservedQuantity(void) {
213 > double NVT::calcConservedQuantity() {
214      double conservedQuantity;
215      double fkBT;
216      double Energy;

Diff Legend

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