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 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC vs.
Revision 600 by gezelter, Mon Jul 14 22:38:13 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 33 | Line 37 | void NPTf::moveA() {
37  
38   void NPTf::moveA() {
39    
40 <  int i,j,k;
37 <  int atomIndex, aMatIndex;
40 >  int i, j, k;
41    DirectionalAtom* dAtom;
42 <  double Tb[3];
43 <  double ji[3];
42 >  double Tb[3], ji[3];
43 >  double A[3][3], I[3][3];
44 >  double angle, mass;
45 >  double vel[3], pos[3], frc[3];
46 >
47    double rj[3];
42 –  double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
43 –  double hm[9];
44 –  double vx, vy, vz;
45 –  double scx, scy, scz;
48    double instaTemp, instaPress, instaVol;
49    double tt2, tb2;
50 <  double angle;
51 <  double press[9];
50 >  double sc[3];
51 >  double eta2ij;
52 >  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
53  
54    tt2 = tauThermostat * tauThermostat;
55    tb2 = tauBarostat * tauBarostat;
# Line 58 | Line 61 | void NPTf::moveA() {
61    // first evolve chi a half step
62    
63    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
64 <  
65 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
66 <    (NkBT*tb2);
67 <  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
68 <  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
69 <  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
70 <  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
71 <    (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/p_convert) /
76 <    (NkBT*tb2);
77 <  
64 >
65 >  for (i = 0; i < 3; i++ ) {
66 >    for (j = 0; j < 3; j++ ) {
67 >      if (i == j) {
68 >        
69 >        eta[i][j] += dt2 * instaVol *
70 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
71 >        
72 >        vScale[i][j] = eta[i][j] + chi;
73 >        
74 >      } else {
75 >        
76 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
77 >
78 >        vScale[i][j] = eta[i][j];
79 >        
80 >      }
81 >    }
82 >  }
83 >
84    for( i=0; i<nAtoms; i++ ){
85 <    atomIndex = i * 3;
86 <    aMatIndex = i * 9;
85 >
86 >    atoms[i]->getVel( vel );
87 >    atoms[i]->getPos( pos );
88 >    atoms[i]->getFrc( frc );
89 >
90 >    mass = atoms[i]->getMass();
91      
92      // velocity half step
93 +        
94 +    info->matVecMul3( vScale, vel, sc );
95      
96 <    vx = vel[atomIndex];
97 <    vy = vel[atomIndex+1];
98 <    vz = vel[atomIndex+2];
99 <    
85 <    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;
88 <    
89 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
90 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
91 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
96 >    for (j = 0; j < 3; j++) {
97 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
98 >      rj[j] = pos[j];
99 >    }
100  
101 <    vel[atomIndex] = vx;
94 <    vel[atomIndex+1] = vy;
95 <    vel[atomIndex+2] = vz;
101 >    atoms[i]->setVel( vel );
102  
103      // position whole step    
104  
99 –    rj[0] = pos[atomIndex];
100 –    rj[1] = pos[atomIndex+1];
101 –    rj[2] = pos[atomIndex+2];
102 –
105      info->wrapVector(rj);
106  
107 <    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];
107 >    info->matVecMul3( eta, rj, sc );
108  
109 <    pos[atomIndex] += dt * (vel[atomIndex] + scx);
110 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
111 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
109 >    for (j = 0; j < 3; j++ )
110 >      pos[j] += dt * (vel[j] + sc[j]);
111    
112      if( atoms[i]->isDirectional() ){
113  
# Line 116 | Line 115 | void NPTf::moveA() {
115            
116        // get and convert the torque to body frame
117        
118 <      Tb[0] = dAtom->getTx();
120 <      Tb[1] = dAtom->getTy();
121 <      Tb[2] = dAtom->getTz();
122 <      
118 >      dAtom->getTrq( Tb );
119        dAtom->lab2Body( Tb );
120        
121        // get the angular momentum, and propagate a half step
122  
123 <      ji[0] = dAtom->getJx();
124 <      ji[1] = dAtom->getJy();
125 <      ji[2] = dAtom->getJz();
123 >      dAtom->getJ( ji );
124 >
125 >      for (j=0; j < 3; j++)
126 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
127        
131 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
132 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
133 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
134 –      
128        // use the angular velocities to propagate the rotation matrix a
129        // full time step
130 <      
130 >
131 >      dAtom->getA(A);
132 >      dAtom->getI(I);
133 >    
134        // rotate about the x-axis      
135 <      angle = dt2 * ji[0] / dAtom->getIxx();
136 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
137 <      
135 >      angle = dt2 * ji[0] / I[0][0];
136 >      this->rotate( 1, 2, angle, ji, A );
137 >
138        // rotate about the y-axis
139 <      angle = dt2 * ji[1] / dAtom->getIyy();
140 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
139 >      angle = dt2 * ji[1] / I[1][1];
140 >      this->rotate( 2, 0, angle, ji, A );
141        
142        // rotate about the z-axis
143 <      angle = dt * ji[2] / dAtom->getIzz();
144 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
143 >      angle = dt * ji[2] / I[2][2];
144 >      this->rotate( 0, 1, angle, ji, A);
145        
146        // rotate about the y-axis
147 <      angle = dt2 * ji[1] / dAtom->getIyy();
148 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
147 >      angle = dt2 * ji[1] / I[1][1];
148 >      this->rotate( 2, 0, angle, ji, A );
149        
150         // rotate about the x-axis
151 <      angle = dt2 * ji[0] / dAtom->getIxx();
152 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
151 >      angle = dt2 * ji[0] / I[0][0];
152 >      this->rotate( 1, 2, angle, ji, A );
153        
154 <      dAtom->setJx( ji[0] );
155 <      dAtom->setJy( ji[1] );
156 <      dAtom->setJz( ji[2] );
161 <    }
162 <    
154 >      dAtom->setJ( ji );
155 >      dAtom->setA( A  );    
156 >    }                    
157    }
158 <
158 >  
159    // Scale the box after all the positions have been moved:
160 <
160 >  
161    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
162    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
163 <
164 <
163 >  
164 >  
165    for(i=0; i<3; i++){
166      for(j=0; j<3; j++){
167 <      ident[i][j] = 0.0;
168 <      eta1[i][j] = eta[3*i+j];
169 <      eta2[i][j] = 0.0;
170 <      for(k=0; k<3; k++){
171 <        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
178 <      }
179 <    }
180 <    ident[i][i] = 1.0;
181 <  }
182 <
183 <  
184 <  info->getBoxM(hm);
185 <
186 <  for(i=0; i<3; i++){
187 <    for(j=0; j<3; j++){      
188 <      hmnew[i][j] = 0.0;
167 >      
168 >      // Calculate the matrix Product of the eta array (we only need
169 >      // the ij element right now):
170 >      
171 >      eta2ij = 0.0;
172        for(k=0; k<3; k++){
173 <        // 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]);
173 >        eta2ij += eta[i][k] * eta[k][j];
174        }
175 +      
176 +      scaleMat[i][j] = 0.0;
177 +      // identity matrix (see above):
178 +      if (i == j) scaleMat[i][j] = 1.0;
179 +      // Taylor expansion for the exponential truncated at second order:
180 +      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
181 +      
182      }
183    }
184    
185 <  for (i = 0; i < 3; i++) {
186 <    for (j = 0; j < 3; j++) {
187 <      // remember that hmat has transpose ordering for Fortran compat:
201 <      hm[3*j + i] = hmnew[i][j];
202 <    }
203 <  }
204 <
205 <  info->setBoxM(hm);
185 >  info->getBoxM(hm);
186 >  info->matMul3(hm, scaleMat, hmnew);
187 >  info->setBoxM(hmnew);
188    
189   }
190  
191   void NPTf::moveB( void ){
192 <  int i,j,k;
193 <  int atomIndex;
192 >
193 >  int i, j;
194    DirectionalAtom* dAtom;
195 <  double Tb[3];
196 <  double ji[3];
197 <  double press[9];
198 <  double instaTemp, instaVol;
195 >  double Tb[3], ji[3];
196 >  double vel[3], frc[3];
197 >  double mass;
198 >
199 >  double instaTemp, instaPress, instaVol;
200    double tt2, tb2;
201 <  double vx, vy, vz;
202 <  double scx, scy, scz;
220 <  const double p_convert = 1.63882576e8;
201 >  double sc[3];
202 >  double press[3][3], vScale[3][3];
203    
204    tt2 = tauThermostat * tauThermostat;
205    tb2 = tauBarostat * tauBarostat;
# Line 230 | Line 212 | void NPTf::moveB( void ){
212    
213    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
214    
215 <  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
216 <    (NkBT*tb2);
217 <  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);
215 >  for (i = 0; i < 3; i++ ) {
216 >    for (j = 0; j < 3; j++ ) {
217 >      if (i == j) {
218  
219 +        eta[i][j] += dt2 * instaVol *
220 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
221 +
222 +        vScale[i][j] = eta[i][j] + chi;
223 +        
224 +      } else {
225 +        
226 +        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
227 +
228 +        vScale[i][j] = eta[i][j];
229 +        
230 +      }
231 +    }
232 +  }
233 +
234    for( i=0; i<nAtoms; i++ ){
247 –    atomIndex = i * 3;
235  
236 +    atoms[i]->getVel( vel );
237 +    atoms[i]->getFrc( frc );
238 +
239 +    mass = atoms[i]->getMass();
240 +    
241      // velocity half step
242 +        
243 +    info->matVecMul3( vScale, vel, sc );
244      
245 <    vx = vel[atomIndex];
246 <    vy = vel[atomIndex+1];
247 <    vz = vel[atomIndex+2];
254 <    
255 <    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;
258 <    
259 <    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
260 <    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
261 <    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
245 >    for (j = 0; j < 3; j++) {
246 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
247 >    }
248  
249 <    vel[atomIndex] = vx;
264 <    vel[atomIndex+1] = vy;
265 <    vel[atomIndex+2] = vz;
249 >    atoms[i]->setVel( vel );
250      
251      if( atoms[i]->isDirectional() ){
252 <      
252 >
253        dAtom = (DirectionalAtom *)atoms[i];
254 <      
254 >          
255        // get and convert the torque to body frame
256        
257 <      Tb[0] = dAtom->getTx();
274 <      Tb[1] = dAtom->getTy();
275 <      Tb[2] = dAtom->getTz();
276 <      
257 >      dAtom->getTrq( Tb );
258        dAtom->lab2Body( Tb );
259        
260 <      // get the angular momentum, and complete the angular momentum
280 <      // half step
260 >      // get the angular momentum, and propagate a half step
261        
262 <      ji[0] = dAtom->getJx();
283 <      ji[1] = dAtom->getJy();
284 <      ji[2] = dAtom->getJz();
262 >      dAtom->getJ( ji );
263        
264 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
265 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
288 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
264 >      for (j=0; j < 3; j++)
265 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
266        
267 <      dAtom->setJx( ji[0] );
268 <      dAtom->setJy( ji[1] );
269 <      dAtom->setJz( ji[2] );
293 <    }
267 >      dAtom->setJ( ji );
268 >
269 >    }                    
270    }
271   }
272  

Diff Legend

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