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 588 by gezelter, Thu Jul 10 17:10:56 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
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, j;
27    chi = 0.0;
28  
29    for(i = 0; i < 3; i++)
30 <    for (j = 0; j < 3; j_++)
30 >    for (j = 0; j < 3; j++)
31        eta[i][j] = 0.0;
32  
33    have_tau_thermostat = 0;
# Line 35 | Line 36 | NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
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;
41 <  int atomIndex, aMatIndex;
41 >  int i, j, k;
42    DirectionalAtom* dAtom;
43 <  double Tb[3];
44 <  double ji[3];
45 <  double ri[3], vi[3], sc[3];
46 <  double instaTemp, instaVol;
47 <  double tt2, tb2, eta2ij;
48 <  double angle;
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 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;
# Line 62 | Line 67 | void NPTf::moveA() {
67    for (i = 0; i < 3; i++ ) {
68      for (j = 0; j < 3; j++ ) {
69        if (i == j) {
70 <
70 >        
71          eta[i][j] += dt2 * instaVol *
72            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
73 <
73 >        
74          vScale[i][j] = eta[i][j] + chi;
75          
76        } else {
# Line 79 | Line 84 | void NPTf::moveA() {
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 <    vi[0] = vel[atomIndex];
99 <    vi[1] = vel[atomIndex+1];
100 <    vi[2] = vel[atomIndex+2];
101 <    
91 <    info->matVecMul3( vScale, vi, sc );
92 <    
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]);
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] = vi[0]
98 <    vel[atomIndex+1] = vi[1];
99 <    vel[atomIndex+2] = vi[2];
103 >    atoms[i]->setVel( vel );
104  
105      // position whole step    
106  
107 <    ri[0] = pos[atomIndex];
104 <    ri[1] = pos[atomIndex+1];
105 <    ri[2] = pos[atomIndex+2];
107 >    info->wrapVector(rj);
108  
109 <    info->wrapVector(ri);
109 >    info->matVecMul3( eta, rj, sc );
110  
111 <    info->matVecMul3( eta, ri, sc );
111 >    for (j = 0; j < 3; j++ )
112 >      pos[j] += dt * (vel[j] + sc[j]);
113  
114 <    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 >    atoms[i]->setPos( pos );
115    
116      if( atoms[i]->isDirectional() ){
117  
# Line 118 | Line 119 | void NPTf::moveA() {
119            
120        // get and convert the torque to body frame
121        
122 <      Tb[0] = dAtom->getTx();
122 <      Tb[1] = dAtom->getTy();
123 <      Tb[2] = dAtom->getTz();
124 <      
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        
133 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
134 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
135 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
136 –      
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] );
163 <    }
164 <    
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 <
174 >      
175        // Calculate the matrix Product of the eta array (we only need
176        // the ij element right now):
177 <
177 >      
178        eta2ij = 0.0;
179        for(k=0; k<3; k++){
180          eta2ij += eta[i][k] * eta[k][j];
# Line 187 | Line 186 | void NPTf::moveA() {
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 +    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
196 +    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
197    }
192 –  
193 –  info->getBoxM(hm);
194 –  info->matMul3(hm, scaleMat, hmnew);
195 –  info->setBoxM(hmnew);
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 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 vi[3], sc[3];
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 sc[3];
242    double press[3][3], vScale[3][3];
243    
244    tt2 = tauThermostat * tauThermostat;
# Line 238 | Line 272 | void NPTf::moveB( void ){
272    }
273  
274    for( i=0; i<nAtoms; i++ ){
241 –    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 <    vi[0] = vel[atomIndex];
286 <    vi[1] = vel[atomIndex+1];
287 <    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]);
285 >    for (j = 0; j < 3; j++) {
286 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 >    }
288  
289 <    vel[atomIndex] = vi[0]
256 <    vel[atomIndex+1] = vi[1];
257 <    vel[atomIndex+2] = vi[2];
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();
266 <      Tb[1] = dAtom->getTy();
267 <      Tb[2] = dAtom->getTz();
268 <      
297 >      dAtom->getTrq( Tb );
298        dAtom->lab2Body( Tb );
299        
300 <      // get the angular momentum, and complete the angular momentum
272 <      // half step
300 >      // get the angular momentum, and propagate a half step
301        
302 <      ji[0] = dAtom->getJx();
275 <      ji[1] = dAtom->getJy();
276 <      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);
280 <      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] );
285 <    }
307 >      dAtom->setJ( ji );
308 >
309 >    }                    
310    }
311   }
312  
313 < int NPTf::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.

Diff Legend

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