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 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 769 by tim, Fri Sep 19 14:22:29 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 +    if (nConstrained){
179 +      constrainA();
180 +    }
181 +  }
182 +    
183 +
184    // Scale the box after all the positions have been moved:
185    
186    scaleFactor = exp(dt*eta);
# Line 139 | Line 194 | template<typename T> void NPTi<T>::moveA() {
194      painCave.isFatal = 1;
195      simError();
196    } else {        
197 <    info->scaleBox(exp(dt*eta));      
198 <  }
197 >    info->scaleBox(scaleFactor);      
198 >  }  
199  
200   }
201  
202   template<typename T> void NPTi<T>::moveB( void ){
203 <
204 <  int i, j;
203 >  
204 >  //new version of NPTi
205 >  int i, j, k;
206    DirectionalAtom* dAtom;
207    double Tb[3], ji[3];
208    double vel[3], frc[3];
209    double mass;
210  
211 <  double instaTemp, instaPress, instaVol;
211 >  double instTemp, instPress, instVol;
212    double tt2, tb2;
213 +  double oldChi, prevChi;
214 +  double oldEta, preEta;
215    
216    tt2 = tauThermostat * tauThermostat;
217    tb2 = tauBarostat * tauBarostat;
218  
161 –  instaTemp = tStats->getTemperature();
162 –  instaPress = tStats->getPressure();
163 –  instaVol = tStats->getVolume();
219  
220 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
221 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
222 <                 (p_convert*NkBT*tb2));
223 <  
220 >  // Set things up for the iteration:
221 >
222 >  oldChi = chi;
223 >  oldEta = eta;
224 >
225    for( i=0; i<nAtoms; i++ ){
226  
227      atoms[i]->getVel( vel );
172 –    atoms[i]->getFrc( frc );
228  
229 <    mass = atoms[i]->getMass();
229 >    for (j=0; j < 3; j++)
230 >      oldVel[3*i + j]  = vel[j];
231  
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 –
232      if( atoms[i]->isDirectional() ){
233  
234        dAtom = (DirectionalAtom *)atoms[i];
235  
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 –
236        dAtom->getJ( ji );
237  
238 <      for (j=0; j < 3; j++)
239 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
238 >      for (j=0; j < 3; j++)
239 >        oldJi[3*i + j] = ji[j];
240  
198 –      dAtom->setJ( ji );
241      }
242    }
243 +
244 +  // do the iteration:
245 +
246 +  instVol = tStats->getVolume();
247 +  
248 +  for (k=0; k < 4; k++) {
249 +    
250 +    instTemp = tStats->getTemperature();
251 +    instPress = tStats->getPressure();
252 +
253 +    // evolve chi another half step using the temperature at t + dt/2
254 +
255 +    prevChi = chi;
256 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
257 +      (tauThermostat*tauThermostat);
258 +
259 +    preEta = eta;
260 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
261 +       (p_convert*NkBT*tb2));
262 +
263 +  
264 +    for( i=0; i<nAtoms; i++ ){
265 +
266 +      atoms[i]->getFrc( frc );
267 +      atoms[i]->getVel(vel);
268 +      
269 +      mass = atoms[i]->getMass();
270 +      
271 +      // velocity half step
272 +      for (j=0; j < 3; j++)
273 +        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
274 +      
275 +      atoms[i]->setVel( vel );
276 +      
277 +      if( atoms[i]->isDirectional() ){
278 +
279 +        dAtom = (DirectionalAtom *)atoms[i];
280 +  
281 +        // get and convert the torque to body frame      
282 +  
283 +        dAtom->getTrq( Tb );
284 +        dAtom->lab2Body( Tb );      
285 +            
286 +        for (j=0; j < 3; j++)
287 +          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
288 +      
289 +          dAtom->setJ( ji );
290 +      }
291 +    }
292 +    
293 +    if (nConstrained){
294 +      constrainB();
295 +    }    
296 +    
297 +    if (fabs(prevChi - chi) <=
298 +        chiTolerance && fabs(preEta -eta) <= etaTolerance)
299 +      break;
300 +  }
301 +
302 +  //calculate integral of chida
303 +  integralOfChidt += dt2*chi;
304 +
305 +
306   }
307  
308 + template<typename T> void NPTi<T>::resetIntegrator() {
309 +  chi = 0.0;
310 +  eta = 0.0;
311 + }
312 +
313   template<typename T> int NPTi<T>::readyCheck() {
314  
315    //check parent's readyCheck() first
# Line 251 | Line 361 | template<typename T> int NPTi<T>::readyCheck() {
361      return -1;
362    }    
363  
364 +  if (!have_chi_tolerance) {
365 +    sprintf( painCave.errMsg,
366 +             "NPTi warning: setting chi tolerance to 1e-6\n");
367 +    chiTolerance = 1e-6;
368 +    have_chi_tolerance = 1;
369 +    painCave.isFatal = 0;
370 +    simError();
371 +  }
372 +
373 +    if (!have_eta_tolerance) {
374 +    sprintf( painCave.errMsg,
375 +             "NPTi warning: setting eta tolerance to 1e-6\n");
376 +    etaTolerance = 1e-6;
377 +    have_eta_tolerance = 1;
378 +    painCave.isFatal = 0;
379 +    simError();
380 +  }
381    // We need NkBT a lot, so just set it here:
382  
383 <  NkBT = (double)info->ndf * kB * targetTemp;
383 >  NkBT = (double)Nparticles * kB * targetTemp;
384 >  fkBT = (double)info->ndf * kB * targetTemp;
385  
386    return 1;
387   }
388 +
389 + template<typename T> double NPTi<T>::getConservedQuantity(void){
390 +
391 +  double conservedQuantity;
392 +  double LkBT;
393 +  double fkBT;
394 +  double f1kBT;
395 +  double f2kBT;
396 +  double NkBT;
397 +  double Energy;
398 +  double thermostat_kinetic;
399 +  double thermostat_potential;
400 +  double barostat_kinetic;
401 +  double barostat_potential;
402 +  double tb2;
403 +  double eta2;  
404 +
405 +  LkBT = (double)(info->getNDF() + 4) * kB * targetTemp;  // 3N + 1
406 +  fkBT = (double)(info->getNDF()    ) * kB * targetTemp;  // 3N - 3
407 +  f1kBT = (double)(info->getNDF()+ 1) * kB * targetTemp;  // 3N - 3 + 1
408 +  NkBT = (double)(info->getNDF() + 3) * kB * targetTemp;  // 3N
409 +  f2kBT = (double)(info->getNDF()+ 2) * kB * targetTemp;  // 3N - 3 + 1
410 +
411 +  Energy = tStats->getTotalE();
412 +
413 +  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
414 +    (2.0 * eConvert);
415 +
416 +  thermostat_potential = fkBT* integralOfChidt / eConvert;
417 +
418 +
419 +  barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta /
420 +    (2.0 * eConvert);
421 +  
422 +  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
423 +    eConvert;
424 +
425 +  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
426 +    barostat_kinetic + barostat_potential;
427 +  
428 +  cout.width(8);
429 +  cout.precision(8);
430 +
431 +  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
432 +      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
433 +      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
434 +
435 +  return conservedQuantity;
436 + }

Diff Legend

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