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 1774 by tim, Tue Nov 23 23:12:23 2004 UTC vs.
Revision 1837 by tim, Thu Dec 2 22:15:31 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 + 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();
# Line 37 | Line 43 | NVT::~NVT() {
43      update();
44   }
45  
40 – NVT::~NVT() {
41 – }
42 –
46   void NVT::update() {
47      oldVel_.resize(info_->getNIntegrableObjects());
48      oldJi_.resize(info_->getNIntegrableObjects());    
49   }
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 56 | Line 59 | void NVT::moveA() {
59      Vector3d pos;
60      Vector3d frc;
61  
62 <    double instTemp;
63 <
62 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
63 >    double chi = currSnapshot->getChi();
64 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
65 >    
66      // We need the temperature at time = t for the chi update below:
67  
68 <    instTemp = tStats->getTemperature();
68 >    double instTemp = thermo.getTemperature();
69  
70      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
71          for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
# Line 73 | Line 78 | void NVT::moveA() {
78          mass = integrableObject->getMass();
79  
80          // velocity half step  (use chi from previous step here):
81 <        //vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
82 <        vel = dt2 *eConvert/mass*frc - dt2*chi*vel;
81 >        //vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
82 >        vel = dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
83          
84          // position whole step
85          //pos[j] += dt * vel[j];
# Line 94 | Line 99 | void NVT::moveA() {
99  
100              ji = integrableObject->getJ();
101  
102 <            //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
103 <            ji += dt2*eConvert*Tb - dt2*ch *ji;
104 <            this->rotationPropagation(integrableObject, ji);
102 >            //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
103 >            ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
104 >            rotAlgo->rotate(integrableObject, ji, dt);
105  
106              integrableObject->setJ(ji);
107          }
# Line 109 | Line 114 | void NVT::moveA() {
114      // Finally, evolve chi a half step (just like a velocity) using
115      // temperature at time t, not time t+dt/2
116  
112 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
113 –    double chi = currSnapshot->getChi();
114 –    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
117      
118      chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
119      integralOfChidt += chi * dt2;
# Line 121 | Line 123 | void NVT::moveB() {
123   }
124  
125   void NVT::moveB() {
126 <    typename SimInfo::MoleculeIterator i;
127 <    typename Molecule::IntegrableObjectIterator  j;
126 >    SimInfo::MoleculeIterator i;
127 >    Molecule::IntegrableObjectIterator  j;
128      Molecule* mol;
129      StuntDouble* integrableObject;
130      
# Line 132 | Line 134 | void NVT::moveB() {
134      Vector3d frc;
135      double mass;
136      double instTemp;
135 –    double oldChi, prevChi;
137      int index;
138      // Set things up for the iteration:
139  
140      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
141      double chi = currSnapshot->getChi();
142 +    double oldChi = chi;
143 +    double  prevChi;
144      double integralOfChidt = currSnapshot->getIntegralOfChiDt();
142 –    oldChi = chi;
145  
146      index = 0;
147      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 155 | Line 157 | void NVT::moveB() {
157  
158      for(int k = 0; k < maxIterNum_; k++) {
159          index = 0;
160 <        instTemp = tStats->getTemperature();
160 >        instTemp = thermo.getTemperature();
161  
162          // evolve chi another half step using the temperature at t + dt/2
163  
# Line 173 | Line 175 | void NVT::moveB() {
175  
176                  // velocity half step
177                  //for(j = 0; j < 3; j++)
178 <                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
179 <                vel = oldVel_ + dt2/mass*eConvert * frc - dt2*chi*oldVel_[index];
178 >                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi);
179 >                vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index];
180              
181                  integrableObject->setVel(vel);
182  
# Line 186 | Line 188 | void NVT::moveB() {
188                      integrableObject->lab2Body(Tb);
189  
190                      //for(j = 0; j < 3; j++)
191 <                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
192 <                    ji += dt2*eConvert*Tb - dt2*ch *oldJi_[index];
191 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
192 >                    ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
193  
194                      integrableObject->setJ(ji);
195                  }
# Line 197 | Line 199 | void NVT::moveB() {
199  
200          //constraintAlgorithm->doConstrainB();
201  
202 <        if (fabs(prevChi - chi) <= chiTolerance)
202 >        if (fabs(prevChi - chi) <= chiTolerance_)
203              break;
204  
205          ++index;
# Line 211 | Line 213 | double NVT::calcConservedQuantity() {
213  
214  
215   double NVT::calcConservedQuantity() {
216 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
217 +    double chi = currSnapshot->getChi();
218 +    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
219      double conservedQuantity;
220      double fkBT;
221      double Energy;
222      double thermostat_kinetic;
223      double thermostat_potential;
224 +    
225 +    fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;
226  
227 <    fkBT = info_->getNdf() *kB *targetTemp_;
227 >    Energy = thermo.getTotalE();
228  
229 <    Energy = tStats->getTotalE();
229 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
230  
231 <    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
231 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
232  
226 –    thermostat_potential = fkBT * integralOfChidt / eConvert;
227 –
233      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
234  
235      return conservedQuantity;

Diff Legend

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