ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/NVT.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/integrators/NVT.cpp (file contents):
Revision 1820 by tim, Tue Nov 23 23:12:23 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), chiTolerance_ (1e-6) {
# Line 37 | Line 39 | NVT::~NVT() {
39      update();
40   }
41  
40 – NVT::~NVT() {
41 – }
42 –
42   void NVT::update() {
43      oldVel_.resize(info_->getNIntegrableObjects());
44      oldJi_.resize(info_->getNIntegrableObjects());    
45   }
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 56 | 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 73 | Line 74 | void NVT::moveA() {
74          mass = integrableObject->getMass();
75  
76          // velocity half step  (use chi from previous step here):
77 <        //vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
78 <        vel = dt2 *eConvert/mass*frc - dt2*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];
# Line 94 | Line 95 | void NVT::moveA() {
95  
96              ji = integrableObject->getJ();
97  
98 <            //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
99 <            ji += dt2*eConvert*Tb - dt2*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 109 | Line 110 | void NVT::moveA() {
110      // Finally, evolve chi a half step (just like a velocity) using
111      // temperature at time t, not time t+dt/2
112  
112 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
113 –    double chi = currSnapshot->getChi();
114 –    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
113      
114      chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
115      integralOfChidt += chi * dt2;
# Line 121 | Line 119 | void NVT::moveB() {
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 132 | Line 130 | void NVT::moveB() {
130      Vector3d frc;
131      double mass;
132      double instTemp;
135 –    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();
142 –    oldChi = chi;
141  
142      index = 0;
143      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
# Line 155 | 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  
# Line 173 | Line 171 | void NVT::moveB() {
171  
172                  // velocity half step
173                  //for(j = 0; j < 3; j++)
174 <                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
175 <                vel = oldVel_ + dt2/mass*eConvert * frc - dt2*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 186 | Line 184 | void NVT::moveB() {
184                      integrableObject->lab2Body(Tb);
185  
186                      //for(j = 0; j < 3; j++)
187 <                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
188 <                    ji += dt2*eConvert*Tb - dt2*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 195 | void NVT::moveB() {
195  
196          //constraintAlgorithm->doConstrainB();
197  
198 <        if (fabs(prevChi - chi) <= chiTolerance)
198 >        if (fabs(prevChi - chi) <= chiTolerance_)
199              break;
200  
201          ++index;
# Line 211 | Line 209 | double NVT::calcConservedQuantity() {
209  
210  
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  
226 –    thermostat_potential = fkBT * integralOfChidt / eConvert;
227 –
229      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
230  
231      return conservedQuantity;

Diff Legend

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