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 577 by gezelter, Wed Jul 9 01:41:11 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 9 | Line 10
10   #include "simError.h"
11  
12  
13 < // Basic isotropic thermostating and barostating via the Melchionna
13 > // Basic non-isotropic thermostating and barostating via the Melchionna
14   // modification of the Hoover algorithm:
15   //
16   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
# Line 19 | Line 20 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
20   //
21   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
22  
23 < NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
24 <  Integrator( theInfo, the_ff )
23 > template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
24 >  T( theInfo, the_ff )
25   {
26 <  int i;
26 >  int i, j;
27    chi = 0.0;
28 <  for(i = 0; i < 9; i++) eta[i] = 0.0;
28 >  integralOfChidt = 0.0;
29 >
30 >  for(i = 0; i < 3; i++)
31 >    for (j = 0; j < 3; j++)
32 >      eta[i][j] = 0.0;
33 >
34    have_tau_thermostat = 0;
35    have_tau_barostat = 0;
36    have_target_temp = 0;
37    have_target_pressure = 0;
38 +
39 +  have_chi_tolerance = 0;
40 +  have_eta_tolerance = 0;
41 +  have_pos_iter_tolerance = 0;
42 +
43 +  oldPos = new double[3*nAtoms];
44 +  oldVel = new double[3*nAtoms];
45 +  oldJi = new double[3*nAtoms];
46 + #ifdef IS_MPI
47 +  Nparticles = mpiSim->getTotAtoms();
48 + #else
49 +  Nparticles = theInfo->n_atoms;
50 + #endif
51   }
52  
53 < void NPTf::moveA() {
53 > template<typename T> NPTf<T>::~NPTf() {
54 >  delete[] oldPos;
55 >  delete[] oldVel;
56 >  delete[] oldJi;
57 > }
58 >
59 > template<typename T> void NPTf<T>::moveA() {
60    
61 <  int i,j,k;
37 <  int atomIndex, aMatIndex;
61 >  int i, j, k;
62    DirectionalAtom* dAtom;
63 <  double Tb[3];
64 <  double ji[3];
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];
67 >
68    double rj[3];
69    double instaTemp, instaPress, instaVol;
70    double tt2, tb2;
71 <  double angle;
72 <  double press[9];
73 <  const double p_convert = 1.63882576e8;
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];
76  
77    tt2 = tauThermostat * tauThermostat;
78    tb2 = tauBarostat * tauBarostat;
79  
80    instaTemp = tStats->getTemperature();
81    tStats->getPressureTensor(press);
53 –
54 –  for (i=0; i < 9; i++) press[i] *= p_convert;
55 –
82    instaVol = tStats->getVolume();
57 –  
58 –  // first evolve chi a half step
83    
84 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
84 >  tStats->getCOM(COM);
85 >
86 >  //calculate scale factor of veloity
87 >  for (i = 0; i < 3; i++ ) {
88 >    for (j = 0; j < 3; j++ ) {
89 >      vScale[i][j] = eta[i][j];
90 >      
91 >      if (i == j) {
92 >        vScale[i][j] += chi;          
93 >      }              
94 >    }
95 >  }
96    
97 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
63 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
64 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
65 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
66 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
67 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
68 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
69 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
70 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
71 <  
97 >  //evolve velocity half step
98    for( i=0; i<nAtoms; i++ ){
73 –    atomIndex = i * 3;
74 –    aMatIndex = i * 9;
75 –    
76 –    // velocity half step
77 –    
78 –    vx = vel[atomIndex];
79 –    vy = vel[atomIndex+1];
80 –    vz = vel[atomIndex+2];
81 –    
82 –    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
83 –    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
84 –    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
85 –    
86 –    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
87 –    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
88 –    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
99  
100 <    vel[atomIndex] = vx;
101 <    vel[atomIndex+1] = vy;
92 <    vel[atomIndex+2] = vz;
100 >    atoms[i]->getVel( vel );
101 >    atoms[i]->getFrc( frc );
102  
103 <    // position whole step    
103 >    mass = atoms[i]->getMass();
104 >    
105 >    info->matVecMul3( vScale, vel, sc );
106  
107 <    rj[0] = pos[atomIndex];
108 <    rj[1] = pos[atomIndex+1];
109 <    rj[2] = pos[atomIndex+2];
107 >    for (j=0; j < 3; j++) {
108 >      // velocity half step  (use chi from previous step here):
109 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
110 >  
111 >    }
112  
113 <    info->wrapVector(rj);
101 <
102 <    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
103 <    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
104 <    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
105 <
106 <    pos[atomIndex] += dt * (vel[atomIndex] + scx);
107 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
108 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
113 >    atoms[i]->setVel( vel );
114    
115      if( atoms[i]->isDirectional() ){
116  
117        dAtom = (DirectionalAtom *)atoms[i];
118 <          
118 >
119        // get and convert the torque to body frame
120        
121 <      Tb[0] = dAtom->getTx();
117 <      Tb[1] = dAtom->getTy();
118 <      Tb[2] = dAtom->getTz();
119 <      
121 >      dAtom->getTrq( Tb );
122        dAtom->lab2Body( Tb );
123        
124        // get the angular momentum, and propagate a half step
125  
126 <      ji[0] = dAtom->getJx();
127 <      ji[1] = dAtom->getJy();
128 <      ji[2] = dAtom->getJz();
126 >      dAtom->getJ( ji );
127 >
128 >      for (j=0; j < 3; j++)
129 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
130        
128 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
129 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
130 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
131 –      
131        // use the angular velocities to propagate the rotation matrix a
132        // full time step
133 <      
133 >
134 >      dAtom->getA(A);
135 >      dAtom->getI(I);
136 >    
137        // rotate about the x-axis      
138 <      angle = dt2 * ji[0] / dAtom->getIxx();
139 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
140 <      
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] / dAtom->getIyy();
143 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
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] / dAtom->getIzz();
147 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
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] / dAtom->getIyy();
151 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
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] / dAtom->getIxx();
155 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
154 >      angle = dt2 * ji[0] / I[0][0];
155 >      this->rotate( 1, 2, angle, ji, A );
156        
157 <      dAtom->setJx( ji[0] );
158 <      dAtom->setJy( ji[1] );
159 <      dAtom->setJz( ji[2] );
157 >      dAtom->setJ( ji );
158 >      dAtom->setA( A  );    
159 >    }    
160 >  }
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    }
161 –
162 –  // Scale the box after all the positions have been moved:
184    
185 +  //the first estimation of r(t+dt) is equal to  r(t)
186 +    
187 +  for(k = 0; k < 4; k ++){
188  
189 +    for(i =0 ; i < nAtoms; i++){
190  
191 <  // Use a taylor expansion for eta products
192 <  
168 <  info->getBoxM(hm);
169 <  
191 >      atoms[i]->getVel(vel);
192 >      atoms[i]->getPos(pos);
193  
194 +      for(j = 0; j < 3; j++)
195 +        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]);
201  
202 +      atoms[i]->setPos( pos );
203  
204 +    }
205  
206 +  }  
207  
208 <   info->scaleBox(exp(dt*eta));
208 >
209 >  // Scale the box after all the positions have been moved:
210 >  
211 >  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
212 >  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
213 >  
214 >  bigScale = 1.0;
215 >  smallScale = 1.0;
216 >  offDiagMax = 0.0;
217 >  
218 >  for(i=0; i<3; i++){
219 >    for(j=0; j<3; j++){
220 >      
221 >      // Calculate the matrix Product of the eta array (we only need
222 >      // the ij element right now):
223 >      
224 >      eta2ij = 0.0;
225 >      for(k=0; k<3; k++){
226 >        eta2ij += eta[i][k] * eta[k][j];
227 >      }
228 >      
229 >      scaleMat[i][j] = 0.0;
230 >      // identity matrix (see above):
231 >      if (i == j) scaleMat[i][j] = 1.0;
232 >      // Taylor expansion for the exponential truncated at second order:
233 >      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
234  
235 +      if (i != j)
236 +        if (fabs(scaleMat[i][j]) > offDiagMax)
237 +          offDiagMax = fabs(scaleMat[i][j]);
238 +    }
239  
240 +    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
241 +    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
242 +  }
243 +  
244 +  if ((bigScale > 1.1) || (smallScale < 0.9)) {
245 +    sprintf( painCave.errMsg,
246 +             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
247 +             " Check your tauBarostat, as it is probably too small!\n\n"
248 +             " scaleMat = [%lf\t%lf\t%lf]\n"
249 +             "            [%lf\t%lf\t%lf]\n"
250 +             "            [%lf\t%lf\t%lf]\n",
251 +             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
252 +             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
253 +             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
254 +    painCave.isFatal = 1;
255 +    simError();
256 +  } else if (offDiagMax > 0.1) {
257 +    sprintf( painCave.errMsg,
258 +             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
259 +             " Check your tauBarostat, as it is probably too small!\n\n"
260 +             " scaleMat = [%lf\t%lf\t%lf]\n"
261 +             "            [%lf\t%lf\t%lf]\n"
262 +             "            [%lf\t%lf\t%lf]\n",
263 +             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
264 +             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
265 +             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
266 +    painCave.isFatal = 1;
267 +    simError();
268 +  } else {
269 +    info->getBoxM(hm);
270 +    info->matMul3(hm, scaleMat, hmnew);
271 +    info->setBoxM(hmnew);
272 +  }
273 +  
274   }
275  
276 < void NPTi::moveB( void ){
277 <  int i,j,k;
278 <  int atomIndex;
276 > template<typename T> void NPTf<T>::moveB( void ){
277 >
278 >  int i, j, k;
279    DirectionalAtom* dAtom;
280 <  double Tb[3];
281 <  double ji[3];
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  
192 –  instaTemp = tStats->getTemperature();
193 –  instaPress = tStats->getPressure();
194 –  instaVol = tStats->getVolume();
294  
295 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
296 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
295 >  // Set things up for the iteration:
296 >
297 >  oldChi = chi;
298    
299 +  for(i = 0; i < 3; i++)
300 +    for(j = 0; j < 3; j++)
301 +      oldEta[i][j] = eta[i][j];
302 +
303    for( i=0; i<nAtoms; i++ ){
304 <    atomIndex = i * 3;
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();
325 >  
326 >  for (k=0; k < 4; k++) {
327      
328 <    // velocity half step
329 <    for( j=atomIndex; j<(atomIndex+3); j++ )
330 <    for( j=atomIndex; j<(atomIndex+3); j++ )
331 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
332 <                       - vel[j]*(chi+eta));
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 <    if( atoms[i]->isDirectional() ){
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 <      dAtom = (DirectionalAtom *)atoms[i];
361 >      mass = atoms[i]->getMass();
362        
363 <      // get and convert the torque to body frame
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 <      Tb[0] = dAtom->getTx();
215 <      Tb[1] = dAtom->getTy();
216 <      Tb[2] = dAtom->getTz();
370 >      atoms[i]->setVel( vel );
371        
372 <      dAtom->lab2Body( Tb );
372 >      if( atoms[i]->isDirectional() ){
373 >
374 >        dAtom = (DirectionalAtom *)atoms[i];
375 >  
376 >        // 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 <      // get the angular momentum, and complete the angular momentum
385 <      // half step
222 <      
223 <      ji[0] = dAtom->getJx();
224 <      ji[1] = dAtom->getJy();
225 <      ji[2] = dAtom->getJz();
226 <      
227 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
228 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
229 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
230 <      
231 <      dAtom->setJx( ji[0] );
232 <      dAtom->setJy( ji[1] );
233 <      dAtom->setJz( ji[2] );
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 +  
401   }
402  
403 < int NPTi::readyCheck() {
403 > template<typename T> void NPTf<T>::resetIntegrator() {
404 >  int i,j;
405 >  
406 >  chi = 0.0;
407 >
408 >  for(i = 0; i < 3; i++)
409 >    for (j = 0; j < 3; j++)
410 >      eta[i][j] = 0.0;
411 >
412 > }
413 >
414 > template<typename T> int NPTf<T>::readyCheck() {
415 >
416 >  //check parent's readyCheck() first
417 >  if (T::readyCheck() == -1)
418 >    return -1;
419  
420    // First check to see if we have a target temperature.
421    // Not having one is fatal.
422    
423    if (!have_target_temp) {
424      sprintf( painCave.errMsg,
425 <             "NPTi error: You can't use the NPTi integrator\n"
425 >             "NPTf error: You can't use the NPTf integrator\n"
426               "   without a targetTemp!\n"
427               );
428      painCave.isFatal = 1;
# Line 252 | Line 432 | int NPTi::readyCheck() {
432  
433    if (!have_target_pressure) {
434      sprintf( painCave.errMsg,
435 <             "NPTi error: You can't use the NPTi integrator\n"
435 >             "NPTf error: You can't use the NPTf integrator\n"
436               "   without a targetPressure!\n"
437               );
438      painCave.isFatal = 1;
# Line 264 | Line 444 | int NPTi::readyCheck() {
444    
445    if (!have_tau_thermostat) {
446      sprintf( painCave.errMsg,
447 <             "NPTi error: If you use the NPTi\n"
447 >             "NPTf error: If you use the NPTf\n"
448               "   integrator, you must set tauThermostat.\n");
449      painCave.isFatal = 1;
450      simError();
# Line 275 | Line 455 | int NPTi::readyCheck() {
455    
456    if (!have_tau_barostat) {
457      sprintf( painCave.errMsg,
458 <             "NPTi error: If you use the NPTi\n"
458 >             "NPTf error: If you use the NPTf\n"
459               "   integrator, you must set tauBarostat.\n");
460      painCave.isFatal = 1;
461      simError();
# Line 284 | Line 464 | int NPTi::readyCheck() {
464  
465    // We need NkBT a lot, so just set it here:
466  
467 <  NkBT = (double)info->ndf * kB * targetTemp;
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);
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;
508 +
509 +  conservedQuantity = U+thermo+integral+PV+baro;
510 +  return conservedQuantity;
511 +  
512 + }

Diff Legend

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