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 763 by tim, Mon Sep 15 16:52:02 2003 UTC vs.
Revision 778 by mmeineke, Fri Sep 19 20:00:27 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 35 | 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];
69 <  double A[3][3], I[3][3];
70 <  double angle, mass;
69 >
70 >  double mass;
71    double vel[3], pos[3], frc[3];
72  
73    double rj[3];
# Line 53 | 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 60 | Line 85 | template<typename T> void NPTf<T>::moveA() {
85    instaTemp = tStats->getTemperature();
86    tStats->getPressureTensor(press);
87    instaVol = tStats->getVolume();
63 –  
64 –  // 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 *
73 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
74 <        
75 <        vScale[i][j] = eta[i][j] + chi;
76 <          
77 <      } else {
78 <        
79 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
80 <
81 <        vScale[i][j] = eta[i][j];
82 <        
83 <      }
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 );
90 –    atoms[i]->getPos( pos );
106      atoms[i]->getFrc( frc );
107  
108      mass = atoms[i]->getMass();
109      
95 –    // velocity half step
96 –        
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]);
101 –      rj[j] = pos[j];
115      }
116  
117      atoms[i]->setVel( vel );
105 –
106 –    // position whole step    
107 –
108 –    info->wrapVector(rj);
109 –
110 –    info->matVecMul3( eta, rj, sc );
111 –
112 –    for (j = 0; j < 3; j++ )
113 –      pos[j] += dt * (vel[j] + sc[j]);
114 –
115 –    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 129 | Line 131 | template<typename T> void NPTf<T>::moveA() {
131  
132        for (j=0; j < 3; j++)
133          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
132 –      
133 –      // use the angular velocities to propagate the rotation matrix a
134 –      // full time step
135 –
136 –      dAtom->getA(A);
137 –      dAtom->getI(I);
138 –    
139 –      // rotate about the x-axis      
140 –      angle = dt2 * ji[0] / I[0][0];
141 –      this->rotate( 1, 2, angle, ji, A );
142 –
143 –      // rotate about the y-axis
144 –      angle = dt2 * ji[1] / I[1][1];
145 –      this->rotate( 2, 0, angle, ji, A );
146 –      
147 –      // rotate about the z-axis
148 –      angle = dt * ji[2] / I[2][2];
149 –      this->rotate( 0, 1, angle, ji, A);
150 –      
151 –      // rotate about the y-axis
152 –      angle = dt2 * ji[1] / I[1][1];
153 –      this->rotate( 2, 0, angle, ji, A );
154 –      
155 –       // rotate about the x-axis
156 –      angle = dt2 * ji[0] / I[0][0];
157 –      this->rotate( 1, 2, angle, ji, A );
134        
135 +      this->rotationPropagation( dAtom, ji );
136 +  
137        dAtom->setJ( ji );
138 <      dAtom->setA( A  );    
161 <    }                    
138 >    }    
139    }
140 +
141 +  // advance chi half step
142 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
143 +
144 +  // calculate the integral of chidt
145 +  integralOfChidt += dt2*chi;
146 +
147 +  // advance eta half step
148 +
149 +  for(i = 0; i < 3; i ++)
150 +    for(j = 0; j < 3; j++){
151 +      if( i == j)
152 +        eta[i][j] += dt2 *  instaVol *
153 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
154 +      else
155 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
156 +    }
157 +    
158 +  //save the old positions
159 +  for(i = 0; i < nAtoms; i++){
160 +    atoms[i]->getPos(pos);
161 +    for(j = 0; j < 3; j++)
162 +      oldPos[i*3 + j] = pos[j];
163 +  }
164    
165 +  //the first estimation of r(t+dt) is equal to  r(t)
166 +    
167 +  for(k = 0; k < 4; k ++){
168 +
169 +    for(i =0 ; i < nAtoms; i++){
170 +
171 +      atoms[i]->getVel(vel);
172 +      atoms[i]->getPos(pos);
173 +
174 +      for(j = 0; j < 3; j++)
175 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
176 +      
177 +      info->matVecMul3( eta, rj, sc );
178 +      
179 +      for(j = 0; j < 3; j++)
180 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
181 +
182 +      atoms[i]->setPos( pos );
183 +
184 +    }
185 +
186 +    if (nConstrained) {
187 +      constrainA();
188 +    }
189 +  }  
190 +
191 +
192    // Scale the box after all the positions have been moved:
193    
194    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
# Line 230 | Line 258 | template<typename T> void NPTf<T>::moveB( void ){
258  
259   template<typename T> void NPTf<T>::moveB( void ){
260  
261 <  int i, j;
261 >  //new version of NPTf
262 >  int i, j, k;
263    DirectionalAtom* dAtom;
264    double Tb[3], ji[3];
265 <  double vel[3], frc[3];
265 >  double vel[3], myVel[3], frc[3];
266    double mass;
267  
268    double instaTemp, instaPress, instaVol;
269    double tt2, tb2;
270    double sc[3];
271    double press[3][3], vScale[3][3];
272 +  double oldChi, prevChi;
273 +  double oldEta[3][3], prevEta[3][3], diffEta;
274    
275    tt2 = tauThermostat * tauThermostat;
276    tb2 = tauBarostat * tauBarostat;
277  
278 <  instaTemp = tStats->getTemperature();
279 <  tStats->getPressureTensor(press);
280 <  instaVol = tStats->getVolume();
250 <  
251 <  // first evolve chi a half step
278 >  // Set things up for the iteration:
279 >
280 >  oldChi = chi;
281    
282 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
283 <  
284 <  for (i = 0; i < 3; i++ ) {
256 <    for (j = 0; j < 3; j++ ) {
257 <      if (i == j) {
282 >  for(i = 0; i < 3; i++)
283 >    for(j = 0; j < 3; j++)
284 >      oldEta[i][j] = eta[i][j];
285  
286 <        eta[i][j] += dt2 * instaVol *
260 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
286 >  for( i=0; i<nAtoms; i++ ){
287  
288 <        vScale[i][j] = eta[i][j] + chi;
263 <        
264 <      } else {
265 <        
266 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
288 >    atoms[i]->getVel( vel );
289  
290 <        vScale[i][j] = eta[i][j];
291 <        
292 <      }
290 >    for (j=0; j < 3; j++)
291 >      oldVel[3*i + j]  = vel[j];
292 >
293 >    if( atoms[i]->isDirectional() ){
294 >
295 >      dAtom = (DirectionalAtom *)atoms[i];
296 >
297 >      dAtom->getJ( ji );
298 >
299 >      for (j=0; j < 3; j++)
300 >        oldJi[3*i + j] = ji[j];
301 >
302      }
303    }
304  
305 <  for( i=0; i<nAtoms; i++ ){
305 >  // do the iteration:
306  
307 <    atoms[i]->getVel( vel );
308 <    atoms[i]->getFrc( frc );
307 >  instaVol = tStats->getVolume();
308 >  
309 >  for (k=0; k < 4; k++) {
310 >    
311 >    instaTemp = tStats->getTemperature();
312 >    tStats->getPressureTensor(press);
313  
314 <    mass = atoms[i]->getMass();
314 >    // evolve chi another half step using the temperature at t + dt/2
315 >
316 >    prevChi = chi;
317 >    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
318      
319 <    // velocity half step
320 <        
321 <    info->matVecMul3( vScale, vel, sc );
284 <    
285 <    for (j = 0; j < 3; j++) {
286 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 <    }
319 >    for(i = 0; i < 3; i++)
320 >      for(j = 0; j < 3; j++)
321 >        prevEta[i][j] = eta[i][j];
322  
323 <    atoms[i]->setVel( vel );
323 >    //advance eta half step and calculate scale factor for velocity
324 >
325 >    for(i = 0; i < 3; i ++)
326 >      for(j = 0; j < 3; j++){
327 >        if( i == j) {
328 >          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
329 >            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
330 >          vScale[i][j] = eta[i][j] + chi;
331 >        } else {
332 >          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
333 >          vScale[i][j] = eta[i][j];
334 >        }
335 >      }  
336      
337 <    if( atoms[i]->isDirectional() ){
337 >    for( i=0; i<nAtoms; i++ ){
338  
339 <      dAtom = (DirectionalAtom *)atoms[i];
340 <          
295 <      // get and convert the torque to body frame
339 >      atoms[i]->getFrc( frc );
340 >      atoms[i]->getVel(vel);
341        
342 <      dAtom->getTrq( Tb );
343 <      dAtom->lab2Body( Tb );
342 >      mass = atoms[i]->getMass();
343 >    
344 >      for (j = 0; j < 3; j++)
345 >        myVel[j] = oldVel[3*i + j];
346        
347 <      // get the angular momentum, and propagate a half step
347 >      info->matVecMul3( vScale, myVel, sc );
348        
349 <      dAtom->getJ( ji );
349 >      // velocity half step
350 >      for (j=0; j < 3; j++) {
351 >        // velocity half step  (use chi from previous step here):
352 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
353 >      }
354        
355 <      for (j=0; j < 3; j++)
305 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
355 >      atoms[i]->setVel( vel );
356        
357 <      dAtom->setJ( ji );
357 >      if( atoms[i]->isDirectional() ){
358  
359 <    }                    
359 >        dAtom = (DirectionalAtom *)atoms[i];
360 >  
361 >        // get and convert the torque to body frame      
362 >  
363 >        dAtom->getTrq( Tb );
364 >        dAtom->lab2Body( Tb );      
365 >            
366 >        for (j=0; j < 3; j++)
367 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
368 >      
369 >          dAtom->setJ( ji );
370 >      }
371 >    }
372 >
373 >    if (nConstrained) {
374 >      constrainB();
375 >    }
376 >    
377 >    diffEta = 0;
378 >    for(i = 0; i < 3; i++)
379 >      diffEta += pow(prevEta[i][i] - eta[i][i], 2);    
380 >    
381 >    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
382 >      break;
383    }
384 +
385 +  //calculate integral of chidt
386 +  integralOfChidt += dt2*chi;
387 +  
388   }
389  
390   template<typename T> void NPTf<T>::resetIntegrator() {
# Line 372 | Line 449 | template<typename T> int NPTf<T>::readyCheck() {
449      return -1;
450    }    
451  
452 <  // We need NkBT a lot, so just set it here:
452 >  
453 >  // We need NkBT a lot, so just set it here: This is the RAW number
454 >  // of particles, so no subtraction or addition of constraints or
455 >  // orientational degrees of freedom:
456 >  
457 >  NkBT = (double)Nparticles * kB * targetTemp;
458 >  
459 >  // fkBT is used because the thermostat operates on more degrees of freedom
460 >  // than the barostat (when there are particles with orientational degrees
461 >  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
462 >  
463 >  fkBT = (double)info->ndf * kB * targetTemp;
464  
377 –  NkBT = (double)info->ndf * kB * targetTemp;
378 –
465    return 1;
466   }
467  
468   template<typename T> double NPTf<T>::getConservedQuantity(void){
469  
470    double conservedQuantity;
471 <  double tb2;
472 <  double eta2[3][3];  
471 >  double Energy;
472 >  double thermostat_kinetic;
473 >  double thermostat_potential;
474 >  double barostat_kinetic;
475 >  double barostat_potential;
476    double trEta;
477 +  double a[3][3], b[3][3];
478  
479 <  //HNVE
390 <  conservedQuantity = tStats->getTotalE();
479 >  Energy = tStats->getTotalE();
480  
481 <  //HNVT
482 <  conservedQuantity += (info->getNDF() * kB * targetTemp *
394 <    (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;
481 >  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
482 >    (2.0 * eConvert);
483  
484 <  //HNPT
397 <  tb2 = tauBarostat *tauBarostat;
484 >  thermostat_potential = fkBT* integralOfChidt / eConvert;
485  
486 <  trEta = info->matTrace3(eta);
486 >  info->transposeMat3(eta, a);
487 >  info->matMul3(a, eta, b);
488 >  trEta = info->matTrace3(b);
489 >
490 >  barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta /
491 >    (2.0 * eConvert);
492    
493 <  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
494 <                        3*NkBT/2 * tb2 * trEta * trEta) / eConvert;
493 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
494 >    eConvert;
495 >
496 >  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
497 >    barostat_kinetic + barostat_potential;
498    
499 +  cout.width(8);
500 +  cout.precision(8);
501 +
502 +  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
503 +      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
504 +      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
505 +
506    return conservedQuantity;
507   }

Diff Legend

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