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 767 by tim, Tue Sep 16 20:02:11 2003 UTC vs.
Revision 772 by gezelter, Fri Sep 19 16:01:07 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   // Basic non-isotropic thermostating and barostating via the Melchionna
17   // modification of the Hoover algorithm:
# Line 48 | Line 51 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
51   #else
52    Nparticles = theInfo->n_atoms;
53   #endif
54 +
55   }
56  
57   template<typename T> NPTf<T>::~NPTf() {
# Line 57 | Line 61 | template<typename T> void NPTf<T>::moveA() {
61   }
62  
63   template<typename T> void NPTf<T>::moveA() {
64 <  
64 >
65 >  // new version of NPTf
66    int i, j, k;
67    DirectionalAtom* dAtom;
68    double Tb[3], ji[3];
# Line 105 | Line 110 | template<typename T> void NPTf<T>::moveA() {
110      info->matVecMul3( vScale, vel, sc );
111  
112      for (j=0; j < 3; j++) {
113 <      // velocity half step  (use chi from previous step here):
113 >      // velocity half step
114        vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
110 –  
115      }
116  
117      atoms[i]->setVel( vel );
# Line 162 | Line 166 | template<typename T> void NPTf<T>::moveA() {
166    // advance chi half step
167    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
168  
169 <  //calculate the integral of chidt
169 >  // calculate the integral of chidt
170    integralOfChidt += dt2*chi;
171  
172 <  //advance eta half step
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);
180 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
181      }
182      
183    //save the old positions
# Line 203 | Line 208 | template<typename T> void NPTf<T>::moveA() {
208  
209      }
210  
211 +    if (nConstrained) {
212 +      constrainA();
213 +    }
214    }  
215  
216  
# Line 275 | Line 283 | template<typename T> void NPTf<T>::moveB( void ){
283  
284   template<typename T> void NPTf<T>::moveB( void ){
285  
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;
# Line 286 | Line 295 | template<typename T> void NPTf<T>::moveB( void ){
295    double sc[3];
296    double press[3][3], vScale[3][3];
297    double oldChi, prevChi;
298 <  double oldEta[3][3], preEta[3][3], diffEta;
298 >  double oldEta[3][3], prevEta[3][3], diffEta;
299    
300    tt2 = tauThermostat * tauThermostat;
301    tb2 = tauBarostat * tauBarostat;
302  
294 –
303    // Set things up for the iteration:
304  
305    oldChi = chi;
# Line 335 | Line 343 | template<typename T> void NPTf<T>::moveB( void ){
343      
344      for(i = 0; i < 3; i++)
345        for(j = 0; j < 3; j++)
346 <        preEta[i][j] = eta[i][j];
346 >        prevEta[i][j] = eta[i][j];
347  
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){
352 >        if( i == j) {
353            eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
354 <            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
354 >            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
355            vScale[i][j] = eta[i][j] + chi;
356 <        }
348 <        else
349 <        {
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 <
355 <    //advance velocity half step
360 >      }  
361 >    
362      for( i=0; i<nAtoms; i++ ){
363  
364        atoms[i]->getFrc( frc );
365        atoms[i]->getVel(vel);
366        
367        mass = atoms[i]->getMass();
368 +    
369 +      for (j = 0; j < 3; j++)
370 +        myVel[j] = oldVel[3*i + j];
371        
372 <      info->matVecMul3( vScale, vel, sc );
373 <
372 >      info->matVecMul3( vScale, myVel, sc );
373 >      
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]);
# Line 385 | Line 395 | template<typename T> void NPTf<T>::moveB( void ){
395        }
396      }
397  
398 +    if (nConstrained) {
399 +      constrainB();
400 +    }
401      
402      diffEta = 0;
403      for(i = 0; i < 3; i++)
404 <      diffEta += pow(preEta[i][i] - eta[i][i], 2);    
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 chida
410 >  //calculate integral of chidt
411    integralOfChidt += dt2*chi;
399 –
412    
413   }
414  
# Line 462 | Line 474 | template<typename T> int NPTf<T>::readyCheck() {
474      return -1;
475    }    
476  
477 <  // We need NkBT a lot, so just set it here:
478 <
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  
490    return 1;
# Line 473 | Line 493 | template<typename T> double NPTf<T>::getConservedQuant
493   template<typename T> double NPTf<T>::getConservedQuantity(void){
494  
495    double conservedQuantity;
496 <  double tb2;
497 <  double trEta;  
498 <  double U;
499 <  double thermo;
500 <  double integral;
501 <  double baro;
502 <  double PV;
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 <  U = tStats->getTotalE();
485 <  thermo = (fkBT * tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
504 >  Energy = tStats->getTotalE();
505  
506 <  tb2 = tauBarostat * tauBarostat;
507 <  trEta = info->matTrace3(eta);
489 <  baro = ((double)info->ndfTrans * kB * targetTemp * tb2 * trEta * trEta / 2.0) / eConvert;
506 >  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
507 >    (2.0 * eConvert);
508  
509 <  integral = ((double)(info->ndf + 1) * kB * targetTemp * integralOfChidt) /eConvert;
509 >  thermostat_potential = fkBT* integralOfChidt / eConvert;
510  
511 <  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
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);
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;
526  
527 <  conservedQuantity = U+thermo+integral+PV+baro;
528 <  return conservedQuantity;
529 <  
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