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 763 by tim, Mon Sep 15 16:52:02 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 9 | Line 10
10   #include "simError.h"
11  
12  
13 < // Basic isotropic thermostating and barostating via the Melchionna
13 > // Basic non-isotropic thermostating and barostating via the Melchionna
14   // modification of the Hoover algorithm:
15   //
16   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
# 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 >  integralOfChidt = 0.0;
29 >
30 >  for(i = 0; i < 3; i++)
31 >    for (j = 0; j < 3; j++)
32 >      eta[i][j] = 0.0;
33 >
34    have_tau_thermostat = 0;
35    have_tau_barostat = 0;
36    have_target_temp = 0;
37    have_target_pressure = 0;
38   }
39  
40 < void NPTf::moveA() {
40 > template<typename T> void NPTf<T>::moveA() {
41    
42 <  int i,j,k;
37 <  int atomIndex, aMatIndex;
42 >  int i, j, k;
43    DirectionalAtom* dAtom;
44 <  double Tb[3];
45 <  double ji[3];
44 >  double Tb[3], ji[3];
45 >  double A[3][3], I[3][3];
46 >  double angle, mass;
47 >  double vel[3], pos[3], frc[3];
48 >
49    double rj[3];
50    double instaTemp, instaPress, instaVol;
51    double tt2, tb2;
52 <  double angle;
53 <  double press[9];
54 <  const double p_convert = 1.63882576e8;
52 >  double sc[3];
53 >  double eta2ij;
54 >  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
55 >  double bigScale, smallScale, offDiagMax;
56  
57    tt2 = tauThermostat * tauThermostat;
58    tb2 = tauBarostat * tauBarostat;
59  
60    instaTemp = tStats->getTemperature();
61    tStats->getPressureTensor(press);
53 –
54 –  for (i=0; i < 9; i++) press[i] *= p_convert;
55 –
62    instaVol = tStats->getVolume();
63    
64    // first evolve chi a half step
65    
66    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
67 <  
68 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
69 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
70 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
71 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
72 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
73 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
74 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
75 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
76 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
77 <  
67 >
68 >  for (i = 0; i < 3; i++ ) {
69 >    for (j = 0; j < 3; j++ ) {
70 >      if (i == j) {
71 >        
72 >        eta[i][j] += dt2 * instaVol *
73 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
74 >        
75 >        vScale[i][j] = eta[i][j] + chi;
76 >          
77 >      } else {
78 >        
79 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
80 >
81 >        vScale[i][j] = eta[i][j];
82 >        
83 >      }
84 >    }
85 >  }
86 >
87    for( i=0; i<nAtoms; i++ ){
88 <    atomIndex = i * 3;
89 <    aMatIndex = i * 9;
88 >
89 >    atoms[i]->getVel( vel );
90 >    atoms[i]->getPos( pos );
91 >    atoms[i]->getFrc( frc );
92 >
93 >    mass = atoms[i]->getMass();
94      
95      // velocity half step
96 +        
97 +    info->matVecMul3( vScale, vel, sc );
98      
99 <    vx = vel[atomIndex];
100 <    vy = vel[atomIndex+1];
101 <    vz = vel[atomIndex+2];
102 <    
82 <    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;
85 <    
86 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
87 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
88 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
99 >    for (j = 0; j < 3; j++) {
100 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
101 >      rj[j] = pos[j];
102 >    }
103  
104 <    vel[atomIndex] = vx;
91 <    vel[atomIndex+1] = vy;
92 <    vel[atomIndex+2] = vz;
104 >    atoms[i]->setVel( vel );
105  
106      // position whole step    
107  
96 –    rj[0] = pos[atomIndex];
97 –    rj[1] = pos[atomIndex+1];
98 –    rj[2] = pos[atomIndex+2];
99 –
108      info->wrapVector(rj);
109  
110 <    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];
110 >    info->matVecMul3( eta, rj, sc );
111  
112 <    pos[atomIndex] += dt * (vel[atomIndex] + scx);
113 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
114 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
112 >    for (j = 0; j < 3; j++ )
113 >      pos[j] += dt * (vel[j] + sc[j]);
114 >
115 >    atoms[i]->setPos( pos );
116    
117      if( atoms[i]->isDirectional() ){
118  
# Line 113 | Line 120 | void NPTf::moveA() {
120            
121        // get and convert the torque to body frame
122        
123 <      Tb[0] = dAtom->getTx();
117 <      Tb[1] = dAtom->getTy();
118 <      Tb[2] = dAtom->getTz();
119 <      
123 >      dAtom->getTrq( Tb );
124        dAtom->lab2Body( Tb );
125        
126        // get the angular momentum, and propagate a half step
127  
128 <      ji[0] = dAtom->getJx();
129 <      ji[1] = dAtom->getJy();
130 <      ji[2] = dAtom->getJz();
128 >      dAtom->getJ( ji );
129 >
130 >      for (j=0; j < 3; j++)
131 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
132        
128 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
129 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
130 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
131 –      
133        // use the angular velocities to propagate the rotation matrix a
134        // full time step
135 <      
135 >
136 >      dAtom->getA(A);
137 >      dAtom->getI(I);
138 >    
139        // rotate about the x-axis      
140 <      angle = dt2 * ji[0] / dAtom->getIxx();
141 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
142 <      
140 >      angle = dt2 * ji[0] / I[0][0];
141 >      this->rotate( 1, 2, angle, ji, A );
142 >
143        // rotate about the y-axis
144 <      angle = dt2 * ji[1] / dAtom->getIyy();
145 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
144 >      angle = dt2 * ji[1] / I[1][1];
145 >      this->rotate( 2, 0, angle, ji, A );
146        
147        // rotate about the z-axis
148 <      angle = dt * ji[2] / dAtom->getIzz();
149 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
148 >      angle = dt * ji[2] / I[2][2];
149 >      this->rotate( 0, 1, angle, ji, A);
150        
151        // rotate about the y-axis
152 <      angle = dt2 * ji[1] / dAtom->getIyy();
153 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
152 >      angle = dt2 * ji[1] / I[1][1];
153 >      this->rotate( 2, 0, angle, ji, A );
154        
155         // rotate about the x-axis
156 <      angle = dt2 * ji[0] / dAtom->getIxx();
157 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
156 >      angle = dt2 * ji[0] / I[0][0];
157 >      this->rotate( 1, 2, angle, ji, A );
158        
159 <      dAtom->setJx( ji[0] );
160 <      dAtom->setJy( ji[1] );
161 <      dAtom->setJz( ji[2] );
158 <    }
159 <    
159 >      dAtom->setJ( ji );
160 >      dAtom->setA( A  );    
161 >    }                    
162    }
163 <
163 >  
164    // Scale the box after all the positions have been moved:
165    
166 <
167 <
166 <  // Use a taylor expansion for eta products
167 <  
168 <  info->getBoxM(hm);
166 >  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
167 >  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
168    
169 +  bigScale = 1.0;
170 +  smallScale = 1.0;
171 +  offDiagMax = 0.0;
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 <
190 >      if (i != j)
191 >        if (fabs(scaleMat[i][j]) > offDiagMax)
192 >          offDiagMax = fabs(scaleMat[i][j]);
193 >    }
194  
195 <
196 <
197 <   info->scaleBox(exp(dt*eta));
198 <
199 <
195 >    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
196 >    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
197 >  }
198 >  
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 >  }
228 >  
229   }
230  
231 < void NPTi::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];
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 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 <  instaPress = tStats->getPressure();
248 >  tStats->getPressureTensor(press);
249    instaVol = tStats->getVolume();
250 <
250 >  
251 >  // first evolve chi a half step
252 >  
253    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
197 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
254    
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++ ){
275 <    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 <    for( j=atomIndex; j<(atomIndex+3); j++ )
283 <    for( j=atomIndex; j<(atomIndex+3); j++ )
205 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
206 <                       - vel[j]*(chi+eta));
282 >        
283 >    info->matVecMul3( vScale, vel, sc );
284      
285 +    for (j = 0; j < 3; j++) {
286 +      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 +    }
288 +
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();
215 <      Tb[1] = dAtom->getTy();
216 <      Tb[2] = dAtom->getTz();
217 <      
297 >      dAtom->getTrq( Tb );
298        dAtom->lab2Body( Tb );
299        
300 <      // get the angular momentum, and complete the angular momentum
221 <      // half step
300 >      // get the angular momentum, and propagate a half step
301        
302 <      ji[0] = dAtom->getJx();
224 <      ji[1] = dAtom->getJy();
225 <      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);
229 <      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] );
234 <    }
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 252 | 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 264 | 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 275 | 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();
# Line 288 | Line 378 | int NPTi::readyCheck() {
378  
379    return 1;
380   }
381 +
382 + template<typename T> double NPTf<T>::getConservedQuantity(void){
383 +
384 +  double conservedQuantity;
385 +  double tb2;
386 +  double eta2[3][3];  
387 +  double trEta;
388 +
389 +  //HNVE
390 +  conservedQuantity = tStats->getTotalE();
391 +
392 +  //HNVT
393 +  conservedQuantity += (info->getNDF() * kB * targetTemp *
394 +    (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;
395 +
396 +  //HNPT
397 +  tb2 = tauBarostat *tauBarostat;
398 +
399 +  trEta = info->matTrace3(eta);
400 +  
401 +  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
402 +                        3*NkBT/2 * tb2 * trEta * trEta) / eConvert;
403 +  
404 +  return conservedQuantity;
405 + }

Diff Legend

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