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 1871 by tim, Thu Dec 9 16:22:42 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 > static IntegratorBuilder<NVT>* NVTCreator = new IntegratorBuilder<NVT>("NVT");
10  
11 <    chiValue = NULL;
11 <    integralOfChidtValue = NULL;
11 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
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 >    Globals* simParams = info_->getSimParams();
14  
15 <    if (theInfo->useInitXSstate) {
15 >    if (simParams->getUseInitXSstate()) {
16 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
17 >        currSnapshot->setChi(0.0);
18 >        currSnapshot->setIntegralOfChiDt(0.0);
19 >    }
20 >    
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_ = simParams->getTargetTemp();
28 >    }
29  
30 <        // retrieve chi and integralOfChidt from simInfo
22 <        data = info->getPropertyByName(CHIVALUE_ID);
30 >    // We must set tauThermostat_.
31  
32 <        if (data) {
33 <            chiValue = dynamic_cast<DoubleGenericData *>(data);
34 <        }
32 >    if (!simParams->haveTauThermostat()) {
33 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
34 >                                     "\tintegrator, you must set tauThermostat_.\n");
35  
36 <        data = info->getPropertyByName(INTEGRALOFCHIDT_ID);
37 <
38 <        if (data) {
39 <            integralOfChidtValue = dynamic_cast<DoubleGenericData *>(data);
40 <        }
33 <
34 <        // chi and integralOfChidt should appear by pair
35 <        if (chiValue && integralOfChidtValue) {
36 <            chi = chiValue->getData();
37 <            integralOfChidt = integralOfChidtValue->getData();
38 <        }
36 >        painCave.severity = OOPSE_ERROR;
37 >        painCave.isFatal = 1;
38 >        simError();
39 >    } else {
40 >        tauThermostat_ = simParams->getTauThermostat();
41      }
42  
43 <    oldVel_.resize(info_->getNIntegrableObjects());
42 <    oldJi_.resize(info_->getNIntegrableObjects());
43 >    update();
44   }
45  
46 < NVT::~NVT() {
46 > void NVT::doUpdate() {
47 >    oldVel_.resize(info_->getNIntegrableObjects());
48 >    oldJi_.resize(info_->getNIntegrableObjects());    
49   }
47 –
50   void NVT::moveA() {
51 <    typename SimInfo::MoleculeIterator i;
52 <    typename Molecule::IntegrableObjectIterator  j;
51 >    SimInfo::MoleculeIterator i;
52 >    Molecule::IntegrableObjectIterator  j;
53      Molecule* mol;
54      StuntDouble* integrableObject;
55      Vector3d Tb;
# Line 57 | Line 59 | void NVT::moveA() {
59      Vector3d pos;
60      Vector3d frc;
61  
62 <    double instTemp;
63 <
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  
67 <    instTemp = tStats->getTemperature();
67 >    double instTemp = thermo.getTemperature();
68  
69      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
70          for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 74 | Line 77 | void NVT::moveA() {
77          mass = integrableObject->getMass();
78  
79          // velocity half step  (use chi from previous step here):
80 <        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
81 <        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
80 >        //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
81 >        vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
82          
83          // position whole step
84          //pos[j] += dt * vel[j];
85 <        pos += fullStep * vel;
85 >        pos += dt * vel;
86  
87          integrableObject->setVel(vel);
88          integrableObject->setPos(pos);
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();
98  
99 <            //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
100 <            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
101 <            this->rotationPropagation(integrableObject, ji);
99 >            //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
100 >            ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
101 >            rotAlgo->rotate(integrableObject, ji, dt);
102  
103              integrableObject->setJ(ji);
104          }
# Line 105 | Line 106 | void NVT::moveA() {
106  
107      }
108      
109 <    if (nConstrained)
109 <        constrainA();
109 >    //constraintAlgorithm->doConstrainA();
110  
111      // Finally, evolve chi a half step (just like a velocity) using
112      // temperature at time t, not time t+dt/2
113  
114 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115 –    double chi = currSnapshot->getChi();
116 –    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
114      
115 <    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
116 <    integralOfChidt += chi * halfStep;
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() {
123 <    typename SimInfo::MoleculeIterator i;
124 <    typename Molecule::IntegrableObjectIterator  j;
123 >    SimInfo::MoleculeIterator i;
124 >    Molecule::IntegrableObjectIterator  j;
125      Molecule* mol;
126      StuntDouble* integrableObject;
127      
# Line 134 | Line 131 | void NVT::moveB() {
131      Vector3d frc;
132      double mass;
133      double instTemp;
137 –    double oldChi, prevChi;
134      int index;
135      // Set things up for the iteration:
136  
137 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
138 <    double chi = currSnapshot->getChi();
139 <    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
140 <    oldChi = chi;
137 >    double chi = currentSnapshot_->getChi();
138 >    double oldChi = chi;
139 >    double  prevChi;
140 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
141  
142      index = 0;
143      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 149 | 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:
155  
156      for(int k = 0; k < maxIterNum_; k++) {
157          index = 0;
158 <        instTemp = tStats->getTemperature();
158 >        instTemp = thermo.getTemperature();
159  
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 ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi);
177 >                vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index];
178              
179                  integrableObject->setVel(vel);
180  
# Line 184 | Line 182 | void NVT::moveB() {
182  
183                      // get and convert the torque to body frame
184  
185 <                    Tb = integrableObject->getTrq();
188 <                    integrableObject->lab2Body(Tb);
185 >                    Tb =  integrableObject->lab2Body(integrableObject->getTrq());
186  
187                      //for(j = 0; j < 3; j++)
188 <                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
189 <                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];
188 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
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      
199  
200 <        if (nConstrained)
201 <            constrainB();
200 >        //constraintAlgorithm->doConstrainB();
201  
202 <        if (fabs(prevChi - chi) <= chiTolerance)
202 >        if (fabs(prevChi - chi) <= chiTolerance_)
203              break;
204  
206 –        ++index;
205      }
206  
207 <    integralOfChidt += halfStep * chi;
207 >    integralOfChidt += dt2 * chi;
208  
209 <    currSnapshot->setChi(chi);
210 <    currSnapshot->setIntegralOfChiDt(integralOfChidt);
209 >    currentSnapshot_->setChi(chi);
210 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
211   }
212  
215 – template <typename T>
216 – int NVT<T>::readyCheck() {
213  
214 <    //check parent's readyCheck() first
219 <    if (T::readyCheck() == -1)
220 <        return -1;
214 > double NVT::calcConservedQuantity() {
215  
216 <    // First check to see if we have a target temperature.
217 <    // 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) {
216 >    double chi = currentSnapshot_->getChi();
217 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
218      double conservedQuantity;
219      double fkBT;
220      double Energy;
221      double thermostat_kinetic;
222      double thermostat_potential;
223 +    
224 +    fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;
225  
226 <    fkBT = info_->getNdf() *kB *targetTemp_;
226 >    Energy = thermo.getTotalE();
227  
228 <    Energy = tStats->getTotalE();
228 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
229  
230 <    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
230 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
231  
271 –    thermostat_potential = fkBT * integralOfChidt / eConvert;
272 –
232      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
233  
234      return conservedQuantity;

Diff Legend

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