ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPTfm.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTfm.cpp (file contents):
Revision 617 by gezelter, Tue Jul 15 19:56:08 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 10 | Line 10
10   #include "Integrator.hpp"
11   #include "simError.h"
12  
13 +
14   // Basic non-isotropic thermostating and barostating via the Melchionna
15   // modification of the Hoover algorithm:
16   //
# Line 23 | Line 24
24   //  The NPTfm variant scales the molecular center-of-mass coordinates
25   //  instead of the atomic coordinates
26  
27 < NPTfm::NPTfm ( SimInfo *theInfo, ForceFields* the_ff):
28 <  Integrator( theInfo, the_ff )
27 > template<typename T> NPTfm<T>::NPTfm ( SimInfo *theInfo, ForceFields* the_ff):
28 >  T( theInfo, the_ff )
29   {
30    int i, j;
31    chi = 0.0;
32 <
32 >  integralOfChidt = 0.0;
33 >  
34    for(i = 0; i < 3; i++)
35      for (j = 0; j < 3; j++)
36        eta[i][j] = 0.0;
# Line 39 | Line 41 | NPTfm::NPTfm ( SimInfo *theInfo, ForceFields* the_ff):
41    have_target_pressure = 0;
42   }
43  
44 < void NPTfm::moveA() {
44 > template<typename T> void NPTfm<T>::moveA() {
45    
46 <  int i, j, k;
47 <  DirectionalAtom* dAtom;
48 <  double Tb[3], ji[3];
49 <  double A[3][3], I[3][3];
50 <  double angle, mass;
51 <  double vel[3], pos[3], frc[3];
46 > //   int i, j, k;
47 > //   DirectionalAtom* dAtom;
48 > //   double Tb[3], ji[3];
49 > //   double A[3][3], I[3][3];
50 > //   double angle, mass;
51 > //   double vel[3], pos[3], frc[3];
52  
53 <  double rj[3];
54 <  double instaTemp, instaPress, instaVol;
55 <  double tt2, tb2;
56 <  double sc[3];
57 <  double eta2ij, smallScale, bigScale, offDiagMax;
58 <  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
53 > //   double rj[3];
54 > //   double instaTemp, instaPress, instaVol;
55 > //   double tt2, tb2;
56 > //   double sc[3];
57 > //   double eta2ij, smallScale, bigScale, offDiagMax;
58 > //   double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
59  
60 <  int nInMol;
61 <  double rc[3];
60 > //   int nInMol;
61 > //   double rc[3];
62  
63 <  nMols = info->n_mol;
64 <  myMolecules = info->molecules;
63 > // /*
64 > //   nMols = info->n_mol;
65 > //   myMolecules = info->molecules;
66  
67 <  tt2 = tauThermostat * tauThermostat;
68 <  tb2 = tauBarostat * tauBarostat;
67 > //   tt2 = tauThermostat * tauThermostat;
68 > //   tb2 = tauBarostat * tauBarostat;
69  
70 <  instaTemp = tStats->getTemperature();
71 <  tStats->getPressureTensor(press);
72 <  instaVol = tStats->getVolume();
70 > //   instaTemp = tStats->getTemperature();
71 > //   tStats->getPressureTensor(press);
72 > //   instaVol = tStats->getVolume();
73    
74 <  // first evolve chi a half step
74 > //   // first evolve chi a half step
75    
76 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
76 > //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
77  
78 <  for (i = 0; i < 3; i++ ) {
79 <    for (j = 0; j < 3; j++ ) {
80 <      if (i == j) {
78 > //   for (i = 0; i < 3; i++ ) {
79 > //     for (j = 0; j < 3; j++ ) {
80 > //       if (i == j) {
81          
82 <        eta[i][j] += dt2 * instaVol *
83 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
82 > //         eta[i][j] += dt2 * instaVol *
83 > //           (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
84          
85 <        vScale[i][j] = eta[i][j] + chi;
85 > //         vScale[i][j] = eta[i][j] + chi;
86          
87 <      } else {
87 > //       } else {
88          
89 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
89 > //         eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
90  
91 <        vScale[i][j] = eta[i][j];
91 > //         vScale[i][j] = eta[i][j];
92          
93 <      }
94 <    }
95 <  }
93 > //       }
94 > //     }
95 > //   }
96  
97  
98 <  for (i = 0; i < nMols; i++) {
98 > //   for (i = 0; i < nMols; i++) {
99  
100 <    myMolecules[i].getCOM(rc);
100 > //     myMolecules[i].getCOM(rc);
101      
102 <    nInMol = myMolecules[i].getNAtoms();
103 <    myAtoms = myMolecules[i].getMyAtoms();
102 > //     nInMol = myMolecules[i].getNAtoms();
103 > //     myAtoms = myMolecules[i].getMyAtoms();
104  
105 <    // find the minimum image coordinates of the molecular centers of mass:
105 > //     // find the minimum image coordinates of the molecular centers of mass:
106  
107 <    info->wrapVector(rc);
107 > //     info->wrapVector(rc);
108  
109 <    for( j=0; j< nInMol; j++ ){
109 > //     for( j=0; j< nInMol; j++ ){
110        
111 <      if(myAtoms[j] != NULL) {
111 > //       if(myAtoms[j] != NULL) {
112        
113 <        myAtoms[j]->getVel( vel );
114 <        myAtoms[j]->getPos( pos );
115 <        myAtoms[j]->getFrc( frc );
113 > //         myAtoms[j]->getVel( vel );
114 > //         myAtoms[j]->getPos( pos );
115 > //         myAtoms[j]->getFrc( frc );
116  
117 <        mass = myAtoms[j]->getMass();
117 > //         mass = myAtoms[j]->getMass();
118      
119 <        // velocity half step
119 > //         // velocity half step
120          
121 <        info->matVecMul3( vScale, vel, sc );
121 > //         info->matVecMul3( vScale, vel, sc );
122      
123 <        for (k = 0; k < 3; k++)
124 <          vel[k] += dt2 * ((frc[k]  / mass) * eConvert - sc[k]);
123 > //         for (k = 0; k < 3; k++)
124 > //           vel[k] += dt2 * ((frc[k]  / mass) * eConvert - sc[k]);
125  
126 <        myAtoms[j]->setVel( vel );
126 > //         myAtoms[j]->setVel( vel );
127  
128 <        // position whole step    
128 > //         // position whole step    
129          
130 <        info->matVecMul3( eta, rc, sc );
130 > //         info->matVecMul3( eta, rc, sc );
131  
132 <        for (k = 0; k < 3; k++ )
133 <          pos[k] += dt * (vel[k] + sc[k]);
132 > //         for (k = 0; k < 3; k++ )
133 > //           pos[k] += dt * (vel[k] + sc[k]);
134  
135 <        myAtoms[j]->setPos( pos );
135 > //         myAtoms[j]->setPos( pos );
136    
137 <        if( myAtoms[j]->isDirectional() ){
137 > //         if( myAtoms[j]->isDirectional() ){
138  
139 <          dAtom = (DirectionalAtom *)myAtoms[j];
139 > //           dAtom = (DirectionalAtom *)myAtoms[j];
140            
141 <          // get and convert the torque to body frame
141 > //           // get and convert the torque to body frame
142            
143 <          dAtom->getTrq( Tb );
144 <          dAtom->lab2Body( Tb );
143 > //           dAtom->getTrq( Tb );
144 > //           dAtom->lab2Body( Tb );
145        
146 <          // get the angular momentum, and propagate a half step
146 > //           // get the angular momentum, and propagate a half step
147            
148 <          dAtom->getJ( ji );
148 > //           dAtom->getJ( ji );
149  
150 <          for (k=0; k < 3; k++)
151 <            ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
150 > //           for (k=0; k < 3; k++)
151 > //             ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
152        
153 <          // use the angular velocities to propagate the rotation matrix a
154 <          // full time step
153 > //           // use the angular velocities to propagate the rotation matrix a
154 > //           // full time step
155            
156 <          dAtom->getA(A);
157 <          dAtom->getI(I);
156 > //           dAtom->getA(A);
157 > //           dAtom->getI(I);
158            
159 <          // rotate about the x-axis      
160 <          angle = dt2 * ji[0] / I[0][0];
161 <          this->rotate( 1, 2, angle, ji, A );
159 > //           // rotate about the x-axis      
160 > //           angle = dt2 * ji[0] / I[0][0];
161 > //           this->rotate( 1, 2, angle, ji, A );
162            
163 <          // rotate about the y-axis
164 <          angle = dt2 * ji[1] / I[1][1];
165 <          this->rotate( 2, 0, angle, ji, A );
163 > //           // rotate about the y-axis
164 > //           angle = dt2 * ji[1] / I[1][1];
165 > //           this->rotate( 2, 0, angle, ji, A );
166            
167 <          // rotate about the z-axis
168 <          angle = dt * ji[2] / I[2][2];
169 <          this->rotate( 0, 1, angle, ji, A);
167 > //           // rotate about the z-axis
168 > //           angle = dt * ji[2] / I[2][2];
169 > //           this->rotate( 0, 1, angle, ji, A);
170            
171 <          // rotate about the y-axis
172 <          angle = dt2 * ji[1] / I[1][1];
173 <          this->rotate( 2, 0, angle, ji, A );
171 > //           // rotate about the y-axis
172 > //           angle = dt2 * ji[1] / I[1][1];
173 > //           this->rotate( 2, 0, angle, ji, A );
174            
175 <          // rotate about the x-axis
176 <          angle = dt2 * ji[0] / I[0][0];
177 <          this->rotate( 1, 2, angle, ji, A );
175 > //           // rotate about the x-axis
176 > //           angle = dt2 * ji[0] / I[0][0];
177 > //           this->rotate( 1, 2, angle, ji, A );
178        
179 <          dAtom->setJ( ji );
180 <          dAtom->setA( A  );    
181 <        }                    
182 <      }
183 <    }
184 <  }
179 > //           dAtom->setJ( ji );
180 > //           dAtom->setA( A  );    
181 > //         }                    
182 > //       }
183 > //     }
184 > //   }
185    
186 <  // Scale the box after all the positions have been moved:
186 > //   // Scale the box after all the positions have been moved:
187    
188 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
189 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
188 > //   // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
189 > //   //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
190    
191    
192 <  bigScale = 1.0;
193 <  smallScale = 1.0;
194 <  offDiagMax = 0.0;
192 > //   bigScale = 1.0;
193 > //   smallScale = 1.0;
194 > //   offDiagMax = 0.0;
195  
196 <  for(i=0; i<3; i++){
197 <    for(j=0; j<3; j++){
196 > //   for(i=0; i<3; i++){
197 > //     for(j=0; j<3; j++){
198        
199 <      // Calculate the matrix Product of the eta array (we only need
200 <      // the ij element right now):
199 > //       // Calculate the matrix Product of the eta array (we only need
200 > //       // the ij element right now):
201        
202 <      eta2ij = 0.0;
203 <      for(k=0; k<3; k++){
204 <        eta2ij += eta[i][k] * eta[k][j];
205 <      }
206 <      
207 <      scaleMat[i][j] = 0.0;
208 <      // identity matrix (see above):
209 <      if (i == j) scaleMat[i][j] = 1.0;
210 <      // Taylor expansion for the exponential truncated at second order:
211 <      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
202 > //       eta2ij = 0.0;
203 > //       for(k=0; k<3; k++){
204 > //         eta2ij += eta[i][k] * eta[k][j];
205 > //       }
206 >      
207 > //       scaleMat[i][j] = 0.0;
208 > //       // identity matrix (see above):
209 > //       if (i == j) scaleMat[i][j] = 1.0;
210 > //       // Taylor expansion for the exponential truncated at second order:
211 > //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
212  
213 <      if (i != j)
214 <        if (fabs(scaleMat[i][j]) > offDiagMax)
215 <          offDiagMax = fabs(scaleMat[i][j]);      
216 <    }
217 <    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
218 <    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
219 <  }
213 > //       if (i != j)
214 > //         if (fabs(scaleMat[i][j]) > offDiagMax)
215 > //           offDiagMax = fabs(scaleMat[i][j]);      
216 > //     }
217 > //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
218 > //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
219 > //   }
220    
221 <  if ((bigScale > 1.1) || (smallScale < 0.9)) {
222 <    sprintf( painCave.errMsg,
223 <             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
224 <             " Check your tauBarostat, as it is probably too small!\n\n"
225 <             " scaleMat = [%lf\t%lf\t%lf]\n"
226 <             "            [%lf\t%lf\t%lf]\n"
227 <             "            [%lf\t%lf\t%lf]\n",
228 <             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
229 <             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
230 <             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
231 <    painCave.isFatal = 1;
232 <    simError();
233 <  } else if (offDiagMax > 0.1) {
234 <    sprintf( painCave.errMsg,
235 <             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
236 <             " Check your tauBarostat, as it is probably too small!\n\n"
237 <             " scaleMat = [%lf\t%lf\t%lf]\n"
238 <             "            [%lf\t%lf\t%lf]\n"
239 <             "            [%lf\t%lf\t%lf]\n",
240 <             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
241 <             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
242 <             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
243 <    painCave.isFatal = 1;
244 <    simError();
245 <  } else {
246 <    info->getBoxM(hm);
247 <    info->matMul3(hm, scaleMat, hmnew);
248 <    info->setBoxM(hmnew);
249 <  }  
250 < }
221 > //   if ((bigScale > 1.1) || (smallScale < 0.9)) {
222 > //     sprintf( painCave.errMsg,
223 > //              "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
224 > //              " Check your tauBarostat, as it is probably too small!\n\n"
225 > //              " scaleMat = [%lf\t%lf\t%lf]\n"
226 > //              "            [%lf\t%lf\t%lf]\n"
227 > //              "            [%lf\t%lf\t%lf]\n",
228 > //              scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
229 > //              scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
230 > //              scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
231 > //     painCave.isFatal = 1;
232 > //     simError();
233 > //   } else if (offDiagMax > 0.1) {
234 > //     sprintf( painCave.errMsg,
235 > //              "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
236 > //              " Check your tauBarostat, as it is probably too small!\n\n"
237 > //              " scaleMat = [%lf\t%lf\t%lf]\n"
238 > //              "            [%lf\t%lf\t%lf]\n"
239 > //              "            [%lf\t%lf\t%lf]\n",
240 > //              scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
241 > //              scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
242 > //              scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
243 > //     painCave.isFatal = 1;
244 > //     simError();
245 > //   } else {
246 > //     info->getBoxM(hm);
247 > //     info->matMul3(hm, scaleMat, hmnew);
248 > //     info->setBoxM(hmnew);
249 > //   }  
250 > // */
251  
252 < void NPTfm::moveB( void ){
252 > //   tt2 = tauThermostat * tauThermostat;
253 > //   tb2 = tauBarostat * tauBarostat;
254  
255 <  int i, j;
256 <  DirectionalAtom* dAtom;
257 <  double Tb[3], ji[3];
254 <  double vel[3], frc[3];
255 <  double mass;
256 <
257 <  double instaTemp, instaPress, instaVol;
258 <  double tt2, tb2;
259 <  double sc[3];
260 <  double press[3][3], vScale[3][3];
255 > //   instaTemp = tStats->getTemperature();
256 > //   tStats->getPressureTensor(press);
257 > //   instaVol = tStats->getVolume();
258    
259 <  tt2 = tauThermostat * tauThermostat;
263 <  tb2 = tauBarostat * tauBarostat;
259 > //   tStats->getCOM(COM);
260  
261 <  instaTemp = tStats->getTemperature();
262 <  tStats->getPressureTensor(press);
263 <  instaVol = tStats->getVolume();
264 <  
265 <  // first evolve chi a half step
266 <  
267 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
268 <  
269 <  for (i = 0; i < 3; i++ ) {
270 <    for (j = 0; j < 3; j++ ) {
275 <      if (i == j) {
261 > //   //calculate scale factor of veloity
262 > //   for (i = 0; i < 3; i++ ) {
263 > //     for (j = 0; j < 3; j++ ) {
264 > //       vScale[i][j] = eta[i][j];
265 >      
266 > //       if (i == j) {
267 > //         vScale[i][j] += chi;          
268 > //       }              
269 > //     }
270 > //   }
271  
277 –        eta[i][j] += dt2 * instaVol *
278 –          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
272  
273 <        vScale[i][j] = eta[i][j] + chi;
281 <        
282 <      } else {
283 <        
284 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
273 > //   for (i = 0; i < nMols; i++) {
274  
275 <        vScale[i][j] = eta[i][j];
276 <        
277 <      }
278 <    }
290 <  }
275 > //     myMolecules[i].getCOM(rc);
276 >    
277 > //     nInMol = myMolecules[i].getNAtoms();
278 > //     myAtoms = myMolecules[i].getMyAtoms();
279  
292 –  for( i=0; i<nAtoms; i++ ){
280  
281 <    atoms[i]->getVel( vel );
282 <    atoms[i]->getFrc( frc );
281 > //     for( j=0; j< nInMol; j++ ){
282 >      
283 > //       if(myAtoms[j] != NULL) {
284 >      
285 > //         myAtoms[j]->getVel( vel );
286 > //         myAtoms[j]->getFrc( frc );
287  
288 <    mass = atoms[i]->getMass();
288 > //         mass = myAtoms[j]->getMass();
289      
290 <    // velocity half step
290 > //         // velocity half step
291          
292 <    info->matVecMul3( vScale, vel, sc );
292 > //         info->matVecMul3( vScale, vel, sc );
293      
294 <    for (j = 0; j < 3; j++) {
295 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
305 <    }
294 > //         for (k = 0; k < 3; k++)
295 > //           vel[k] += dt2 * ((frc[k]  / mass) * eConvert - sc[k]);
296  
297 <    atoms[i]->setVel( vel );
298 <    
299 <    if( atoms[i]->isDirectional() ){
297 > //         myAtoms[j]->setVel( vel );
298 >    
299 > //         if( myAtoms[j]->isDirectional() ){
300  
301 <      dAtom = (DirectionalAtom *)atoms[i];
301 > //           dAtom = (DirectionalAtom *)myAtoms[j];
302            
303 <      // get and convert the torque to body frame
303 > //           // get and convert the torque to body frame
304 >          
305 > //           dAtom->getTrq( Tb );
306 > //           dAtom->lab2Body( Tb );
307        
308 <      dAtom->getTrq( Tb );
309 <      dAtom->lab2Body( Tb );
310 <      
318 <      // get the angular momentum, and propagate a half step
319 <      
320 <      dAtom->getJ( ji );
321 <      
322 <      for (j=0; j < 3; j++)
323 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
324 <      
325 <      dAtom->setJ( ji );
308 > //           // get the angular momentum, and propagate a half step
309 >          
310 > //           dAtom->getJ( ji );
311  
312 <    }                    
313 <  }
314 < }
312 > //           for (k=0; k < 3; k++)
313 > //             ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
314 >      
315 > //           // use the angular velocities to propagate the rotation matrix a
316 > //           // full time step
317 >          
318 > //           dAtom->getA(A);
319 > //           dAtom->getI(I);
320 >          
321 > //           // rotate about the x-axis      
322 > //           angle = dt2 * ji[0] / I[0][0];
323 > //           this->rotate( 1, 2, angle, ji, A );
324 >          
325 > //           // rotate about the y-axis
326 > //           angle = dt2 * ji[1] / I[1][1];
327 > //           this->rotate( 2, 0, angle, ji, A );
328 >          
329 > //           // rotate about the z-axis
330 > //           angle = dt * ji[2] / I[2][2];
331 > //           this->rotate( 0, 1, angle, ji, A);
332 >          
333 > //           // rotate about the y-axis
334 > //           angle = dt2 * ji[1] / I[1][1];
335 > //           this->rotate( 2, 0, angle, ji, A );
336 >          
337 > //           // rotate about the x-axis
338 > //           angle = dt2 * ji[0] / I[0][0];
339 > //           this->rotate( 1, 2, angle, ji, A );
340 >      
341 > //           dAtom->setJ( ji );
342 > //           dAtom->setA( A  );    
343 > //         }                    
344 > //       }
345 > //     }
346 > //   }
347  
331 – int NPTfm::readyCheck() {
348  
349 + //   // advance chi half step
350 + //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
351 +
352 + //   //calculate the integral of chidt
353 + //   integralOfChidt += dt2*chi;
354 +
355 + //   //advance eta half step
356 + //   for(i = 0; i < 3; i ++)
357 + //     for(j = 0; j < 3; j++){
358 + //       if( i == j)
359 + //         eta[i][j] += dt2 *  instaVol *
360 + //        (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
361 + //       else
362 + //         eta[i][j] += dt2 * instaVol * press[i][j] / ( NkBT*tb2);
363 + //     }
364 +    
365 + //   //save the old positions
366 + //   for(i = 0; i < nAtoms; i++){
367 + //     atoms[i]->getPos(pos);
368 + //     for(j = 0; j < 3; j++)
369 + //       oldPos[i*3 + j] = pos[j];
370 + //   }
371 +  
372 + //   //the first estimation of r(t+dt) is equal to  r(t)
373 +    
374 + //   for(k = 0; k < 4; k ++){
375 +
376 + //     for(i =0 ; i < nAtoms; i++){
377 +
378 + //       atoms[i]->getVel(vel);
379 + //       atoms[i]->getPos(pos);
380 +
381 + //       for(j = 0; j < 3; j++)
382 + //         rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
383 +      
384 + //       info->matVecMul3( eta, rj, sc );
385 +      
386 + //       for(j = 0; j < 3; j++)
387 + //         pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
388 +
389 + //       atoms[i]->setPos( pos );
390 +
391 + //     }
392 +
393 + //   }  
394 +
395 +
396 + //   // Scale the box after all the positions have been moved:
397 +  
398 + //   // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
399 + //   //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
400 +  
401 + //   bigScale = 1.0;
402 + //   smallScale = 1.0;
403 + //   offDiagMax = 0.0;
404 +  
405 + //   for(i=0; i<3; i++){
406 + //     for(j=0; j<3; j++){
407 +      
408 + //       // Calculate the matrix Product of the eta array (we only need
409 + //       // the ij element right now):
410 +      
411 + //       eta2ij = 0.0;
412 + //       for(k=0; k<3; k++){
413 + //         eta2ij += eta[i][k] * eta[k][j];
414 + //       }
415 +      
416 + //       scaleMat[i][j] = 0.0;
417 + //       // identity matrix (see above):
418 + //       if (i == j) scaleMat[i][j] = 1.0;
419 + //       // Taylor expansion for the exponential truncated at second order:
420 + //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
421 +
422 + //       if (i != j)
423 + //         if (fabs(scaleMat[i][j]) > offDiagMax)
424 + //           offDiagMax = fabs(scaleMat[i][j]);
425 + //     }
426 +
427 + //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
428 + //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
429 + //   }
430 +  
431 + //   if ((bigScale > 1.1) || (smallScale < 0.9)) {
432 + //     sprintf( painCave.errMsg,
433 + //              "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
434 + //              " Check your tauBarostat, as it is probably too small!\n\n"
435 + //              " scaleMat = [%lf\t%lf\t%lf]\n"
436 + //              "            [%lf\t%lf\t%lf]\n"
437 + //              "            [%lf\t%lf\t%lf]\n",
438 + //              scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
439 + //              scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
440 + //              scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
441 + //     painCave.isFatal = 1;
442 + //     simError();
443 + //   } else if (offDiagMax > 0.1) {
444 + //     sprintf( painCave.errMsg,
445 + //              "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
446 + //              " Check your tauBarostat, as it is probably too small!\n\n"
447 + //              " scaleMat = [%lf\t%lf\t%lf]\n"
448 + //              "            [%lf\t%lf\t%lf]\n"
449 + //              "            [%lf\t%lf\t%lf]\n",
450 + //              scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
451 + //              scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
452 + //              scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
453 + //     painCave.isFatal = 1;
454 + //     simError();
455 + //   } else {
456 + //     info->getBoxM(hm);
457 + //     info->matMul3(hm, scaleMat, hmnew);
458 + //     info->setBoxM(hmnew);
459 + //   }
460 +  
461 +
462 +
463 + }
464 +
465 + template<typename T> void NPTfm<T>::moveB( void ){
466 +
467 + //   int i, j;
468 + //   DirectionalAtom* dAtom;
469 + //   double Tb[3], ji[3];
470 + //   double vel[3], frc[3];
471 + //   double mass;
472 +
473 + //   double instaTemp, instaPress, instaVol;
474 + //   double tt2, tb2;
475 + //   double sc[3];
476 + //   double press[3][3], vScale[3][3];
477 + //   double oldChi, prevChi;
478 + //   double oldEta[3][3], preEta[3][3], diffEta;  
479 +
480 + // /*  
481 + //   tt2 = tauThermostat * tauThermostat;
482 + //   tb2 = tauBarostat * tauBarostat;
483 +
484 + //   instaTemp = tStats->getTemperature();
485 + //   tStats->getPressureTensor(press);
486 + //   instaVol = tStats->getVolume();
487 +  
488 + //   // first evolve chi a half step
489 +  
490 + //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
491 +  
492 + //   for (i = 0; i < 3; i++ ) {
493 + //     for (j = 0; j < 3; j++ ) {
494 + //       if (i == j) {
495 +
496 + //         eta[i][j] += dt2 * instaVol *
497 + //           (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
498 +
499 + //         vScale[i][j] = eta[i][j] + chi;
500 +        
501 + //       } else {
502 +        
503 + //         eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
504 +
505 + //         vScale[i][j] = eta[i][j];
506 +        
507 + //       }
508 + //     }
509 + //   }
510 +
511 + //   for( i=0; i<nAtoms; i++ ){
512 +
513 + //     atoms[i]->getVel( vel );
514 + //     atoms[i]->getFrc( frc );
515 +
516 + //     mass = atoms[i]->getMass();
517 +    
518 + //     // velocity half step
519 +        
520 + //     info->matVecMul3( vScale, vel, sc );
521 +    
522 + //     for (j = 0; j < 3; j++) {
523 + //       vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
524 + //     }
525 +
526 + //     atoms[i]->setVel( vel );
527 +    
528 + //     if( atoms[i]->isDirectional() ){
529 +
530 + //       dAtom = (DirectionalAtom *)atoms[i];
531 +          
532 + //       // get and convert the torque to body frame
533 +      
534 + //       dAtom->getTrq( Tb );
535 + //       dAtom->lab2Body( Tb );
536 +      
537 + //       // get the angular momentum, and propagate a half step
538 +      
539 + //       dAtom->getJ( ji );
540 +      
541 + //       for (j=0; j < 3; j++)
542 + //         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
543 +      
544 + //       dAtom->setJ( ji );
545 +
546 + //     }                    
547 + //   }
548 + // */
549 +
550 + //   tt2 = tauThermostat * tauThermostat;
551 + //   tb2 = tauBarostat * tauBarostat;
552 +
553 +
554 + //   // Set things up for the iteration:
555 +
556 + //   oldChi = chi;
557 +  
558 + //   for(i = 0; i < 3; i++)
559 + //     for(j = 0; j < 3; j++)
560 + //       oldEta[i][j] = eta[i][j];
561 +
562 + //   for( i=0; i<nAtoms; i++ ){
563 +
564 + //     atoms[i]->getVel( vel );
565 +
566 + //     for (j=0; j < 3; j++)
567 + //       oldVel[3*i + j]  = vel[j];
568 +
569 + //     if( atoms[i]->isDirectional() ){
570 +
571 + //       dAtom = (DirectionalAtom *)atoms[i];
572 +
573 + //       dAtom->getJ( ji );
574 +
575 + //       for (j=0; j < 3; j++)
576 + //         oldJi[3*i + j] = ji[j];
577 +
578 + //     }
579 + //   }
580 +
581 + //   // do the iteration:
582 +
583 + //   instaVol = tStats->getVolume();
584 +  
585 + //   for (k=0; k < 4; k++) {
586 +    
587 + //     instaTemp = tStats->getTemperature();
588 + //     tStats->getPressureTensor(press);
589 +
590 + //     // evolve chi another half step using the temperature at t + dt/2
591 +
592 + //     prevChi = chi;
593 + //     chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
594 +    
595 + //     for(i = 0; i < 3; i++)
596 + //       for(j = 0; j < 3; j++)
597 + //         preEta[i][j] = eta[i][j];
598 +
599 + //     //advance eta half step and calculate scale factor for velocity
600 + //     for(i = 0; i < 3; i ++)
601 + //       for(j = 0; j < 3; j++){
602 + //         if( i == j){
603 + //           eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
604 + //          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
605 + //           vScale[i][j] = eta[i][j] + chi;
606 + //         }
607 + //         else
608 + //         {
609 + //           eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
610 + //           vScale[i][j] = eta[i][j];
611 + //         }
612 + //     }      
613 +
614 + //     //advance velocity half step
615 + //     for( i=0; i<nAtoms; i++ ){
616 +
617 + //       atoms[i]->getFrc( frc );
618 + //       atoms[i]->getVel(vel);
619 +      
620 + //       mass = atoms[i]->getMass();
621 +      
622 + //       info->matVecMul3( vScale, vel, sc );
623 +
624 + //       for (j=0; j < 3; j++) {
625 + //         // velocity half step  (use chi from previous step here):
626 + //         vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
627 + //       }
628 +      
629 + //       atoms[i]->setVel( vel );
630 +      
631 + //       if( atoms[i]->isDirectional() ){
632 +
633 + //         dAtom = (DirectionalAtom *)atoms[i];
634 +  
635 + //         // get and convert the torque to body frame      
636 +  
637 + //         dAtom->getTrq( Tb );
638 + //         dAtom->lab2Body( Tb );      
639 +            
640 + //         for (j=0; j < 3; j++)
641 + //           ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
642 +      
643 + //           dAtom->setJ( ji );
644 + //       }
645 + //     }
646 +
647 +    
648 + //     diffEta = 0;
649 + //     for(i = 0; i < 3; i++)
650 + //       diffEta += pow(preEta[i][i] - eta[i][i], 2);    
651 +    
652 + //     if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
653 + //       break;
654 + //   }
655 +
656 + //   //calculate integral of chida
657 + //   integralOfChidt += dt2*chi;
658 + }
659 +
660 + template<typename T> void NPTfm<T>::resetIntegrator() {
661 +  int i,j;
662 +  
663 +  chi = 0.0;
664 +
665 +  for(i = 0; i < 3; i++)
666 +    for (j = 0; j < 3; j++)
667 +      eta[i][j] = 0.0;
668 + }
669 +
670 + template<typename T> int NPTfm<T>::readyCheck() {
671 +
672 +  //check parent's readyCheck() first
673 +  if (T::readyCheck() == -1)
674 +    return -1;
675 +  
676    // First check to see if we have a target temperature.
677    // Not having one is fatal.
678    
# Line 381 | Line 724 | int NPTfm::readyCheck() {
724  
725    return 1;
726   }
727 +
728 + template<typename T> double NPTfm<T>::getConservedQuantity(void){
729 +
730 +
731 + //   double conservedQuantity;
732 + //   double tb2;
733 + //   double trEta;  
734 + //   double E_NPT;
735 + //   double U;
736 + //   double TS;
737 + //   double PV;
738 + //   double extra;
739 +
740 + //   U = tStats->getTotalE();
741 +
742 + //   TS = fkBT *
743 + //     (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
744 +
745 + //   PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
746 +
747 + //   tb2 = tauBarostat * tauBarostat;
748 +
749 + //   trEta = info->matTrace3(eta);
750 +
751 + //   extra = (fkBT * tb2 * trEta * trEta / 2.0 ) / eConvert;
752 +
753 + //   cout.width(8);
754 + //   cout.precision(8);
755 +  
756 + //   cout << info->getTime() << "\t"
757 + //        << chi << "\t"
758 + //        << trEta << "\t"
759 + //        << U << "\t"
760 + //        << TS << "\t"
761 + //        << PV << "\t"
762 + //        << extra << "\t"
763 + //        << U+TS+PV+extra << endl;
764 +
765 + //   conservedQuantity = U+TS+PV+extra;
766 + //   return conservedQuantity;
767 + }

Diff Legend

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