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 1883 by tim, Mon Dec 13 22:30:27 2004 UTC

# Line 1 | Line 1
1 < #inlcude "integrators/NVT.hpp"
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){
6 <    GenericData *data;
7 <    DoubleGenericData *chiValue;
8 <    DoubleGenericData *integralOfChidtValue;
9 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
10  
11 <    chiValue = NULL;
11 <    integralOfChidtValue = NULL;
11 >    Globals* simParams = info_->getSimParams();
12  
13 <    chi = 0.0;
14 <    have_tau_thermostat = 0;
15 <    have_target_temp = 0;
16 <    have_chi_tolerance = 0;
17 <    integralOfChidt = 0.0;
13 >    if (simParams->getUseInitXSstate()) {
14 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
15 >        currSnapshot->setChi(0.0);
16 >        currSnapshot->setIntegralOfChiDt(0.0);
17 >    }
18 >    
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_ = simParams->getTargetTemp();
26 >    }
27  
28 <    if (theInfo->useInitXSstate) {
28 >    // We must set tauThermostat_.
29  
30 <        // retrieve chi and integralOfChidt from simInfo
31 <        data = info->getPropertyByName(CHIVALUE_ID);
30 >    if (!simParams->haveTauThermostat()) {
31 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
32 >                                     "\tintegrator, you must set tauThermostat_.\n");
33  
34 <        if (data) {
35 <            chiValue = dynamic_cast<DoubleGenericData *>(data);
36 <        }
37 <
38 <        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 <        }
34 >        painCave.severity = OOPSE_ERROR;
35 >        painCave.isFatal = 1;
36 >        simError();
37 >    } else {
38 >        tauThermostat_ = simParams->getTauThermostat();
39      }
40  
41 <    oldVel_.resize(info_->getNIntegrableObjects());
42 <    oldJi_.resize(info_->getNIntegrableObjects());
41 >    update();
42   }
43  
44 < NVT::~NVT() {
44 > void NVT::doUpdate() {
45 >    oldVel_.resize(info_->getNIntegrableObjects());
46 >    oldJi_.resize(info_->getNIntegrableObjects());    
47   }
47 –
48   void NVT::moveA() {
49 <    typename SimInfo::MoleculeIterator i;
50 <    typename Molecule::IntegrableObjectIterator  j;
49 >    SimInfo::MoleculeIterator i;
50 >    Molecule::IntegrableObjectIterator  j;
51      Molecule* mol;
52      StuntDouble* integrableObject;
53      Vector3d Tb;
# Line 57 | Line 57 | void NVT::moveA() {
57      Vector3d pos;
58      Vector3d frc;
59  
60 <    double instTemp;
61 <
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  
65 <    instTemp = tStats->getTemperature();
65 >    double instTemp = thermo.getTemperature();
66  
67      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
68          for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 74 | Line 75 | void NVT::moveA() {
75          mass = integrableObject->getMass();
76  
77          // velocity half step  (use chi from previous step here):
78 <        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
79 <        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
78 >        //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
79 >        vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
80          
81          // position whole step
82          //pos[j] += dt * vel[j];
83 <        pos += fullStep * vel;
83 >        pos += dt * vel;
84  
85          integrableObject->setVel(vel);
86          integrableObject->setPos(pos);
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();
96  
97 <            //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
98 <            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
99 <            this->rotationPropagation(integrableObject, ji);
97 >            //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
98 >            ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
99 >            rotAlgo->rotate(integrableObject, ji, dt);
100  
101              integrableObject->setJ(ji);
102          }
# 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
111  
114 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115 –    double chi = currSnapshot->getChi();
116 –    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
112      
113 <    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
114 <    integralOfChidt += chi * halfStep;
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() {
121 <    typename SimInfo::MoleculeIterator i;
122 <    typename Molecule::IntegrableObjectIterator  j;
121 >    SimInfo::MoleculeIterator i;
122 >    Molecule::IntegrableObjectIterator  j;
123      Molecule* mol;
124      StuntDouble* integrableObject;
125      
# Line 134 | Line 129 | void NVT::moveB() {
129      Vector3d frc;
130      double mass;
131      double instTemp;
137 –    double oldChi, prevChi;
132      int index;
133      // Set things up for the iteration:
134  
135 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
136 <    double chi = currSnapshot->getChi();
137 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
138 <    oldChi = chi;
135 >    double chi = currentSnapshot_->getChi();
136 >    double oldChi = chi;
137 >    double  prevChi;
138 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
139  
140      index = 0;
141      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 149 | 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:
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 184 | Line 180 | void NVT::moveB() {
180  
181                      // get and convert the torque to body frame
182  
183 <                    Tb = integrableObject->getTrq();
188 <                    integrableObject->lab2Body(Tb);
183 >                    Tb =  integrableObject->lab2Body(integrableObject->getTrq());
184  
185                      //for(j = 0; j < 3; j++)
186 <                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
187 <                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];
186 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
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 <        if (nConstrained)
201 <            constrainB();
198 >        //constraintAlgorithm->doConstrainB();
199  
200 <        if (fabs(prevChi - chi) <= chiTolerance)
200 >        if (fabs(prevChi - chi) <= chiTolerance_)
201              break;
202  
206 –        ++index;
203      }
204  
205 <    integralOfChidt += halfStep * chi;
205 >    integralOfChidt += dt2 * chi;
206  
207 <    currSnapshot->setChi(chi);
208 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
207 >    currentSnapshot_->setChi(chi);
208 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
209   }
210  
215 – template <typename T>
216 – int NVT<T>::readyCheck() {
211  
212 <    //check parent's readyCheck() first
219 <    if (T::readyCheck() == -1)
220 <        return -1;
212 > double NVT::calcConservedQuantity() {
213  
214 <    // First check to see if we have a target temperature.
215 <    // 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) {
214 >    double chi = currentSnapshot_->getChi();
215 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
216      double conservedQuantity;
217      double fkBT;
218      double Energy;
219      double thermostat_kinetic;
220      double thermostat_potential;
221 +    
222 +    fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;
223  
224 <    fkBT = info_->getNdf() *kB *targetTemp_;
224 >    Energy = thermo.getTotalE();
225  
226 <    Energy = tStats->getTotalE();
226 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
227  
228 <    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
228 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
229  
271 –    thermostat_potential = fkBT * integralOfChidt / eConvert;
272 –
230      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
231  
232      return conservedQuantity;

Diff Legend

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