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

Comparing:
branches/mmeineke/OOPSE/libmdtools/Symplectic.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Symplectic.cpp (file contents), Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

# Line 1 | Line 1
1   #include <iostream>
2   #include <cstdlib>
3  
4 + #ifdef IS_MPI
5 + #include "mpiSimulation.hpp"
6 + #include <unistd.h>
7 + #endif //is_mpi
8 +
9   #include "Integrator.hpp"
5 – #include "Thermo.hpp"
6 – #include "ReadWrite.hpp"
7 – #include "ForceFields.hpp"
10   #include "simError.h"
11  
10 – extern "C"{
11 –  
12 –  void v_constrain_a_( double &dt, int &n_atoms, double* mass,
13 –                       double* Rx, double* Ry, double* Rz,
14 –                       double* Vx, double* Vy, double* Vz,
15 –                       double* Fx, double* Fy, double* Fz,
16 –                       int &n_constrained, double *constr_sqr,
17 –                       int* constr_i, int* constr_j,
18 –                       double &box_x, double &box_y, double &box_z );
12  
13 <  void v_constrain_b_( double &dt, int &n_atoms, double* mass,
14 <                       double* Rx, double* Ry, double* Rz,
15 <                       double* Vx, double* Vy, double* Vz,
23 <                       double* Fx, double* Fy, double* Fz,
24 <                       double &Kinetic,
25 <                       int &n_constrained, double *constr_sqr,
26 <                       int* constr_i, int* constr_j,
27 <                       double &box_x, double &box_y, double &box_z );
28 < }
29 <
30 <
31 <
32 <
33 < Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff ){
34 <  entry_plug = the_entry_plug;
13 > Symplectic::Symplectic( SimInfo* theInfo, ForceFields* the_ff ){
14 >  
15 >  info = theInfo;
16    myFF = the_ff;
17    isFirst = 1;
18  
19 +  molecules = info->molecules;
20 +  nMols = info->n_mol;
21  
39 –  srInteractions = entry_plug->sr_interactions;
40 –  nSRI           =           entry_plug->n_SRI;
41 –
22    // give a little love back to the SimInfo object
23    
24 <  if( entry_plug->the_integrator != NULL ) delete entry_plug->the_integrator;
25 <  entry_plug->the_integrator = this;
24 >  if( info->the_integrator != NULL ) delete info->the_integrator;
25 >  info->the_integrator = this;
26  
27 <  // grab the masses
27 >  // check for constraints
28 >  
29 >  constrainedI = NULL;
30 >  constrainedJ = NULL;
31 >  constrainedDsqr = NULL;
32 >  nConstrained = 0;
33  
34 <  mass = new double[entry_plug->n_atoms];
35 <  for(int i = 0; i < entry_plug->n_atoms; i++){
36 <    mass[i] = entry_plug->atoms[i]->getMass();
34 >  checkConstraints();
35 > }
36 >
37 > Symplectic::~Symplectic() {
38 >  
39 >  if( nConstrained ){
40 >    delete[] constrainedI;
41 >    delete[] constrainedJ;
42 >    delete[] constrainedDsqr;
43    }
44 +  
45 + }
46  
47 <
55 <  // check for constraints
47 > void Symplectic::checkConstraints( void ){
48  
57 –  is_constrained = 0;
49  
50 +  isConstrained = 0;
51 +
52    Constraint *temp_con;
53    Constraint *dummy_plug;
54 <  temp_con = new Constraint[nSRI];
55 <  n_constrained = 0;
54 >  temp_con = new Constraint[info->n_SRI];
55 >  nConstrained = 0;
56    int constrained = 0;
57    
58 <  for(int i = 0; i < nSRI; i++){
58 >  SRI** theArray;
59 >  for(int i = 0; i < nMols; i++){
60      
61 <    constrained = srInteractions[i]->is_constrained();
62 <    
69 <    if(constrained){
61 >    theArray = (SRI**) molecules[i].getMyBonds();
62 >    for(int j=0; j<molecules[i].getNBonds(); j++){
63        
64 <      dummy_plug = srInteractions[i]->get_constraint();
65 <      temp_con[n_constrained].set_a( dummy_plug->get_a() );
66 <      temp_con[n_constrained].set_b( dummy_plug->get_b() );
67 <      temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
64 >      constrained = theArray[j]->is_constrained();
65 >      
66 >      if(constrained){
67 >        
68 >        dummy_plug = theArray[j]->get_constraint();
69 >        temp_con[nConstrained].set_a( dummy_plug->get_a() );
70 >        temp_con[nConstrained].set_b( dummy_plug->get_b() );
71 >        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
72 >        
73 >        nConstrained++;
74 >        constrained = 0;
75 >      }
76 >    }
77  
78 <      n_constrained++;
79 <      constrained = 0;
78 >    theArray = (SRI**) molecules[i].getMyBends();
79 >    for(int j=0; j<molecules[i].getNBends(); j++){
80 >      
81 >      constrained = theArray[j]->is_constrained();
82 >      
83 >      if(constrained){
84 >        
85 >        dummy_plug = theArray[j]->get_constraint();
86 >        temp_con[nConstrained].set_a( dummy_plug->get_a() );
87 >        temp_con[nConstrained].set_b( dummy_plug->get_b() );
88 >        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
89 >        
90 >        nConstrained++;
91 >        constrained = 0;
92 >      }
93      }
94 +
95 +    theArray = (SRI**) molecules[i].getMyTorsions();
96 +    for(int j=0; j<molecules[i].getNTorsions(); j++){
97 +      
98 +      constrained = theArray[j]->is_constrained();
99 +      
100 +      if(constrained){
101 +        
102 +        dummy_plug = theArray[j]->get_constraint();
103 +        temp_con[nConstrained].set_a( dummy_plug->get_a() );
104 +        temp_con[nConstrained].set_b( dummy_plug->get_b() );
105 +        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
106 +        
107 +        nConstrained++;
108 +        constrained = 0;
109 +      }
110 +    }
111    }
112  
113 <  if(n_constrained > 0){
113 >  if(nConstrained > 0){
114      
115 <    is_constrained = 1;
116 <    constrained_i = new int[n_constrained];
117 <    constrained_j = new int[n_constrained];
118 <    constrained_dsqr = new double[n_constrained];
115 >    isConstrained = 1;
116 >
117 >    if(constrainedI != NULL )    delete[] constrainedI;
118 >    if(constrainedJ != NULL )    delete[] constrainedJ;
119 >    if(constrainedDsqr != NULL ) delete[] constrainedDsqr;
120 >
121 >    constrainedI =    new int[nConstrained];
122 >    constrainedJ =    new int[nConstrained];
123 >    constrainedDsqr = new double[nConstrained];
124      
125 <    for( int i = 0; i < n_constrained; i++){
125 >    for( int i = 0; i < nConstrained; i++){
126        
127 <      /* add 1 to the index for the fortran arrays. */
128 <
129 <      constrained_i[i] = temp_con[i].get_a() + 1;
93 <      constrained_j[i] = temp_con[i].get_b() + 1;
94 <      constrained_dsqr[i] = temp_con[i].get_dsqr();
127 >      constrainedI[i] = temp_con[i].get_a();
128 >      constrainedJ[i] = temp_con[i].get_b();
129 >      constrainedDsqr[i] = temp_con[i].get_dsqr();
130      }
131    }
132    
133    delete[] temp_con;
134   }
135  
101 – Symplectic::~Symplectic() {
102 –  
103 –  if( n_constrained ){
104 –    delete[] constrained_i;
105 –    delete[] constrained_j;
106 –    delete[] constrained_dsqr;
107 –  }
108 –  
109 – }
136  
111 –
137   void Symplectic::integrate( void ){
138  
114 –  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
115 –
139    int i, j;                         // loop counters
140 <  int nAtoms = entry_plug->n_atoms; // the number of atoms
140 >  int nAtoms = info->n_atoms; // the number of atoms
141    double kE = 0.0;                  // the kinetic energy  
142    double rot_kE;
143    double trans_kE;
# Line 122 | Line 145 | void Symplectic::integrate( void ){
145    double dt2;                       // half the dt
146  
147    double vx, vy, vz;    // the velocities
148 < //  double vx2, vy2, vz2; // the square of the velocities
148 >  double vx2, vy2, vz2; // the square of the velocities
149    double rx, ry, rz;    // the postitions
150    
151    double ji[3];   // the body frame angular momentum
# Line 130 | Line 153 | void Symplectic::integrate( void ){
153    double Tb[3];   // torque in the body frame
154    double angle;   // the angle through which to rotate the rotation matrix
155    double A[3][3]; // the rotation matrix
156 +  double press[9];
157  
158    int time;
159  
160 <  double dt          = entry_plug->dt;
161 <  double runTime     = entry_plug->run_time;
162 <  double sampleTime  = entry_plug->sampleTime;
163 <  double statusTime  = entry_plug->statusTime;
164 <  double thermalTime = entry_plug->thermalTime;
160 >  double dt          = info->dt;
161 >  double runTime     = info->run_time;
162 >  double sampleTime  = info->sampleTime;
163 >  double statusTime  = info->statusTime;
164 >  double thermalTime = info->thermalTime;
165  
166    int n_loops  = (int)( runTime / dt );
167    int sample_n = (int)( sampleTime / dt );
168    int status_n = (int)( statusTime / dt );
169    int vel_n    = (int)( thermalTime / dt );
170  
171 <  int calcPot;
171 >  int calcPot, calcStress;
172 >  int isError;
173  
174 <   Thermo *tStats = new Thermo( entry_plug );
174 >  tStats   = new Thermo( info );
175 >  e_out    = new StatWriter( info );
176 >  dump_out = new DumpWriter( info );
177  
178 <  StatWriter*  e_out    = new StatWriter( entry_plug );
152 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
153 <
154 <  Atom** atoms = entry_plug->atoms;
178 >  Atom** atoms = info->atoms;
179    DirectionalAtom* dAtom;
180    dt2 = 0.5 * dt;
181  
182 <  // initialize the forces the before the first step
182 >  // initialize the forces before the first step
183  
184 <  myFF->doForces(1,0);
184 >  myFF->doForces(1,1);
185    
186 <  if( entry_plug->setTemp ){
186 >  if( info->setTemp ){
187      
188      tStats->velocitize();
189    }
# Line 169 | Line 193 | void Symplectic::integrate( void ){
193    
194    calcPot = 0;
195  
196 <  if( n_constrained ){
196 >  for( tl=0; tl<nLoops; tl++){
197  
198 <    double *Rx = new double[nAtoms];
199 <    double *Ry = new double[nAtoms];
200 <    double *Rz = new double[nAtoms];
198 >    integrateStep( calcPot, calcStress );
199 >      
200 >    time = tl + 1;
201      
202 <    double *Vx = new double[nAtoms];
203 <    double *Vy = new double[nAtoms];
204 <    double *Vz = new double[nAtoms];
202 >    if( info->setTemp ){
203 >      if( !(time % vel_n) ) tStats->velocitize();
204 >    }
205 >    if( !(time % sample_n) ) dump_out->writeDump( time * dt );
206 >    if( !((time+1) % status_n) ) {
207 >      calcPot = 1;
208 >      calcStress = 1;
209 >    }
210 >    if( !(time % status_n) ){
211 >      e_out->writeStat( time * dt );
212 >      calcPot = 0;
213 >      if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
214 >      else calcStress = 0;
215 >    }    
216 >
217 >  
218 >  }
219 >
220 >  dump_out->writeFinal();
221 >
222 >  delete dump_out;
223 >  delete e_out;
224 > }
225 >
226 >
227 > void Symplectic::moveA( void ){
228 >  
229 >  int i,j,k;
230 >  int atomIndex, aMatIndex;
231 >  DirectionalAtom* dAtom;
232 >  double Tb[3];
233 >  double ji[3];
234 >
235 >  for( i=0; i<nAtoms; i++ ){
236 >    atomIndex = i * 3;
237 >    aMatIndex = i * 9;
238      
239 <    double *Fx = new double[nAtoms];
240 <    double *Fy = new double[nAtoms];
241 <    double *Fz = new double[nAtoms];
185 <    
239 >    // velocity half step
240 >    for( j=atomIndex; j<(atomIndex+3); j++ )
241 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
242  
243 <    for( tl=0; tl < n_loops; tl++ ){
243 >    // position whole step    
244 >    for( j=atomIndex; j<(atomIndex+3); j++ )
245 >      pos[j] += dt * vel[j];
246 >
247 >  
248 >    if( atoms[i]->isDirectional() ){
249 >
250 >      dAtom = (DirectionalAtom *)atoms[i];
251 >          
252 >      // get and convert the torque to body frame
253        
254 <      for( j=0; j<nAtoms; j++ ){
254 >      Tb[0] = dAtom->getTx();
255 >      Tb[1] = dAtom->getTy();
256 >      Tb[2] = dAtom->getTz();
257 >      
258 >      dAtom->lab2Body( Tb );
259 >      
260 >      // get the angular momentum, and propagate a half step
261 >      
262 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
263 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
264 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
265 >      
266 >      // use the angular velocities to propagate the rotation matrix a
267 >      // full time step
268 >      
269 >      // rotate about the x-axis      
270 >      angle = dt2 * ji[0] / dAtom->getIxx();
271 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
272 >      
273 >      // rotate about the y-axis
274 >      angle = dt2 * ji[1] / dAtom->getIyy();
275 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
276 >      
277 >      // rotate about the z-axis
278 >      angle = dt * ji[2] / dAtom->getIzz();
279 >      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
280 >      
281 >      // rotate about the y-axis
282 >      angle = dt2 * ji[1] / dAtom->getIyy();
283 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
284 >      
285 >       // rotate about the x-axis
286 >      angle = dt2 * ji[0] / dAtom->getIxx();
287 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
288 >      
289 >      dAtom->setJx( ji[0] );
290 >      dAtom->setJy( ji[1] );
291 >      dAtom->setJz( ji[2] );
292 >    }
293 >    
294 >  }
295 > }
296  
191 –        Rx[j] = atoms[j]->getX();
192 –        Ry[j] = atoms[j]->getY();
193 –        Rz[j] = atoms[j]->getZ();
297  
298 <        Vx[j] = atoms[j]->get_vx();
299 <        Vy[j] = atoms[j]->get_vy();
300 <        Vz[j] = atoms[j]->get_vz();
298 > void Integrator::moveB( void ){
299 >  int i,j,k;
300 >  int atomIndex;
301 >  DirectionalAtom* dAtom;
302 >  double Tb[3];
303 >  double ji[3];
304  
305 <        Fx[j] = atoms[j]->getFx();
306 <        Fy[j] = atoms[j]->getFy();
201 <        Fz[j] = atoms[j]->getFz();
305 >  for( i=0; i<nAtoms; i++ ){
306 >    atomIndex = i * 3;
307  
308 <      }
309 <        
310 <      v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
206 <                      Fx, Fy, Fz,
207 <                      n_constrained, constrained_dsqr,
208 <                      constrained_i, constrained_j,
209 <                      entry_plug->box_x,
210 <                      entry_plug->box_y,
211 <                      entry_plug->box_z );
212 <      
213 <      for( j=0; j<nAtoms; j++ ){
308 >    // velocity half step
309 >    for( j=atomIndex; j<(atomIndex+3); j++ )
310 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
311  
312 <        atoms[j]->setX(Rx[j]);
216 <        atoms[j]->setY(Ry[j]);
217 <        atoms[j]->setZ(Rz[j]);
218 <        
219 <        atoms[j]->set_vx(Vx[j]);
220 <        atoms[j]->set_vy(Vy[j]);
221 <        atoms[j]->set_vz(Vz[j]);
222 <      }
223 <
224 <
225 <      for( i=0; i<nAtoms; i++ ){
226 <        if( atoms[i]->isDirectional() ){
227 <                  
228 <          dAtom = (DirectionalAtom *)atoms[i];
229 <          
230 <          // get and convert the torque to body frame
231 <          
232 <          Tb[0] = dAtom->getTx();
233 <          Tb[1] = dAtom->getTy();
234 <          Tb[2] = dAtom->getTz();
235 <          
236 <          dAtom->lab2Body( Tb );
237 <          
238 <          // get the angular momentum, and propagate a half step
239 <          
240 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
241 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
242 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
243 <          
244 <          // get the atom's rotation matrix
245 <          
246 <          A[0][0] = dAtom->getAxx();
247 <          A[0][1] = dAtom->getAxy();
248 <          A[0][2] = dAtom->getAxz();
249 <          
250 <          A[1][0] = dAtom->getAyx();
251 <          A[1][1] = dAtom->getAyy();
252 <          A[1][2] = dAtom->getAyz();
253 <          
254 <          A[2][0] = dAtom->getAzx();
255 <          A[2][1] = dAtom->getAzy();
256 <          A[2][2] = dAtom->getAzz();
257 <          
258 <          
259 <          // use the angular velocities to propagate the rotation matrix a
260 <          // full time step
261 <          
262 <          
263 <          angle = dt2 * ji[0] / dAtom->getIxx();
264 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
265 <          
266 <          angle = dt2 * ji[1] / dAtom->getIyy();
267 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
268 <          
269 <          angle = dt * ji[2] / dAtom->getIzz();
270 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
271 <          
272 <          angle = dt2 * ji[1] / dAtom->getIyy();
273 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
274 <          
275 <          angle = dt2 * ji[0] / dAtom->getIxx();
276 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
277 <          
278 <          
279 <          dAtom->setA( A );
280 <          dAtom->setJx( ji[0] );
281 <          dAtom->setJy( ji[1] );
282 <          dAtom->setJz( ji[2] );
283 <        }
284 <      }
312 >    if( atoms[i]->isDirectional() ){
313        
314 <      // calculate the forces
314 >      dAtom = (DirectionalAtom *)atoms[i];
315        
316 <      myFF->doForces(calcPot, 0);
316 >      // get and convert the torque to body frame
317        
318 <      // move b
318 >      Tb[0] = dAtom->getTx();
319 >      Tb[1] = dAtom->getTy();
320 >      Tb[2] = dAtom->getTz();
321 >      
322 >      dAtom->lab2Body( Tb );
323 >      
324 >      // get the angular momentum, and complete the angular momentum
325 >      // half step
326 >      
327 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
328 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
329 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
330 >      
331 >      jx2 = ji[0] * ji[0];
332 >      jy2 = ji[1] * ji[1];
333 >      jz2 = ji[2] * ji[2];
334 >      
335 >      dAtom->setJx( ji[0] );
336 >      dAtom->setJy( ji[1] );
337 >      dAtom->setJz( ji[2] );
338 >    }
339 >  }
340  
341 <      for( j=0; j<nAtoms; j++ ){
341 > }
342  
294 –        Rx[j] = atoms[j]->getX();
295 –        Ry[j] = atoms[j]->getY();
296 –        Rz[j] = atoms[j]->getZ();
343  
344 <        Vx[j] = atoms[j]->get_vx();
299 <        Vy[j] = atoms[j]->get_vy();
300 <        Vz[j] = atoms[j]->get_vz();
344 > void Integrator::constrainA(){
345  
346 <        Fx[j] = atoms[j]->getFx();
303 <        Fy[j] = atoms[j]->getFy();
304 <        Fz[j] = atoms[j]->getFz();
305 <      }
306 <        
307 <      v_constrain_b_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
308 <                      Fx, Fy, Fz,
309 <                      kE, n_constrained, constrained_dsqr,
310 <                      constrained_i, constrained_j,
311 <                      entry_plug->box_x,
312 <                      entry_plug->box_y,
313 <                      entry_plug->box_z );
314 <      
315 <      for( j=0; j<nAtoms; j++ ){
346 >  
347  
317 –        atoms[j]->setX(Rx[j]);
318 –        atoms[j]->setY(Ry[j]);
319 –        atoms[j]->setZ(Rz[j]);
348  
349 <        atoms[j]->set_vx(Vx[j]);
322 <        atoms[j]->set_vy(Vy[j]);
323 <        atoms[j]->set_vz(Vz[j]);
324 <      }
325 <      
326 <      for( i=0; i< nAtoms; i++ ){
349 > }
350  
328 –        if( atoms[i]->isDirectional() ){
351  
330 –          dAtom = (DirectionalAtom *)atoms[i];
331 –          
332 –          // get and convert the torque to body frame
333 –          
334 –          Tb[0] = dAtom->getTx();
335 –          Tb[1] = dAtom->getTy();
336 –          Tb[2] = dAtom->getTz();
337 –          
338 –          dAtom->lab2Body( Tb );
339 –          
340 –          // get the angular momentum, and complete the angular momentum
341 –          // half step
342 –          
343 –          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
344 –          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
345 –          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
346 –          
347 –          dAtom->setJx( ji[0] );
348 –          dAtom->setJy( ji[1] );
349 –          dAtom->setJz( ji[2] );
350 –        }
351 –      }
352 –    
353 –      time = tl + 1;
354 –      
355 –      if( entry_plug->setTemp ){
356 –        if( !(time % vel_n) ) tStats->velocitize();
357 –      }
358 –      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
359 –      if( !((time+1) % status_n) ) calcPot = 1;
360 –      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
361 –    }
362 –  }
363 –  else{
352  
365 –    for( tl=0; tl<n_loops; tl++ ){
366 –      
367 –      kE = 0.0;
368 –      rot_kE= 0.0;
369 –      trans_kE = 0.0;
370 –      
371 –      for( i=0; i<nAtoms; i++ ){
372 –        
373 –        // velocity half step
374 –        
375 –        vx = atoms[i]->get_vx() +
376 –          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
377 –        vy = atoms[i]->get_vy() +
378 –          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
379 –        vz = atoms[i]->get_vz() +
380 –          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
381 –        
382 –        // position whole step
383 –        
384 –        rx = atoms[i]->getX() + dt * vx;
385 –        ry = atoms[i]->getY() + dt * vy;
386 –        rz = atoms[i]->getZ() + dt * vz;
387 –        
388 –        atoms[i]->setX( rx );
389 –        atoms[i]->setY( ry );
390 –        atoms[i]->setZ( rz );
391 –        
392 –        atoms[i]->set_vx( vx );
393 –        atoms[i]->set_vy( vy );
394 –        atoms[i]->set_vz( vz );
395 –        
396 –        if( atoms[i]->isDirectional() ){
353  
398 –          dAtom = (DirectionalAtom *)atoms[i];
399 –          
400 –          // get and convert the torque to body frame
401 –          
402 –          Tb[0] = dAtom->getTx();
403 –          Tb[1] = dAtom->getTy();
404 –          Tb[2] = dAtom->getTz();
405 –          
406 –          dAtom->lab2Body( Tb );
407 –          
408 –          // get the angular momentum, and propagate a half step
409 –          
410 –          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
411 –          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
412 –          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
413 –          
414 –          // get the atom's rotation matrix
415 –          
416 –          A[0][0] = dAtom->getAxx();
417 –          A[0][1] = dAtom->getAxy();
418 –          A[0][2] = dAtom->getAxz();
419 –          
420 –          A[1][0] = dAtom->getAyx();
421 –          A[1][1] = dAtom->getAyy();
422 –          A[1][2] = dAtom->getAyz();
423 –          
424 –          A[2][0] = dAtom->getAzx();
425 –          A[2][1] = dAtom->getAzy();
426 –          A[2][2] = dAtom->getAzz();
427 –          
428 –          
429 –          // use the angular velocities to propagate the rotation matrix a
430 –          // full time step
431 –          
432 –          
433 –          angle = dt2 * ji[0] / dAtom->getIxx();
434 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
435 –          
436 –          angle = dt2 * ji[1] / dAtom->getIyy();
437 –          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
438 –          
439 –          angle = dt * ji[2] / dAtom->getIzz();
440 –          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
441 –          
442 –          angle = dt2 * ji[1] / dAtom->getIyy();
443 –          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
444 –          
445 –          angle = dt2 * ji[0] / dAtom->getIxx();
446 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
447 –          
448 –          
449 –          dAtom->setA( A );
450 –          dAtom->setJx( ji[0] );
451 –          dAtom->setJy( ji[1] );
452 –          dAtom->setJz( ji[2] );
453 –        }
454 –      }
455 –      
456 –      // calculate the forces
457 –      
458 –      myFF->doForces(calcPot,0);
459 –      
460 –      for( i=0; i< nAtoms; i++ ){
461 –        
462 –        // complete the velocity half step
463 –        
464 –        vx = atoms[i]->get_vx() +
465 –          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
466 –        vy = atoms[i]->get_vy() +
467 –          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
468 –        vz = atoms[i]->get_vz() +
469 –          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
470 –        
471 –        atoms[i]->set_vx( vx );
472 –        atoms[i]->set_vy( vy );
473 –        atoms[i]->set_vz( vz );
474 –        
475 – //      vx2 = vx * vx;
476 – //      vy2 = vy * vy;
477 – //      vz2 = vz * vz;
478 –        
479 –        if( atoms[i]->isDirectional() ){
354  
481 –          dAtom = (DirectionalAtom *)atoms[i];
482 –          
483 –          // get and convert the torque to body frame
484 –          
485 –          Tb[0] = dAtom->getTx();
486 –          Tb[1] = dAtom->getTy();
487 –          Tb[2] = dAtom->getTz();
488 –          
489 –          dAtom->lab2Body( Tb );
490 –          
491 –          // get the angular momentum, and complete the angular momentum
492 –          // half step
493 –          
494 –          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
495 –          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
496 –          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
497 –          
498 –          jx2 = ji[0] * ji[0];
499 –          jy2 = ji[1] * ji[1];
500 –          jz2 = ji[2] * ji[2];
501 –          
502 –          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
503 –            + (jz2 / dAtom->getIzz());
504 –          
505 –          dAtom->setJx( ji[0] );
506 –          dAtom->setJy( ji[1] );
507 –          dAtom->setJz( ji[2] );
508 –        }
509 –      }
510 –      
511 –      time = tl + 1;
512 –      
513 –      if( entry_plug->setTemp ){
514 –        if( !(time % vel_n) ) tStats->velocitize();
515 –      }
516 –      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
517 –      if( !((time+1) % status_n) ) calcPot = 1;
518 –      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
519 –    }
520 –  }
355  
522 –  dump_out->writeFinal();
356  
524 –  delete dump_out;
525 –  delete e_out;
526 – }
357  
358 +
359   void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
360                           double A[3][3] ){
361  
# Line 542 | Line 373 | void Symplectic::rotate( int axes1, int axes2, double
373  
374    for(i=0; i<3; i++){
375      for(j=0; j<3; j++){
376 <      tempA[i][j] = A[i][j];
376 >      tempA[j][i] = A[i][j];
377      }
378    }
379  
# Line 601 | Line 432 | void Symplectic::rotate( int axes1, int axes2, double
432      for(j=0; j<3; j++){
433        A[j][i] = 0.0;
434        for(k=0; k<3; k++){
435 <        A[j][i] += tempA[k][i] * rot[j][k];
435 >        A[j][i] += tempA[i][k] * rot[j][k];
436        }
437      }
438    }

Diff Legend

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