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 590 by mmeineke, Thu Jul 10 22:15:53 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 22 | Line 22 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
22   NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
23    Integrator( theInfo, the_ff )
24   {
25 <  int i;
25 >  int i, j;
26    chi = 0.0;
27 <  for(i = 0; i < 9; i++) eta[i] = 0.0;
27 >
28 >  for(i = 0; i < 3; i++)
29 >    for (j = 0; j < 3; j++)
30 >      eta[i][j] = 0.0;
31 >
32    have_tau_thermostat = 0;
33    have_tau_barostat = 0;
34    have_target_temp = 0;
# Line 38 | Line 42 | void NPTf::moveA() {
42    DirectionalAtom* dAtom;
43    double Tb[3];
44    double ji[3];
45 <  double rj[3];
46 <  double instaTemp, instaPress, instaVol;
47 <  double tt2, tb2;
45 >  double ri[3], vi[3], sc[3];
46 >  double instaTemp, instaVol;
47 >  double tt2, tb2, eta2ij;
48    double angle;
49 <  double press[9];
46 <  const double p_convert = 1.63882576e8;
49 >  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
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);
63 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
64 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
65 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
66 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
67 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
68 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
69 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
70 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
71 <  
61 >
62 >  for (i = 0; i < 3; i++ ) {
63 >    for (j = 0; j < 3; j++ ) {
64 >      if (i == j) {
65 >
66 >        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 >      }
78 >    }
79 >  }
80 >
81    for( i=0; i<nAtoms; i++ ){
82      atomIndex = i * 3;
83      aMatIndex = i * 9;
84      
85      // velocity half step
86      
87 <    vx = vel[atomIndex];
88 <    vy = vel[atomIndex+1];
89 <    vz = vel[atomIndex+2];
87 >    vi[0] = vel[atomIndex];
88 >    vi[1] = vel[atomIndex+1];
89 >    vi[2] = vel[atomIndex+2];
90      
91 <    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
83 <    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
84 <    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
91 >    info->matVecMul3( vScale, vi, sc );
92      
93 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
94 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
95 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
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]);
96  
97 <    vel[atomIndex] = vx;
98 <    vel[atomIndex+1] = vy;
99 <    vel[atomIndex+2] = vz;
97 >    vel[atomIndex]   = vi[0];
98 >    vel[atomIndex+1] = vi[1];
99 >    vel[atomIndex+2] = vi[2];
100  
101      // position whole step    
102  
103 <    rj[0] = pos[atomIndex];
104 <    rj[1] = pos[atomIndex+1];
105 <    rj[2] = pos[atomIndex+2];
103 >    ri[0] = pos[atomIndex];
104 >    ri[1] = pos[atomIndex+1];
105 >    ri[2] = pos[atomIndex+2];
106  
107 <    info->wrapVector(rj);
107 >    info->wrapVector(ri);
108  
109 <    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
103 <    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
104 <    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
109 >    info->matVecMul3( eta, ri, sc );
110  
111 <    pos[atomIndex] += dt * (vel[atomIndex] + scx);
112 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
113 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
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]);
114    
115      if( atoms[i]->isDirectional() ){
116  
# Line 160 | Line 165 | void NPTf::moveA() {
165    }
166  
167    // Scale the box after all the positions have been moved:
163 –  
168  
169 <
170 <  // Use a taylor expansion for eta products
167 <  
168 <  info->getBoxM(hm);
169 <  
169 >  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
170 >  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
171  
172  
173 +  for(i=0; i<3; i++){
174 +    for(j=0; j<3; j++){
175  
176 +      // Calculate the matrix Product of the eta array (we only need
177 +      // the ij element right now):
178  
179 +      eta2ij = 0.0;
180 +      for(k=0; k<3; k++){
181 +        eta2ij += eta[i][k] * eta[k][j];
182 +      }
183 +      
184 +      scaleMat[i][j] = 0.0;
185 +      // identity matrix (see above):
186 +      if (i == j) scaleMat[i][j] = 1.0;
187 +      // Taylor expansion for the exponential truncated at second order:
188 +      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
189  
190 <   info->scaleBox(exp(dt*eta));
191 <
192 <
190 >    }
191 >  }
192 >  
193 >  info->getBoxM(hm);
194 >  info->matMul3(hm, scaleMat, hmnew);
195 >  info->setBoxM(hmnew);
196 >  
197   }
198  
199 < void NPTi::moveB( void ){
200 <  int i,j,k;
199 > void NPTf::moveB( void ){
200 >  int i,j, k;
201    int atomIndex;
202    DirectionalAtom* dAtom;
203    double Tb[3];
204    double ji[3];
205 <  double instaTemp, instaPress, instaVol;
205 >  double vi[3], sc[3];
206 >  double instaTemp, instaVol;
207    double tt2, tb2;
208 +  double press[3][3], vScale[3][3];
209    
210    tt2 = tauThermostat * tauThermostat;
211    tb2 = tauBarostat * tauBarostat;
212  
213    instaTemp = tStats->getTemperature();
214 <  instaPress = tStats->getPressure();
214 >  tStats->getPressureTensor(press);
215    instaVol = tStats->getVolume();
216 <
216 >  
217 >  // first evolve chi a half step
218 >  
219    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
197 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
220    
221 +  for (i = 0; i < 3; i++ ) {
222 +    for (j = 0; j < 3; j++ ) {
223 +      if (i == j) {
224 +
225 +        eta[i][j] += dt2 * instaVol *
226 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
227 +
228 +        vScale[i][j] = eta[i][j] + chi;
229 +        
230 +      } else {
231 +        
232 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
233 +
234 +        vScale[i][j] = eta[i][j];
235 +        
236 +      }
237 +    }
238 +  }
239 +
240    for( i=0; i<nAtoms; i++ ){
241      atomIndex = i * 3;
242 <    
242 >
243      // 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));
244      
245 +    vi[0] = vel[atomIndex];
246 +    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]);
254 +
255 +    vel[atomIndex]   = vi[0];
256 +    vel[atomIndex+1] = vi[1];
257 +    vel[atomIndex+2] = vi[2];
258 +    
259      if( atoms[i]->isDirectional() ){
260        
261        dAtom = (DirectionalAtom *)atoms[i];
# Line 235 | Line 286 | int NPTi::readyCheck() {
286    }
287   }
288  
289 < int NPTi::readyCheck() {
289 > int NPTf::readyCheck() {
290  
291    // First check to see if we have a target temperature.
292    // Not having one is fatal.
293    
294    if (!have_target_temp) {
295      sprintf( painCave.errMsg,
296 <             "NPTi error: You can't use the NPTi integrator\n"
296 >             "NPTf error: You can't use the NPTf integrator\n"
297               "   without a targetTemp!\n"
298               );
299      painCave.isFatal = 1;
# Line 252 | Line 303 | int NPTi::readyCheck() {
303  
304    if (!have_target_pressure) {
305      sprintf( painCave.errMsg,
306 <             "NPTi error: You can't use the NPTi integrator\n"
306 >             "NPTf error: You can't use the NPTf integrator\n"
307               "   without a targetPressure!\n"
308               );
309      painCave.isFatal = 1;
# Line 264 | Line 315 | int NPTi::readyCheck() {
315    
316    if (!have_tau_thermostat) {
317      sprintf( painCave.errMsg,
318 <             "NPTi error: If you use the NPTi\n"
318 >             "NPTf error: If you use the NPTf\n"
319               "   integrator, you must set tauThermostat.\n");
320      painCave.isFatal = 1;
321      simError();
# Line 275 | Line 326 | int NPTi::readyCheck() {
326    
327    if (!have_tau_barostat) {
328      sprintf( painCave.errMsg,
329 <             "NPTi error: If you use the NPTi\n"
329 >             "NPTf error: If you use the NPTf\n"
330               "   integrator, you must set tauBarostat.\n");
331      painCave.isFatal = 1;
332      simError();

Diff Legend

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