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 1765 by tim, Mon Nov 22 20:55:52 2004 UTC

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

Diff Legend

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