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

Comparing trunk/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 600 by gezelter, Mon Jul 14 22:38:13 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 19 | Line 20 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
20   //
21   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
22  
23 < NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
24 <  Integrator( theInfo, the_ff )
23 > template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
24 >  T( theInfo, the_ff )
25   {
26    int i, j;
27    chi = 0.0;
28 +  integralOfChidt = 0.0;
29  
30    for(i = 0; i < 3; i++)
31      for (j = 0; j < 3; j++)
# Line 33 | Line 35 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
35    have_tau_barostat = 0;
36    have_target_temp = 0;
37    have_target_pressure = 0;
38 +
39 +  have_chi_tolerance = 0;
40 +  have_eta_tolerance = 0;
41 +  have_pos_iter_tolerance = 0;
42 +
43 +  oldPos = new double[3*nAtoms];
44 +  oldVel = new double[3*nAtoms];
45 +  oldJi = new double[3*nAtoms];
46 + #ifdef IS_MPI
47 +  Nparticles = mpiSim->getTotAtoms();
48 + #else
49 +  Nparticles = theInfo->n_atoms;
50 + #endif
51   }
52  
53 < void NPTf::moveA() {
53 > template<typename T> NPTf<T>::~NPTf() {
54 >  delete[] oldPos;
55 >  delete[] oldVel;
56 >  delete[] oldJi;
57 > }
58 >
59 > template<typename T> void NPTf<T>::moveA() {
60    
61    int i, j, k;
62    DirectionalAtom* dAtom;
# Line 50 | Line 71 | void NPTf::moveA() {
71    double sc[3];
72    double eta2ij;
73    double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
74 +  double bigScale, smallScale, offDiagMax;
75 +  double COM[3];
76  
77    tt2 = tauThermostat * tauThermostat;
78    tb2 = tauBarostat * tauBarostat;
# Line 57 | Line 80 | void NPTf::moveA() {
80    instaTemp = tStats->getTemperature();
81    tStats->getPressureTensor(press);
82    instaVol = tStats->getVolume();
60 –  
61 –  // first evolve chi a half step
83    
84 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
84 >  tStats->getCOM(COM);
85  
86 +  //calculate scale factor of veloity
87    for (i = 0; i < 3; i++ ) {
88      for (j = 0; j < 3; j++ ) {
89 +      vScale[i][j] = eta[i][j];
90 +      
91        if (i == j) {
92 <        
93 <        eta[i][j] += dt2 * instaVol *
70 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
71 <        
72 <        vScale[i][j] = eta[i][j] + chi;
73 <        
74 <      } else {
75 <        
76 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
77 <
78 <        vScale[i][j] = eta[i][j];
79 <        
80 <      }
92 >        vScale[i][j] += chi;          
93 >      }              
94      }
95    }
96 <
96 >  
97 >  //evolve velocity half step
98    for( i=0; i<nAtoms; i++ ){
99  
100      atoms[i]->getVel( vel );
87 –    atoms[i]->getPos( pos );
101      atoms[i]->getFrc( frc );
102  
103      mass = atoms[i]->getMass();
104      
92 –    // velocity half step
93 –        
105      info->matVecMul3( vScale, vel, sc );
106 <    
107 <    for (j = 0; j < 3; j++) {
106 >
107 >    for (j=0; j < 3; j++) {
108 >      // velocity half step  (use chi from previous step here):
109        vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
110 <      rj[j] = pos[j];
110 >  
111      }
112  
113      atoms[i]->setVel( vel );
102 –
103 –    // position whole step    
104 –
105 –    info->wrapVector(rj);
106 –
107 –    info->matVecMul3( eta, rj, sc );
108 –
109 –    for (j = 0; j < 3; j++ )
110 –      pos[j] += dt * (vel[j] + sc[j]);
114    
115      if( atoms[i]->isDirectional() ){
116  
117        dAtom = (DirectionalAtom *)atoms[i];
118 <          
118 >
119        // get and convert the torque to body frame
120        
121        dAtom->getTrq( Tb );
# Line 153 | Line 156 | void NPTf::moveA() {
156        
157        dAtom->setJ( ji );
158        dAtom->setA( A  );    
159 <    }                    
159 >    }    
160    }
161 +
162 +  // advance chi half step
163 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
164 +
165 +  //calculate the integral of chidt
166 +  integralOfChidt += dt2*chi;
167 +
168 +  //advance eta half step
169 +  for(i = 0; i < 3; i ++)
170 +    for(j = 0; j < 3; j++){
171 +      if( i == j)
172 +        eta[i][j] += dt2 *  instaVol *
173 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
174 +      else
175 +        eta[i][j] += dt2 * instaVol * press[i][j] / ( NkBT*tb2);
176 +    }
177 +    
178 +  //save the old positions
179 +  for(i = 0; i < nAtoms; i++){
180 +    atoms[i]->getPos(pos);
181 +    for(j = 0; j < 3; j++)
182 +      oldPos[i*3 + j] = pos[j];
183 +  }
184    
185 +  //the first estimation of r(t+dt) is equal to  r(t)
186 +    
187 +  for(k = 0; k < 4; k ++){
188 +
189 +    for(i =0 ; i < nAtoms; i++){
190 +
191 +      atoms[i]->getVel(vel);
192 +      atoms[i]->getPos(pos);
193 +
194 +      for(j = 0; j < 3; j++)
195 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
196 +      
197 +      info->matVecMul3( eta, rj, sc );
198 +      
199 +      for(j = 0; j < 3; j++)
200 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
201 +
202 +      atoms[i]->setPos( pos );
203 +
204 +    }
205 +
206 +  }  
207 +
208 +
209    // Scale the box after all the positions have been moved:
210    
211    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
212    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
213    
214 +  bigScale = 1.0;
215 +  smallScale = 1.0;
216 +  offDiagMax = 0.0;
217    
218    for(i=0; i<3; i++){
219      for(j=0; j<3; j++){
# Line 178 | Line 231 | void NPTf::moveA() {
231        if (i == j) scaleMat[i][j] = 1.0;
232        // Taylor expansion for the exponential truncated at second order:
233        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
234 <      
234 >
235 >      if (i != j)
236 >        if (fabs(scaleMat[i][j]) > offDiagMax)
237 >          offDiagMax = fabs(scaleMat[i][j]);
238      }
239 +
240 +    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
241 +    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
242    }
243    
244 <  info->getBoxM(hm);
245 <  info->matMul3(hm, scaleMat, hmnew);
246 <  info->setBoxM(hmnew);
244 >  if ((bigScale > 1.1) || (smallScale < 0.9)) {
245 >    sprintf( painCave.errMsg,
246 >             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
247 >             " Check your tauBarostat, as it is probably too small!\n\n"
248 >             " scaleMat = [%lf\t%lf\t%lf]\n"
249 >             "            [%lf\t%lf\t%lf]\n"
250 >             "            [%lf\t%lf\t%lf]\n",
251 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
252 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
253 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
254 >    painCave.isFatal = 1;
255 >    simError();
256 >  } else if (offDiagMax > 0.1) {
257 >    sprintf( painCave.errMsg,
258 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
259 >             " Check your tauBarostat, as it is probably too small!\n\n"
260 >             " scaleMat = [%lf\t%lf\t%lf]\n"
261 >             "            [%lf\t%lf\t%lf]\n"
262 >             "            [%lf\t%lf\t%lf]\n",
263 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
264 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
265 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
266 >    painCave.isFatal = 1;
267 >    simError();
268 >  } else {
269 >    info->getBoxM(hm);
270 >    info->matMul3(hm, scaleMat, hmnew);
271 >    info->setBoxM(hmnew);
272 >  }
273    
274   }
275  
276 < void NPTf::moveB( void ){
276 > template<typename T> void NPTf<T>::moveB( void ){
277  
278 <  int i, j;
278 >  int i, j, k;
279    DirectionalAtom* dAtom;
280    double Tb[3], ji[3];
281    double vel[3], frc[3];
# Line 200 | Line 285 | void NPTf::moveB( void ){
285    double tt2, tb2;
286    double sc[3];
287    double press[3][3], vScale[3][3];
288 +  double oldChi, prevChi;
289 +  double oldEta[3][3], preEta[3][3], diffEta;
290    
291    tt2 = tauThermostat * tauThermostat;
292    tb2 = tauBarostat * tauBarostat;
293  
207 –  instaTemp = tStats->getTemperature();
208 –  tStats->getPressureTensor(press);
209 –  instaVol = tStats->getVolume();
210 –  
211 –  // first evolve chi a half step
212 –  
213 –  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
214 –  
215 –  for (i = 0; i < 3; i++ ) {
216 –    for (j = 0; j < 3; j++ ) {
217 –      if (i == j) {
294  
295 <        eta[i][j] += dt2 * instaVol *
220 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
295 >  // Set things up for the iteration:
296  
297 <        vScale[i][j] = eta[i][j] + chi;
298 <        
299 <      } else {
300 <        
301 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
297 >  oldChi = chi;
298 >  
299 >  for(i = 0; i < 3; i++)
300 >    for(j = 0; j < 3; j++)
301 >      oldEta[i][j] = eta[i][j];
302  
228 –        vScale[i][j] = eta[i][j];
229 –        
230 –      }
231 –    }
232 –  }
233 –
303    for( i=0; i<nAtoms; i++ ){
304  
305      atoms[i]->getVel( vel );
237 –    atoms[i]->getFrc( frc );
306  
307 <    mass = atoms[i]->getMass();
308 <    
241 <    // velocity half step
242 <        
243 <    info->matVecMul3( vScale, vel, sc );
244 <    
245 <    for (j = 0; j < 3; j++) {
246 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
247 <    }
307 >    for (j=0; j < 3; j++)
308 >      oldVel[3*i + j]  = vel[j];
309  
249 –    atoms[i]->setVel( vel );
250 –    
310      if( atoms[i]->isDirectional() ){
311  
312        dAtom = (DirectionalAtom *)atoms[i];
313 <          
314 <      // get and convert the torque to body frame
313 >
314 >      dAtom->getJ( ji );
315 >
316 >      for (j=0; j < 3; j++)
317 >        oldJi[3*i + j] = ji[j];
318 >
319 >    }
320 >  }
321 >
322 >  // do the iteration:
323 >
324 >  instaVol = tStats->getVolume();
325 >  
326 >  for (k=0; k < 4; k++) {
327 >    
328 >    instaTemp = tStats->getTemperature();
329 >    tStats->getPressureTensor(press);
330 >
331 >    // evolve chi another half step using the temperature at t + dt/2
332 >
333 >    prevChi = chi;
334 >    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
335 >    
336 >    for(i = 0; i < 3; i++)
337 >      for(j = 0; j < 3; j++)
338 >        preEta[i][j] = eta[i][j];
339 >
340 >    //advance eta half step and calculate scale factor for velocity
341 >    for(i = 0; i < 3; i ++)
342 >      for(j = 0; j < 3; j++){
343 >        if( i == j){
344 >          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
345 >            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
346 >          vScale[i][j] = eta[i][j] + chi;
347 >        }
348 >        else
349 >        {
350 >          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
351 >          vScale[i][j] = eta[i][j];
352 >        }
353 >    }      
354 >
355 >    //advance velocity half step
356 >    for( i=0; i<nAtoms; i++ ){
357 >
358 >      atoms[i]->getFrc( frc );
359 >      atoms[i]->getVel(vel);
360        
361 <      dAtom->getTrq( Tb );
258 <      dAtom->lab2Body( Tb );
361 >      mass = atoms[i]->getMass();
362        
363 <      // get the angular momentum, and propagate a half step
363 >      info->matVecMul3( vScale, vel, sc );
364 >
365 >      for (j=0; j < 3; j++) {
366 >        // velocity half step  (use chi from previous step here):
367 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
368 >      }
369        
370 <      dAtom->getJ( ji );
370 >      atoms[i]->setVel( vel );
371        
372 <      for (j=0; j < 3; j++)
373 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
372 >      if( atoms[i]->isDirectional() ){
373 >
374 >        dAtom = (DirectionalAtom *)atoms[i];
375 >  
376 >        // get and convert the torque to body frame      
377 >  
378 >        dAtom->getTrq( Tb );
379 >        dAtom->lab2Body( Tb );      
380 >            
381 >        for (j=0; j < 3; j++)
382 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
383        
384 <      dAtom->setJ( ji );
384 >          dAtom->setJ( ji );
385 >      }
386 >    }
387  
388 <    }                    
388 >    
389 >    diffEta = 0;
390 >    for(i = 0; i < 3; i++)
391 >      diffEta += pow(preEta[i][i] - eta[i][i], 2);    
392 >    
393 >    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
394 >      break;
395    }
396 +
397 +  //calculate integral of chida
398 +  integralOfChidt += dt2*chi;
399 +
400 +  
401   }
402  
403 < int NPTf::readyCheck() {
403 > template<typename T> void NPTf<T>::resetIntegrator() {
404 >  int i,j;
405 >  
406 >  chi = 0.0;
407 >
408 >  for(i = 0; i < 3; i++)
409 >    for (j = 0; j < 3; j++)
410 >      eta[i][j] = 0.0;
411 >
412 > }
413 >
414 > template<typename T> int NPTf<T>::readyCheck() {
415 >
416 >  //check parent's readyCheck() first
417 >  if (T::readyCheck() == -1)
418 >    return -1;
419  
420    // First check to see if we have a target temperature.
421    // Not having one is fatal.
# Line 319 | Line 464 | int NPTf::readyCheck() {
464  
465    // We need NkBT a lot, so just set it here:
466  
467 <  NkBT = (double)info->ndf * kB * targetTemp;
467 >  NkBT = (double)Nparticles * kB * targetTemp;
468 >  fkBT = (double)info->ndf * kB * targetTemp;
469  
470    return 1;
471   }
472 +
473 + template<typename T> double NPTf<T>::getConservedQuantity(void){
474 +
475 +  double conservedQuantity;
476 +  double tb2;
477 +  double trEta;  
478 +  double U;
479 +  double thermo;
480 +  double integral;
481 +  double baro;
482 +  double PV;
483 +
484 +  U = tStats->getTotalE();
485 +  thermo = (fkBT * tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
486 +
487 +  tb2 = tauBarostat * tauBarostat;
488 +  trEta = info->matTrace3(eta);
489 +  baro = ((double)info->ndfTrans * kB * targetTemp * tb2 * trEta * trEta / 2.0) / eConvert;
490 +
491 +  integral = ((double)(info->ndf + 1) * kB * targetTemp * integralOfChidt) /eConvert;
492 +
493 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
494 +
495 +
496 +  cout.width(8);
497 +  cout.precision(8);
498 +  
499 +  cout << info->getTime() << "\t"
500 +       << chi << "\t"
501 +       << trEta << "\t"
502 +       << U << "\t"
503 +       << thermo << "\t"
504 +       << baro << "\t"
505 +       << integral << "\t"
506 +       << PV << "\t"
507 +       << U+thermo+integral+PV+baro << endl;
508 +
509 +  conservedQuantity = U+thermo+integral+PV+baro;
510 +  return conservedQuantity;
511 +  
512 + }

Diff Legend

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