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 658 by tim, Thu Jul 31 15:35:07 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;
# Line 34 | Line 35 | NPTim::NPTim ( SimInfo *theInfo, ForceFields* the_ff):
35    have_target_pressure = 0;
36   }
37  
38 < void NPTim::moveA() {
38 > template<typename T> void NPTim<T>::moveA() {
39    
40 <  int i,j,k;
40 <  int nInMol, aMatIndex;
40 >  int i, j, k;
41    DirectionalAtom* dAtom;
42 <  Atom** theAtoms;
43 <  Molecule** myMols;
44 <  double Tb[3];
45 <  double ji[3];
46 <  double rc[3];
47 <  double mass;
48 <  double rx, ry, rz, vx, vy, vz, fx, fy, fz;
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;
51 <  double angle;
49 >  double tt2, tb2, scaleFactor;
50  
51 +  int nInMol;
52 +  double rc[3];
53 +
54    nMols = info->n_mol;
55 <  myMols = info->molecules;
55 >  myMolecules = info->molecules;
56  
57    tt2 = tauThermostat * tauThermostat;
58    tb2 = tauBarostat * tauBarostat;
# Line 68 | Line 69 | void NPTim::moveA() {
69  
70    for( i = 0; i < nMols; i++) {
71  
72 <    myMols[i].getCOM(rc);
72 >    myMolecules[i].getCOM(rc);
73  
74 <    nInMol = myMols[i]->getNAtoms();
75 <    theAtoms = myMols[i]->getMyAtoms();
74 >    nInMol = myMolecules[i].getNAtoms();
75 >    myAtoms = myMolecules[i].getMyAtoms();
76  
77      // find the minimum image coordinates of the molecular centers of mass:
78  
# Line 79 | Line 80 | void NPTim::moveA() {
80  
81      for (j = 0; j < nInMol; j++) {
82  
83 <      if(theAtoms[j] != NULL) {
83 >      if(myAtoms[j] != NULL) {
84  
85 <        aMatIndex = 9 * theAtoms[j]->getIndex();
86 <        
87 <        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();
85 >        myAtoms[j]->getVel( vel );
86 >        myAtoms[j]->getPos( pos );
87 >        myAtoms[j]->getFrc( frc );
88  
89 <        rx = theAtoms[j]->getX();
97 <        ry = theAtoms[j]->getY();
98 <        rz = theAtoms[j]->getZ();
89 >        mass = myAtoms[j]->getMass();
90  
91 <        // velocity half step
91 >        for (k=0; k < 3; k++)
92 >          vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));
93  
94 <        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));
94 >        myAtoms[j]->setVel( vel );
95  
96 <        // position whole step    
96 >        for (k = 0; k < 3; k++)
97 >          pos[k] += dt * (vel[k] + eta*rc[k]);
98  
99 <        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);
99 >        myAtoms[j]->setPos( pos );
100  
101 <        theAtoms[j]->setX(rx);        
117 <        theAtoms[j]->setY(ry);
118 <        theAtoms[j]->setZ(rz);
101 >        if( myAtoms[j]->isDirectional() ){
102  
103 <        if( theAtoms[j]->isDirectional() ){
121 <
122 <          dAtom = (DirectionalAtom *)theAtoms[j];
103 >          dAtom = (DirectionalAtom *)myAtoms[j];
104            
105            // get and convert the torque to body frame
106 <          
107 <          Tb[0] = dAtom->getTx();
127 <          Tb[1] = dAtom->getTy();
128 <          Tb[2] = dAtom->getTz();
129 <          
106 >      
107 >          dAtom->getTrq( Tb );
108            dAtom->lab2Body( Tb );
109 <          
109 >      
110            // get the angular momentum, and propagate a half step
111            
112 <          ji[0] = dAtom->getJx();
113 <          ji[1] = dAtom->getJy();
114 <          ji[2] = dAtom->getJz();
112 >          dAtom->getJ( ji );
113 >
114 >          for (k=0; k < 3; k++)
115 >            ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
116            
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 –          
117            // use the angular velocities to propagate the rotation matrix a
118            // full time step
119            
120 +          dAtom->getA(A);
121 +          dAtom->getI(I);
122 +          
123            // rotate about the x-axis      
124 <          angle = dt2 * ji[0] / dAtom->getIxx();
125 <          this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
124 >          angle = dt2 * ji[0] / I[0][0];
125 >          this->rotate( 1, 2, angle, ji, A );
126            
127            // rotate about the y-axis
128 <          angle = dt2 * ji[1] / dAtom->getIyy();
129 <          this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
128 >          angle = dt2 * ji[1] / I[1][1];
129 >          this->rotate( 2, 0, angle, ji, A );
130            
131            // rotate about the z-axis
132 <          angle = dt * ji[2] / dAtom->getIzz();
133 <          this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
132 >          angle = dt * ji[2] / I[2][2];
133 >          this->rotate( 0, 1, angle, ji, A);
134            
135            // rotate about the y-axis
136 <          angle = dt2 * ji[1] / dAtom->getIyy();
137 <          this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
136 >          angle = dt2 * ji[1] / I[1][1];
137 >          this->rotate( 2, 0, angle, ji, A );
138            
139            // rotate about the x-axis
140 <          angle = dt2 * ji[0] / dAtom->getIxx();
141 <          this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
140 >          angle = dt2 * ji[0] / I[0][0];
141 >          this->rotate( 1, 2, angle, ji, A );
142            
143 <          dAtom->setJx( ji[0] );
144 <          dAtom->setJy( ji[1] );
145 <          dAtom->setJz( ji[2] );
168 <        }
169 <        
143 >          dAtom->setJ( ji );
144 >          dAtom->setA( A  );    
145 >        }                        
146        }
147      }
148    }
149 +
150    // Scale the box after all the positions have been moved:
151    
152 <  cerr << "eta = " << eta
153 <       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";
154 <  
155 <  info->scaleBox(exp(dt*eta));
152 >  scaleFactor = exp(dt*eta);
153 >
154 >  if (scaleFactor > 1.1 || scaleFactor < 0.9) {
155 >    sprintf( painCave.errMsg,
156 >             "NPTi error: Attempting a Box scaling of more than 10 percent"
157 >             " check your tauBarostat, as it is probably too small!\n"
158 >             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
159 >             );
160 >    painCave.isFatal = 1;
161 >    simError();
162 >  } else {        
163 >    info->scaleBox(exp(dt*eta));      
164 >  }
165   }
166  
167 < void NPTi::moveB( void ){
168 <  int i,j,k;
183 <  int atomIndex;
167 > template<typename T> void NPTim<T>::moveB( void ){
168 >  int i, j;
169    DirectionalAtom* dAtom;
170 <  double Tb[3];
171 <  double ji[3];
170 >  double Tb[3], ji[3];
171 >  double vel[3], frc[3];
172 >  double mass;
173 >
174    double instaTemp, instaPress, instaVol;
175    double tt2, tb2;
176 <  
176 >
177    tt2 = tauThermostat * tauThermostat;
178    tb2 = tauBarostat * tauBarostat;
179  
# Line 199 | Line 186 | void NPTi::moveB( void ){
186                   (p_convert*NkBT*tb2));
187    
188    for( i=0; i<nAtoms; i++ ){
202 –    atomIndex = i * 3;
189      
190 +    atoms[i]->getVel( vel );
191 +    atoms[i]->getFrc( frc );
192 +    
193 +    mass = atoms[i]->getMass();
194 +    
195      // velocity half step
196 <    for( j=atomIndex; j<(atomIndex+3); j++ )
197 <    for( j=atomIndex; j<(atomIndex+3); j++ )
207 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
208 <                       - vel[j]*(chi+eta));
196 >    for (j=0; j < 3; j++)
197 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
198      
199 +    atoms[i]->setVel( vel );
200 +    
201      if( atoms[i]->isDirectional() ){
202        
203        dAtom = (DirectionalAtom *)atoms[i];
204        
205 <      // get and convert the torque to body frame
205 >      // get and convert the torque to body frame      
206        
207 <      Tb[0] = dAtom->getTx();
217 <      Tb[1] = dAtom->getTy();
218 <      Tb[2] = dAtom->getTz();
219 <      
207 >      dAtom->getTrq( Tb );
208        dAtom->lab2Body( Tb );
209        
210 <      // get the angular momentum, and complete the angular momentum
223 <      // half step
210 >      // get the angular momentum, and propagate a half step
211        
212 <      ji[0] = dAtom->getJx();
226 <      ji[1] = dAtom->getJy();
227 <      ji[2] = dAtom->getJz();
212 >      dAtom->getJ( ji );
213        
214 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
215 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
231 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
214 >      for (j=0; j < 3; j++)
215 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
216        
217 <      dAtom->setJx( ji[0] );
234 <      dAtom->setJy( ji[1] );
235 <      dAtom->setJz( ji[2] );
217 >      dAtom->setJ( ji );
218      }
219    }
220   }
221  
222 < int NPTi::readyCheck() {
222 > template<typename T> int NPTim<T>::readyCheck() {
223 >
224 >  //check parent's readyCheck() first
225 >  if (T::readyCheck() == -1)
226 >    return -1;
227  
228    // First check to see if we have a target temperature.
229    // Not having one is fatal.
230    
231    if (!have_target_temp) {
232      sprintf( painCave.errMsg,
233 <             "NPTi error: You can't use the NPTi integrator\n"
233 >             "NPTim error: You can't use the NPTim integrator\n"
234               "   without a targetTemp!\n"
235               );
236      painCave.isFatal = 1;
# Line 254 | Line 240 | int NPTi::readyCheck() {
240  
241    if (!have_target_pressure) {
242      sprintf( painCave.errMsg,
243 <             "NPTi error: You can't use the NPTi integrator\n"
243 >             "NPTim error: You can't use the NPTim integrator\n"
244               "   without a targetPressure!\n"
245               );
246      painCave.isFatal = 1;
# Line 266 | Line 252 | int NPTi::readyCheck() {
252    
253    if (!have_tau_thermostat) {
254      sprintf( painCave.errMsg,
255 <             "NPTi error: If you use the NPTi\n"
255 >             "NPTim error: If you use the NPTim\n"
256               "   integrator, you must set tauThermostat.\n");
257      painCave.isFatal = 1;
258      simError();
# Line 277 | Line 263 | int NPTi::readyCheck() {
263    
264    if (!have_tau_barostat) {
265      sprintf( painCave.errMsg,
266 <             "NPTi error: If you use the NPTi\n"
266 >             "NPTim error: If you use the NPTim\n"
267               "   integrator, you must set tauBarostat.\n");
268      painCave.isFatal = 1;
269      simError();

Diff Legend

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