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 746 by mmeineke, Thu Sep 4 21:48:35 2003 UTC vs.
Revision 763 by tim, Mon Sep 15 16:52:02 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 +
61 + //   int i, j;
62 + //   DirectionalAtom* dAtom;
63 + //   double Tb[3], ji[3];
64 + //   double A[3][3], I[3][3];
65 + //   double angle, mass;
66 + //   double vel[3], pos[3], frc[3];
67 +
68 + //   double rj[3];
69 + //   double instaTemp, instaPress, instaVol;
70 + //   double tt2, tb2, scaleFactor;
71 +
72 + //   tt2 = tauThermostat * tauThermostat;
73 + //   tb2 = tauBarostat * tauBarostat;
74 +
75 + //   instaTemp = tStats->getTemperature();
76 + //   instaPress = tStats->getPressure();
77 + //   instaVol = tStats->getVolume();
78 +  
79 + //    // first evolve chi a half step
80    
81 <  int i, j;
81 > //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
82 > //   eta += dt2 * ( instaVol * (instaPress - targetPressure) /
83 > //               (p_convert*NkBT*tb2));
84 >
85 > //   integralOfChidt += dt2* chi;
86 >
87 > //   for( i=0; i<nAtoms; i++ ){
88 > //     atoms[i]->getVel( vel );
89 > //     atoms[i]->getPos( pos );
90 > //     atoms[i]->getFrc( frc );
91 >
92 > //     mass = atoms[i]->getMass();
93 >
94 > //     for (j=0; j < 3; j++) {
95 > //       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
96 > //       rj[j] = pos[j];
97 > //     }
98 >
99 > //     atoms[i]->setVel( vel );
100 >
101 > //     info->wrapVector(rj);
102 >
103 > //     for (j = 0; j < 3; j++)
104 > //       pos[j] += dt * (vel[j] + eta*rj[j]);
105 >
106 > //     atoms[i]->setPos( pos );
107 >
108 > //     if( atoms[i]->isDirectional() ){
109 >
110 > //       dAtom = (DirectionalAtom *)atoms[i];
111 >          
112 > //       // get and convert the torque to body frame
113 >      
114 > //       dAtom->getTrq( Tb );
115 > //       dAtom->lab2Body( Tb );
116 >      
117 > //       // get the angular momentum, and propagate a half step
118 >
119 > //       dAtom->getJ( ji );
120 >
121 > //       for (j=0; j < 3; j++)
122 > //         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
123 >      
124 > //       // use the angular velocities to propagate the rotation matrix a
125 > //       // full time step
126 >
127 > //       dAtom->getA(A);
128 > //       dAtom->getI(I);
129 >    
130 > //       // rotate about the x-axis      
131 > //       angle = dt2 * ji[0] / I[0][0];
132 > //       this->rotate( 1, 2, angle, ji, A );
133 >
134 > //       // rotate about the y-axis
135 > //       angle = dt2 * ji[1] / I[1][1];
136 > //       this->rotate( 2, 0, angle, ji, A );
137 >      
138 > //       // rotate about the z-axis
139 > //       angle = dt * ji[2] / I[2][2];
140 > //       this->rotate( 0, 1, angle, ji, A);
141 >      
142 > //       // rotate about the y-axis
143 > //       angle = dt2 * ji[1] / I[1][1];
144 > //       this->rotate( 2, 0, angle, ji, A );
145 >      
146 > //        // rotate about the x-axis
147 > //       angle = dt2 * ji[0] / I[0][0];
148 > //       this->rotate( 1, 2, angle, ji, A );
149 >      
150 > //       dAtom->setJ( ji );
151 > //       dAtom->setA( A  );    
152 > //     }                
153 >
154 > //   }
155 >
156 > //   // Scale the box after all the positions have been moved:
157 >  
158 > //   scaleFactor = exp(dt*eta);
159 >
160 > //   if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
161 > //     sprintf( painCave.errMsg,
162 > //              "NPTi error: Attempting a Box scaling of more than 10 percent"
163 > //              " check your tauBarostat, as it is probably too small!\n"
164 > //              " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
165 > //              );
166 > //     painCave.isFatal = 1;
167 > //     simError();
168 > //   } else {        
169 > //     info->scaleBox(exp(dt*eta));      
170 > //   }
171 >  
172 >
173 >  //new version of NPTi
174 >  int i, j, k;
175    DirectionalAtom* dAtom;
176    double Tb[3], ji[3];
177    double A[3][3], I[3][3];
# Line 43 | Line 181 | template<typename T> void NPTi<T>::moveA() {
181    double rj[3];
182    double instaTemp, instaPress, instaVol;
183    double tt2, tb2, scaleFactor;
184 +  double COM[3];
185  
186    tt2 = tauThermostat * tauThermostat;
187    tb2 = tauBarostat * tauBarostat;
# Line 50 | Line 189 | template<typename T> void NPTi<T>::moveA() {
189    instaTemp = tStats->getTemperature();
190    instaPress = tStats->getPressure();
191    instaVol = tStats->getVolume();
53 –  
54 –   // first evolve chi a half step
192    
193 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
194 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
195 <                 (p_convert*NkBT*tb2));
59 <
193 >  tStats->getCOM(COM);
194 >  
195 >  //evolve velocity half step
196    for( i=0; i<nAtoms; i++ ){
197 +
198      atoms[i]->getVel( vel );
62 –    atoms[i]->getPos( pos );
199      atoms[i]->getFrc( frc );
200  
201      mass = atoms[i]->getMass();
202  
203      for (j=0; j < 3; j++) {
204 <      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
205 <      rj[j] = pos[j];
204 >      // velocity half step  (use chi from previous step here):
205 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
206 >  
207      }
208  
209      atoms[i]->setVel( vel );
210 <
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 <
210 >  
211      if( atoms[i]->isDirectional() ){
212  
213        dAtom = (DirectionalAtom *)atoms[i];
214 <          
214 >
215        // get and convert the torque to body frame
216        
217        dAtom->getTrq( Tb );
# Line 122 | Line 252 | template<typename T> void NPTi<T>::moveA() {
252        
253        dAtom->setJ( ji );
254        dAtom->setA( A  );    
255 <    }                
255 >    }    
256 >  }
257  
258 +  // evolve chi and eta  half step
259 +  
260 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
261 +  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));
262 +
263 +  //calculate the integral of chidt
264 +  integralOfChidt += dt2*chi;
265 +
266 +  //save the old positions
267 +  for(i = 0; i < nAtoms; i++){
268 +    atoms[i]->getPos(pos);
269 +    for(j = 0; j < 3; j++)
270 +      oldPos[i*3 + j] = pos[j];
271    }
272 +  
273 +  //the first estimation of r(t+dt) is equal to  r(t)
274 +    
275 +  for(k = 0; k < 4; k ++){
276  
277 +    for(i =0 ; i < nAtoms; i++){
278 +
279 +      atoms[i]->getVel(vel);
280 +      atoms[i]->getPos(pos);
281 +
282 +      for(j = 0; j < 3; j++)
283 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
284 +
285 +      
286 +      //wrapVector(r(t)) = r(t)-R0
287 +      //info->wrapVector(rj);
288 +      
289 +      for(j = 0; j < 3; j++)
290 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
291 +
292 +      atoms[i]->setPos( pos );
293 +
294 +    }
295 +
296 +  }
297 +    
298 +
299    // Scale the box after all the positions have been moved:
300    
301    scaleFactor = exp(dt*eta);
# Line 139 | Line 309 | template<typename T> void NPTi<T>::moveA() {
309      painCave.isFatal = 1;
310      simError();
311    } else {        
312 <    info->scaleBox(exp(dt*eta));      
313 <  }
312 >    info->scaleBox(scaleFactor);      
313 >  }  
314  
315 +  //advance volume;
316 +  volume = volume * exp(dt*eta);
317   }
318  
319   template<typename T> void NPTi<T>::moveB( void ){
320  
321 + /*
322    int i, j;
323    DirectionalAtom* dAtom;
324    double Tb[3], ji[3];
# Line 165 | Line 338 | template<typename T> void NPTi<T>::moveB( void ){
338    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
339    eta += dt2 * ( instaVol * (instaPress - targetPressure) /
340                   (p_convert*NkBT*tb2));
341 +  integralOfChidt += dt2*chi;
342    
343    for( i=0; i<nAtoms; i++ ){
344  
# Line 198 | Line 372 | template<typename T> void NPTi<T>::moveB( void ){
372        dAtom->setJ( ji );
373      }
374    }
375 +
376 + */
377 +  
378 +  //new version of NPTi
379 +  int i, j, k;
380 +  DirectionalAtom* dAtom;
381 +  double Tb[3], ji[3];
382 +  double vel[3], frc[3];
383 +  double mass;
384 +
385 +  double instTemp, instPress, instVol;
386 +  double tt2, tb2;
387 +  double oldChi, prevChi;
388 +  double oldEta, preEta;
389 +  
390 +  tt2 = tauThermostat * tauThermostat;
391 +  tb2 = tauBarostat * tauBarostat;
392 +
393 +
394 +  // Set things up for the iteration:
395 +
396 +  oldChi = chi;
397 +  oldEta = eta;
398 +
399 +  for( i=0; i<nAtoms; i++ ){
400 +
401 +    atoms[i]->getVel( vel );
402 +
403 +    for (j=0; j < 3; j++)
404 +      oldVel[3*i + j]  = vel[j];
405 +
406 +    if( atoms[i]->isDirectional() ){
407 +
408 +      dAtom = (DirectionalAtom *)atoms[i];
409 +
410 +      dAtom->getJ( ji );
411 +
412 +      for (j=0; j < 3; j++)
413 +        oldJi[3*i + j] = ji[j];
414 +
415 +    }
416 +  }
417 +
418 +  // do the iteration:
419 +
420 +  instVol = tStats->getVolume();
421 +  
422 +  for (k=0; k < 4; k++) {
423 +    
424 +    instTemp = tStats->getTemperature();
425 +    instPress = tStats->getPressure();
426 +
427 +    // evolve chi another half step using the temperature at t + dt/2
428 +
429 +    prevChi = chi;
430 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
431 +      (tauThermostat*tauThermostat);
432 +
433 +    preEta = eta;
434 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
435 +       (p_convert*NkBT*tb2));
436 +
437 +  
438 +    for( i=0; i<nAtoms; i++ ){
439 +
440 +      atoms[i]->getFrc( frc );
441 +      atoms[i]->getVel(vel);
442 +      
443 +      mass = atoms[i]->getMass();
444 +      
445 +      // velocity half step
446 +      for (j=0; j < 3; j++)
447 +        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
448 +      
449 +      atoms[i]->setVel( vel );
450 +      
451 +      if( atoms[i]->isDirectional() ){
452 +
453 +        dAtom = (DirectionalAtom *)atoms[i];
454 +  
455 +        // get and convert the torque to body frame      
456 +  
457 +        dAtom->getTrq( Tb );
458 +        dAtom->lab2Body( Tb );      
459 +            
460 +        for (j=0; j < 3; j++)
461 +          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
462 +      
463 +          dAtom->setJ( ji );
464 +      }
465 +    }
466 +
467 +    if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance)
468 +      break;
469 +  }
470 +
471 +  //calculate integral of chida
472 +  integralOfChidt += dt2*chi;
473 +
474 +
475   }
476  
477   template<typename T> void NPTi<T>::resetIntegrator() {
# Line 256 | Line 530 | template<typename T> int NPTi<T>::readyCheck() {
530      return -1;
531    }    
532  
533 +  if (!have_chi_tolerance) {
534 +    sprintf( painCave.errMsg,
535 +             "NPTi warning: setting chi tolerance to 1e-6\n");
536 +    chiTolerance = 1e-6;
537 +    have_chi_tolerance = 1;
538 +    painCave.isFatal = 0;
539 +    simError();
540 +  }
541 +
542 +    if (!have_eta_tolerance) {
543 +    sprintf( painCave.errMsg,
544 +             "NPTi warning: setting eta tolerance to 1e-6\n");
545 +    etaTolerance = 1e-6;
546 +    have_eta_tolerance = 1;
547 +    painCave.isFatal = 0;
548 +    simError();
549 +  }
550    // We need NkBT a lot, so just set it here:
551  
552 <  NkBT = (double)info->ndf * kB * targetTemp;
552 >  NkBT = (double)Nparticles * kB * targetTemp;
553 >  fkBT = (double)info->ndf * kB * targetTemp;
554  
555    return 1;
556   }
557 +
558 + template<typename T> double NPTi<T>::getConservedQuantity(void){
559 +
560 +  double conservedQuantity;
561 +  double tb2;
562 +  double eta2;  
563 +  double E_NPT;
564 +  double U;
565 +  double TS;
566 +  double PV;
567 +  double extra;
568 +
569 +  static double pre_U;
570 +  static double pre_TS;
571 +  static double pre_PV;
572 +  static double pre_extra;
573 +  static int hackCount = 0;
574 +
575 +  double delta_U;
576 +  double delta_TS;
577 +  double delta_PV;
578 +  double delta_extra;
579 +
580 +  U = tStats->getTotalE();
581 +
582 +  TS = fkBT *
583 +    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
584 +
585 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
586 +
587 +  tb2 = tauBarostat * tauBarostat;
588 +  eta2 = eta * eta;
589 +
590 +  extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
591 +  /*
592 +  if(hackCount == 0){
593 +    pre_U = U;
594 +    pre_TS =TS;
595 +    pre_PV = PV;
596 +    pre_extra =extra;
597 +    hackCount ++;
598 +  }
599 +
600 +  delta_U = U - pre_U;
601 +  delta_TS = TS - pre_TS;
602 +  delta_PV = PV - pre_PV;
603 +  delta_extra = extra - pre_extra;
604 + */
605 +  cout.width(8);
606 +  cout.precision(8);
607 +
608 +  
609 +  cout << info->getTime() << "\t"
610 +       << chi << "\t"
611 +       << eta << "\t"
612 +       << U << "\t"
613 +       << TS << "\t"
614 +       << PV << "\t"
615 +       << extra << "\t"
616 +       << U+TS+PV+extra << endl;
617 +
618 + /*
619 +    pre_U = U;
620 +    pre_TS =TS;
621 +    pre_PV = PV;
622 +    pre_extra =extra;
623 +
624 +
625 +  cout << info->getTime() << "\t"
626 +       << U << "\t"
627 +       << U+TS << "\t"
628 +       << U+TS+PV << "\t"
629 +       << U+TS+PV+extra << endl;
630 + */
631 +  conservedQuantity = U+TS+PV+extra;
632 +  return conservedQuantity;
633 + }

Diff Legend

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