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 586 by mmeineke, Wed Jul 9 22:14:06 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;
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];
42 –  double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
43 –  double hm[9];
44 –  double vx, vy, vz;
45 –  double scx, scy, scz;
69    double instaTemp, instaPress, instaVol;
70    double tt2, tb2;
71 <  double angle;
72 <  double press[9];
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 >        vScale[i][j] += chi;          
93 >      }              
94 >    }
95 >  }
96    
97 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
63 <    (NkBT*tb2);
64 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
65 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
66 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
67 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
68 <    (NkBT*tb2);
69 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
70 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
71 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
72 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
73 <    (NkBT*tb2);
74 <  
97 >  //evolve velocity half step
98    for( i=0; i<nAtoms; i++ ){
76 –    atomIndex = i * 3;
77 –    aMatIndex = i * 9;
78 –    
79 –    // velocity half step
80 –    
81 –    vx = vel[atomIndex];
82 –    vy = vel[atomIndex+1];
83 –    vz = vel[atomIndex+2];
84 –    
85 –    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
86 –    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
87 –    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
88 –    
89 –    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
90 –    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
91 –    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
99  
100 <    vel[atomIndex] = vx;
101 <    vel[atomIndex+1] = vy;
95 <    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);
104 <
105 <    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
106 <    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
107 <    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
108 <
109 <    pos[atomIndex] += dt * (vel[atomIndex] + scx);
110 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
111 <    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();
120 <      Tb[1] = dAtom->getTy();
121 <      Tb[2] = dAtom->getTz();
122 <      
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        
131 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
132 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
133 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
134 –      
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 <  for(i=0; i<3; i++){
203 <    for(j=0; j<3; j++){
173 <      ident[i][j] = 0.0;
174 <      eta1[i][j] = eta[3*i+j];
175 <      eta2[i][j] = 0.0;
176 <      for(k=0; k<3; k++){
177 <        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
178 <      }
202 >      atoms[i]->setPos( pos );
203 >
204      }
180 –    ident[i][i] = 1.0;
181 –  }
205  
206 <  
207 <  info->getBoxM(hm);
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 <      hmnew[i][j] = 0.0;
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 <        // remember that hmat has transpose ordering for Fortran compat:
191 <        hmnew[i][j] += hm[3*k+i] * (ident[k][j]
192 <                                    + dt * eta1[k][j]
193 <                                    + 0.5 * dt * dt * eta2[k][j]);
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 <  for (i = 0; i < 3; i++) {
245 <    for (j = 0; j < 3; j++) {
246 <      // remember that hmat has transpose ordering for Fortran compat:
247 <      hm[3*j + i] = hmnew[i][j];
248 <    }
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    }
204 –
205 –  info->setBoxM(hm);
273    
274   }
275  
276 < void NPTf::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];
282 <  double press[9];
283 <  double instaTemp, instaVol;
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 vx, vy, vz;
287 <  double scx, scy, scz;
288 <  const double p_convert = 1.63882576e8;
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 <  instaTemp = tStats->getTemperature();
295 <  tStats->getPressureTensor(press);
296 <  instaVol = tStats->getVolume();
297 <  
229 <  // first evolve chi a half step
294 >
295 >  // Set things up for the iteration:
296 >
297 >  oldChi = chi;
298    
299 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
300 <  
301 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
234 <    (NkBT*tb2);
235 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
236 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
237 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
238 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
239 <    (NkBT*tb2);
240 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
241 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
242 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
243 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
244 <    (NkBT*tb2);
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++ ){
247 –    atomIndex = i * 3;
304  
305 <    // velocity half step
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 <    vx = vel[atomIndex];
329 <    vy = vel[atomIndex+1];
253 <    vz = vel[atomIndex+2];
254 <    
255 <    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
256 <    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
257 <    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
258 <    
259 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
260 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
261 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
328 >    instaTemp = tStats->getTemperature();
329 >    tStats->getPressureTensor(press);
330  
331 <    vel[atomIndex] = vx;
332 <    vel[atomIndex+1] = vy;
333 <    vel[atomIndex+2] = vz;
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();
274 <      Tb[1] = dAtom->getTy();
275 <      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
281 <      
282 <      ji[0] = dAtom->getJx();
283 <      ji[1] = dAtom->getJy();
284 <      ji[2] = dAtom->getJz();
285 <      
286 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
287 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
288 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
289 <      
290 <      dAtom->setJx( ji[0] );
291 <      dAtom->setJy( ji[1] );
292 <      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 343 | 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