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 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC

# Line 9 | Line 9
9   #include "simError.h"
10  
11  
12 < // Basic isotropic thermostating and barostating via the Melchionna
12 > // Basic non-isotropic thermostating and barostating via the Melchionna
13   // modification of the Hoover algorithm:
14   //
15   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
# Line 39 | Line 39 | void NPTf::moveA() {
39    double Tb[3];
40    double ji[3];
41    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;
46    double instaTemp, instaPress, instaVol;
47    double tt2, tb2;
48    double angle;
49    double press[9];
46 –  const double p_convert = 1.63882576e8;
50  
51    tt2 = tauThermostat * tauThermostat;
52    tb2 = tauBarostat * tauBarostat;
53  
54    instaTemp = tStats->getTemperature();
55    tStats->getPressureTensor(press);
53 –
54 –  for (i=0; i < 9; i++) press[i] *= p_convert;
55 –
56    instaVol = tStats->getVolume();
57    
58    // first evolve chi a half step
59    
60    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
61    
62 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
62 >  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) / (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) / (NkBT*tb2);
72 >  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
73 >    (NkBT*tb2);
74    
75    for( i=0; i<nAtoms; i++ ){
76      atomIndex = i * 3;
# Line 160 | Line 163 | void NPTf::moveA() {
163    }
164  
165    // Scale the box after all the positions have been moved:
163 –  
166  
167 +  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
168 +  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
169  
170 <  // Use a taylor expansion for eta products
170 >
171 >  for(i=0; i<3; i++){
172 >    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 >      }
179 >    }
180 >    ident[i][i] = 1.0;
181 >  }
182 >
183    
184    info->getBoxM(hm);
185 +
186 +  for(i=0; i<3; i++){
187 +    for(j=0; j<3; j++){      
188 +      hmnew[i][j] = 0.0;
189 +      for(k=0; k<3; k++){
190 +        // 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]);
194 +      }
195 +    }
196 +  }
197    
198 +  for (i = 0; i < 3; i++) {
199 +    for (j = 0; j < 3; j++) {
200 +      // remember that hmat has transpose ordering for Fortran compat:
201 +      hm[3*j + i] = hmnew[i][j];
202 +    }
203 +  }
204  
205 <
206 <
173 <
174 <
175 <   info->scaleBox(exp(dt*eta));
176 <
177 <
205 >  info->setBoxM(hm);
206 >  
207   }
208  
209 < void NPTi::moveB( void ){
209 > void NPTf::moveB( void ){
210    int i,j,k;
211    int atomIndex;
212    DirectionalAtom* dAtom;
213    double Tb[3];
214    double ji[3];
215 <  double instaTemp, instaPress, instaVol;
215 >  double press[9];
216 >  double instaTemp, instaVol;
217    double tt2, tb2;
218 +  double vx, vy, vz;
219 +  double scx, scy, scz;
220 +  const double p_convert = 1.63882576e8;
221    
222    tt2 = tauThermostat * tauThermostat;
223    tb2 = tauBarostat * tauBarostat;
224  
225    instaTemp = tStats->getTemperature();
226 <  instaPress = tStats->getPressure();
226 >  tStats->getPressureTensor(press);
227    instaVol = tStats->getVolume();
228 <
228 >  
229 >  // first evolve chi a half step
230 >  
231    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
197 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
232    
233 +  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);
245 +
246    for( i=0; i<nAtoms; i++ ){
247      atomIndex = i * 3;
248 <    
248 >
249      // velocity half step
203 –    for( j=atomIndex; j<(atomIndex+3); j++ )
204 –    for( j=atomIndex; j<(atomIndex+3); j++ )
205 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
206 –                       - vel[j]*(chi+eta));
250      
251 +    vx = vel[atomIndex];
252 +    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);
262 +
263 +    vel[atomIndex] = vx;
264 +    vel[atomIndex+1] = vy;
265 +    vel[atomIndex+2] = vz;
266 +    
267      if( atoms[i]->isDirectional() ){
268        
269        dAtom = (DirectionalAtom *)atoms[i];
# Line 235 | Line 294 | int NPTi::readyCheck() {
294    }
295   }
296  
297 < int NPTi::readyCheck() {
297 > int NPTf::readyCheck() {
298  
299    // First check to see if we have a target temperature.
300    // Not having one is fatal.
301    
302    if (!have_target_temp) {
303      sprintf( painCave.errMsg,
304 <             "NPTi error: You can't use the NPTi integrator\n"
304 >             "NPTf error: You can't use the NPTf integrator\n"
305               "   without a targetTemp!\n"
306               );
307      painCave.isFatal = 1;
# Line 252 | Line 311 | int NPTi::readyCheck() {
311  
312    if (!have_target_pressure) {
313      sprintf( painCave.errMsg,
314 <             "NPTi error: You can't use the NPTi integrator\n"
314 >             "NPTf error: You can't use the NPTf integrator\n"
315               "   without a targetPressure!\n"
316               );
317      painCave.isFatal = 1;
# Line 264 | Line 323 | int NPTi::readyCheck() {
323    
324    if (!have_tau_thermostat) {
325      sprintf( painCave.errMsg,
326 <             "NPTi error: If you use the NPTi\n"
326 >             "NPTf error: If you use the NPTf\n"
327               "   integrator, you must set tauThermostat.\n");
328      painCave.isFatal = 1;
329      simError();
# Line 275 | Line 334 | int NPTi::readyCheck() {
334    
335    if (!have_tau_barostat) {
336      sprintf( painCave.errMsg,
337 <             "NPTi error: If you use the NPTi\n"
337 >             "NPTf error: If you use the NPTf\n"
338               "   integrator, you must set tauBarostat.\n");
339      painCave.isFatal = 1;
340      simError();

Diff Legend

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