ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPTi.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 645 by tim, Tue Jul 22 19:54:52 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 9 | Line 9
9   #include "Integrator.hpp"
10   #include "simError.h"
11  
12 + #ifdef IS_MPI
13 + #include "mpiSimulation.hpp"
14 + #endif
15  
16 +
17   // Basic isotropic thermostating and barostating via the Melchionna
18   // modification of the Hoover algorithm:
19   //
# Line 25 | Line 29 | template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo,
29   {
30    chi = 0.0;
31    eta = 0.0;
32 +  integralOfChidt = 0.0;
33    have_tau_thermostat = 0;
34    have_tau_barostat = 0;
35    have_target_temp = 0;
36    have_target_pressure = 0;
37 +  have_chi_tolerance = 0;
38 +  have_eta_tolerance = 0;
39 +  have_pos_iter_tolerance = 0;
40 +
41 +  oldPos = new double[3*nAtoms];
42 +  oldVel = new double[3*nAtoms];
43 +  oldJi = new double[3*nAtoms];
44 + #ifdef IS_MPI
45 +  Nparticles = mpiSim->getTotAtoms();
46 + #else
47 +  Nparticles = theInfo->n_atoms;
48 + #endif
49 +
50   }
51  
52 + template<typename T> NPTi<T>::~NPTi() {
53 +  delete[] oldPos;
54 +  delete[] oldVel;
55 +  delete[] oldJi;
56 + }
57 +
58   template<typename T> void NPTi<T>::moveA() {
59 <  
60 <  int i, j;
59 >
60 >  //new version of NPTi
61 >  int i, j, k;
62    DirectionalAtom* dAtom;
63    double Tb[3], ji[3];
64    double A[3][3], I[3][3];
# Line 43 | Line 68 | template<typename T> void NPTi<T>::moveA() {
68    double rj[3];
69    double instaTemp, instaPress, instaVol;
70    double tt2, tb2, scaleFactor;
71 +  double COM[3];
72  
73    tt2 = tauThermostat * tauThermostat;
74    tb2 = tauBarostat * tauBarostat;
# Line 50 | Line 76 | template<typename T> void NPTi<T>::moveA() {
76    instaTemp = tStats->getTemperature();
77    instaPress = tStats->getPressure();
78    instaVol = tStats->getVolume();
53 –  
54 –   // first evolve chi a half step
79    
80 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
81 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
82 <                 (p_convert*NkBT*tb2));
59 <
80 >  tStats->getCOM(COM);
81 >  
82 >  //evolve velocity half step
83    for( i=0; i<nAtoms; i++ ){
84 +
85      atoms[i]->getVel( vel );
62 –    atoms[i]->getPos( pos );
86      atoms[i]->getFrc( frc );
87  
88      mass = atoms[i]->getMass();
89  
90      for (j=0; j < 3; j++) {
91 <      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
92 <      rj[j] = pos[j];
91 >      // velocity half step  (use chi from previous step here):
92 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
93 >  
94      }
95  
96      atoms[i]->setVel( vel );
97 <
74 <    info->wrapVector(rj);
75 <
76 <    for (j = 0; j < 3; j++)
77 <      pos[j] += dt * (vel[j] + eta*rj[j]);
78 <
79 <    atoms[i]->setPos( pos );
80 <
97 >  
98      if( atoms[i]->isDirectional() ){
99  
100        dAtom = (DirectionalAtom *)atoms[i];
101 <          
101 >
102        // get and convert the torque to body frame
103        
104        dAtom->getTrq( Tb );
# Line 122 | Line 139 | template<typename T> void NPTi<T>::moveA() {
139        
140        dAtom->setJ( ji );
141        dAtom->setA( A  );    
142 <    }                
142 >    }    
143 >  }
144  
145 +  // evolve chi and eta  half step
146 +  
147 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
148 +  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));
149 +
150 +  //calculate the integral of chidt
151 +  integralOfChidt += dt2*chi;
152 +
153 +  //save the old positions
154 +  for(i = 0; i < nAtoms; i++){
155 +    atoms[i]->getPos(pos);
156 +    for(j = 0; j < 3; j++)
157 +      oldPos[i*3 + j] = pos[j];
158    }
159 +  
160 +  //the first estimation of r(t+dt) is equal to  r(t)
161 +    
162 +  for(k = 0; k < 4; k ++){
163  
164 +    for(i =0 ; i < nAtoms; i++){
165 +
166 +      atoms[i]->getVel(vel);
167 +      atoms[i]->getPos(pos);
168 +
169 +      for(j = 0; j < 3; j++)
170 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];    
171 +      
172 +      for(j = 0; j < 3; j++)
173 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
174 +
175 +      atoms[i]->setPos( pos );
176 +
177 +    }
178 +
179 +  }
180 +    
181 +
182    // Scale the box after all the positions have been moved:
183    
184    scaleFactor = exp(dt*eta);
# Line 139 | Line 192 | template<typename T> void NPTi<T>::moveA() {
192      painCave.isFatal = 1;
193      simError();
194    } else {        
195 <    info->scaleBox(exp(dt*eta));      
196 <  }
195 >    info->scaleBox(scaleFactor);      
196 >  }  
197  
198   }
199  
200   template<typename T> void NPTi<T>::moveB( void ){
201 <
202 <  int i, j;
201 >  
202 >  //new version of NPTi
203 >  int i, j, k;
204    DirectionalAtom* dAtom;
205    double Tb[3], ji[3];
206    double vel[3], frc[3];
207    double mass;
208  
209 <  double instaTemp, instaPress, instaVol;
209 >  double instTemp, instPress, instVol;
210    double tt2, tb2;
211 +  double oldChi, prevChi;
212 +  double oldEta, preEta;
213    
214    tt2 = tauThermostat * tauThermostat;
215    tb2 = tauBarostat * tauBarostat;
216  
161 –  instaTemp = tStats->getTemperature();
162 –  instaPress = tStats->getPressure();
163 –  instaVol = tStats->getVolume();
217  
218 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
219 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
220 <                 (p_convert*NkBT*tb2));
221 <  
218 >  // Set things up for the iteration:
219 >
220 >  oldChi = chi;
221 >  oldEta = eta;
222 >
223    for( i=0; i<nAtoms; i++ ){
224  
225      atoms[i]->getVel( vel );
172 –    atoms[i]->getFrc( frc );
226  
227 <    mass = atoms[i]->getMass();
227 >    for (j=0; j < 3; j++)
228 >      oldVel[3*i + j]  = vel[j];
229  
176 –    // velocity half step
177 –    for (j=0; j < 3; j++)
178 –      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
179 –    
180 –    atoms[i]->setVel( vel );
181 –
230      if( atoms[i]->isDirectional() ){
231  
232        dAtom = (DirectionalAtom *)atoms[i];
233  
186 –      // get and convert the torque to body frame      
187 –
188 –      dAtom->getTrq( Tb );
189 –      dAtom->lab2Body( Tb );
190 –
191 –      // get the angular momentum, and propagate a half step
192 –
234        dAtom->getJ( ji );
235  
236 <      for (j=0; j < 3; j++)
237 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
236 >      for (j=0; j < 3; j++)
237 >        oldJi[3*i + j] = ji[j];
238  
198 –      dAtom->setJ( ji );
239      }
240    }
241 +
242 +  // do the iteration:
243 +
244 +  instVol = tStats->getVolume();
245 +  
246 +  for (k=0; k < 4; k++) {
247 +    
248 +    instTemp = tStats->getTemperature();
249 +    instPress = tStats->getPressure();
250 +
251 +    // evolve chi another half step using the temperature at t + dt/2
252 +
253 +    prevChi = chi;
254 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
255 +      (tauThermostat*tauThermostat);
256 +
257 +    preEta = eta;
258 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
259 +       (p_convert*NkBT*tb2));
260 +
261 +  
262 +    for( i=0; i<nAtoms; i++ ){
263 +
264 +      atoms[i]->getFrc( frc );
265 +      atoms[i]->getVel(vel);
266 +      
267 +      mass = atoms[i]->getMass();
268 +      
269 +      // velocity half step
270 +      for (j=0; j < 3; j++)
271 +        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
272 +      
273 +      atoms[i]->setVel( vel );
274 +      
275 +      if( atoms[i]->isDirectional() ){
276 +
277 +        dAtom = (DirectionalAtom *)atoms[i];
278 +  
279 +        // get and convert the torque to body frame      
280 +  
281 +        dAtom->getTrq( Tb );
282 +        dAtom->lab2Body( Tb );      
283 +            
284 +        for (j=0; j < 3; j++)
285 +          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
286 +      
287 +          dAtom->setJ( ji );
288 +      }
289 +    }
290 +
291 +    if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance)
292 +      break;
293 +  }
294 +
295 +  //calculate integral of chida
296 +  integralOfChidt += dt2*chi;
297 +
298 +
299   }
300  
301 + template<typename T> void NPTi<T>::resetIntegrator() {
302 +  chi = 0.0;
303 +  eta = 0.0;
304 + }
305 +
306   template<typename T> int NPTi<T>::readyCheck() {
307 +
308 +  //check parent's readyCheck() first
309 +  if (T::readyCheck() == -1)
310 +    return -1;
311  
312    // First check to see if we have a target temperature.
313    // Not having one is fatal.
# Line 247 | Line 354 | template<typename T> int NPTi<T>::readyCheck() {
354      return -1;
355    }    
356  
357 +  if (!have_chi_tolerance) {
358 +    sprintf( painCave.errMsg,
359 +             "NPTi warning: setting chi tolerance to 1e-6\n");
360 +    chiTolerance = 1e-6;
361 +    have_chi_tolerance = 1;
362 +    painCave.isFatal = 0;
363 +    simError();
364 +  }
365 +
366 +    if (!have_eta_tolerance) {
367 +    sprintf( painCave.errMsg,
368 +             "NPTi warning: setting eta tolerance to 1e-6\n");
369 +    etaTolerance = 1e-6;
370 +    have_eta_tolerance = 1;
371 +    painCave.isFatal = 0;
372 +    simError();
373 +  }
374    // We need NkBT a lot, so just set it here:
375  
376 <  NkBT = (double)info->ndf * kB * targetTemp;
376 >  NkBT = (double)Nparticles * kB * targetTemp;
377 >  fkBT = (double)info->ndf * kB * targetTemp;
378  
379    return 1;
380   }
381 +
382 + template<typename T> double NPTi<T>::getConservedQuantity(void){
383 +
384 +  double conservedQuantity;
385 +  double tb2;
386 +  double eta2;  
387 +  double E_NPT;
388 +  double U;
389 +  double TS;
390 +  double PV;
391 +  double extra;
392 +
393 +  U = tStats->getTotalE();
394 +
395 +  TS = fkBT *
396 +    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
397 +
398 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
399 +
400 +  tb2 = tauBarostat * tauBarostat;
401 +  eta2 = eta * eta;
402 +
403 +
404 +  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;
405 +
406 +  cout.width(8);
407 +  cout.precision(8);
408 +
409 +  
410 +  cout << info->getTime() << "\t"
411 +       << chi << "\t"
412 +       << eta << "\t"
413 +       << U << "\t"
414 +       << TS << "\t"
415 +       << PV << "\t"
416 +       << extra << "\t"
417 +       << U+TS+PV+extra << endl;
418 +
419 +  conservedQuantity = U+TS+PV+extra;
420 +  return conservedQuantity;
421 + }

Diff Legend

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