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 590 by mmeineke, Thu Jul 10 22:15:53 2003 UTC vs.
Revision 600 by gezelter, Mon Jul 14 22:38:13 2003 UTC

# Line 37 | Line 37 | void NPTf::moveA() {
37  
38   void NPTf::moveA() {
39    
40 <  int i,j,k;
41 <  int atomIndex, aMatIndex;
40 >  int i, j, k;
41    DirectionalAtom* dAtom;
42 <  double Tb[3];
43 <  double ji[3];
44 <  double ri[3], vi[3], sc[3];
45 <  double instaTemp, instaVol;
46 <  double tt2, tb2, eta2ij;
47 <  double angle;
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];
48 >  double instaTemp, instaPress, instaVol;
49 >  double tt2, tb2;
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;
# Line 62 | Line 65 | void NPTf::moveA() {
65    for (i = 0; i < 3; i++ ) {
66      for (j = 0; j < 3; j++ ) {
67        if (i == j) {
68 <
68 >        
69          eta[i][j] += dt2 * instaVol *
70            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
71 <
71 >        
72          vScale[i][j] = eta[i][j] + chi;
73          
74        } else {
# Line 79 | Line 82 | void NPTf::moveA() {
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 <    vi[0] = vel[atomIndex];
97 <    vi[1] = vel[atomIndex+1];
98 <    vi[2] = vel[atomIndex+2];
99 <    
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 >    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]   = vi[0];
98 <    vel[atomIndex+1] = vi[1];
99 <    vel[atomIndex+2] = vi[2];
101 >    atoms[i]->setVel( vel );
102  
103      // position whole step    
104  
105 <    ri[0] = pos[atomIndex];
104 <    ri[1] = pos[atomIndex+1];
105 <    ri[2] = pos[atomIndex+2];
105 >    info->wrapVector(rj);
106  
107 <    info->wrapVector(ri);
107 >    info->matVecMul3( eta, rj, sc );
108  
109 <    info->matVecMul3( eta, ri, sc );
110 <
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]);
109 >    for (j = 0; j < 3; j++ )
110 >      pos[j] += dt * (vel[j] + sc[j]);
111    
112      if( atoms[i]->isDirectional() ){
113  
# Line 118 | Line 115 | void NPTf::moveA() {
115            
116        // get and convert the torque to body frame
117        
118 <      Tb[0] = dAtom->getTx();
122 <      Tb[1] = dAtom->getTy();
123 <      Tb[2] = dAtom->getTz();
124 <      
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        
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 –      
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] );
163 <    }
164 <    
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 <
167 >      
168        // Calculate the matrix Product of the eta array (we only need
169        // the ij element right now):
170 <
170 >      
171        eta2ij = 0.0;
172        for(k=0; k<3; k++){
173          eta2ij += eta[i][k] * eta[k][j];
# Line 186 | Line 178 | void NPTf::moveA() {
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 <
181 >      
182      }
183    }
184 <  
184 >  
185    info->getBoxM(hm);
186    info->matMul3(hm, scaleMat, hmnew);
187    info->setBoxM(hmnew);
# Line 197 | Line 189 | void NPTf::moveB( void ){
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 vi[3], sc[3];
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 sc[3];
202    double press[3][3], vScale[3][3];
203    
204    tt2 = tauThermostat * tauThermostat;
# Line 238 | Line 232 | void NPTf::moveB( void ){
232    }
233  
234    for( i=0; i<nAtoms; i++ ){
241 –    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 <    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]);
245 >    for (j = 0; j < 3; j++) {
246 >      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
247 >    }
248  
249 <    vel[atomIndex]   = vi[0];
256 <    vel[atomIndex+1] = vi[1];
257 <    vel[atomIndex+2] = vi[2];
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();
266 <      Tb[1] = dAtom->getTy();
267 <      Tb[2] = dAtom->getTz();
268 <      
257 >      dAtom->getTrq( Tb );
258        dAtom->lab2Body( Tb );
259        
260 <      // get the angular momentum, and complete the angular momentum
272 <      // half step
260 >      // get the angular momentum, and propagate a half step
261        
262 <      ji[0] = dAtom->getJx();
275 <      ji[1] = dAtom->getJy();
276 <      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);
280 <      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] );
285 <    }
267 >      dAtom->setJ( ji );
268 >
269 >    }                    
270    }
271   }
272  

Diff Legend

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