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 604 by gezelter, Tue Jul 15 03:08:00 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 36 | Line 37 | void NPTim::moveA() {
37  
38   void NPTim::moveA() {
39    
40 <  int i,j,k;
40 <  int nInMol, aMatIndex;
40 >  int i, j;
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;
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[i]->getVel( vel );
86 >        myAtoms[i]->getPos( pos );
87 >        myAtoms[i]->getFrc( frc );
88  
89 <        rx = theAtoms[j]->getX();
97 <        ry = theAtoms[j]->getY();
98 <        rz = theAtoms[j]->getZ();
89 >        mass = myAtoms[i]->getMass();
90  
91 <        // velocity half step
91 >        for (j=0; j < 3; j++)
92 >          vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(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[i]->setVel( vel );
95  
96 <        // position whole step    
96 >        for (j = 0; j < 3; j++)
97 >          pos[j] += dt * (vel[j] + eta*rc[j]);
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 >        atoms[i]->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 (j=0; j < 3; j++)
115 >            ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*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    }
# Line 178 | Line 154 | void NPTi::moveB( void ){
154    info->scaleBox(exp(dt*eta));
155   }
156  
157 < void NPTi::moveB( void ){
158 <  int i,j,k;
183 <  int atomIndex;
157 > void NPTim::moveB( void ){
158 >  int i, j;
159    DirectionalAtom* dAtom;
160 <  double Tb[3];
161 <  double ji[3];
160 >  double Tb[3], ji[3];
161 >  double vel[3], frc[3];
162 >  double mass;
163 >
164    double instaTemp, instaPress, instaVol;
165    double tt2, tb2;
166 <  
166 >
167    tt2 = tauThermostat * tauThermostat;
168    tb2 = tauBarostat * tauBarostat;
169  
# Line 199 | Line 176 | void NPTi::moveB( void ){
176                   (p_convert*NkBT*tb2));
177    
178    for( i=0; i<nAtoms; i++ ){
202 –    atomIndex = i * 3;
179      
180 +    atoms[i]->getVel( vel );
181 +    atoms[i]->getFrc( frc );
182 +    
183 +    mass = atoms[i]->getMass();
184 +    
185      // velocity half step
186 <    for( j=atomIndex; j<(atomIndex+3); j++ )
187 <    for( j=atomIndex; j<(atomIndex+3); j++ )
207 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
208 <                       - vel[j]*(chi+eta));
186 >    for (j=0; j < 3; j++)
187 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
188      
189 +    atoms[i]->setVel( vel );
190 +    
191      if( atoms[i]->isDirectional() ){
192        
193        dAtom = (DirectionalAtom *)atoms[i];
194        
195 <      // get and convert the torque to body frame
195 >      // get and convert the torque to body frame      
196        
197 <      Tb[0] = dAtom->getTx();
217 <      Tb[1] = dAtom->getTy();
218 <      Tb[2] = dAtom->getTz();
219 <      
197 >      dAtom->getTrq( Tb );
198        dAtom->lab2Body( Tb );
199        
200 <      // get the angular momentum, and complete the angular momentum
223 <      // half step
200 >      // get the angular momentum, and propagate a half step
201        
202 <      ji[0] = dAtom->getJx();
226 <      ji[1] = dAtom->getJy();
227 <      ji[2] = dAtom->getJz();
202 >      dAtom->getJ( ji );
203        
204 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
205 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
231 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
204 >      for (j=0; j < 3; j++)
205 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
206        
207 <      dAtom->setJx( ji[0] );
234 <      dAtom->setJy( ji[1] );
235 <      dAtom->setJz( ji[2] );
207 >      dAtom->setJ( ji );
208      }
209    }
210   }
211  
212 < int NPTi::readyCheck() {
212 > int NPTim::readyCheck() {
213  
214    // First check to see if we have a target temperature.
215    // Not having one is fatal.
216    
217    if (!have_target_temp) {
218      sprintf( painCave.errMsg,
219 <             "NPTi error: You can't use the NPTi integrator\n"
219 >             "NPTim error: You can't use the NPTim integrator\n"
220               "   without a targetTemp!\n"
221               );
222      painCave.isFatal = 1;
# Line 254 | Line 226 | int NPTi::readyCheck() {
226  
227    if (!have_target_pressure) {
228      sprintf( painCave.errMsg,
229 <             "NPTi error: You can't use the NPTi integrator\n"
229 >             "NPTim error: You can't use the NPTim integrator\n"
230               "   without a targetPressure!\n"
231               );
232      painCave.isFatal = 1;
# Line 266 | Line 238 | int NPTi::readyCheck() {
238    
239    if (!have_tau_thermostat) {
240      sprintf( painCave.errMsg,
241 <             "NPTi error: If you use the NPTi\n"
241 >             "NPTim error: If you use the NPTim\n"
242               "   integrator, you must set tauThermostat.\n");
243      painCave.isFatal = 1;
244      simError();
# Line 277 | Line 249 | int NPTi::readyCheck() {
249    
250    if (!have_tau_barostat) {
251      sprintf( painCave.errMsg,
252 <             "NPTi error: If you use the NPTi\n"
252 >             "NPTim error: If you use the NPTim\n"
253               "   integrator, you must set tauBarostat.\n");
254      painCave.isFatal = 1;
255      simError();

Diff Legend

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