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 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# 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 23 | Line 26 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
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++)
29 >  
30 >  int i,j;
31 >  
32 >  for(i = 0; i < 3; i++){
33 >    for (j = 0; j < 3; j++){
34 >      
35        eta[i][j] = 0.0;
36 +      oldEta[i][j] = 0.0;
37 +    }
38 +  }
39 + }
40  
41 <  have_tau_thermostat = 0;
35 <  have_tau_barostat = 0;
36 <  have_target_temp = 0;
37 <  have_target_pressure = 0;
41 > template<typename T> NPTf<T>::~NPTf() {
42  
43 <  have_chi_tolerance = 0;
44 <  have_eta_tolerance = 0;
41 <  have_pos_iter_tolerance = 0;
43 >  // empty for now
44 > }
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
46 > template<typename T> void NPTf<T>::resetIntegrator() {
47 >  
48 >  int i, j;
49 >  
50 >  for(i = 0; i < 3; i++)
51 >    for (j = 0; j < 3; j++)
52 >      eta[i][j] = 0.0;
53 >  
54 >  T::resetIntegrator();
55   }
56  
57 < template<typename T> NPTf<T>::~NPTf() {
58 <  delete[] oldPos;
59 <  delete[] oldVel;
60 <  delete[] oldJi;
57 > template<typename T> void NPTf<T>::evolveEtaA() {
58 >  
59 >  int i, j;
60 >  
61 >  for(i = 0; i < 3; i ++){
62 >    for(j = 0; j < 3; j++){
63 >      if( i == j)
64 >        eta[i][j] += dt2 *  instaVol *
65 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
66 >      else
67 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
68 >    }
69 >  }
70 >  
71 >  for(i = 0; i < 3; i++)
72 >    for (j = 0; j < 3; j++)
73 >      oldEta[i][j] = eta[i][j];
74   }
75  
76 < template<typename T> void NPTf<T>::moveA() {
76 > template<typename T> void NPTf<T>::evolveEtaB() {
77    
78 <  int i, j, k;
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];
78 >  int i,j;
79  
80 <  double rj[3];
81 <  double instaTemp, instaPress, instaVol;
82 <  double tt2, tb2;
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];
80 >  for(i = 0; i < 3; i++)
81 >    for (j = 0; j < 3; j++)
82 >      prevEta[i][j] = eta[i][j];
83  
84 <  tt2 = tauThermostat * tauThermostat;
85 <  tb2 = tauBarostat * tauBarostat;
84 >  for(i = 0; i < 3; i ++){
85 >    for(j = 0; j < 3; j++){
86 >      if( i == j) {
87 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
88 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
89 >      } else {
90 >        eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
91 >      }
92 >    }
93 >  }
94 > }
95  
96 <  instaTemp = tStats->getTemperature();
97 <  tStats->getPressureTensor(press);
98 <  instaVol = tStats->getVolume();
83 <  
84 <  tStats->getCOM(COM);
96 > template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
97 >  int i,j;
98 >  double vScale[3][3];
99  
86 –  //calculate scale factor of veloity
100    for (i = 0; i < 3; i++ ) {
101      for (j = 0; j < 3; j++ ) {
102        vScale[i][j] = eta[i][j];
# Line 94 | Line 107 | template<typename T> void NPTf<T>::moveA() {
107      }
108    }
109    
110 <  //evolve velocity half step
111 <  for( i=0; i<nAtoms; i++ ){
110 >  info->matVecMul3( vScale, vel, sc );
111 > }
112  
113 <    atoms[i]->getVel( vel );
114 <    atoms[i]->getFrc( frc );
113 > template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
114 >  int i,j;
115 >  double myVel[3];
116 >  double vScale[3][3];
117  
118 <    mass = atoms[i]->getMass();
119 <    
120 <    info->matVecMul3( vScale, vel, sc );
121 <
122 <    for (j=0; j < 3; j++) {
123 <      // velocity half step  (use chi from previous step here):
124 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
110 <  
118 >  for (i = 0; i < 3; i++ ) {
119 >    for (j = 0; j < 3; j++ ) {
120 >      vScale[i][j] = eta[i][j];
121 >      
122 >      if (i == j) {
123 >        vScale[i][j] += chi;          
124 >      }              
125      }
112 –
113 –    atoms[i]->setVel( vel );
114 –  
115 –    if( atoms[i]->isDirectional() ){
116 –
117 –      dAtom = (DirectionalAtom *)atoms[i];
118 –
119 –      // get and convert the torque to body frame
120 –      
121 –      dAtom->getTrq( Tb );
122 –      dAtom->lab2Body( Tb );
123 –      
124 –      // get the angular momentum, and propagate a half step
125 –
126 –      dAtom->getJ( ji );
127 –
128 –      for (j=0; j < 3; j++)
129 –        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
130 –      
131 –      // use the angular velocities to propagate the rotation matrix a
132 –      // full time step
133 –
134 –      dAtom->getA(A);
135 –      dAtom->getI(I);
136 –    
137 –      // rotate about the x-axis      
138 –      angle = dt2 * ji[0] / I[0][0];
139 –      this->rotate( 1, 2, angle, ji, A );
140 –
141 –      // rotate about the y-axis
142 –      angle = dt2 * ji[1] / I[1][1];
143 –      this->rotate( 2, 0, angle, ji, A );
144 –      
145 –      // rotate about the z-axis
146 –      angle = dt * ji[2] / I[2][2];
147 –      this->rotate( 0, 1, angle, ji, A);
148 –      
149 –      // rotate about the y-axis
150 –      angle = dt2 * ji[1] / I[1][1];
151 –      this->rotate( 2, 0, angle, ji, A );
152 –      
153 –       // rotate about the x-axis
154 –      angle = dt2 * ji[0] / I[0][0];
155 –      this->rotate( 1, 2, angle, ji, A );
156 –      
157 –      dAtom->setJ( ji );
158 –      dAtom->setA( A  );    
159 –    }    
126    }
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 –  }
127    
128 <  //the first estimation of r(t+dt) is equal to  r(t)
129 <    
187 <  for(k = 0; k < 4; k ++){
128 >  for (j = 0; j < 3; j++)
129 >    myVel[j] = oldVel[3*index + j];
130  
131 <    for(i =0 ; i < nAtoms; i++){
131 >  info->matVecMul3( vScale, myVel, sc );
132 > }
133  
134 <      atoms[i]->getVel(vel);
135 <      atoms[i]->getPos(pos);
134 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
135 >                                               int index, double sc[3]){
136 >  int j;
137 >  double rj[3];
138  
139 <      for(j = 0; j < 3; j++)
140 <        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]);
139 >  for(j=0; j<3; j++)
140 >    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
141  
142 <      atoms[i]->setPos( pos );
142 >  info->matVecMul3( eta, rj, sc );
143 > }
144  
145 <    }
145 > template<typename T> void NPTf<T>::scaleSimBox( void ){
146  
147 <  }  
147 >  int i,j,k;
148 >  double scaleMat[3][3];
149 >  double eta2ij;
150 >  double bigScale, smallScale, offDiagMax;
151 >  double hm[3][3], hmnew[3][3];
152 >  
153  
154 <
154 >
155    // Scale the box after all the positions have been moved:
156    
157    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
# Line 270 | Line 216 | template<typename T> void NPTf<T>::moveA() {
216      info->matMul3(hm, scaleMat, hmnew);
217      info->setBoxM(hmnew);
218    }
273 –  
219   }
220  
221 < template<typename T> void NPTf<T>::moveB( void ){
221 > template<typename T> bool NPTf<T>::etaConverged() {
222 >  int i;
223 >  double diffEta, sumEta;
224  
225 <  int i, j, k;
279 <  DirectionalAtom* dAtom;
280 <  double Tb[3], ji[3];
281 <  double vel[3], frc[3];
282 <  double mass;
283 <
284 <  double instaTemp, instaPress, instaVol;
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 <
294 <
295 <  // Set things up for the iteration:
296 <
297 <  oldChi = chi;
298 <  
225 >  sumEta = 0;
226    for(i = 0; i < 3; i++)
227 <    for(j = 0; j < 3; j++)
301 <      oldEta[i][j] = eta[i][j];
302 <
303 <  for( i=0; i<nAtoms; i++ ){
304 <
305 <    atoms[i]->getVel( vel );
306 <
307 <    for (j=0; j < 3; j++)
308 <      oldVel[3*i + j]  = vel[j];
309 <
310 <    if( atoms[i]->isDirectional() ){
311 <
312 <      dAtom = (DirectionalAtom *)atoms[i];
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();
227 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
228    
229 <  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 <      mass = atoms[i]->getMass();
362 <      
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 <      atoms[i]->setVel( vel );
371 <      
372 <      if( atoms[i]->isDirectional() ){
373 <
374 <        dAtom = (DirectionalAtom *)atoms[i];
229 >  diffEta = sqrt( sumEta / 3.0 );
230    
231 <        // 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 );
385 <      }
386 <    }
387 <
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 <  
231 >  return ( diffEta <= etaTolerance );
232   }
233  
234 < template<typename T> void NPTf<T>::resetIntegrator() {
404 <  int i,j;
234 > template<typename T> double NPTf<T>::getConservedQuantity(void){
235    
236 <  chi = 0.0;
236 >  double conservedQuantity;
237 >  double totalEnergy;
238 >  double thermostat_kinetic;
239 >  double thermostat_potential;
240 >  double barostat_kinetic;
241 >  double barostat_potential;
242 >  double trEta;
243 >  double a[3][3], b[3][3];
244  
245 <  for(i = 0; i < 3; i++)
409 <    for (j = 0; j < 3; j++)
410 <      eta[i][j] = 0.0;
245 >  totalEnergy = tStats->getTotalE();
246  
247 < }
247 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
248 >    (2.0 * eConvert);
249  
250 < template<typename T> int NPTf<T>::readyCheck() {
250 >  thermostat_potential = fkBT* integralOfChidt / eConvert;
251  
252 <  //check parent's readyCheck() first
253 <  if (T::readyCheck() == -1)
254 <    return -1;
255 <
256 <  // First check to see if we have a target temperature.
257 <  // Not having one is fatal.
252 >  info->transposeMat3(eta, a);
253 >  info->matMul3(a, eta, b);
254 >  trEta = info->matTrace3(b);
255 >
256 >  barostat_kinetic = NkBT * tb2 * trEta /
257 >    (2.0 * eConvert);
258    
259 <  if (!have_target_temp) {
260 <    sprintf( painCave.errMsg,
425 <             "NPTf error: You can't use the NPTf integrator\n"
426 <             "   without a targetTemp!\n"
427 <             );
428 <    painCave.isFatal = 1;
429 <    simError();
430 <    return -1;
431 <  }
259 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
260 >    eConvert;
261  
262 <  if (!have_target_pressure) {
263 <    sprintf( painCave.errMsg,
435 <             "NPTf error: You can't use the NPTf integrator\n"
436 <             "   without a targetPressure!\n"
437 <             );
438 <    painCave.isFatal = 1;
439 <    simError();
440 <    return -1;
441 <  }
262 >  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
263 >    barostat_kinetic + barostat_potential;
264    
265 <  // We must set tauThermostat.
266 <  
445 <  if (!have_tau_thermostat) {
446 <    sprintf( painCave.errMsg,
447 <             "NPTf error: If you use the NPTf\n"
448 <             "   integrator, you must set tauThermostat.\n");
449 <    painCave.isFatal = 1;
450 <    simError();
451 <    return -1;
452 <  }    
265 > //   cout.width(8);
266 > //   cout.precision(8);
267  
268 <  // We must set tauBarostat.
269 <  
270 <  if (!have_tau_barostat) {
457 <    sprintf( painCave.errMsg,
458 <             "NPTf error: If you use the NPTf\n"
459 <             "   integrator, you must set tauBarostat.\n");
460 <    painCave.isFatal = 1;
461 <    simError();
462 <    return -1;
463 <  }    
268 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270 > //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
271  
272 <  // We need NkBT a lot, so just set it here:
466 <
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);
272 >  return conservedQuantity;
273    
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 –  
274   }

Diff Legend

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