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 588 by gezelter, Thu Jul 10 17:10:56 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_++)
34 >    for (j = 0; j < 3; j++)
35        eta[i][j] = 0.0;
36  
37    have_tau_thermostat = 0;
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 <  
59 <  int i,j,k;
60 <  int atomIndex, aMatIndex;
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];
69 <  double ji[3];
70 <  double ri[3], vi[3], sc[3];
71 <  double instaTemp, instaVol;
72 <  double tt2, tb2, eta2ij;
73 <  double angle;
68 >  double Tb[3], ji[3];
69 >  double A[3][3], I[3][3];
70 >  double angle, mass;
71 >  double vel[3], pos[3], frc[3];
72 >
73 >  double rj[3];
74 >  double instaTemp, instaPress, instaVol;
75 >  double tt2, tb2;
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 54 | Line 85 | void NPTf::moveA() {
85    instaTemp = tStats->getTemperature();
86    tStats->getPressureTensor(press);
87    instaVol = tStats->getVolume();
57 –  
58 –  // 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 *
67 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
68 <
69 <        vScale[i][j] = eta[i][j] + chi;
70 <        
71 <      } else {
72 <        
73 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
74 <
75 <        vScale[i][j] = eta[i][j];
76 <        
77 <      }
97 >        vScale[i][j] += chi;          
98 >      }              
99      }
100    }
101 <
101 >  
102 >  //evolve velocity half step
103    for( i=0; i<nAtoms; i++ ){
82 –    atomIndex = i * 3;
83 –    aMatIndex = i * 9;
84 –    
85 –    // velocity half step
86 –    
87 –    vi[0] = vel[atomIndex];
88 –    vi[1] = vel[atomIndex+1];
89 –    vi[2] = vel[atomIndex+2];
90 –    
91 –    info->matVecMul3( vScale, vi, sc );
92 –    
93 –    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
94 –    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
95 –    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);
104  
105 <    vel[atomIndex] = vi[0]
106 <    vel[atomIndex+1] = vi[1];
99 <    vel[atomIndex+2] = vi[2];
105 >    atoms[i]->getVel( vel );
106 >    atoms[i]->getFrc( frc );
107  
108 <    // position whole step    
108 >    mass = atoms[i]->getMass();
109 >    
110 >    info->matVecMul3( vScale, vel, sc );
111  
112 <    ri[0] = pos[atomIndex];
113 <    ri[1] = pos[atomIndex+1];
114 <    ri[2] = pos[atomIndex+2];
112 >    for (j=0; j < 3; j++) {
113 >      // velocity half step
114 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
115 >    }
116  
117 <    info->wrapVector(ri);
108 <
109 <    info->matVecMul3( eta, ri, sc );
110 <
111 <    pos[atomIndex] += dt * (vel[atomIndex] + sc[0]);
112 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + sc[1]);
113 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + sc[2]);
117 >    atoms[i]->setVel( vel );
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 <      Tb[0] = dAtom->getTx();
122 <      Tb[1] = dAtom->getTy();
123 <      Tb[2] = dAtom->getTz();
124 <      
125 >      dAtom->getTrq( Tb );
126        dAtom->lab2Body( Tb );
127        
128        // get the angular momentum, and propagate a half step
129  
130 <      ji[0] = dAtom->getJx();
131 <      ji[1] = dAtom->getJy();
132 <      ji[2] = dAtom->getJz();
130 >      dAtom->getJ( ji );
131 >
132 >      for (j=0; j < 3; j++)
133 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
134        
133 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
134 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
135 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
136 –      
135        // use the angular velocities to propagate the rotation matrix a
136        // full time step
137 <      
137 >
138 >      dAtom->getA(A);
139 >      dAtom->getI(I);
140 >    
141        // rotate about the x-axis      
142 <      angle = dt2 * ji[0] / dAtom->getIxx();
143 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
144 <      
142 >      angle = dt2 * ji[0] / I[0][0];
143 >      this->rotate( 1, 2, angle, ji, A );
144 >
145        // rotate about the y-axis
146 <      angle = dt2 * ji[1] / dAtom->getIyy();
147 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
146 >      angle = dt2 * ji[1] / I[1][1];
147 >      this->rotate( 2, 0, angle, ji, A );
148        
149        // rotate about the z-axis
150 <      angle = dt * ji[2] / dAtom->getIzz();
151 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
150 >      angle = dt * ji[2] / I[2][2];
151 >      this->rotate( 0, 1, angle, ji, A);
152        
153        // rotate about the y-axis
154 <      angle = dt2 * ji[1] / dAtom->getIyy();
155 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
154 >      angle = dt2 * ji[1] / I[1][1];
155 >      this->rotate( 2, 0, angle, ji, A );
156        
157         // rotate about the x-axis
158 <      angle = dt2 * ji[0] / dAtom->getIxx();
159 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
158 >      angle = dt2 * ji[0] / I[0][0];
159 >      this->rotate( 1, 2, angle, ji, A );
160        
161 <      dAtom->setJx( ji[0] );
162 <      dAtom->setJy( ji[1] );
163 <      dAtom->setJz( ji[2] );
161 >      dAtom->setJ( ji );
162 >      dAtom->setA( A  );    
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 <  // Scale the box after all the positions have been moved:
194 >    for(i =0 ; i < nAtoms; i++){
195  
196 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
197 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
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++){
228 <
228 >      
229        // Calculate the matrix Product of the eta array (we only need
230        // the ij element right now):
231 <
231 >      
232        eta2ij = 0.0;
233        for(k=0; k<3; k++){
234          eta2ij += eta[i][k] * eta[k][j];
# Line 187 | Line 240 | void NPTf::moveA() {
240        // Taylor expansion for the exponential truncated at second order:
241        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
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    }
192 –  
193 –  info->getBoxM(hm);
194 –  info->matMul3(hm, scaleMat, hmnew);
195 –  info->setBoxM(hmnew);
251    
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 ){
285 <  int i,j, k;
286 <  int atomIndex;
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];
290 <  double ji[3];
291 <  double vi[3], sc[3];
292 <  double instaTemp, instaVol;
289 >  double Tb[3], ji[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();
216 <  
217 <  // 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++ ) {
222 <    for (j = 0; j < 3; j++ ) {
223 <      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 *
226 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
311 >  for( i=0; i<nAtoms; i++ ){
312  
313 <        vScale[i][j] = eta[i][j] + chi;
229 <        
230 <      } else {
231 <        
232 <        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++ ){
241 <    atomIndex = i * 3;
330 >  // do the iteration:
331  
332 <    // velocity half step
332 >  instaVol = tStats->getVolume();
333 >  
334 >  for (k=0; k < 4; k++) {
335      
336 <    vi[0] = vel[atomIndex];
337 <    vi[1] = vel[atomIndex+1];
338 <    vi[2] = vel[atomIndex+2];
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 <    info->matVecMul3( vScale, vi, sc );
345 <    
346 <    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
252 <    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
253 <    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);
344 >    for(i = 0; i < 3; i++)
345 >      for(j = 0; j < 3; j++)
346 >        prevEta[i][j] = eta[i][j];
347  
348 <    vel[atomIndex] = vi[0]
349 <    vel[atomIndex+1] = vi[1];
350 <    vel[atomIndex+2] = vi[2];
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 >      atoms[i]->getFrc( frc );
365 >      atoms[i]->getVel(vel);
366        
367 <      dAtom = (DirectionalAtom *)atoms[i];
367 >      mass = atoms[i]->getMass();
368 >    
369 >      for (j = 0; j < 3; j++)
370 >        myVel[j] = oldVel[3*i + j];
371        
372 <      // get and convert the torque to body frame
372 >      info->matVecMul3( vScale, myVel, sc );
373        
374 <      Tb[0] = dAtom->getTx();
375 <      Tb[1] = dAtom->getTy();
376 <      Tb[2] = dAtom->getTz();
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 <      dAtom->lab2Body( Tb );
380 >      atoms[i]->setVel( vel );
381        
382 <      // get the angular momentum, and complete the angular momentum
383 <      // half step
382 >      if( atoms[i]->isDirectional() ){
383 >
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 <      ji[0] = dAtom->getJx();
395 <      ji[1] = dAtom->getJy();
276 <      ji[2] = dAtom->getJz();
277 <      
278 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
279 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
280 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
281 <      
282 <      dAtom->setJx( ji[0] );
283 <      dAtom->setJy( ji[1] );
284 <      dAtom->setJz( ji[2] );
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 333 | 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  
338 –  NkBT = (double)info->ndf * kB * targetTemp;
339 –
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