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 1765 by tim, Mon Nov 22 20:55:52 2004 UTC vs.
Revision 1821 by tim, Thu Dec 2 00:26:23 2004 UTC

# Line 1 | Line 1
1 < #inlcude "integrators/NVT.hpp"
2 <
1 > #include "integrators/NVT.hpp"
2 > #include "primitives/Molecule.hpp"
3 > #include "utils/simError.h"
4 > #include "utils/OOPSEConstant.hpp"
5   namespace oopse {
6  
7 < NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){
6 <    GenericData *data;
7 <    DoubleGenericData *chiValue;
8 <    DoubleGenericData *integralOfChidtValue;
7 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
8  
9 <    chiValue = NULL;
11 <    integralOfChidtValue = NULL;
9 >    Globals* globals = info_->getGlobals();
10  
11 <    chi = 0.0;
12 <    have_tau_thermostat = 0;
13 <    have_target_temp = 0;
14 <    have_chi_tolerance = 0;
15 <    integralOfChidt = 0.0;
11 >    if (globals->getUseInitXSstate()) {
12 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
13 >        currSnapshot->setChi(0.0);
14 >        currSnapshot->setIntegralOfChiDt(0.0);
15 >    }
16 >    
17 >    if (!globals->haveTargetTemp()) {
18 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
19 >        painCave.isFatal = 1;
20 >        painCave.severity = OOPSE_ERROR;
21 >        simError();
22 >    } else {
23 >        targetTemp_ = globals->getTargetTemp();
24 >    }
25  
26 <    if (theInfo->useInitXSstate) {
26 >    // We must set tauThermostat_.
27  
28 <        // retrieve chi and integralOfChidt from simInfo
29 <        data = info->getPropertyByName(CHIVALUE_ID);
28 >    if (!globals->haveTauThermostat()) {
29 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
30 >                                     "\tintegrator, you must set tauThermostat_.\n");
31  
32 <        if (data) {
33 <            chiValue = dynamic_cast<DoubleGenericData *>(data);
34 <        }
35 <
36 <        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 <        }
32 >        painCave.severity = OOPSE_ERROR;
33 >        painCave.isFatal = 1;
34 >        simError();
35 >    } else {
36 >        tauThermostat_ = globals->getTauThermostat();
37      }
38  
39 <    oldVel_.resize(info_->getNIntegrableObjects());
42 <    oldJi_.resize(info_->getNIntegrableObjects());
39 >    update();
40   }
41  
42 < NVT::~NVT() {
42 > void NVT::update() {
43 >    oldVel_.resize(info_->getNIntegrableObjects());
44 >    oldJi_.resize(info_->getNIntegrableObjects());    
45   }
47 –
46   void NVT::moveA() {
47 <    typename SimInfo::MoleculeIterator i;
48 <    typename Molecule::IntegrableObjectIterator  j;
47 >    SimInfo::MoleculeIterator i;
48 >    Molecule::IntegrableObjectIterator  j;
49      Molecule* mol;
50      StuntDouble* integrableObject;
51      Vector3d Tb;
# Line 57 | Line 55 | void NVT::moveA() {
55      Vector3d pos;
56      Vector3d frc;
57  
58 <    double instTemp;
59 <
58 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
59 >    double chi = currSnapshot->getChi();
60 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
61 >    
62      // We need the temperature at time = t for the chi update below:
63  
64 <    instTemp = tStats->getTemperature();
64 >    double instTemp = thermo.getTemperature();
65  
66      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
67          for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 74 | Line 74 | void NVT::moveA() {
74          mass = integrableObject->getMass();
75  
76          // velocity half step  (use chi from previous step here):
77 <        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
78 <        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
77 >        //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
78 >        vel = dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
79          
80          // position whole step
81          //pos[j] += dt * vel[j];
82 <        pos += fullStep * vel;
82 >        pos += dt * vel;
83  
84          integrableObject->setVel(vel);
85          integrableObject->setPos(pos);
# Line 95 | Line 95 | void NVT::moveA() {
95  
96              ji = integrableObject->getJ();
97  
98 <            //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
99 <            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
100 <            this->rotationPropagation(integrableObject, ji);
98 >            //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
99 >            ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
100 >            rotAlgo->rotate(integrableObject, ji, dt);
101  
102              integrableObject->setJ(ji);
103          }
# Line 105 | Line 105 | void NVT::moveA() {
105  
106      }
107      
108 <    if (nConstrained)
109 <        constrainA();
108 >    //constraintAlgorithm->doConstrainA();
109  
110      // Finally, evolve chi a half step (just like a velocity) using
111      // temperature at time t, not time t+dt/2
112  
113 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
114 <    double chi = currSnapshot->getChi();
115 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
117 <    
118 <    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
119 <    integralOfChidt += chi * halfStep;
113 >    
114 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
115 >    integralOfChidt += chi * dt2;
116  
117      currSnapshot->setChi(chi);
118      currSnapshot->setIntegralOfChiDt(integralOfChidt);
119   }
120  
121   void NVT::moveB() {
122 <    typename SimInfo::MoleculeIterator i;
123 <    typename Molecule::IntegrableObjectIterator  j;
122 >    SimInfo::MoleculeIterator i;
123 >    Molecule::IntegrableObjectIterator  j;
124      Molecule* mol;
125      StuntDouble* integrableObject;
126      
# Line 134 | Line 130 | void NVT::moveB() {
130      Vector3d frc;
131      double mass;
132      double instTemp;
137 –    double oldChi, prevChi;
133      int index;
134      // Set things up for the iteration:
135  
136      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
137      double chi = currSnapshot->getChi();
138 +    double oldChi = chi;
139 +    double  prevChi;
140      double integralOfChidt = currSnapshot->getIntegralOfChiDt();
144 –    oldChi = chi;
141  
142      index = 0;
143      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 157 | Line 153 | void NVT::moveB() {
153  
154      for(int k = 0; k < maxIterNum_; k++) {
155          index = 0;
156 <        instTemp = tStats->getTemperature();
156 >        instTemp = thermo.getTemperature();
157  
158          // evolve chi another half step using the temperature at t + dt/2
159  
160          prevChi = chi;
161 <        chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
161 >        chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
162  
163          for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164              for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 175 | Line 171 | void NVT::moveB() {
171  
172                  // velocity half step
173                  //for(j = 0; j < 3; j++)
174 <                //    vel[j] = oldVel_[3*i+j] + halfStep * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
175 <                vel = oldVel_ + halfStep/mass*eConvert * frc - halfStep*chi*oldVel_[index];
174 >                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi);
175 >                vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index];
176              
177                  integrableObject->setVel(vel);
178  
# Line 188 | Line 184 | void NVT::moveB() {
184                      integrableObject->lab2Body(Tb);
185  
186                      //for(j = 0; j < 3; j++)
187 <                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
188 <                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];
187 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
188 >                    ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
189  
190                      integrableObject->setJ(ji);
191                  }
# Line 197 | Line 193 | void NVT::moveB() {
193          }
194      
195  
196 <        if (nConstrained)
201 <            constrainB();
196 >        //constraintAlgorithm->doConstrainB();
197  
198 <        if (fabs(prevChi - chi) <= chiTolerance)
198 >        if (fabs(prevChi - chi) <= chiTolerance_)
199              break;
200  
201          ++index;
202      }
203  
204 <    integralOfChidt += halfStep * chi;
204 >    integralOfChidt += dt2 * chi;
205  
206      currSnapshot->setChi(chi);
207      currSnapshot->setIntegralOfChiDt(integralOfChidt);
208   }
209  
215 – template <typename T>
216 – int NVT<T>::readyCheck() {
210  
211 <    //check parent's readyCheck() first
212 <    if (T::readyCheck() == -1)
213 <        return -1;
214 <
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) {
211 > double NVT::calcConservedQuantity() {
212 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
213 >    double chi = currSnapshot->getChi();
214 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
215      double conservedQuantity;
216      double fkBT;
217      double Energy;
218      double thermostat_kinetic;
219      double thermostat_potential;
220 +    
221 +    fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;
222  
223 <    fkBT = info_->getNdf() *kB *targetTemp_;
223 >    Energy = thermo.getTotalE();
224  
225 <    Energy = tStats->getTotalE();
225 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
226  
227 <    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
227 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
228  
271 –    thermostat_potential = fkBT * integralOfChidt / eConvert;
272 –
229      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
230  
231      return conservedQuantity;

Diff Legend

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