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:
trunk/OOPSE-4/src/integrators/NVT.cpp (file contents), Revision 1625 by tim, Thu Oct 21 16:22:01 2004 UTC vs.
branches/new_design/OOPSE-4/src/integrators/NVT.cpp (file contents), Revision 1867 by tim, Tue Dec 7 23:08:14 2004 UTC

# Line 1 | Line 1
1 < #include <math.h>
2 <
3 < #include "primitives/Atom.hpp"
4 < #include "primitives/SRI.hpp"
5 < #include "primitives/AbstractClasses.hpp"
6 < #include "brains/SimInfo.hpp"
7 < #include "UseTheForce/ForceFields.hpp"
8 < #include "brains/Thermo.hpp"
9 < #include "io/ReadWrite.hpp"
10 < #include "integrators/Integrator.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 < // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
9 > static IntegratorBuilder<NVT>* NVTCreator = new IntegratorBuilder<NVT>("NVT");
10  
11 < template<typename T> NVT<T>::NVT ( SimInfo *theInfo, ForceFields* the_ff):
17 <  T( theInfo, the_ff )
18 < {
19 <  GenericData* data;
20 <  DoubleGenericData * chiValue;
21 <  DoubleGenericData * integralOfChidtValue;
11 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
12  
13 <  chiValue = NULL;
24 <  integralOfChidtValue = NULL;
13 >    Globals* simParams = info_->getSimParams();
14  
15 <  chi = 0.0;
16 <  have_tau_thermostat = 0;
17 <  have_target_temp = 0;
18 <  have_chi_tolerance = 0;
19 <  integralOfChidt = 0.0;
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 +    // We must set tauThermostat_.
31  
32 <  if( theInfo->useInitXSstate ){
32 >    if (!simParams->haveTauThermostat()) {
33 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
34 >                                     "\tintegrator, you must set tauThermostat_.\n");
35  
36 <    // retrieve chi and integralOfChidt from simInfo
37 <    data = info->getProperty(CHIVALUE_ID);
38 <    if(data){
39 <      chiValue = dynamic_cast<DoubleGenericData*>(data);
36 >        painCave.severity = OOPSE_ERROR;
37 >        painCave.isFatal = 1;
38 >        simError();
39 >    } else {
40 >        tauThermostat_ = simParams->getTauThermostat();
41      }
40 –    
41 –    data = info->getProperty(INTEGRALOFCHIDT_ID);
42 –    if(data){
43 –      integralOfChidtValue = dynamic_cast<DoubleGenericData*>(data);
44 –    }
45 –    
46 –    // chi and integralOfChidt should appear by pair
47 –    if(chiValue && integralOfChidtValue){
48 –      chi = chiValue->getData();
49 –      integralOfChidt = integralOfChidtValue->getData();
50 –    }
51 –  }
42  
43 <  oldVel = new double[3*integrableObjects.size()];
54 <  oldJi = new double[3*integrableObjects.size()];
43 >    update();
44   }
45  
46 < template<typename T> NVT<T>::~NVT() {
47 <  delete[] oldVel;
48 <  delete[] oldJi;
46 > void NVT::doUpdate() {
47 >    oldVel_.resize(info_->getNIntegrableObjects());
48 >    oldJi_.resize(info_->getNIntegrableObjects());    
49   }
50 + void NVT::moveA() {
51 +    SimInfo::MoleculeIterator i;
52 +    Molecule::IntegrableObjectIterator  j;
53 +    Molecule* mol;
54 +    StuntDouble* integrableObject;
55 +    Vector3d Tb;
56 +    Vector3d ji;
57 +    double mass;
58 +    Vector3d vel;
59 +    Vector3d pos;
60 +    Vector3d frc;
61  
62 < template<typename T> void NVT<T>::moveA() {
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 <  int i, j;
65 <  DirectionalAtom* dAtom;
66 <  double Tb[3], ji[3];
67 <  double mass;
68 <  double vel[3], pos[3], frc[3];
67 >    double instTemp = thermo.getTemperature();
68  
69 <  double instTemp;
69 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
70 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
71 >               integrableObject = mol->nextIntegrableObject(j)) {
72  
73 <  // We need the temperature at time = t for the chi update below:
73 >        vel = integrableObject->getVel();
74 >        pos = integrableObject->getPos();
75 >        frc = integrableObject->getFrc();
76  
77 <  instTemp = tStats->getTemperature();
77 >        mass = integrableObject->getMass();
78  
79 <  for( i=0; i < integrableObjects.size(); i++ ){
79 >        // velocity half step  (use chi from previous step here):
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 += dt * vel;
86  
87 <    integrableObjects[i]->getVel( vel );
88 <    integrableObjects[i]->getPos( pos );
80 <    integrableObjects[i]->getFrc( frc );
87 >        integrableObject->setVel(vel);
88 >        integrableObject->setPos(pos);
89  
90 <    mass = integrableObjects[i]->getMass();
90 >        if (integrableObject->isDirectional()) {
91  
92 <    for (j=0; j < 3; j++) {
85 <      // velocity half step  (use chi from previous step here):
86 <      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
87 <      // position whole step
88 <      pos[j] += dt * vel[j];
89 <    }
92 >            // get and convert the torque to body frame
93  
94 <    integrableObjects[i]->setVel( vel );
95 <    integrableObjects[i]->setPos( pos );
94 >            Tb = integrableObject->getTrq();
95 >            integrableObject->lab2Body(Tb);
96  
97 <    if( integrableObjects[i]->isDirectional() ){
97 >            // get the angular momentum, and propagate a half step
98  
99 <      // get and convert the torque to body frame
99 >            ji = integrableObject->getJ();
100  
101 <      integrableObjects[i]->getTrq( Tb );
102 <      integrableObjects[i]->lab2Body( Tb );
101 >            //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
102 >            ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
103 >            rotAlgo->rotate(integrableObject, ji, dt);
104  
105 <      // get the angular momentum, and propagate a half step
105 >            integrableObject->setJ(ji);
106 >        }
107 >    }
108  
103 –      integrableObjects[i]->getJ( ji );
104 –
105 –      for (j=0; j < 3; j++)
106 –        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
107 –
108 –      this->rotationPropagation( integrableObjects[i], ji );
109 –
110 –      integrableObjects[i]->setJ( ji );
109      }
110 <  }
111 <  
114 <  if(nConstrained)
115 <    constrainA();
110 >    
111 >    //constraintAlgorithm->doConstrainA();
112  
113 <  // Finally, evolve chi a half step (just like a velocity) using
114 <  // temperature at time t, not time t+dt/2
113 >    // Finally, evolve chi a half step (just like a velocity) using
114 >    // temperature at time t, not time t+dt/2
115  
116 <  //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
117 <  
118 <  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
123 <  integralOfChidt += chi*dt2;
116 >    
117 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
118 >    integralOfChidt += chi * dt2;
119  
120 +    currentSnapshot_->setChi(chi);
121 +    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
122   }
123  
124 < template<typename T> void NVT<T>::moveB( void ){
125 <  int i, j, k;
126 <  double Tb[3], ji[3];
127 <  double vel[3], frc[3];
128 <  double mass;
129 <  double instTemp;
130 <  double oldChi, prevChi;
124 > void NVT::moveB() {
125 >    SimInfo::MoleculeIterator i;
126 >    Molecule::IntegrableObjectIterator  j;
127 >    Molecule* mol;
128 >    StuntDouble* integrableObject;
129 >    
130 >    Vector3d Tb;
131 >    Vector3d ji;    
132 >    Vector3d vel;
133 >    Vector3d frc;
134 >    double mass;
135 >    double instTemp;
136 >    int index;
137 >    // Set things up for the iteration:
138  
139 <  // Set things up for the iteration:
139 >    double chi = currentSnapshot_->getChi();
140 >    double oldChi = chi;
141 >    double  prevChi;
142 >    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
143  
144 <  oldChi = chi;
144 >    index = 0;
145 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
146 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
147 >               integrableObject = mol->nextIntegrableObject(j)) {
148 >                oldVel_[index] = integrableObject->getVel();
149 >                oldJi_[index] = integrableObject->getJ();                
150 >        }
151 >        ++index;              
152 >    }
153  
154 <  for( i=0; i < integrableObjects.size(); i++ ){
154 >    // do the iteration:
155  
156 <    integrableObjects[i]->getVel( vel );
156 >    for(int k = 0; k < maxIterNum_; k++) {
157 >        index = 0;
158 >        instTemp = thermo.getTemperature();
159  
160 <    for (j=0; j < 3; j++)
144 <      oldVel[3*i + j]  = vel[j];
160 >        // evolve chi another half step using the temperature at t + dt/2
161  
162 <    if( integrableObjects[i]->isDirectional() ){
162 >        prevChi = chi;
163 >        chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
164  
165 <      integrableObjects[i]->getJ( ji );
165 >        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
166 >            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
167 >                   integrableObject = mol->nextIntegrableObject(j)) {
168  
169 <      for (j=0; j < 3; j++)
170 <        oldJi[3*i + j] = ji[j];
169 >                frc = integrableObject->getFrc();
170 >                vel = integrableObject->getVel();
171  
172 <    }
154 <  }
172 >                mass = integrableObject->getMass();
173  
174 <  // do the iteration:
174 >                // velocity half step
175 >                //for(j = 0; j < 3; j++)
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  
181 <  for (k=0; k < 4; k++) {
181 >                if (integrableObject->isDirectional()) {
182  
183 <    instTemp = tStats->getTemperature();
183 >                    // get and convert the torque to body frame
184  
185 <    // evolve chi another half step using the temperature at t + dt/2
185 >                    Tb = integrableObject->getTrq();
186 >                    integrableObject->lab2Body(Tb);
187  
188 <    prevChi = chi;
189 <    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
190 <      (tauThermostat*tauThermostat);
188 >                    //for(j = 0; j < 3; j++)
189 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
190 >                    ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
191  
192 <    for( i=0; i < integrableObjects.size(); i++ ){
192 >                    integrableObject->setJ(ji);
193 >                }
194  
170 –      integrableObjects[i]->getFrc( frc );
171 –      integrableObjects[i]->getVel(vel);
195  
196 <      mass = integrableObjects[i]->getMass();
196 >                ++index;
197 >            }
198 >        }
199 >    
200  
201 <      // velocity half step
176 <      for (j=0; j < 3; j++)
177 <        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*chi);
201 >        //constraintAlgorithm->doConstrainB();
202  
203 <      integrableObjects[i]->setVel( vel );
203 >        if (fabs(prevChi - chi) <= chiTolerance_)
204 >            break;
205  
181 –      if( integrableObjects[i]->isDirectional() ){
182 –
183 –        // get and convert the torque to body frame
184 –
185 –        integrableObjects[i]->getTrq( Tb );
186 –        integrableObjects[i]->lab2Body( Tb );
187 –
188 –        for (j=0; j < 3; j++)
189 –          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
190 –
191 –        integrableObjects[i]->setJ( ji );
192 –      }
206      }
194 –    
195 –    if(nConstrained)
196 –      constrainB();
207  
208 <    if (fabs(prevChi - chi) <= chiTolerance) break;
199 <  }
208 >    integralOfChidt += dt2 * chi;
209  
210 <  integralOfChidt += dt2*chi;
210 >    currentSnapshot_->setChi(chi);
211 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
212   }
213  
204 – template<typename T> void NVT<T>::resetIntegrator( void ){
214  
215 <  chi = 0.0;
207 <  integralOfChidt = 0.0;
208 < }
215 > double NVT::calcConservedQuantity() {
216  
217 < template<typename T> int NVT<T>::readyCheck() {
217 >    double chi = currentSnapshot_->getChi();
218 >    double integralOfChidt = currentSnapshot_->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 <  //check parent's readyCheck() first
213 <  if (T::readyCheck() == -1)
214 <    return -1;
227 >    Energy = thermo.getTotalE();
228  
229 <  // First check to see if we have a target temperature.
217 <  // Not having one is fatal.
229 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
230  
231 <  if (!have_target_temp) {
220 <    sprintf( painCave.errMsg,
221 <             "You can't use the NVT integrator without a targetTemp!\n"
222 <             );
223 <    painCave.isFatal = 1;
224 <    painCave.severity = OOPSE_ERROR;
225 <    simError();
226 <    return -1;
227 <  }
231 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
232  
233 <  // We must set tauThermostat.
233 >    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
234  
235 <  if (!have_tau_thermostat) {
232 <    sprintf( painCave.errMsg,
233 <             "If you use the constant temperature\n"
234 <             "\tintegrator, you must set tauThermostat.\n");
235 <    painCave.severity = OOPSE_ERROR;
236 <    painCave.isFatal = 1;
237 <    simError();
238 <    return -1;
239 <  }
240 <
241 <  if (!have_chi_tolerance) {
242 <    sprintf( painCave.errMsg,
243 <             "In NVT integrator: setting chi tolerance to 1e-6\n");
244 <    chiTolerance = 1e-6;
245 <    have_chi_tolerance = 1;
246 <    painCave.severity = OOPSE_INFO;
247 <    painCave.isFatal = 0;
248 <    simError();
249 <  }
250 <
251 <  return 1;
252 <
235 >    return conservedQuantity;
236   }
237  
255 – template<typename T> double NVT<T>::getConservedQuantity(void){
238  
239 <  double conservedQuantity;
258 <  double fkBT;
259 <  double Energy;
260 <  double thermostat_kinetic;
261 <  double thermostat_potential;
262 <
263 <  fkBT = (double)(info->ndf) * kB * targetTemp;
264 <
265 <  Energy = tStats->getTotalE();
266 <
267 <  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
268 <    (2.0 * eConvert);
269 <
270 <  thermostat_potential = fkBT * integralOfChidt / eConvert;
271 <
272 <  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
273 <
274 <  return conservedQuantity;
275 < }
276 <
277 < template<typename T> string NVT<T>::getAdditionalParameters(void){
278 <  string parameters;
279 <  const int BUFFERSIZE = 2000; // size of the read buffer
280 <  char buffer[BUFFERSIZE];
281 <
282 <  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
283 <  parameters += buffer;
284 <
285 <  return parameters;
286 < }
239 > }//end namespace oopse

Diff Legend

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