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

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 477 by gezelter, Tue Apr 8 14:34:30 2003 UTC vs.
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"
8 – #include "ExtendedSystem.hpp"
10   #include "simError.h"
11  
11 – extern "C"{
12 –  
13 –  void v_constrain_a_( double &dt, int &n_atoms, double* mass,
14 –                       double* Rx, double* Ry, double* Rz,
15 –                       double* Vx, double* Vy, double* Vz,
16 –                       double* Fx, double* Fy, double* Fz,
17 –                       int &n_constrained, double *constr_sqr,
18 –                       int* constr_i, int* constr_j,
19 –                       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,
24 <                       double* Fx, double* Fy, double* Fz,
25 <                       double &Kinetic,
26 <                       int &n_constrained, double *constr_sqr,
27 <                       int* constr_i, int* constr_j,
28 <                       double &box_x, double &box_y, double &box_z );
29 < }
30 <
31 <
32 <
33 <
34 < Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff,
35 <                        ExtendedSystem* the_es ){
36 <  entry_plug = the_entry_plug;
13 > Symplectic::Symplectic( SimInfo* theInfo, ForceFields* the_ff ){
14 >  
15 >  info = theInfo;
16    myFF = the_ff;
38 –  myES = the_es;
17    isFirst = 1;
18  
19 <  std::cerr<< "calling symplectic constructor\n";
19 >  molecules = info->molecules;
20 >  nMols = info->n_mol;
21  
43 –  molecules = entry_plug->molecules;
44 –  nMols = entry_plug->n_mol;
45 –
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  
51 –  // grab the masses
52 –
53 –  mass = new double[entry_plug->n_atoms];
54 –  for(int i = 0; i < entry_plug->n_atoms; i++){
55 –    mass[i] = entry_plug->atoms[i]->getMass();
56 –  }  
57 –
27    // check for constraints
28 +  
29 +  constrainedI = NULL;
30 +  constrainedJ = NULL;
31 +  constrainedDsqr = NULL;
32 +  nConstrained = 0;
33  
34 <  is_constrained = 0;
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 + void Symplectic::checkConstraints( void ){
48 +
49 +
50 +  isConstrained = 0;
51 +
52    Constraint *temp_con;
53    Constraint *dummy_plug;
54 <  temp_con = new Constraint[entry_plug->n_SRI];
55 <  n_constrained = 0;
54 >  temp_con = new Constraint[info->n_SRI];
55 >  nConstrained = 0;
56    int constrained = 0;
57    
58    SRI** theArray;
# Line 76 | Line 66 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
66        if(constrained){
67          
68          dummy_plug = theArray[j]->get_constraint();
69 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
70 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
71 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        n_constrained++;
73 >        nConstrained++;
74          constrained = 0;
75        }
76      }
# Line 93 | Line 83 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
83        if(constrained){
84          
85          dummy_plug = theArray[j]->get_constraint();
86 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
87 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
88 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        n_constrained++;
90 >        nConstrained++;
91          constrained = 0;
92        }
93      }
# Line 110 | Line 100 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
100        if(constrained){
101          
102          dummy_plug = theArray[j]->get_constraint();
103 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
104 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
105 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
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 <        n_constrained++;
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;
135 <      constrained_j[i] = temp_con[i].get_b() + 1;
136 <      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  
143 – Symplectic::~Symplectic() {
144 –  
145 –  if( n_constrained ){
146 –    delete[] constrained_i;
147 –    delete[] constrained_j;
148 –    delete[] constrained_dsqr;
149 –  }
150 –  
151 – }
136  
153 –
137   void Symplectic::integrate( void ){
138  
156 –  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
157 –
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 164 | 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 172 | 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 );
# Line 187 | Line 169 | void Symplectic::integrate( void ){
169    int vel_n    = (int)( thermalTime / dt );
170  
171    int calcPot, calcStress;
172 +  int isError;
173  
174 <  Thermo *tStats;
175 <  StatWriter*  e_out;
176 <  DumpWriter*  dump_out;
174 >  tStats   = new Thermo( info );
175 >  e_out    = new StatWriter( info );
176 >  dump_out = new DumpWriter( info );
177  
178 <  std::cerr << "about to call new thermo\n";
196 <
197 <  tStats   = new Thermo( entry_plug );
198 <  e_out    = new StatWriter( entry_plug );
199 <
200 <  std::cerr << "calling dumpWriter \n";
201 <  dump_out = new DumpWriter( entry_plug );
202 <  std::cerr << "called dumpWriter \n";
203 <
204 <  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,1);
185    
186 <  if( entry_plug->setTemp ){
186 >  if( info->setTemp ){
187      
188      tStats->velocitize();
189    }
# Line 219 | Line 193 | void Symplectic::integrate( void ){
193    
194    calcPot = 0;
195  
196 <  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
197 <    calcStress = 1;
198 <  } else {
199 <    calcStress = 0;
200 <  }
227 <
228 <  if( n_constrained ){
229 <
230 <    double *Rx = new double[nAtoms];
231 <    double *Ry = new double[nAtoms];
232 <    double *Rz = new double[nAtoms];
196 >  for( tl=0; tl<nLoops; tl++){
197 >
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];
205 <    
206 <    double *Fx = new double[nAtoms];
207 <    double *Fy = new double[nAtoms];
208 <    double *Fz = new double[nAtoms];
209 <    
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 <    for( tl=0; tl < n_loops; tl++ ){
217 >  
218 >  }
219  
220 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
246 <        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
247 <      
248 <      for( j=0; j<nAtoms; j++ ){
220 >  dump_out->writeFinal();
221  
222 <        Rx[j] = atoms[j]->getX();
223 <        Ry[j] = atoms[j]->getY();
224 <        Rz[j] = atoms[j]->getZ();
222 >  delete dump_out;
223 >  delete e_out;
224 > }
225  
254 –        Vx[j] = atoms[j]->get_vx();
255 –        Vy[j] = atoms[j]->get_vy();
256 –        Vz[j] = atoms[j]->get_vz();
226  
227 <        Fx[j] = atoms[j]->getFx();
228 <        Fy[j] = atoms[j]->getFy();
229 <        Fz[j] = atoms[j]->getFz();
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 <      }
236 <        
237 <      v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
238 <                      Fx, Fy, Fz,
239 <                      n_constrained, constrained_dsqr,
240 <                      constrained_i, constrained_j,
241 <                      entry_plug->box_x,
269 <                      entry_plug->box_y,
270 <                      entry_plug->box_z );
271 <      
272 <      for( j=0; j<nAtoms; j++ ){
235 >  for( i=0; i<nAtoms; i++ ){
236 >    atomIndex = i * 3;
237 >    aMatIndex = i * 9;
238 >    
239 >    // velocity half step
240 >    for( j=atomIndex; j<(atomIndex+3); j++ )
241 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
242  
243 <        atoms[j]->setX(Rx[j]);
244 <        atoms[j]->setY(Ry[j]);
245 <        atoms[j]->setZ(Rz[j]);
277 <        
278 <        atoms[j]->set_vx(Vx[j]);
279 <        atoms[j]->set_vy(Vy[j]);
280 <        atoms[j]->set_vz(Vz[j]);
281 <      }
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 < //       for( i=0; i<nAtoms; i++ ){
285 < // //   if( atoms[i]->isDirectional() ){
286 <                  
287 < // //     dAtom = (DirectionalAtom *)atoms[i];
250 >      dAtom = (DirectionalAtom *)atoms[i];
251            
252 < // //     // get and convert the torque to body frame
290 <          
291 < // //     Tb[0] = dAtom->getTx();
292 < // //     Tb[1] = dAtom->getTy();
293 < // //     Tb[2] = dAtom->getTz();
294 <          
295 < // //     dAtom->lab2Body( Tb );
296 <          
297 < // //     // get the angular momentum, and propagate a half step
298 <          
299 < // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300 < // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301 < // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
302 <          
303 < // //     // get the atom's rotation matrix
304 <          
305 < // //     A[0][0] = dAtom->getAxx();
306 < // //     A[0][1] = dAtom->getAxy();
307 < // //     A[0][2] = dAtom->getAxz();
308 <          
309 < // //     A[1][0] = dAtom->getAyx();
310 < // //     A[1][1] = dAtom->getAyy();
311 < // //     A[1][2] = dAtom->getAyz();
312 <          
313 < // //     A[2][0] = dAtom->getAzx();
314 < // //     A[2][1] = dAtom->getAzy();
315 < // //     A[2][2] = dAtom->getAzz();
316 <          
317 <          
318 < // //     // use the angular velocities to propagate the rotation matrix a
319 < // //     // full time step
320 <          
321 <          
322 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
323 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
324 <          
325 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
326 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
327 <          
328 < // //     angle = dt * ji[2] / dAtom->getIzz();
329 < // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
330 <          
331 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
332 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
333 <          
334 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
335 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
336 <          
337 <          
338 < // //     dAtom->setA( A );
339 < // //     dAtom->setJx( ji[0] );
340 < // //     dAtom->setJy( ji[1] );
341 < // //     dAtom->setJz( ji[2] );
342 < // //   }
343 < //       }
252 >      // get and convert the torque to body frame
253        
254 <      // calculate the forces
254 >      Tb[0] = dAtom->getTx();
255 >      Tb[1] = dAtom->getTy();
256 >      Tb[2] = dAtom->getTz();
257        
258 <      myFF->doForces(calcPot, calcStress);
258 >      dAtom->lab2Body( Tb );
259        
260 <      // move b
350 <
351 <      for( j=0; j<nAtoms; j++ ){
352 <
353 <        Rx[j] = atoms[j]->getX();
354 <        Ry[j] = atoms[j]->getY();
355 <        Rz[j] = atoms[j]->getZ();
356 <
357 <        Vx[j] = atoms[j]->get_vx();
358 <        Vy[j] = atoms[j]->get_vy();
359 <        Vz[j] = atoms[j]->get_vz();
360 <
361 <        Fx[j] = atoms[j]->getFx();
362 <        Fy[j] = atoms[j]->getFy();
363 <        Fz[j] = atoms[j]->getFz();
364 <      }
365 <        
366 <      v_constrain_b_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
367 <                      Fx, Fy, Fz,
368 <                      kE, n_constrained, constrained_dsqr,
369 <                      constrained_i, constrained_j,
370 <                      entry_plug->box_x,
371 <                      entry_plug->box_y,
372 <                      entry_plug->box_z );
260 >      // get the angular momentum, and propagate a half step
261        
262 <      for( j=0; j<nAtoms; j++ ){
263 <
264 <        atoms[j]->setX(Rx[j]);
377 <        atoms[j]->setY(Ry[j]);
378 <        atoms[j]->setZ(Rz[j]);
379 <
380 <        atoms[j]->set_vx(Vx[j]);
381 <        atoms[j]->set_vy(Vy[j]);
382 <        atoms[j]->set_vz(Vz[j]);
383 <      }
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 < //       for( i=0; i< nAtoms; i++ ){
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  
387 – //      if( atoms[i]->isDirectional() ){
297  
298 < //        dAtom = (DirectionalAtom *)atoms[i];
299 <          
300 < //        // get and convert the torque to body frame
301 <          
302 < //        Tb[0] = dAtom->getTx();
303 < //        Tb[1] = dAtom->getTy();
395 < //        Tb[2] = dAtom->getTz();
396 <          
397 < //        dAtom->lab2Body( Tb );
398 <          
399 < //        // get the angular momentum, and complete the angular momentum
400 < //        // half step
401 <          
402 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
405 <          
406 < //        dAtom->setJx( ji[0] );
407 < //        dAtom->setJy( ji[1] );
408 < //        dAtom->setJz( ji[2] );
409 < //      }
410 < //       }
411 <    
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 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
306 <        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
305 >  for( i=0; i<nAtoms; i++ ){
306 >    atomIndex = i * 3;
307  
308 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
309 <        myES->NoseHooverAndersonNPT( dt,
310 <                                     tStats->getKinetic(),
419 <                                     tStats->getPressure());
308 >    // velocity half step
309 >    for( j=atomIndex; j<(atomIndex+3); j++ )
310 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
311  
312 <      time = tl + 1;
312 >    if( atoms[i]->isDirectional() ){
313        
314 <      if( entry_plug->setTemp ){
315 <        if( !(time % vel_n) ) tStats->velocitize();
316 <      }
317 <      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
318 <      if( !((time+1) % status_n) ) {
319 <        calcPot = 1;
320 <        // bitwise masking in case we need it for NPT
321 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
322 <      }
323 <      if( !(time % status_n) ){
324 <        e_out->writeStat( time * dt );
325 <        calcPot = 0;
326 <        // bitwise masking in case we need it for NPT
327 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
328 <      }
314 >      dAtom = (DirectionalAtom *)atoms[i];
315 >      
316 >      // get and convert the torque to body frame
317 >      
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    }
440 –  else{
340  
341 <    for( tl=0; tl<n_loops; tl++ ){
443 <      
444 <      kE = 0.0;
445 <      rot_kE= 0.0;
446 <      trans_kE = 0.0;
341 > }
342  
448 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
449 –        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
450 –      
451 –      for( i=0; i<nAtoms; i++ ){
452 –        
453 –        // velocity half step
454 –        
455 –        vx = atoms[i]->get_vx() +
456 –          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
457 –        vy = atoms[i]->get_vy() +
458 –          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
459 –        vz = atoms[i]->get_vz() +
460 –          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
461 –        
462 –        // position whole step
463 –        
464 –        rx = atoms[i]->getX() + dt * vx;
465 –        ry = atoms[i]->getY() + dt * vy;
466 –        rz = atoms[i]->getZ() + dt * vz;
467 –        
468 –        atoms[i]->setX( rx );
469 –        atoms[i]->setY( ry );
470 –        atoms[i]->setZ( rz );
471 –        
472 –        atoms[i]->set_vx( vx );
473 –        atoms[i]->set_vy( vy );
474 –        atoms[i]->set_vz( vz );
475 –        
476 – //      if( atoms[i]->isDirectional() ){
343  
344 < //        dAtom = (DirectionalAtom *)atoms[i];
479 <          
480 < //        // get and convert the torque to body frame
481 <          
482 < //        Tb[0] = dAtom->getTx();
483 < //        Tb[1] = dAtom->getTy();
484 < //        Tb[2] = dAtom->getTz();
485 <          
486 < //        dAtom->lab2Body( Tb );
487 <          
488 < //        // get the angular momentum, and propagate a half step
489 <          
490 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
493 <          
494 < //        // get the atom's rotation matrix
495 <          
496 < //        A[0][0] = dAtom->getAxx();
497 < //        A[0][1] = dAtom->getAxy();
498 < //        A[0][2] = dAtom->getAxz();
499 <          
500 < //        A[1][0] = dAtom->getAyx();
501 < //        A[1][1] = dAtom->getAyy();
502 < //        A[1][2] = dAtom->getAyz();
503 <          
504 < //        A[2][0] = dAtom->getAzx();
505 < //        A[2][1] = dAtom->getAzy();
506 < //        A[2][2] = dAtom->getAzz();
507 <          
508 <          
509 < //        // use the angular velocities to propagate the rotation matrix a
510 < //        // full time step
511 <          
512 <          
513 < //        angle = dt2 * ji[0] / dAtom->getIxx();
514 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
515 <          
516 < //        angle = dt2 * ji[1] / dAtom->getIyy();
517 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
518 <          
519 < //        angle = dt * ji[2] / dAtom->getIzz();
520 < //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
521 <          
522 < //        angle = dt2 * ji[1] / dAtom->getIyy();
523 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
524 <          
525 < //        angle = dt2 * ji[0] / dAtom->getIxx();
526 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
527 <          
528 <          
529 < //        dAtom->setA( A );
530 < //        dAtom->setJx( ji[0] );
531 < //        dAtom->setJy( ji[1] );
532 < //        dAtom->setJz( ji[2] );
533 < //      }
534 <      }
535 <      
536 <      // calculate the forces
537 <      
538 <      myFF->doForces(calcPot,calcStress);
539 <      
540 <      for( i=0; i< nAtoms; i++ ){
541 <        
542 <        // complete the velocity half step
543 <        
544 <        vx = atoms[i]->get_vx() +
545 <          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
546 <        vy = atoms[i]->get_vy() +
547 <          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
548 <        vz = atoms[i]->get_vz() +
549 <          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
550 <        
551 <        atoms[i]->set_vx( vx );
552 <        atoms[i]->set_vy( vy );
553 <        atoms[i]->set_vz( vz );
554 <        
555 < //      vx2 = vx * vx;
556 < //      vy2 = vy * vy;
557 < //      vz2 = vz * vz;
558 <        
559 < //      if( atoms[i]->isDirectional() ){
344 > void Integrator::constrainA(){
345  
561 – //        dAtom = (DirectionalAtom *)atoms[i];
562 –          
563 – //        // get and convert the torque to body frame
564 –          
565 – //        Tb[0] = dAtom->getTx();
566 – //        Tb[1] = dAtom->getTy();
567 – //        Tb[2] = dAtom->getTz();
568 –          
569 – //        dAtom->lab2Body( Tb );
570 –          
571 – //        // get the angular momentum, and complete the angular momentum
572 – //        // half step
573 –          
574 – //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575 – //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576 – //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
577 –          
578 – //        jx2 = ji[0] * ji[0];
579 – //        jy2 = ji[1] * ji[1];
580 – //        jz2 = ji[2] * ji[2];
581 –          
582 – //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583 – //          + (jz2 / dAtom->getIzz());
584 –          
585 – //        dAtom->setJx( ji[0] );
586 – //        dAtom->setJy( ji[1] );
587 – //        dAtom->setJz( ji[2] );
588 – //      }
589 –      }
590 –    
591 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
592 –        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
593 –
594 –      if (!strcasecmp( entry_plug->ensemble, "NPT") )
595 –        myES->NoseHooverAndersonNPT( dt,
596 –                                     tStats->getKinetic(),
597 –                                     tStats->getPressure());
346    
599 –      time = tl + 1;
600 –      
601 –      if( entry_plug->setTemp ){
602 –        if( !(time % vel_n) ) tStats->velocitize();
603 –      }
604 –      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605 –      if( !((time+1) % status_n) ) {
606 –        calcPot = 1;
607 –        // bitwise masking in case we need it for NPT
608 –        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
609 –      }
610 –      if( !(time % status_n) ){
611 –        e_out->writeStat( time * dt );
612 –        calcPot = 0;
613 –        // bitwise masking in case we need it for NPT
614 –        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
615 –      }      
616 –    }
617 –  }
347  
619 –  dump_out->writeFinal();
348  
621 –  delete dump_out;
622 –  delete e_out;
349   }
350  
351 +
352 +
353 +
354 +
355 +
356 +
357 +
358 +
359   void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
360                           double A[3][3] ){
361  

Diff Legend

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