ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/OOPSE/libmdtools/NPTf.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 576 by gezelter, Tue Jul 8 21:10:16 2003 UTC vs.
Revision 588 by gezelter, Thu Jul 10 17:10:56 2003 UTC

# Line 9 | Line 9
9   #include "simError.h"
10  
11  
12 < // Basic isotropic thermostating and barostating via the Melchionna
12 > // Basic non-isotropic thermostating and barostating via the Melchionna
13   // modification of the Hoover algorithm:
14   //
15   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
# Line 19 | Line 19
19   //
20   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
21  
22 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
22 > NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
23    Integrator( theInfo, the_ff )
24   {
25 <  int i;
25 >  int i, j;
26    chi = 0.0;
27 <  for(i = 0; i < 9; i++) eta[i] = 0.0;
27 >
28 >  for(i = 0; i < 3; i++)
29 >    for (j = 0; j < 3; j_++)
30 >      eta[i][j] = 0.0;
31 >
32    have_tau_thermostat = 0;
33    have_tau_barostat = 0;
34    have_target_temp = 0;
35    have_target_pressure = 0;
36   }
37  
38 < void NPTi::moveA() {
38 > void NPTf::moveA() {
39    
40    int i,j,k;
41    int atomIndex, aMatIndex;
42    DirectionalAtom* dAtom;
43    double Tb[3];
44    double ji[3];
45 <  double rj[3];
46 <  double instaTemp, instaPress, instaVol;
47 <  double tt2, tb2;
45 >  double ri[3], vi[3], sc[3];
46 >  double instaTemp, instaVol;
47 >  double tt2, tb2, eta2ij;
48    double angle;
49 +  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
50  
51    tt2 = tauThermostat * tauThermostat;
52    tb2 = tauBarostat * tauBarostat;
53  
54    instaTemp = tStats->getTemperature();
55 <  instaPress = tStats->getPressure();
55 >  tStats->getPressureTensor(press);
56    instaVol = tStats->getVolume();
57    
58    // first evolve chi a half step
59    
60    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
56 –  
57 –  for (i = 0; i < 9; i++) {
58 –    eta[i] += dt2 * ( instaVol * (sigma[i] - targetPressure*identMat[i]))
59 –      / (NkBT*tb2));
60 – }
61  
62 +  for (i = 0; i < 3; i++ ) {
63 +    for (j = 0; j < 3; j++ ) {
64 +      if (i == j) {
65 +
66 +        eta[i][j] += dt2 * instaVol *
67 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
68 +
69 +        vScale[i][j] = eta[i][j] + chi;
70 +        
71 +      } else {
72 +        
73 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
74 +
75 +        vScale[i][j] = eta[i][j];
76 +        
77 +      }
78 +    }
79 +  }
80 +
81    for( i=0; i<nAtoms; i++ ){
82      atomIndex = i * 3;
83      aMatIndex = i * 9;
84      
85      // velocity half step
86 <    for( j=atomIndex; j<(atomIndex+3); j++ )
87 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
88 <                       - vel[j]*(chi+eta));
86 >    
87 >    vi[0] = vel[atomIndex];
88 >    vi[1] = vel[atomIndex+1];
89 >    vi[2] = vel[atomIndex+2];
90 >    
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]);
96  
97 +    vel[atomIndex] = vi[0]
98 +    vel[atomIndex+1] = vi[1];
99 +    vel[atomIndex+2] = vi[2];
100 +
101      // position whole step    
102  
103 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
104 <      rj[0] = pos[j];
105 <      rj[1] = pos[j+1];
76 <      rj[2] = pos[j+2];
103 >    ri[0] = pos[atomIndex];
104 >    ri[1] = pos[atomIndex+1];
105 >    ri[2] = pos[atomIndex+2];
106  
107 <      info->wrapVector(rj);
107 >    info->wrapVector(ri);
108  
109 <      pos[j] += dt * (vel[j] + eta*rj[0]);
81 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
82 <      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
83 <    }
109 >    info->matVecMul3( eta, ri, sc );
110  
111 <    // Scale the box after all the positions have been moved:
112 <
113 <    info->scaleBox(exp(dt*eta));
111 >    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    
115      if( atoms[i]->isDirectional() ){
116  
# Line 137 | Line 163 | void NPTi::moveA() {
163      }
164      
165    }
166 +
167 +  // Scale the box after all the positions have been moved:
168 +
169 +  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
170 +  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
171 +
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 +    }
191 +  }
192 +  
193 +  info->getBoxM(hm);
194 +  info->matMul3(hm, scaleMat, hmnew);
195 +  info->setBoxM(hmnew);
196 +  
197   }
198  
199 < void NPTi::moveB( void ){
200 <  int i,j,k;
199 > void NPTf::moveB( void ){
200 >  int i,j, k;
201    int atomIndex;
202    DirectionalAtom* dAtom;
203    double Tb[3];
204    double ji[3];
205 <  double instaTemp, instaPress, instaVol;
205 >  double vi[3], sc[3];
206 >  double instaTemp, instaVol;
207    double tt2, tb2;
208 +  double press[3][3], vScale[3][3];
209    
210    tt2 = tauThermostat * tauThermostat;
211    tb2 = tauBarostat * tauBarostat;
212  
213    instaTemp = tStats->getTemperature();
214 <  instaPress = tStats->getPressure();
214 >  tStats->getPressureTensor(press);
215    instaVol = tStats->getVolume();
216 <
216 >  
217 >  // first evolve chi a half step
218 >  
219    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
159 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
220    
221 +  for (i = 0; i < 3; i++ ) {
222 +    for (j = 0; j < 3; j++ ) {
223 +      if (i == j) {
224 +
225 +        eta[i][j] += dt2 * instaVol *
226 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
227 +
228 +        vScale[i][j] = eta[i][j] + chi;
229 +        
230 +      } else {
231 +        
232 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
233 +
234 +        vScale[i][j] = eta[i][j];
235 +        
236 +      }
237 +    }
238 +  }
239 +
240    for( i=0; i<nAtoms; i++ ){
241      atomIndex = i * 3;
242 <    
242 >
243      // velocity half step
165 –    for( j=atomIndex; j<(atomIndex+3); j++ )
166 –    for( j=atomIndex; j<(atomIndex+3); j++ )
167 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168 –                       - vel[j]*(chi+eta));
244      
245 +    vi[0] = vel[atomIndex];
246 +    vi[1] = vel[atomIndex+1];
247 +    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]);
254 +
255 +    vel[atomIndex] = vi[0]
256 +    vel[atomIndex+1] = vi[1];
257 +    vel[atomIndex+2] = vi[2];
258 +    
259      if( atoms[i]->isDirectional() ){
260        
261        dAtom = (DirectionalAtom *)atoms[i];
# Line 197 | Line 286 | void NPTi::moveB( void ){
286    }
287   }
288  
289 < int NPTi::readyCheck() {
289 > int NPTf::readyCheck() {
290  
291    // First check to see if we have a target temperature.
292    // Not having one is fatal.
293    
294    if (!have_target_temp) {
295      sprintf( painCave.errMsg,
296 <             "NPTi error: You can't use the NPTi integrator\n"
296 >             "NPTf error: You can't use the NPTf integrator\n"
297               "   without a targetTemp!\n"
298               );
299      painCave.isFatal = 1;
# Line 214 | Line 303 | int NPTi::readyCheck() {
303  
304    if (!have_target_pressure) {
305      sprintf( painCave.errMsg,
306 <             "NPTi error: You can't use the NPTi integrator\n"
306 >             "NPTf error: You can't use the NPTf integrator\n"
307               "   without a targetPressure!\n"
308               );
309      painCave.isFatal = 1;
# Line 226 | Line 315 | int NPTi::readyCheck() {
315    
316    if (!have_tau_thermostat) {
317      sprintf( painCave.errMsg,
318 <             "NPTi error: If you use the NPTi\n"
318 >             "NPTf error: If you use the NPTf\n"
319               "   integrator, you must set tauThermostat.\n");
320      painCave.isFatal = 1;
321      simError();
# Line 237 | Line 326 | int NPTi::readyCheck() {
326    
327    if (!have_tau_barostat) {
328      sprintf( painCave.errMsg,
329 <             "NPTi error: If you use the NPTi\n"
329 >             "NPTf error: If you use the NPTf\n"
330               "   integrator, you must set tauBarostat.\n");
331      painCave.isFatal = 1;
332      simError();

Diff Legend

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