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 578 by gezelter, Wed Jul 9 02:15:29 2003 UTC vs.
Revision 746 by mmeineke, Thu Sep 4 21:48:35 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 >
29 >  for(i = 0; i < 3; i++)
30 >    for (j = 0; j < 3; j++)
31 >      eta[i][j] = 0.0;
32 >
33    have_tau_thermostat = 0;
34    have_tau_barostat = 0;
35    have_target_temp = 0;
36    have_target_pressure = 0;
37   }
38  
39 < void NPTf::moveA() {
39 > template<typename T> void NPTf<T>::moveA() {
40    
41 <  int i,j,k;
37 <  int atomIndex, aMatIndex;
41 >  int i, j, k;
42    DirectionalAtom* dAtom;
43 <  double Tb[3];
44 <  double ji[3];
43 >  double Tb[3], ji[3];
44 >  double A[3][3], I[3][3];
45 >  double angle, mass;
46 >  double vel[3], pos[3], frc[3];
47 >
48    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;
49    double instaTemp, instaPress, instaVol;
50    double tt2, tb2;
51 <  double angle;
52 <  double press[9];
53 <  const double p_convert = 1.63882576e8;
51 >  double sc[3];
52 >  double eta2ij;
53 >  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
54 >  double bigScale, smallScale, offDiagMax;
55  
56    tt2 = tauThermostat * tauThermostat;
57    tb2 = tauBarostat * tauBarostat;
58  
59    instaTemp = tStats->getTemperature();
60    tStats->getPressureTensor(press);
57 –
58 –  for (i=0; i < 9; i++) press[i] *= p_convert;
59 –
61    instaVol = tStats->getVolume();
62    
63    // first evolve chi a half step
64    
65    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
66 <  
67 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
68 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
69 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
70 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
71 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
72 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
73 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
74 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
75 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
76 <  
66 >
67 >  for (i = 0; i < 3; i++ ) {
68 >    for (j = 0; j < 3; j++ ) {
69 >      if (i == j) {
70 >        
71 >        eta[i][j] += dt2 * instaVol *
72 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
73 >        
74 >        vScale[i][j] = eta[i][j] + chi;
75 >        
76 >      } else {
77 >        
78 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
79 >
80 >        vScale[i][j] = eta[i][j];
81 >        
82 >      }
83 >    }
84 >  }
85 >
86    for( i=0; i<nAtoms; i++ ){
87 <    atomIndex = i * 3;
88 <    aMatIndex = i * 9;
87 >
88 >    atoms[i]->getVel( vel );
89 >    atoms[i]->getPos( pos );
90 >    atoms[i]->getFrc( frc );
91 >
92 >    mass = atoms[i]->getMass();
93      
94      // velocity half step
95 +        
96 +    info->matVecMul3( vScale, vel, sc );
97      
98 <    vx = vel[atomIndex];
99 <    vy = vel[atomIndex+1];
100 <    vz = vel[atomIndex+2];
101 <    
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);
98 >    for (j = 0; j < 3; j++) {
99 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
100 >      rj[j] = pos[j];
101 >    }
102  
103 <    vel[atomIndex] = vx;
95 <    vel[atomIndex+1] = vy;
96 <    vel[atomIndex+2] = vz;
103 >    atoms[i]->setVel( vel );
104  
105      // position whole step    
106  
100 –    rj[0] = pos[atomIndex];
101 –    rj[1] = pos[atomIndex+1];
102 –    rj[2] = pos[atomIndex+2];
103 –
107      info->wrapVector(rj);
108  
109 <    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 >    info->matVecMul3( eta, rj, 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 >    for (j = 0; j < 3; j++ )
112 >      pos[j] += dt * (vel[j] + sc[j]);
113 >
114 >    atoms[i]->setPos( pos );
115    
116      if( atoms[i]->isDirectional() ){
117  
# Line 117 | Line 119 | void NPTf::moveA() {
119            
120        // get and convert the torque to body frame
121        
122 <      Tb[0] = dAtom->getTx();
121 <      Tb[1] = dAtom->getTy();
122 <      Tb[2] = dAtom->getTz();
123 <      
122 >      dAtom->getTrq( Tb );
123        dAtom->lab2Body( Tb );
124        
125        // get the angular momentum, and propagate a half step
126  
127 <      ji[0] = dAtom->getJx();
128 <      ji[1] = dAtom->getJy();
129 <      ji[2] = dAtom->getJz();
127 >      dAtom->getJ( ji );
128 >
129 >      for (j=0; j < 3; j++)
130 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
131        
132 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
133 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
134 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
135 –      
132        // use the angular velocities to propagate the rotation matrix a
133        // full time step
134 <      
134 >
135 >      dAtom->getA(A);
136 >      dAtom->getI(I);
137 >    
138        // rotate about the x-axis      
139 <      angle = dt2 * ji[0] / dAtom->getIxx();
140 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
141 <      
139 >      angle = dt2 * ji[0] / I[0][0];
140 >      this->rotate( 1, 2, angle, ji, A );
141 >
142        // rotate about the y-axis
143 <      angle = dt2 * ji[1] / dAtom->getIyy();
144 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
143 >      angle = dt2 * ji[1] / I[1][1];
144 >      this->rotate( 2, 0, angle, ji, A );
145        
146        // rotate about the z-axis
147 <      angle = dt * ji[2] / dAtom->getIzz();
148 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
147 >      angle = dt * ji[2] / I[2][2];
148 >      this->rotate( 0, 1, angle, ji, A);
149        
150        // rotate about the y-axis
151 <      angle = dt2 * ji[1] / dAtom->getIyy();
152 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
151 >      angle = dt2 * ji[1] / I[1][1];
152 >      this->rotate( 2, 0, angle, ji, A );
153        
154         // rotate about the x-axis
155 <      angle = dt2 * ji[0] / dAtom->getIxx();
156 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
155 >      angle = dt2 * ji[0] / I[0][0];
156 >      this->rotate( 1, 2, angle, ji, A );
157        
158 <      dAtom->setJx( ji[0] );
159 <      dAtom->setJy( ji[1] );
160 <      dAtom->setJz( ji[2] );
162 <    }
163 <    
158 >      dAtom->setJ( ji );
159 >      dAtom->setA( A  );    
160 >    }                    
161    }
162 <
162 >  
163    // Scale the box after all the positions have been moved:
164 <
164 >  
165    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
166    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
167 <
168 <
167 >  
168 >  bigScale = 1.0;
169 >  smallScale = 1.0;
170 >  offDiagMax = 0.0;
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;
174 >      
175 >      // Calculate the matrix Product of the eta array (we only need
176 >      // the ij element right now):
177 >      
178 >      eta2ij = 0.0;
179        for(k=0; k<3; k++){
180 <        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
180 >        eta2ij += eta[i][k] * eta[k][j];
181        }
182 +      
183 +      scaleMat[i][j] = 0.0;
184 +      // identity matrix (see above):
185 +      if (i == j) scaleMat[i][j] = 1.0;
186 +      // Taylor expansion for the exponential truncated at second order:
187 +      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
188 +
189 +      if (i != j)
190 +        if (fabs(scaleMat[i][j]) > offDiagMax)
191 +          offDiagMax = fabs(scaleMat[i][j]);
192 +      
193      }
181 –    ident[i][i] = 1.0;
182 –  }
194  
195 <  
196 <  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 <    }
195 >    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
196 >    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
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 <    }
199 >  if ((bigScale > 1.1) || (smallScale < 0.9)) {
200 >    sprintf( painCave.errMsg,
201 >             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
202 >             " Check your tauBarostat, as it is probably too small!\n\n"
203 >             " scaleMat = [%lf\t%lf\t%lf]\n"
204 >             "            [%lf\t%lf\t%lf]\n"
205 >             "            [%lf\t%lf\t%lf]\n",
206 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
207 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
208 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
209 >    painCave.isFatal = 1;
210 >    simError();
211 >  } else if (offDiagMax > 0.1) {
212 >    sprintf( painCave.errMsg,
213 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
214 >             " Check your tauBarostat, as it is probably too small!\n\n"
215 >             " scaleMat = [%lf\t%lf\t%lf]\n"
216 >             "            [%lf\t%lf\t%lf]\n"
217 >             "            [%lf\t%lf\t%lf]\n",
218 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
219 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
220 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
221 >    painCave.isFatal = 1;
222 >    simError();
223 >  } else {
224 >    info->getBoxM(hm);
225 >    info->matMul3(hm, scaleMat, hmnew);
226 >    info->setBoxM(hmnew);
227    }
205 –
206 –  info->setBoxM(hm);
228    
229   }
230  
231 < void NPTf::moveB( void ){
232 <  int i,j,k;
233 <  int atomIndex;
231 > template<typename T> void NPTf<T>::moveB( void ){
232 >
233 >  int i, j;
234    DirectionalAtom* dAtom;
235 <  double Tb[3];
236 <  double ji[3];
237 <  double press[9];
238 <  double instaTemp, instaVol;
235 >  double Tb[3], ji[3];
236 >  double vel[3], frc[3];
237 >  double mass;
238 >
239 >  double instaTemp, instaPress, instaVol;
240    double tt2, tb2;
241 <  double vx, vy, vz;
242 <  double scx, scy, scz;
221 <  const double p_convert = 1.63882576e8;
241 >  double sc[3];
242 >  double press[3][3], vScale[3][3];
243    
244    tt2 = tauThermostat * tauThermostat;
245    tb2 = tauBarostat * tauBarostat;
246  
247    instaTemp = tStats->getTemperature();
248    tStats->getPressureTensor(press);
228 –
229 –  for (i=0; i < 9; i++) press[i] *= p_convert;
230 –
249    instaVol = tStats->getVolume();
250    
251    // first evolve chi a half step
252    
253    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
254    
255 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
256 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
257 <  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);
255 >  for (i = 0; i < 3; i++ ) {
256 >    for (j = 0; j < 3; j++ ) {
257 >      if (i == j) {
258  
259 +        eta[i][j] += dt2 * instaVol *
260 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
261 +
262 +        vScale[i][j] = eta[i][j] + chi;
263 +        
264 +      } else {
265 +        
266 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
267 +
268 +        vScale[i][j] = eta[i][j];
269 +        
270 +      }
271 +    }
272 +  }
273 +
274    for( i=0; i<nAtoms; i++ ){
248 –    atomIndex = i * 3;
275  
276 +    atoms[i]->getVel( vel );
277 +    atoms[i]->getFrc( frc );
278 +
279 +    mass = atoms[i]->getMass();
280 +    
281      // velocity half step
282 +        
283 +    info->matVecMul3( vScale, vel, sc );
284      
285 <    vx = vel[atomIndex];
286 <    vy = vel[atomIndex+1];
287 <    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);
285 >    for (j = 0; j < 3; j++) {
286 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 >    }
288  
289 <    vel[atomIndex] = vx;
265 <    vel[atomIndex+1] = vy;
266 <    vel[atomIndex+2] = vz;
289 >    atoms[i]->setVel( vel );
290      
291      if( atoms[i]->isDirectional() ){
292 <      
292 >
293        dAtom = (DirectionalAtom *)atoms[i];
294 <      
294 >          
295        // get and convert the torque to body frame
296        
297 <      Tb[0] = dAtom->getTx();
275 <      Tb[1] = dAtom->getTy();
276 <      Tb[2] = dAtom->getTz();
277 <      
297 >      dAtom->getTrq( Tb );
298        dAtom->lab2Body( Tb );
299        
300 <      // get the angular momentum, and complete the angular momentum
281 <      // half step
300 >      // get the angular momentum, and propagate a half step
301        
302 <      ji[0] = dAtom->getJx();
284 <      ji[1] = dAtom->getJy();
285 <      ji[2] = dAtom->getJz();
302 >      dAtom->getJ( ji );
303        
304 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
305 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
289 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
304 >      for (j=0; j < 3; j++)
305 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
306        
307 <      dAtom->setJx( ji[0] );
308 <      dAtom->setJy( ji[1] );
309 <      dAtom->setJz( ji[2] );
294 <    }
307 >      dAtom->setJ( ji );
308 >
309 >    }                    
310    }
311   }
312  
313 < int NPTi::readyCheck() {
313 > template<typename T> void NPTf<T>::resetIntegrator() {
314 >  int i,j;
315 >  
316 >  chi = 0.0;
317 >
318 >  for(i = 0; i < 3; i++)
319 >    for (j = 0; j < 3; j++)
320 >      eta[i][j] = 0.0;
321 >
322 > }
323 >
324 > template<typename T> int NPTf<T>::readyCheck() {
325 >
326 >  //check parent's readyCheck() first
327 >  if (T::readyCheck() == -1)
328 >    return -1;
329  
330    // First check to see if we have a target temperature.
331    // Not having one is fatal.
332    
333    if (!have_target_temp) {
334      sprintf( painCave.errMsg,
335 <             "NPTi error: You can't use the NPTi integrator\n"
335 >             "NPTf error: You can't use the NPTf integrator\n"
336               "   without a targetTemp!\n"
337               );
338      painCave.isFatal = 1;
# Line 312 | Line 342 | int NPTi::readyCheck() {
342  
343    if (!have_target_pressure) {
344      sprintf( painCave.errMsg,
345 <             "NPTi error: You can't use the NPTi integrator\n"
345 >             "NPTf error: You can't use the NPTf integrator\n"
346               "   without a targetPressure!\n"
347               );
348      painCave.isFatal = 1;
# Line 324 | Line 354 | int NPTi::readyCheck() {
354    
355    if (!have_tau_thermostat) {
356      sprintf( painCave.errMsg,
357 <             "NPTi error: If you use the NPTi\n"
357 >             "NPTf error: If you use the NPTf\n"
358               "   integrator, you must set tauThermostat.\n");
359      painCave.isFatal = 1;
360      simError();
# Line 335 | Line 365 | int NPTi::readyCheck() {
365    
366    if (!have_tau_barostat) {
367      sprintf( painCave.errMsg,
368 <             "NPTi error: If you use the NPTi\n"
368 >             "NPTf error: If you use the NPTf\n"
369               "   integrator, you must set tauBarostat.\n");
370      painCave.isFatal = 1;
371      simError();

Diff Legend

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