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 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/new_design/OOPSE-4/src/integrators/NVT.cpp (file contents), Revision 1821 by tim, Thu Dec 2 00:26:23 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/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) {
8  
9 < // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
9 >    Globals* globals = info_->getGlobals();
10  
11 < template<typename T> NVT<T>::NVT ( SimInfo *theInfo, ForceFields* the_ff):
12 <  T( theInfo, the_ff )
13 < {
14 <  GenericData* data;
15 <  DoubleData * chiValue;
16 <  DoubleData * integralOfChidtValue;
11 >    if (globals->getUseInitXSstate()) {
12 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
13 >        currSnapshot->setChi(0.0);
14 >        currSnapshot->setIntegralOfChiDt(0.0);
15 >    }
16 >    
17 >    if (!globals->haveTargetTemp()) {
18 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
19 >        painCave.isFatal = 1;
20 >        painCave.severity = OOPSE_ERROR;
21 >        simError();
22 >    } else {
23 >        targetTemp_ = globals->getTargetTemp();
24 >    }
25  
26 <  chiValue = NULL;
24 <  integralOfChidtValue = NULL;
26 >    // We must set tauThermostat_.
27  
28 <  chi = 0.0;
29 <  have_tau_thermostat = 0;
30 <  have_target_temp = 0;
29 <  have_chi_tolerance = 0;
30 <  integralOfChidt = 0.0;
28 >    if (!globals->haveTauThermostat()) {
29 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
30 >                                     "\tintegrator, you must set tauThermostat_.\n");
31  
32 <
33 <  if( theInfo->useInitXSstate ){
34 <
35 <    // retrieve chi and integralOfChidt from simInfo
36 <    data = info->getProperty(CHIVALUE_ID);
37 <    if(data){
38 <      chiValue = dynamic_cast<DoubleData*>(data);
32 >        painCave.severity = OOPSE_ERROR;
33 >        painCave.isFatal = 1;
34 >        simError();
35 >    } else {
36 >        tauThermostat_ = globals->getTauThermostat();
37      }
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 –  }
38  
39 <  oldVel = new double[3*integrableObjects.size()];
54 <  oldJi = new double[3*integrableObjects.size()];
39 >    update();
40   }
41  
42 < template<typename T> NVT<T>::~NVT() {
43 <  delete[] oldVel;
44 <  delete[] oldJi;
42 > void NVT::update() {
43 >    oldVel_.resize(info_->getNIntegrableObjects());
44 >    oldJi_.resize(info_->getNIntegrableObjects());    
45   }
46 + void NVT::moveA() {
47 +    SimInfo::MoleculeIterator i;
48 +    Molecule::IntegrableObjectIterator  j;
49 +    Molecule* mol;
50 +    StuntDouble* integrableObject;
51 +    Vector3d Tb;
52 +    Vector3d ji;
53 +    double mass;
54 +    Vector3d vel;
55 +    Vector3d pos;
56 +    Vector3d frc;
57  
58 < template<typename T> void NVT<T>::moveA() {
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 <  int i, j;
65 <  DirectionalAtom* dAtom;
66 <  double Tb[3], ji[3];
67 <  double mass;
68 <  double vel[3], pos[3], frc[3];
64 >    double instTemp = thermo.getTemperature();
65  
66 <  double instTemp;
66 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
67 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
68 >               integrableObject = mol->nextIntegrableObject(j)) {
69  
70 <  // We need the temperature at time = t for the chi update below:
70 >        vel = integrableObject->getVel();
71 >        pos = integrableObject->getPos();
72 >        frc = integrableObject->getFrc();
73  
74 <  instTemp = tStats->getTemperature();
74 >        mass = integrableObject->getMass();
75  
76 <  for( i=0; i < integrableObjects.size(); i++ ){
76 >        // velocity half step  (use chi from previous step here):
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];
82 >        pos += dt * vel;
83  
84 <    integrableObjects[i]->getVel( vel );
85 <    integrableObjects[i]->getPos( pos );
80 <    integrableObjects[i]->getFrc( frc );
84 >        integrableObject->setVel(vel);
85 >        integrableObject->setPos(pos);
86  
87 <    mass = integrableObjects[i]->getMass();
87 >        if (integrableObject->isDirectional()) {
88  
89 <    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 <    }
89 >            // get and convert the torque to body frame
90  
91 <    integrableObjects[i]->setVel( vel );
92 <    integrableObjects[i]->setPos( pos );
91 >            Tb = integrableObject->getTrq();
92 >            integrableObject->lab2Body(Tb);
93  
94 <    if( integrableObjects[i]->isDirectional() ){
94 >            // get the angular momentum, and propagate a half step
95  
96 <      // get and convert the torque to body frame
96 >            ji = integrableObject->getJ();
97  
98 <      integrableObjects[i]->getTrq( Tb );
99 <      integrableObjects[i]->lab2Body( Tb );
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 <      // get the angular momentum, and propagate a half step
102 >            integrableObject->setJ(ji);
103 >        }
104 >    }
105  
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 );
106      }
107 <  }
108 <  
114 <  if(nConstrained)
115 <    constrainA();
107 >    
108 >    //constraintAlgorithm->doConstrainA();
109  
110 <  // Finally, evolve chi a half step (just like a velocity) using
111 <  // temperature at time t, not time t+dt/2
110 >    // Finally, evolve chi a half step (just like a velocity) using
111 >    // temperature at time t, not time t+dt/2
112  
113 <  //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
114 <  
115 <  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
123 <  integralOfChidt += chi*dt2;
113 >    
114 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
115 >    integralOfChidt += chi * dt2;
116  
117 +    currSnapshot->setChi(chi);
118 +    currSnapshot->setIntegralOfChiDt(integralOfChidt);
119   }
120  
121 < template<typename T> void NVT<T>::moveB( void ){
122 <  int i, j, k;
123 <  double Tb[3], ji[3];
124 <  double vel[3], frc[3];
125 <  double mass;
126 <  double instTemp;
127 <  double oldChi, prevChi;
121 > void NVT::moveB() {
122 >    SimInfo::MoleculeIterator i;
123 >    Molecule::IntegrableObjectIterator  j;
124 >    Molecule* mol;
125 >    StuntDouble* integrableObject;
126 >    
127 >    Vector3d Tb;
128 >    Vector3d ji;    
129 >    Vector3d vel;
130 >    Vector3d frc;
131 >    double mass;
132 >    double instTemp;
133 >    int index;
134 >    // Set things up for the iteration:
135  
136 <  // Set things up for the iteration:
136 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
137 >    double chi = currSnapshot->getChi();
138 >    double oldChi = chi;
139 >    double  prevChi;
140 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
141  
142 <  oldChi = chi;
142 >    index = 0;
143 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
144 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
145 >               integrableObject = mol->nextIntegrableObject(j)) {
146 >                oldVel_[index] = integrableObject->getVel();
147 >                oldJi_[index] = integrableObject->getJ();                
148 >        }
149 >        ++index;              
150 >    }
151  
152 <  for( i=0; i < integrableObjects.size(); i++ ){
152 >    // do the iteration:
153  
154 <    integrableObjects[i]->getVel( vel );
155 <
156 <    for (j=0; j < 3; j++)
144 <      oldVel[3*i + j]  = vel[j];
154 >    for(int k = 0; k < maxIterNum_; k++) {
155 >        index = 0;
156 >        instTemp = thermo.getTemperature();
157  
158 <    if( integrableObjects[i]->isDirectional() ){
158 >        // evolve chi another half step using the temperature at t + dt/2
159  
160 <      integrableObjects[i]->getJ( ji );
160 >        prevChi = chi;
161 >        chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
162  
163 <      for (j=0; j < 3; j++)
164 <        oldJi[3*i + j] = ji[j];
163 >        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164 >            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165 >                   integrableObject = mol->nextIntegrableObject(j)) {
166  
167 <    }
168 <  }
167 >                frc = integrableObject->getFrc();
168 >                vel = integrableObject->getVel();
169  
170 <  // do the iteration:
170 >                mass = integrableObject->getMass();
171  
172 <  for (k=0; k < 4; k++) {
172 >                // velocity half step
173 >                //for(j = 0; j < 3; j++)
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  
179 <    instTemp = tStats->getTemperature();
179 >                if (integrableObject->isDirectional()) {
180  
181 <    // evolve chi another half step using the temperature at t + dt/2
181 >                    // get and convert the torque to body frame
182  
183 <    prevChi = chi;
184 <    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
166 <      (tauThermostat*tauThermostat);
183 >                    Tb = integrableObject->getTrq();
184 >                    integrableObject->lab2Body(Tb);
185  
186 <    for( i=0; i < integrableObjects.size(); i++ ){
186 >                    //for(j = 0; j < 3; j++)
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 <      integrableObjects[i]->getFrc( frc );
191 <      integrableObjects[i]->getVel(vel);
190 >                    integrableObject->setJ(ji);
191 >                }
192 >            }
193 >        }
194 >    
195  
196 <      mass = integrableObjects[i]->getMass();
196 >        //constraintAlgorithm->doConstrainB();
197  
198 <      // velocity half step
199 <      for (j=0; j < 3; j++)
177 <        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*chi);
198 >        if (fabs(prevChi - chi) <= chiTolerance_)
199 >            break;
200  
201 <      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 <      }
201 >        ++index;
202      }
194 –    
195 –    if(nConstrained)
196 –      constrainB();
203  
204 <    if (fabs(prevChi - chi) <= chiTolerance) break;
199 <  }
204 >    integralOfChidt += dt2 * chi;
205  
206 <  integralOfChidt += dt2*chi;
206 >    currSnapshot->setChi(chi);
207 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);
208   }
209  
204 – template<typename T> void NVT<T>::resetIntegrator( void ){
210  
211 <  chi = 0.0;
212 <  integralOfChidt = 0.0;
213 < }
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 < template<typename T> int NVT<T>::readyCheck() {
223 >    Energy = thermo.getTotalE();
224  
225 <  //check parent's readyCheck() first
213 <  if (T::readyCheck() == -1)
214 <    return -1;
225 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
226  
227 <  // First check to see if we have a target temperature.
217 <  // Not having one is fatal.
227 >    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
228  
229 <  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 <  }
229 >    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
230  
231 <  // We must set tauThermostat.
230 <
231 <  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 <
231 >    return conservedQuantity;
232   }
233  
255 – template<typename T> double NVT<T>::getConservedQuantity(void){
234  
235 <  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 < }
235 > }//end namespace oopse

Diff Legend

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