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

Comparing trunk/OOPSE/libmdtools/NPTim.cpp (file contents):
Revision 596 by gezelter, Mon Jul 14 15:04:55 2003 UTC vs.
Revision 763 by tim, Mon Sep 15 16:52:02 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2   #include "Atom.hpp"
3 + #include "Molecule.hpp"
4   #include "SRI.hpp"
5   #include "AbstractClasses.hpp"
6   #include "SimInfo.hpp"
# Line 23 | Line 24
24   //  The NPTim variant scales the molecular center-of-mass coordinates
25   //  instead of the atomic coordinates
26  
27 < NPTim::NPTim ( SimInfo *theInfo, ForceFields* the_ff):
28 <  Integrator( theInfo, the_ff )
27 > template<typename T> NPTim<T>::NPTim ( SimInfo *theInfo, ForceFields* the_ff):
28 >  T( theInfo, the_ff )
29   {
30    chi = 0.0;
31    eta = 0.0;
32 +  integralOfChidt = 0.0;
33    have_tau_thermostat = 0;
34    have_tau_barostat = 0;
35    have_target_temp = 0;
36    have_target_pressure = 0;
37   }
38  
39 < void NPTim::moveA() {
39 > template<typename T> void NPTim<T>::moveA() {
40    
41 <  int i,j,k;
40 <  int nInMol, aMatIndex;
41 >  int i, j, k;
42    DirectionalAtom* dAtom;
43 <  Atom** theAtoms;
44 <  Molecule** myMols;
45 <  double Tb[3];
46 <  double ji[3];
47 <  double rc[3];
48 <  double mass;
48 <  double rx, ry, rz, vx, vy, vz, fx, fy, fz;
43 >  double Tb[3], ji[3];
44 >  double A[3][3], I[3][3];
45 >  double angle, mass;
46 >  double vel[3], pos[3], frc[3];
47 >
48 >  double rj[3];
49    double instaTemp, instaPress, instaVol;
50 <  double tt2, tb2;
51 <  double angle;
50 >  double tt2, tb2, scaleFactor;
51  
52 +  int nInMol;
53 +  double rc[3];
54 +
55    nMols = info->n_mol;
56 <  myMols = info->molecules;
56 >  myMolecules = info->molecules;
57  
58    tt2 = tauThermostat * tauThermostat;
59    tb2 = tauBarostat * tauBarostat;
# Line 68 | Line 70 | void NPTim::moveA() {
70  
71    for( i = 0; i < nMols; i++) {
72  
73 <    myMols[i].getCOM(rc);
73 >    myMolecules[i].getCOM(rc);
74  
75 <    nInMol = myMols[i]->getNAtoms();
76 <    theAtoms = myMols[i]->getMyAtoms();
75 >    nInMol = myMolecules[i].getNAtoms();
76 >    myAtoms = myMolecules[i].getMyAtoms();
77  
78      // find the minimum image coordinates of the molecular centers of mass:
79  
# Line 79 | Line 81 | void NPTim::moveA() {
81  
82      for (j = 0; j < nInMol; j++) {
83  
84 <      if(theAtoms[j] != NULL) {
84 >      if(myAtoms[j] != NULL) {
85  
86 <        aMatIndex = 9 * theAtoms[j]->getIndex();
87 <        
88 <        mass = theAtoms[j]->getMass();
87 <        
88 <        vx = theAtoms[j]->get_vx();
89 <        vy = theAtoms[j]->get_vy();
90 <        vz = theAtoms[j]->get_vz();
91 <        
92 <        fx = theAtoms[j]->getFx();
93 <        fy = theAtoms[j]->getFy();
94 <        fz = theAtoms[j]->getFz();
86 >        myAtoms[j]->getVel( vel );
87 >        myAtoms[j]->getPos( pos );
88 >        myAtoms[j]->getFrc( frc );
89  
90 <        rx = theAtoms[j]->getX();
97 <        ry = theAtoms[j]->getY();
98 <        rz = theAtoms[j]->getZ();
90 >        mass = myAtoms[j]->getMass();
91  
92 <        // velocity half step
92 >        for (k=0; k < 3; k++)
93 >          vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));
94  
95 <        vx += dt2 * ((fx / mass)*eConvert - vx*(chi+eta));
103 <        vy += dt2 * ((fy / mass)*eConvert - vy*(chi+eta));
104 <        vz += dt2 * ((fz / mass)*eConvert - vz*(chi+eta));
95 >        myAtoms[j]->setVel( vel );
96  
97 <        // position whole step    
97 >        for (k = 0; k < 3; k++)
98 >          pos[k] += dt * (vel[k] + eta*rc[k]);
99  
100 <        rx += dt*(vx + eta*rc[0]);
109 <        ry += dt*(vy + eta*rc[1]);
110 <        rz += dt*(vz + eta*rc[2]);
111 <        
112 <        theAtoms[j]->set_vx(vx);        
113 <        theAtoms[j]->set_vy(vy);
114 <        theAtoms[j]->set_vz(vz);
100 >        myAtoms[j]->setPos( pos );
101  
102 <        theAtoms[j]->setX(rx);        
117 <        theAtoms[j]->setY(ry);
118 <        theAtoms[j]->setZ(rz);
102 >        if( myAtoms[j]->isDirectional() ){
103  
104 <        if( theAtoms[j]->isDirectional() ){
121 <
122 <          dAtom = (DirectionalAtom *)theAtoms[j];
104 >          dAtom = (DirectionalAtom *)myAtoms[j];
105            
106            // get and convert the torque to body frame
107 <          
108 <          Tb[0] = dAtom->getTx();
127 <          Tb[1] = dAtom->getTy();
128 <          Tb[2] = dAtom->getTz();
129 <          
107 >      
108 >          dAtom->getTrq( Tb );
109            dAtom->lab2Body( Tb );
110 <          
110 >      
111            // get the angular momentum, and propagate a half step
112            
113 <          ji[0] = dAtom->getJx();
114 <          ji[1] = dAtom->getJy();
115 <          ji[2] = dAtom->getJz();
113 >          dAtom->getJ( ji );
114 >
115 >          for (k=0; k < 3; k++)
116 >            ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
117            
138 –          ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
139 –          ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
140 –          ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
141 –          
118            // use the angular velocities to propagate the rotation matrix a
119            // full time step
120            
121 +          dAtom->getA(A);
122 +          dAtom->getI(I);
123 +          
124            // rotate about the x-axis      
125 <          angle = dt2 * ji[0] / dAtom->getIxx();
126 <          this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
125 >          angle = dt2 * ji[0] / I[0][0];
126 >          this->rotate( 1, 2, angle, ji, A );
127            
128            // rotate about the y-axis
129 <          angle = dt2 * ji[1] / dAtom->getIyy();
130 <          this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
129 >          angle = dt2 * ji[1] / I[1][1];
130 >          this->rotate( 2, 0, angle, ji, A );
131            
132            // rotate about the z-axis
133 <          angle = dt * ji[2] / dAtom->getIzz();
134 <          this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
133 >          angle = dt * ji[2] / I[2][2];
134 >          this->rotate( 0, 1, angle, ji, A);
135            
136            // rotate about the y-axis
137 <          angle = dt2 * ji[1] / dAtom->getIyy();
138 <          this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
137 >          angle = dt2 * ji[1] / I[1][1];
138 >          this->rotate( 2, 0, angle, ji, A );
139            
140            // rotate about the x-axis
141 <          angle = dt2 * ji[0] / dAtom->getIxx();
142 <          this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
141 >          angle = dt2 * ji[0] / I[0][0];
142 >          this->rotate( 1, 2, angle, ji, A );
143            
144 <          dAtom->setJx( ji[0] );
145 <          dAtom->setJy( ji[1] );
146 <          dAtom->setJz( ji[2] );
168 <        }
169 <        
144 >          dAtom->setJ( ji );
145 >          dAtom->setA( A  );    
146 >        }                        
147        }
148      }
149    }
150 +
151    // Scale the box after all the positions have been moved:
152    
153 <  cerr << "eta = " << eta
154 <       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";
155 <  
156 <  info->scaleBox(exp(dt*eta));
153 >  scaleFactor = exp(dt*eta);
154 >
155 >  if (scaleFactor > 1.1 || scaleFactor < 0.9) {
156 >    sprintf( painCave.errMsg,
157 >             "NPTi error: Attempting a Box scaling of more than 10 percent"
158 >             " check your tauBarostat, as it is probably too small!\n"
159 >             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
160 >             );
161 >    painCave.isFatal = 1;
162 >    simError();
163 >  } else {        
164 >    info->scaleBox(exp(dt*eta));      
165 >  }
166   }
167  
168 < void NPTi::moveB( void ){
169 <  int i,j,k;
183 <  int atomIndex;
168 > template<typename T> void NPTim<T>::moveB( void ){
169 >  int i, j;
170    DirectionalAtom* dAtom;
171 <  double Tb[3];
172 <  double ji[3];
171 >  double Tb[3], ji[3];
172 >  double vel[3], frc[3];
173 >  double mass;
174 >
175    double instaTemp, instaPress, instaVol;
176    double tt2, tb2;
177 <  
177 >
178    tt2 = tauThermostat * tauThermostat;
179    tb2 = tauBarostat * tauBarostat;
180  
# Line 199 | Line 187 | void NPTi::moveB( void ){
187                   (p_convert*NkBT*tb2));
188    
189    for( i=0; i<nAtoms; i++ ){
202 –    atomIndex = i * 3;
190      
191 +    atoms[i]->getVel( vel );
192 +    atoms[i]->getFrc( frc );
193 +    
194 +    mass = atoms[i]->getMass();
195 +    
196      // velocity half step
197 <    for( j=atomIndex; j<(atomIndex+3); j++ )
198 <    for( j=atomIndex; j<(atomIndex+3); j++ )
207 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
208 <                       - vel[j]*(chi+eta));
197 >    for (j=0; j < 3; j++)
198 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
199      
200 +    atoms[i]->setVel( vel );
201 +    
202      if( atoms[i]->isDirectional() ){
203        
204        dAtom = (DirectionalAtom *)atoms[i];
205        
206 <      // get and convert the torque to body frame
206 >      // get and convert the torque to body frame      
207        
208 <      Tb[0] = dAtom->getTx();
217 <      Tb[1] = dAtom->getTy();
218 <      Tb[2] = dAtom->getTz();
219 <      
208 >      dAtom->getTrq( Tb );
209        dAtom->lab2Body( Tb );
210        
211 <      // get the angular momentum, and complete the angular momentum
223 <      // half step
211 >      // get the angular momentum, and propagate a half step
212        
213 <      ji[0] = dAtom->getJx();
226 <      ji[1] = dAtom->getJy();
227 <      ji[2] = dAtom->getJz();
213 >      dAtom->getJ( ji );
214        
215 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
216 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
231 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
215 >      for (j=0; j < 3; j++)
216 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
217        
218 <      dAtom->setJx( ji[0] );
234 <      dAtom->setJy( ji[1] );
235 <      dAtom->setJz( ji[2] );
218 >      dAtom->setJ( ji );
219      }
220    }
221   }
222  
223 < int NPTi::readyCheck() {
223 > template<typename T> void NPTim<T>::resetIntegrator() {
224 >  chi = 0.0;
225 >  eta = 0.0;
226 > }
227 >
228 > template<typename T> int NPTim<T>::readyCheck() {
229 >
230 >  //check parent's readyCheck() first
231 >  if (T::readyCheck() == -1)
232 >    return -1;
233  
234    // First check to see if we have a target temperature.
235    // Not having one is fatal.
236    
237    if (!have_target_temp) {
238      sprintf( painCave.errMsg,
239 <             "NPTi error: You can't use the NPTi integrator\n"
239 >             "NPTim error: You can't use the NPTim integrator\n"
240               "   without a targetTemp!\n"
241               );
242      painCave.isFatal = 1;
# Line 254 | Line 246 | int NPTi::readyCheck() {
246  
247    if (!have_target_pressure) {
248      sprintf( painCave.errMsg,
249 <             "NPTi error: You can't use the NPTi integrator\n"
249 >             "NPTim error: You can't use the NPTim integrator\n"
250               "   without a targetPressure!\n"
251               );
252      painCave.isFatal = 1;
# Line 266 | Line 258 | int NPTi::readyCheck() {
258    
259    if (!have_tau_thermostat) {
260      sprintf( painCave.errMsg,
261 <             "NPTi error: If you use the NPTi\n"
261 >             "NPTim error: If you use the NPTim\n"
262               "   integrator, you must set tauThermostat.\n");
263      painCave.isFatal = 1;
264      simError();
# Line 277 | Line 269 | int NPTi::readyCheck() {
269    
270    if (!have_tau_barostat) {
271      sprintf( painCave.errMsg,
272 <             "NPTi error: If you use the NPTi\n"
272 >             "NPTim error: If you use the NPTim\n"
273               "   integrator, you must set tauBarostat.\n");
274      painCave.isFatal = 1;
275      simError();
# Line 290 | Line 282 | int NPTi::readyCheck() {
282  
283    return 1;
284   }
285 +
286 + template<typename T> double NPTim<T>::getConservedQuantity(void){
287 +
288 +  double conservedQuantity;
289 +  double tb2;
290 +  double eta2;  
291 +
292 +
293 +  //HNVE
294 +  conservedQuantity = tStats->getTotalE();
295 +
296 +  //HNVT
297 +  conservedQuantity += (info->getNDF() * kB * targetTemp *
298 +                        (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2))/ eConvert ;
299 +
300 +  //HNPT
301 +  tb2 = tauBarostat *tauBarostat;
302 +  eta2 = eta * eta;
303 +  
304 +  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
305 +                        3*NkBT/2 * tb2 * eta2) / eConvert;
306 +  
307 +  return conservedQuantity;
308 + }

Diff Legend

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