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 772 by gezelter, Fri Sep 19 16:01:07 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 8 | Line 9
9   #include "Integrator.hpp"
10   #include "simError.h"
11  
12 + #ifdef IS_MPI
13 + #include "mpiSimulation.hpp"
14 + #endif
15  
16   // Basic non-isotropic thermostating and barostating via the Melchionna
17   // modification of the Hoover algorithm:
# Line 19 | Line 23 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
23   //
24   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
25  
26 < NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
27 <  Integrator( theInfo, the_ff )
26 > template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
27 >  T( theInfo, the_ff )
28   {
29    int i, j;
30    chi = 0.0;
31 +  integralOfChidt = 0.0;
32  
33    for(i = 0; i < 3; i++)
34      for (j = 0; j < 3; j++)
# Line 33 | Line 38 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
38    have_tau_barostat = 0;
39    have_target_temp = 0;
40    have_target_pressure = 0;
41 +
42 +  have_chi_tolerance = 0;
43 +  have_eta_tolerance = 0;
44 +  have_pos_iter_tolerance = 0;
45 +
46 +  oldPos = new double[3*nAtoms];
47 +  oldVel = new double[3*nAtoms];
48 +  oldJi = new double[3*nAtoms];
49 + #ifdef IS_MPI
50 +  Nparticles = mpiSim->getTotAtoms();
51 + #else
52 +  Nparticles = theInfo->n_atoms;
53 + #endif
54 +
55   }
56  
57 < void NPTf::moveA() {
58 <  
57 > template<typename T> NPTf<T>::~NPTf() {
58 >  delete[] oldPos;
59 >  delete[] oldVel;
60 >  delete[] oldJi;
61 > }
62 >
63 > template<typename T> void NPTf<T>::moveA() {
64 >
65 >  // new version of NPTf
66    int i, j, k;
67    DirectionalAtom* dAtom;
68    double Tb[3], ji[3];
# Line 50 | Line 76 | void NPTf::moveA() {
76    double sc[3];
77    double eta2ij;
78    double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
79 +  double bigScale, smallScale, offDiagMax;
80 +  double COM[3];
81  
82    tt2 = tauThermostat * tauThermostat;
83    tb2 = tauBarostat * tauBarostat;
# Line 57 | Line 85 | void NPTf::moveA() {
85    instaTemp = tStats->getTemperature();
86    tStats->getPressureTensor(press);
87    instaVol = tStats->getVolume();
60 –  
61 –  // first evolve chi a half step
88    
89 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
89 >  tStats->getCOM(COM);
90  
91 +  //calculate scale factor of veloity
92    for (i = 0; i < 3; i++ ) {
93      for (j = 0; j < 3; j++ ) {
94 +      vScale[i][j] = eta[i][j];
95 +      
96        if (i == j) {
97 <        
98 <        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 <      }
97 >        vScale[i][j] += chi;          
98 >      }              
99      }
100    }
101 <
101 >  
102 >  //evolve velocity half step
103    for( i=0; i<nAtoms; i++ ){
104  
105      atoms[i]->getVel( vel );
87 –    atoms[i]->getPos( pos );
106      atoms[i]->getFrc( frc );
107  
108      mass = atoms[i]->getMass();
109      
92 –    // velocity half step
93 –        
110      info->matVecMul3( vScale, vel, sc );
111 <    
112 <    for (j = 0; j < 3; j++) {
111 >
112 >    for (j=0; j < 3; j++) {
113 >      // velocity half step
114        vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
98 –      rj[j] = pos[j];
115      }
116  
117      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]);
118    
119      if( atoms[i]->isDirectional() ){
120  
121        dAtom = (DirectionalAtom *)atoms[i];
122 <          
122 >
123        // get and convert the torque to body frame
124        
125        dAtom->getTrq( Tb );
# Line 153 | Line 160 | void NPTf::moveA() {
160        
161        dAtom->setJ( ji );
162        dAtom->setA( A  );    
163 <    }                    
163 >    }    
164    }
165 +
166 +  // advance chi half step
167 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
168 +
169 +  // calculate the integral of chidt
170 +  integralOfChidt += dt2*chi;
171 +
172 +  // advance eta half step
173 +
174 +  for(i = 0; i < 3; i ++)
175 +    for(j = 0; j < 3; j++){
176 +      if( i == j)
177 +        eta[i][j] += dt2 *  instaVol *
178 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
179 +      else
180 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
181 +    }
182 +    
183 +  //save the old positions
184 +  for(i = 0; i < nAtoms; i++){
185 +    atoms[i]->getPos(pos);
186 +    for(j = 0; j < 3; j++)
187 +      oldPos[i*3 + j] = pos[j];
188 +  }
189    
190 +  //the first estimation of r(t+dt) is equal to  r(t)
191 +    
192 +  for(k = 0; k < 4; k ++){
193 +
194 +    for(i =0 ; i < nAtoms; i++){
195 +
196 +      atoms[i]->getVel(vel);
197 +      atoms[i]->getPos(pos);
198 +
199 +      for(j = 0; j < 3; j++)
200 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
201 +      
202 +      info->matVecMul3( eta, rj, sc );
203 +      
204 +      for(j = 0; j < 3; j++)
205 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
206 +
207 +      atoms[i]->setPos( pos );
208 +
209 +    }
210 +
211 +    if (nConstrained) {
212 +      constrainA();
213 +    }
214 +  }  
215 +
216 +
217    // Scale the box after all the positions have been moved:
218    
219    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
220    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
221    
222 +  bigScale = 1.0;
223 +  smallScale = 1.0;
224 +  offDiagMax = 0.0;
225    
226    for(i=0; i<3; i++){
227      for(j=0; j<3; j++){
# Line 178 | Line 239 | void NPTf::moveA() {
239        if (i == j) scaleMat[i][j] = 1.0;
240        // Taylor expansion for the exponential truncated at second order:
241        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
242 <      
242 >
243 >      if (i != j)
244 >        if (fabs(scaleMat[i][j]) > offDiagMax)
245 >          offDiagMax = fabs(scaleMat[i][j]);
246      }
247 +
248 +    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
249 +    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
250    }
251    
252 <  info->getBoxM(hm);
253 <  info->matMul3(hm, scaleMat, hmnew);
254 <  info->setBoxM(hmnew);
252 >  if ((bigScale > 1.1) || (smallScale < 0.9)) {
253 >    sprintf( painCave.errMsg,
254 >             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
255 >             " Check your tauBarostat, as it is probably too small!\n\n"
256 >             " scaleMat = [%lf\t%lf\t%lf]\n"
257 >             "            [%lf\t%lf\t%lf]\n"
258 >             "            [%lf\t%lf\t%lf]\n",
259 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
260 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
261 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
262 >    painCave.isFatal = 1;
263 >    simError();
264 >  } else if (offDiagMax > 0.1) {
265 >    sprintf( painCave.errMsg,
266 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
267 >             " Check your tauBarostat, as it is probably too small!\n\n"
268 >             " scaleMat = [%lf\t%lf\t%lf]\n"
269 >             "            [%lf\t%lf\t%lf]\n"
270 >             "            [%lf\t%lf\t%lf]\n",
271 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
272 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
273 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
274 >    painCave.isFatal = 1;
275 >    simError();
276 >  } else {
277 >    info->getBoxM(hm);
278 >    info->matMul3(hm, scaleMat, hmnew);
279 >    info->setBoxM(hmnew);
280 >  }
281    
282   }
283  
284 < void NPTf::moveB( void ){
284 > template<typename T> void NPTf<T>::moveB( void ){
285  
286 <  int i, j;
286 >  //new version of NPTf
287 >  int i, j, k;
288    DirectionalAtom* dAtom;
289    double Tb[3], ji[3];
290 <  double vel[3], frc[3];
290 >  double vel[3], myVel[3], frc[3];
291    double mass;
292  
293    double instaTemp, instaPress, instaVol;
294    double tt2, tb2;
295    double sc[3];
296    double press[3][3], vScale[3][3];
297 +  double oldChi, prevChi;
298 +  double oldEta[3][3], prevEta[3][3], diffEta;
299    
300    tt2 = tauThermostat * tauThermostat;
301    tb2 = tauBarostat * tauBarostat;
302  
303 <  instaTemp = tStats->getTemperature();
304 <  tStats->getPressureTensor(press);
305 <  instaVol = tStats->getVolume();
210 <  
211 <  // first evolve chi a half step
303 >  // Set things up for the iteration:
304 >
305 >  oldChi = chi;
306    
307 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
308 <  
309 <  for (i = 0; i < 3; i++ ) {
216 <    for (j = 0; j < 3; j++ ) {
217 <      if (i == j) {
307 >  for(i = 0; i < 3; i++)
308 >    for(j = 0; j < 3; j++)
309 >      oldEta[i][j] = eta[i][j];
310  
311 <        eta[i][j] += dt2 * instaVol *
220 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
311 >  for( i=0; i<nAtoms; i++ ){
312  
313 <        vScale[i][j] = eta[i][j] + chi;
223 <        
224 <      } else {
225 <        
226 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
313 >    atoms[i]->getVel( vel );
314  
315 <        vScale[i][j] = eta[i][j];
316 <        
317 <      }
315 >    for (j=0; j < 3; j++)
316 >      oldVel[3*i + j]  = vel[j];
317 >
318 >    if( atoms[i]->isDirectional() ){
319 >
320 >      dAtom = (DirectionalAtom *)atoms[i];
321 >
322 >      dAtom->getJ( ji );
323 >
324 >      for (j=0; j < 3; j++)
325 >        oldJi[3*i + j] = ji[j];
326 >
327      }
328    }
329  
330 <  for( i=0; i<nAtoms; i++ ){
330 >  // do the iteration:
331  
332 <    atoms[i]->getVel( vel );
333 <    atoms[i]->getFrc( frc );
332 >  instaVol = tStats->getVolume();
333 >  
334 >  for (k=0; k < 4; k++) {
335 >    
336 >    instaTemp = tStats->getTemperature();
337 >    tStats->getPressureTensor(press);
338  
339 <    mass = atoms[i]->getMass();
339 >    // evolve chi another half step using the temperature at t + dt/2
340 >
341 >    prevChi = chi;
342 >    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
343      
344 <    // velocity half step
345 <        
346 <    info->matVecMul3( vScale, vel, sc );
244 <    
245 <    for (j = 0; j < 3; j++) {
246 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
247 <    }
344 >    for(i = 0; i < 3; i++)
345 >      for(j = 0; j < 3; j++)
346 >        prevEta[i][j] = eta[i][j];
347  
348 <    atoms[i]->setVel( vel );
348 >    //advance eta half step and calculate scale factor for velocity
349 >
350 >    for(i = 0; i < 3; i ++)
351 >      for(j = 0; j < 3; j++){
352 >        if( i == j) {
353 >          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
354 >            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
355 >          vScale[i][j] = eta[i][j] + chi;
356 >        } else {
357 >          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
358 >          vScale[i][j] = eta[i][j];
359 >        }
360 >      }  
361      
362 <    if( atoms[i]->isDirectional() ){
362 >    for( i=0; i<nAtoms; i++ ){
363  
364 <      dAtom = (DirectionalAtom *)atoms[i];
365 <          
255 <      // get and convert the torque to body frame
364 >      atoms[i]->getFrc( frc );
365 >      atoms[i]->getVel(vel);
366        
367 <      dAtom->getTrq( Tb );
368 <      dAtom->lab2Body( Tb );
367 >      mass = atoms[i]->getMass();
368 >    
369 >      for (j = 0; j < 3; j++)
370 >        myVel[j] = oldVel[3*i + j];
371        
372 <      // get the angular momentum, and propagate a half step
372 >      info->matVecMul3( vScale, myVel, sc );
373        
374 <      dAtom->getJ( ji );
374 >      // velocity half step
375 >      for (j=0; j < 3; j++) {
376 >        // velocity half step  (use chi from previous step here):
377 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
378 >      }
379        
380 <      for (j=0; j < 3; j++)
265 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
380 >      atoms[i]->setVel( vel );
381        
382 <      dAtom->setJ( ji );
382 >      if( atoms[i]->isDirectional() ){
383  
384 <    }                    
384 >        dAtom = (DirectionalAtom *)atoms[i];
385 >  
386 >        // get and convert the torque to body frame      
387 >  
388 >        dAtom->getTrq( Tb );
389 >        dAtom->lab2Body( Tb );      
390 >            
391 >        for (j=0; j < 3; j++)
392 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
393 >      
394 >          dAtom->setJ( ji );
395 >      }
396 >    }
397 >
398 >    if (nConstrained) {
399 >      constrainB();
400 >    }
401 >    
402 >    diffEta = 0;
403 >    for(i = 0; i < 3; i++)
404 >      diffEta += pow(prevEta[i][i] - eta[i][i], 2);    
405 >    
406 >    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
407 >      break;
408    }
409 +
410 +  //calculate integral of chidt
411 +  integralOfChidt += dt2*chi;
412 +  
413   }
414  
415 < int NPTf::readyCheck() {
415 > template<typename T> void NPTf<T>::resetIntegrator() {
416 >  int i,j;
417 >  
418 >  chi = 0.0;
419 >
420 >  for(i = 0; i < 3; i++)
421 >    for (j = 0; j < 3; j++)
422 >      eta[i][j] = 0.0;
423 >
424 > }
425 >
426 > template<typename T> int NPTf<T>::readyCheck() {
427 >
428 >  //check parent's readyCheck() first
429 >  if (T::readyCheck() == -1)
430 >    return -1;
431  
432    // First check to see if we have a target temperature.
433    // Not having one is fatal.
# Line 317 | Line 474 | int NPTf::readyCheck() {
474      return -1;
475    }    
476  
477 <  // We need NkBT a lot, so just set it here:
477 >  
478 >  // We need NkBT a lot, so just set it here: This is the RAW number
479 >  // of particles, so no subtraction or addition of constraints or
480 >  // orientational degrees of freedom:
481 >  
482 >  NkBT = (double)Nparticles * kB * targetTemp;
483 >  
484 >  // fkBT is used because the thermostat operates on more degrees of freedom
485 >  // than the barostat (when there are particles with orientational degrees
486 >  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
487 >  
488 >  fkBT = (double)info->ndf * kB * targetTemp;
489  
322 –  NkBT = (double)info->ndf * kB * targetTemp;
323 –
490    return 1;
491   }
492 +
493 + template<typename T> double NPTf<T>::getConservedQuantity(void){
494 +
495 +  double conservedQuantity;
496 +  double Energy;
497 +  double thermostat_kinetic;
498 +  double thermostat_potential;
499 +  double barostat_kinetic;
500 +  double barostat_potential;
501 +  double trEta;
502 +  double a[3][3], b[3][3];
503 +
504 +  Energy = tStats->getTotalE();
505 +
506 +  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
507 +    (2.0 * eConvert);
508 +
509 +  thermostat_potential = fkBT* integralOfChidt / eConvert;
510 +
511 +  info->transposeMat3(eta, a);
512 +  info->matMul3(a, eta, b);
513 +  trEta = info->matTrace3(b);
514 +
515 +  barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta /
516 +    (2.0 * eConvert);
517 +  
518 +  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
519 +    eConvert;
520 +
521 +  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
522 +    barostat_kinetic + barostat_potential;
523 +  
524 +  cout.width(8);
525 +  cout.precision(8);
526 +
527 +  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
528 +      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
529 +      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
530 +
531 +  return conservedQuantity;
532 + }

Diff Legend

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