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 746 by mmeineke, Thu Sep 4 21:48:35 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 25 | Line 28 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
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 34 | Line 38 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
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 + 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 <  
64 >
65 >  // new version of NPTf
66    int i, j, k;
67    DirectionalAtom* dAtom;
68    double Tb[3], ji[3];
# Line 52 | Line 77 | template<typename T> void NPTf<T>::moveA() {
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 59 | Line 85 | template<typename T> void NPTf<T>::moveA() {
85    instaTemp = tStats->getTemperature();
86    tStats->getPressureTensor(press);
87    instaVol = tStats->getVolume();
62 –  
63 –  // 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 *
72 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
73 <        
74 <        vScale[i][j] = eta[i][j] + chi;
75 <        
76 <      } else {
77 <        
78 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
79 <
80 <        vScale[i][j] = eta[i][j];
81 <        
82 <      }
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 );
89 –    atoms[i]->getPos( pos );
106      atoms[i]->getFrc( frc );
107  
108      mass = atoms[i]->getMass();
109      
94 –    // velocity half step
95 –        
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]);
100 –      rj[j] = pos[j];
115      }
116  
117      atoms[i]->setVel( vel );
104 –
105 –    // position whole step    
106 –
107 –    info->wrapVector(rj);
108 –
109 –    info->matVecMul3( eta, rj, sc );
110 –
111 –    for (j = 0; j < 3; j++ )
112 –      pos[j] += dt * (vel[j] + sc[j]);
113 –
114 –    atoms[i]->setPos( pos );
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 157 | Line 160 | template<typename T> void NPTf<T>::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)
# Line 189 | Line 243 | template<typename T> void NPTf<T>::moveA() {
243        if (i != j)
244          if (fabs(scaleMat[i][j]) > offDiagMax)
245            offDiagMax = fabs(scaleMat[i][j]);
192 –      
246      }
247  
248      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 230 | Line 283 | template<typename T> void NPTf<T>::moveB( void ){
283  
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();
250 <  
251 <  // 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++ ) {
256 <    for (j = 0; j < 3; j++ ) {
257 <      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  
259 –        eta[i][j] += dt2 * instaVol *
260 –          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
261 –
262 –        vScale[i][j] = eta[i][j] + chi;
263 –        
264 –      } else {
265 –        
266 –        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
267 –
268 –        vScale[i][j] = eta[i][j];
269 –        
270 –      }
271 –    }
272 –  }
273 –
311    for( i=0; i<nAtoms; i++ ){
312  
313      atoms[i]->getVel( vel );
277 –    atoms[i]->getFrc( frc );
314  
315 <    mass = atoms[i]->getMass();
316 <    
281 <    // velocity half step
282 <        
283 <    info->matVecMul3( vScale, vel, sc );
284 <    
285 <    for (j = 0; j < 3; j++) {
286 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 <    }
315 >    for (j=0; j < 3; j++)
316 >      oldVel[3*i + j]  = vel[j];
317  
289 –    atoms[i]->setVel( vel );
290 –    
318      if( atoms[i]->isDirectional() ){
319  
320        dAtom = (DirectionalAtom *)atoms[i];
321 <          
322 <      // get and convert the torque to body frame
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 >  // do the iteration:
331 >
332 >  instaVol = tStats->getVolume();
333 >  
334 >  for (k=0; k < 4; k++) {
335 >    
336 >    instaTemp = tStats->getTemperature();
337 >    tStats->getPressureTensor(press);
338 >
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 >    for(i = 0; i < 3; i++)
345 >      for(j = 0; j < 3; 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) {
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 >    for( i=0; i<nAtoms; i++ ){
363 >
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++)
305 <        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   template<typename T> void NPTf<T>::resetIntegrator() {
# Line 372 | 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:
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  
377 –  NkBT = (double)info->ndf * kB * targetTemp;
378 –
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