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 645 by tim, Tue Jul 22 19:54:52 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 >
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];
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);
53 –
54 –  for (i=0; i < 9; i++) press[i] *= p_convert;
55 –
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 <    
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);
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;
91 <    vel[atomIndex+1] = vy;
92 <    vel[atomIndex+2] = vz;
103 >    atoms[i]->setVel( vel );
104  
105      // position whole step    
106  
96 –    rj[0] = pos[atomIndex];
97 –    rj[1] = pos[atomIndex+1];
98 –    rj[2] = pos[atomIndex+2];
99 –
107      info->wrapVector(rj);
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, 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 113 | Line 119 | void NPTf::moveA() {
119            
120        // get and convert the torque to body frame
121        
122 <      Tb[0] = dAtom->getTx();
117 <      Tb[1] = dAtom->getTy();
118 <      Tb[2] = dAtom->getTz();
119 <      
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        
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 –      
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] );
158 <    }
159 <    
158 >      dAtom->setJ( ji );
159 >      dAtom->setA( A  );    
160 >    }                    
161    }
162 <
162 >  
163    // Scale the box after all the positions have been moved:
164    
165 <
166 <
166 <  // Use a taylor expansion for eta products
167 <  
168 <  info->getBoxM(hm);
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 <
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 >      
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 >        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 +    }
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> int NPTf<T>::readyCheck() {
314  
315    // First check to see if we have a target temperature.
316    // Not having one is fatal.
317    
318    if (!have_target_temp) {
319      sprintf( painCave.errMsg,
320 <             "NPTi error: You can't use the NPTi integrator\n"
320 >             "NPTf error: You can't use the NPTf integrator\n"
321               "   without a targetTemp!\n"
322               );
323      painCave.isFatal = 1;
# Line 252 | Line 327 | int NPTi::readyCheck() {
327  
328    if (!have_target_pressure) {
329      sprintf( painCave.errMsg,
330 <             "NPTi error: You can't use the NPTi integrator\n"
330 >             "NPTf error: You can't use the NPTf integrator\n"
331               "   without a targetPressure!\n"
332               );
333      painCave.isFatal = 1;
# Line 264 | Line 339 | int NPTi::readyCheck() {
339    
340    if (!have_tau_thermostat) {
341      sprintf( painCave.errMsg,
342 <             "NPTi error: If you use the NPTi\n"
342 >             "NPTf error: If you use the NPTf\n"
343               "   integrator, you must set tauThermostat.\n");
344      painCave.isFatal = 1;
345      simError();
# Line 275 | Line 350 | int NPTi::readyCheck() {
350    
351    if (!have_tau_barostat) {
352      sprintf( painCave.errMsg,
353 <             "NPTi error: If you use the NPTi\n"
353 >             "NPTf error: If you use the NPTf\n"
354               "   integrator, you must set tauBarostat.\n");
355      painCave.isFatal = 1;
356      simError();

Diff Legend

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