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 590 by mmeineke, Thu Jul 10 22:15:53 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 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, j;
27    chi = 0.0;
28 +  integralOfChidt = 0.0;
29  
30    for(i = 0; i < 3; i++)
31      for (j = 0; j < 3; j++)
# Line 33 | Line 35 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
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;
41 <  int atomIndex, aMatIndex;
61 >  int i, j, k;
62    DirectionalAtom* dAtom;
63 <  double Tb[3];
64 <  double ji[3];
65 <  double ri[3], vi[3], sc[3];
66 <  double instaTemp, instaVol;
67 <  double tt2, tb2, eta2ij;
68 <  double angle;
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 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;
# Line 54 | Line 80 | void NPTf::moveA() {
80    instaTemp = tStats->getTemperature();
81    tStats->getPressureTensor(press);
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 <
93 <        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 <      }
92 >        vScale[i][j] += chi;          
93 >      }              
94      }
95    }
96 <
96 >  
97 >  //evolve velocity half step
98    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]);
99  
100 <    vel[atomIndex]   = vi[0];
101 <    vel[atomIndex+1] = vi[1];
99 <    vel[atomIndex+2] = vi[2];
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 <    ri[0] = pos[atomIndex];
108 <    ri[1] = pos[atomIndex+1];
109 <    ri[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(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]);
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();
122 <      Tb[1] = dAtom->getTy();
123 <      Tb[2] = dAtom->getTz();
124 <      
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        
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 –      
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    }
184 +  
185 +  //the first estimation of r(t+dt) is equal to  r(t)
186 +    
187 +  for(k = 0; k < 4; k ++){
188  
189 <  // Scale the box after all the positions have been moved:
189 >    for(i =0 ; i < nAtoms; i++){
190  
191 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
192 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
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 +
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 <
220 >      
221        // Calculate the matrix Product of the eta array (we only need
222        // the ij element right now):
223 <
223 >      
224        eta2ij = 0.0;
225        for(k=0; k<3; k++){
226          eta2ij += eta[i][k] * eta[k][j];
# Line 187 | Line 232 | void NPTf::moveA() {
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    }
192 –  
193 –  info->getBoxM(hm);
194 –  info->matMul3(hm, scaleMat, hmnew);
195 –  info->setBoxM(hmnew);
243    
244 < }
245 <
246 < void NPTf::moveB( void ){
247 <  int i,j, k;
248 <  int atomIndex;
249 <  DirectionalAtom* dAtom;
250 <  double Tb[3];
251 <  double ji[3];
252 <  double vi[3], sc[3];
253 <  double instaTemp, instaVol;
254 <  double tt2, tb2;
255 <  double press[3][3], vScale[3][3];
256 <  
257 <  tt2 = tauThermostat * tauThermostat;
258 <  tb2 = tauBarostat * tauBarostat;
259 <
260 <  instaTemp = tStats->getTemperature();
261 <  tStats->getPressureTensor(press);
262 <  instaVol = tStats->getVolume();
263 <  
264 <  // first evolve chi a half step
265 <  
266 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
267 <  
268 <  for (i = 0; i < 3; i++ ) {
269 <    for (j = 0; j < 3; j++ ) {
270 <      if (i == j) {
271 <
272 <        eta[i][j] += dt2 * instaVol *
273 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
274 <
228 <        vScale[i][j] = eta[i][j] + chi;
229 <        
230 <      } else {
231 <        
232 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
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 <        vScale[i][j] = eta[i][j];
277 <        
278 <      }
276 > template<typename T> void NPTf<T>::moveB( void ){
277 >
278 >  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 >  
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 >
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 <  for( i=0; i<nAtoms; i++ ){
241 <    atomIndex = i * 3;
322 >  // do the iteration:
323  
324 <    // velocity half step
324 >  instaVol = tStats->getVolume();
325 >  
326 >  for (k=0; k < 4; k++) {
327      
328 <    vi[0] = vel[atomIndex];
329 <    vi[1] = vel[atomIndex+1];
247 <    vi[2] = vel[atomIndex+2];
248 <    
249 <    info->matVecMul3( vScale, vi, sc );
250 <    
251 <    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]);
328 >    instaTemp = tStats->getTemperature();
329 >    tStats->getPressureTensor(press);
330  
331 <    vel[atomIndex]   = vi[0];
332 <    vel[atomIndex+1] = vi[1];
333 <    vel[atomIndex+2] = vi[2];
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();
266 <      Tb[1] = dAtom->getTy();
267 <      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
273 <      
274 <      ji[0] = dAtom->getJx();
275 <      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] );
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 NPTf::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.
# Line 335 | Line 464 | int NPTf::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