ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 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 19 | Line 19
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 | NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
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  
51    tt2 = tauThermostat * tauThermostat;
52    tb2 = tauBarostat * tauBarostat;
53  
54    instaTemp = tStats->getTemperature();
55 <  instaPress = tStats->getPressure();
55 >  tStats->getPressureTensor(press);
56    instaVol = tStats->getVolume();
57    
58    // first evolve chi a half step
59    
60    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
61    
62 <  for (i = 0; i < 9; i++) {
63 <    eta[i] += dt2 * ( instaVol * (sigma[i] - targetPressure*identMat[i]))
64 <      / (NkBT*tb2));
65 < }
66 <
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/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 >  
75    for( i=0; i<nAtoms; i++ ){
76      atomIndex = i * 3;
77      aMatIndex = i * 9;
78      
79      // velocity half step
80 <    for( j=atomIndex; j<(atomIndex+3); j++ )
81 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
82 <                       - vel[j]*(chi+eta));
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);
92  
93 +    vel[atomIndex] = vx;
94 +    vel[atomIndex+1] = vy;
95 +    vel[atomIndex+2] = vz;
96 +
97      // position whole step    
98  
99 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
100 <      rj[0] = pos[j];
101 <      rj[1] = pos[j+1];
76 <      rj[2] = pos[j+2];
99 >    rj[0] = pos[atomIndex];
100 >    rj[1] = pos[atomIndex+1];
101 >    rj[2] = pos[atomIndex+2];
102  
103 <      info->wrapVector(rj);
103 >    info->wrapVector(rj);
104  
105 <      pos[j] += dt * (vel[j] + eta*rj[0]);
106 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
107 <      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
83 <    }
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 <    // Scale the box after all the positions have been moved:
110 <
111 <    info->scaleBox(exp(dt*eta));
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);
112    
113      if( atoms[i]->isDirectional() ){
114  
# Line 137 | Line 161 | void NPTi::moveA() {
161      }
162      
163    }
164 +
165 +  // Scale the box after all the positions have been moved:
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 +
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 +  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;
159 –  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
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));
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 197 | Line 294 | void NPTi::moveB( void ){
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 214 | 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 226 | 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 237 | 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