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 576 by gezelter, Tue Jul 8 21:10:16 2003 UTC vs.
Revision 580 by gezelter, Wed Jul 9 13:56:36 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 19 | Line 19 | NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
19   //
20   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
21  
22 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
22 > NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
23    Integrator( theInfo, the_ff )
24   {
25    int i;
# Line 31 | Line 31 | void NPTi::moveA() {
31    have_target_pressure = 0;
32   }
33  
34 < void NPTi::moveA() {
34 > void NPTf::moveA() {
35    
36    int i,j,k;
37    int atomIndex, aMatIndex;
# Line 39 | Line 39 | void NPTi::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];
50 +  const double p_convert = 1.63882576e8;
51  
52    tt2 = tauThermostat * tauThermostat;
53    tb2 = tauBarostat * tauBarostat;
54  
55    instaTemp = tStats->getTemperature();
56 <  instaPress = tStats->getPressure();
56 >  tStats->getPressureTensor(press);
57 >
58 >  for (i=0; i < 9; i++) press[i] *= p_convert;
59 >
60    instaVol = tStats->getVolume();
61    
62    // first evolve chi a half step
63    
64    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
65    
66 <  for (i = 0; i < 9; i++) {
67 <    eta[i] += dt2 * ( instaVol * (sigma[i] - targetPressure*identMat[i]))
68 <      / (NkBT*tb2));
69 < }
70 <
66 >  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
67 >  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
68 >  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
69 >  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
70 >  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
71 >  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
72 >  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
73 >  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
74 >  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
75 >  
76    for( i=0; i<nAtoms; i++ ){
77      atomIndex = i * 3;
78      aMatIndex = i * 9;
79      
80      // velocity half step
81 <    for( j=atomIndex; j<(atomIndex+3); j++ )
82 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
83 <                       - vel[j]*(chi+eta));
81 >    
82 >    vx = vel[atomIndex];
83 >    vy = vel[atomIndex+1];
84 >    vz = vel[atomIndex+2];
85 >    
86 >    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
87 >    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
88 >    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
89 >    
90 >    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
91 >    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
92 >    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
93  
94 +    vel[atomIndex] = vx;
95 +    vel[atomIndex+1] = vy;
96 +    vel[atomIndex+2] = vz;
97 +
98      // position whole step    
99  
100 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
101 <      rj[0] = pos[j];
102 <      rj[1] = pos[j+1];
76 <      rj[2] = pos[j+2];
100 >    rj[0] = pos[atomIndex];
101 >    rj[1] = pos[atomIndex+1];
102 >    rj[2] = pos[atomIndex+2];
103  
104 <      info->wrapVector(rj);
104 >    info->wrapVector(rj);
105  
106 <      pos[j] += dt * (vel[j] + eta*rj[0]);
107 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
108 <      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
83 <    }
106 >    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
107 >    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
108 >    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
109  
110 <    // Scale the box after all the positions have been moved:
111 <
112 <    info->scaleBox(exp(dt*eta));
110 >    pos[atomIndex] += dt * (vel[atomIndex] + scx);
111 >    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
112 >    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
113    
114      if( atoms[i]->isDirectional() ){
115  
# Line 136 | Line 161 | void NPTi::moveA() {
161        dAtom->setJz( ji[2] );
162      }
163      
164 +  }
165 +
166 +  // Scale the box after all the positions have been moved:
167 +
168 +  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
169 +  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
170 +
171 +
172 +  for(i=0; i<3; i++){
173 +    for(j=0; j<3; j++){
174 +      ident[i][j] = 0.0;
175 +      eta1[i][j] = eta[3*i+j];
176 +      eta2[i][j] = 0.0;
177 +      for(k=0; k<3; k++){
178 +        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
179 +      }
180 +    }
181 +    ident[i][i] = 1.0;
182 +  }
183 +
184 +  
185 +  info->getBoxM(hm);
186 +
187 +  for(i=0; i<3; i++){
188 +    for(j=0; j<3; j++){      
189 +      hmnew[i][j] = 0.0;
190 +      for(k=0; k<3; k++){
191 +        // remember that hmat has transpose ordering for Fortran compat:
192 +        hmnew[i][j] += hm[3*k+i] * (ident[k][j]
193 +                                    + dt * eta1[k][j]
194 +                                    + 0.5 * dt * dt * eta2[k][j]);
195 +      }
196 +    }
197    }
198 +  
199 +  for (i = 0; i < 3; i++) {
200 +    for (j = 0; j < 3; j++) {
201 +      // remember that hmat has transpose ordering for Fortran compat:
202 +      hm[3*j + 1] = hmnew[i][j];
203 +    }
204 +  }
205 +
206 +  info->setBoxM(hm);
207 +  
208   }
209  
210 < void NPTi::moveB( void ){
210 > void NPTf::moveB( void ){
211    int i,j,k;
212    int atomIndex;
213    DirectionalAtom* dAtom;
214    double Tb[3];
215    double ji[3];
216 <  double instaTemp, instaPress, instaVol;
216 >  double press[9];
217 >  double instaTemp, instaVol;
218    double tt2, tb2;
219 +  double vx, vy, vz;
220 +  double scx, scy, scz;
221 +  const double p_convert = 1.63882576e8;
222    
223    tt2 = tauThermostat * tauThermostat;
224    tb2 = tauBarostat * tauBarostat;
225  
226    instaTemp = tStats->getTemperature();
227 <  instaPress = tStats->getPressure();
156 <  instaVol = tStats->getVolume();
227 >  tStats->getPressureTensor(press);
228  
229 +  for (i=0; i < 9; i++) press[i] *= p_convert;
230 +
231 +  instaVol = tStats->getVolume();
232 +  
233 +  // first evolve chi a half step
234 +  
235    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
159 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
236    
237 +  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
238 +  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
239 +  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
240 +  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
241 +  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
242 +  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
243 +  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
244 +  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
245 +  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
246 +
247    for( i=0; i<nAtoms; i++ ){
248      atomIndex = i * 3;
249 <    
249 >
250      // velocity half step
165 –    for( j=atomIndex; j<(atomIndex+3); j++ )
166 –    for( j=atomIndex; j<(atomIndex+3); j++ )
167 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168 –                       - vel[j]*(chi+eta));
251      
252 +    vx = vel[atomIndex];
253 +    vy = vel[atomIndex+1];
254 +    vz = vel[atomIndex+2];
255 +    
256 +    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
257 +    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
258 +    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
259 +    
260 +    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
261 +    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
262 +    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
263 +
264 +    vel[atomIndex] = vx;
265 +    vel[atomIndex+1] = vy;
266 +    vel[atomIndex+2] = vz;
267 +    
268      if( atoms[i]->isDirectional() ){
269        
270        dAtom = (DirectionalAtom *)atoms[i];
# Line 197 | Line 295 | int NPTi::readyCheck() {
295    }
296   }
297  
298 < int NPTi::readyCheck() {
298 > int NPTf::readyCheck() {
299  
300    // First check to see if we have a target temperature.
301    // Not having one is fatal.
302    
303    if (!have_target_temp) {
304      sprintf( painCave.errMsg,
305 <             "NPTi error: You can't use the NPTi integrator\n"
305 >             "NPTf error: You can't use the NPTf integrator\n"
306               "   without a targetTemp!\n"
307               );
308      painCave.isFatal = 1;
# Line 214 | Line 312 | int NPTi::readyCheck() {
312  
313    if (!have_target_pressure) {
314      sprintf( painCave.errMsg,
315 <             "NPTi error: You can't use the NPTi integrator\n"
315 >             "NPTf error: You can't use the NPTf integrator\n"
316               "   without a targetPressure!\n"
317               );
318      painCave.isFatal = 1;
# Line 226 | Line 324 | int NPTi::readyCheck() {
324    
325    if (!have_tau_thermostat) {
326      sprintf( painCave.errMsg,
327 <             "NPTi error: If you use the NPTi\n"
327 >             "NPTf error: If you use the NPTf\n"
328               "   integrator, you must set tauThermostat.\n");
329      painCave.isFatal = 1;
330      simError();
# Line 237 | Line 335 | int NPTi::readyCheck() {
335    
336    if (!have_tau_barostat) {
337      sprintf( painCave.errMsg,
338 <             "NPTi error: If you use the NPTi\n"
338 >             "NPTf error: If you use the NPTf\n"
339               "   integrator, you must set tauBarostat.\n");
340      painCave.isFatal = 1;
341      simError();

Diff Legend

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