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:
trunk/OOPSE-3.0/src/integrators/NVT.cpp (file contents), Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/integrators/NVT.cpp (file contents), Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC

# Line 1 | Line 1
1 < #include <math.h>
1 > #inlcude "integrators/NVT.hpp"
2  
3 < #include "Atom.hpp"
4 < #include "SRI.hpp"
5 < #include "AbstractClasses.hpp"
6 < #include "SimInfo.hpp"
7 < #include "ForceFields.hpp"
8 < #include "Thermo.hpp"
9 < #include "ReadWrite.hpp"
10 < #include "Integrator.hpp"
11 < #include "simError.h"
3 > namespace oopse {
4  
5 + NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
6  
7 < // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
7 >    Globals* globals = info_->getGlobals();
8  
9 < template<typename T> NVT<T>::NVT ( SimInfo *theInfo, ForceFields* the_ff):
10 <  T( theInfo, the_ff )
11 < {
12 <  GenericData* data;
13 <  DoubleData * chiValue;
14 <  DoubleData * integralOfChidtValue;
9 >    if (globals->getUseInitXSstate()) {
10 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
11 >        currSnapshot->setChi(0.0);
12 >        currSnapshot->setIntegralOfChiDt(0.0);
13 >    }
14 >    
15 >    if (!globals->haveTargetTemp()) {
16 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
17 >        painCave.isFatal = 1;
18 >        painCave.severity = OOPSE_ERROR;
19 >        simError();
20 >    } else {
21 >        targetTemp_ = globals->getTargetTemp();
22 >    }
23  
24 <  chiValue = NULL;
24 <  integralOfChidtValue = NULL;
24 >    // We must set tauThermostat_.
25  
26 <  chi = 0.0;
27 <  have_tau_thermostat = 0;
28 <  have_target_temp = 0;
29 <  have_chi_tolerance = 0;
30 <  integralOfChidt = 0.0;
26 >    if (!globals->haveTauThermostat()) {
27 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
28 >                                     "\tintegrator, you must set tauThermostat_.\n");
29  
30 <
31 <  if( theInfo->useInitXSstate ){
32 <
33 <    // retrieve chi and integralOfChidt from simInfo
34 <    data = info->getProperty(CHIVALUE_ID);
37 <    if(data){
38 <      chiValue = dynamic_cast<DoubleData*>(data);
30 >        painCave.severity = OOPSE_ERROR;
31 >        painCave.isFatal = 1;
32 >        simError();
33 >    } else {
34 >        tauThermostat_ = globals->getTauThermostat();
35      }
40 –    
41 –    data = info->getProperty(INTEGRALOFCHIDT_ID);
42 –    if(data){
43 –      integralOfChidtValue = dynamic_cast<DoubleData*>(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 –  }
36  
37 <  oldVel = new double[3*integrableObjects.size()];
54 <  oldJi = new double[3*integrableObjects.size()];
37 >    update();
38   }
39  
40 < template<typename T> NVT<T>::~NVT() {
58 <  delete[] oldVel;
59 <  delete[] oldJi;
40 > NVT::~NVT() {
41   }
42  
43 < template<typename T> void NVT<T>::moveA() {
43 > void NVT::update() {
44 >    oldVel_.resize(info_->getNIntegrableObjects());
45 >    oldJi_.resize(info_->getNIntegrableObjects());    
46 > }
47 > void NVT::moveA() {
48 >    typename SimInfo::MoleculeIterator i;
49 >    typename Molecule::IntegrableObjectIterator  j;
50 >    Molecule* mol;
51 >    StuntDouble* integrableObject;
52 >    Vector3d Tb;
53 >    Vector3d ji;
54 >    double mass;
55 >    Vector3d vel;
56 >    Vector3d pos;
57 >    Vector3d frc;
58  
59 <  int i, j;
65 <  DirectionalAtom* dAtom;
66 <  double Tb[3], ji[3];
67 <  double mass;
68 <  double vel[3], pos[3], frc[3];
59 >    double instTemp;
60  
61 <  double instTemp;
61 >    // We need the temperature at time = t for the chi update below:
62  
63 <  // We need the temperature at time = t for the chi update below:
63 >    instTemp = tStats->getTemperature();
64  
65 <  instTemp = tStats->getTemperature();
65 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
66 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
67 >               integrableObject = mol->nextIntegrableObject(j)) {
68  
69 <  for( i=0; i < integrableObjects.size(); i++ ){
69 >        vel = integrableObject->getVel();
70 >        pos = integrableObject->getPos();
71 >        frc = integrableObject->getFrc();
72  
73 <    integrableObjects[i]->getVel( vel );
79 <    integrableObjects[i]->getPos( pos );
80 <    integrableObjects[i]->getFrc( frc );
73 >        mass = integrableObject->getMass();
74  
75 <    mass = integrableObjects[i]->getMass();
75 >        // velocity half step  (use chi from previous step here):
76 >        //vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
77 >        vel = dt2 *eConvert/mass*frc - dt2*chi*vel;
78 >        
79 >        // position whole step
80 >        //pos[j] += dt * vel[j];
81 >        pos += dt * vel;
82  
83 <    for (j=0; j < 3; j++) {
84 <      // 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 <    }
83 >        integrableObject->setVel(vel);
84 >        integrableObject->setPos(pos);
85  
86 <    integrableObjects[i]->setVel( vel );
92 <    integrableObjects[i]->setPos( pos );
86 >        if (integrableObject->isDirectional()) {
87  
88 <    if( integrableObjects[i]->isDirectional() ){
88 >            // get and convert the torque to body frame
89  
90 <      // get and convert the torque to body frame
90 >            Tb = integrableObject->getTrq();
91 >            integrableObject->lab2Body(Tb);
92  
93 <      integrableObjects[i]->getTrq( Tb );
99 <      integrableObjects[i]->lab2Body( Tb );
93 >            // get the angular momentum, and propagate a half step
94  
95 <      // get the angular momentum, and propagate a half step
95 >            ji = integrableObject->getJ();
96  
97 <      integrableObjects[i]->getJ( ji );
97 >            //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
98 >            ji += dt2*eConvert*Tb - dt2*ch *ji;
99 >            this->rotationPropagation(integrableObject, ji);
100  
101 <      for (j=0; j < 3; j++)
102 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
101 >            integrableObject->setJ(ji);
102 >        }
103 >    }
104  
108 –      this->rotationPropagation( integrableObjects[i], ji );
109 –
110 –      integrableObjects[i]->setJ( ji );
105      }
106 <  }
107 <  
114 <  if(nConstrained)
115 <    constrainA();
106 >    
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
109 >    // Finally, evolve chi a half step (just like a velocity) using
110 >    // temperature at time t, not time t+dt/2
111  
112 <  //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
113 <  
114 <  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
115 <  integralOfChidt += chi*dt2;
112 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
113 >    double chi = currSnapshot->getChi();
114 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
115 >    
116 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
117 >    integralOfChidt += chi * dt2;
118  
119 +    currSnapshot->setChi(chi);
120 +    currSnapshot->setIntegralOfChiDt(integralOfChidt);
121   }
122  
123 < template<typename T> void NVT<T>::moveB( void ){
124 <  int i, j, k;
125 <  double Tb[3], ji[3];
126 <  double vel[3], frc[3];
127 <  double mass;
128 <  double instTemp;
129 <  double oldChi, prevChi;
123 > void NVT::moveB() {
124 >    typename SimInfo::MoleculeIterator i;
125 >    typename Molecule::IntegrableObjectIterator  j;
126 >    Molecule* mol;
127 >    StuntDouble* integrableObject;
128 >    
129 >    Vector3d Tb;
130 >    Vector3d ji;    
131 >    Vector3d vel;
132 >    Vector3d frc;
133 >    double mass;
134 >    double instTemp;
135 >    double oldChi, prevChi;
136 >    int index;
137 >    // Set things up for the iteration:
138  
139 <  // Set things up for the iteration:
139 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
140 >    double chi = currSnapshot->getChi();
141 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
142 >    oldChi = chi;
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 = tStats->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 ) * eConvert - oldVel_[3*i + j]*chi);
177 >                vel = oldVel_ + dt2/mass*eConvert * 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] * eConvert - oldJi_[3*i+j]*chi);
190 >                    ji += dt2*eConvert*Tb - dt2*ch *oldJi_[index];
191  
192 <    for( i=0; i < integrableObjects.size(); i++ ){
192 >                    integrableObject->setJ(ji);
193 >                }
194 >            }
195 >        }
196 >    
197  
198 <      integrableObjects[i]->getFrc( frc );
171 <      integrableObjects[i]->getVel(vel);
198 >        //constraintAlgorithm->doConstrainB();
199  
200 <      mass = integrableObjects[i]->getMass();
200 >        if (fabs(prevChi - chi) <= chiTolerance)
201 >            break;
202  
203 <      // 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);
178 <
179 <      integrableObjects[i]->setVel( vel );
180 <
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 <      }
203 >        ++index;
204      }
194 –    
195 –    if(nConstrained)
196 –      constrainB();
205  
206 <    if (fabs(prevChi - chi) <= chiTolerance) break;
199 <  }
206 >    integralOfChidt += dt2 * chi;
207  
208 <  integralOfChidt += dt2*chi;
208 >    currSnapshot->setChi(chi);
209 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);
210   }
211  
204 – template<typename T> void NVT<T>::resetIntegrator( void ){
212  
213 <  chi = 0.0;
214 <  integralOfChidt = 0.0;
215 < }
213 > double NVT::calcConservedQuantity() {
214 >    double conservedQuantity;
215 >    double fkBT;
216 >    double Energy;
217 >    double thermostat_kinetic;
218 >    double thermostat_potential;
219  
220 < template<typename T> int NVT<T>::readyCheck() {
220 >    fkBT = info_->getNdf() *kB *targetTemp_;
221  
222 <  //check parent's readyCheck() first
213 <  if (T::readyCheck() == -1)
214 <    return -1;
222 >    Energy = tStats->getTotalE();
223  
224 <  // First check to see if we have a target temperature.
217 <  // Not having one is fatal.
224 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
225  
226 <  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 <  }
226 >    thermostat_potential = fkBT * integralOfChidt / eConvert;
227  
228 <  // We must set tauThermostat.
228 >    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
229  
230 <  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 <
230 >    return conservedQuantity;
231   }
232  
255 – template<typename T> double NVT<T>::getConservedQuantity(void){
233  
234 <  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 < }
234 > }//end namespace oopse

Diff Legend

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