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 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC vs.
Revision 588 by gezelter, Thu Jul 10 17:10:56 2003 UTC

# Line 22 | Line 22 | NPTf::NPTf ( 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;
# Line 38 | Line 42 | void NPTf::moveA() {
42    DirectionalAtom* dAtom;
43    double Tb[3];
44    double ji[3];
45 <  double rj[3];
46 <  double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
47 <  double hm[9];
44 <  double vx, vy, vz;
45 <  double scx, scy, scz;
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[9];
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;
# Line 58 | Line 58 | void NPTf::moveA() {
58    // first evolve chi a half step
59    
60    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
61 <  
62 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
63 <    (NkBT*tb2);
64 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
65 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
66 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
67 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
68 <    (NkBT*tb2);
69 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
70 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
71 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
72 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
73 <    (NkBT*tb2);
74 <  
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      
87 <    vx = vel[atomIndex];
88 <    vy = vel[atomIndex+1];
89 <    vz = vel[atomIndex+2];
87 >    vi[0] = vel[atomIndex];
88 >    vi[1] = vel[atomIndex+1];
89 >    vi[2] = vel[atomIndex+2];
90      
91 <    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
86 <    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
87 <    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
91 >    info->matVecMul3( vScale, vi, sc );
92      
93 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
94 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
95 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
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] = vx;
98 <    vel[atomIndex+1] = vy;
99 <    vel[atomIndex+2] = vz;
97 >    vel[atomIndex] = vi[0]
98 >    vel[atomIndex+1] = vi[1];
99 >    vel[atomIndex+2] = vi[2];
100  
101      // position whole step    
102  
103 <    rj[0] = pos[atomIndex];
104 <    rj[1] = pos[atomIndex+1];
105 <    rj[2] = pos[atomIndex+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 <    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
106 <    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
107 <    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
109 >    info->matVecMul3( eta, ri, 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 >    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 170 | Line 172 | void NPTf::moveA() {
172  
173    for(i=0; i<3; i++){
174      for(j=0; j<3; j++){
175 <      ident[i][j] = 0.0;
176 <      eta1[i][j] = eta[3*i+j];
177 <      eta2[i][j] = 0.0;
178 <      for(k=0; k<3; k++){
179 <        eta2[i][j] += eta[3*i+k] * eta[3*k+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      }
180 –    ident[i][i] = 1.0;
191    }
182 –
192    
193    info->getBoxM(hm);
194 <
195 <  for(i=0; i<3; i++){
187 <    for(j=0; j<3; j++){      
188 <      hmnew[i][j] = 0.0;
189 <      for(k=0; k<3; k++){
190 <        // remember that hmat has transpose ordering for Fortran compat:
191 <        hmnew[i][j] += hm[3*k+i] * (ident[k][j]
192 <                                    + dt * eta1[k][j]
193 <                                    + 0.5 * dt * dt * eta2[k][j]);
194 <      }
195 <    }
196 <  }
194 >  info->matMul3(hm, scaleMat, hmnew);
195 >  info->setBoxM(hmnew);
196    
198 –  for (i = 0; i < 3; i++) {
199 –    for (j = 0; j < 3; j++) {
200 –      // remember that hmat has transpose ordering for Fortran compat:
201 –      hm[3*j + i] = hmnew[i][j];
202 –    }
203 –  }
204 –
205 –  info->setBoxM(hm);
206 –  
197   }
198  
199   void NPTf::moveB( void ){
200 <  int i,j,k;
200 >  int i,j, k;
201    int atomIndex;
202    DirectionalAtom* dAtom;
203    double Tb[3];
204    double ji[3];
205 <  double press[9];
205 >  double vi[3], sc[3];
206    double instaTemp, instaVol;
207    double tt2, tb2;
208 <  double vx, vy, vz;
219 <  double scx, scy, scz;
220 <  const double p_convert = 1.63882576e8;
208 >  double press[3][3], vScale[3][3];
209    
210    tt2 = tauThermostat * tauThermostat;
211    tb2 = tauBarostat * tauBarostat;
# Line 230 | Line 218 | void NPTf::moveB( void ){
218    
219    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
220    
221 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
222 <    (NkBT*tb2);
223 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
236 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
237 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
238 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
239 <    (NkBT*tb2);
240 <  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
241 <  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
242 <  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
243 <  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
244 <    (NkBT*tb2);
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  
243      // velocity half step
244      
245 <    vx = vel[atomIndex];
246 <    vy = vel[atomIndex+1];
247 <    vz = vel[atomIndex+2];
245 >    vi[0] = vel[atomIndex];
246 >    vi[1] = vel[atomIndex+1];
247 >    vi[2] = vel[atomIndex+2];
248      
249 <    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
256 <    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
257 <    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
249 >    info->matVecMul3( vScale, vi, sc );
250      
251 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
252 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
253 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
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] = vx;
256 <    vel[atomIndex+1] = vy;
257 <    vel[atomIndex+2] = vz;
255 >    vel[atomIndex] = vi[0]
256 >    vel[atomIndex+1] = vi[1];
257 >    vel[atomIndex+2] = vi[2];
258      
259      if( atoms[i]->isDirectional() ){
260        

Diff Legend

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