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 475 by gezelter, Tue Apr 8 12:44:18 2003 UTC vs.
Revision 548 by mmeineke, Wed Jun 4 21:06:45 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 > Integrator::Integrator( 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  
27 <  // grab the masses
27 >  nAtoms = info->n_atoms;
28  
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 –
29    // check for constraints
30 +  
31 +  constrainedA    = NULL;
32 +  constrainedB    = NULL;
33 +  constrainedDsqr = NULL;
34 +  moving          = NULL;
35 +  moved           = NULL;
36 +  prePos          = NULL;
37 +  
38 +  nConstrained = 0;
39  
40 <  is_constrained = 0;
40 >  checkConstraints();
41 > }
42  
43 + Integrator::~Integrator() {
44 +  
45 +  if( nConstrained ){
46 +    delete[] constrainedA;
47 +    delete[] constrainedB;
48 +    delete[] constrainedDsqr;
49 +    delete[] moving;
50 +    delete[] moved;
51 +    delete[] prePos;
52 + k
53 +  }
54 +  
55 + }
56 +
57 + void Integrator::checkConstraints( void ){
58 +
59 +
60 +  isConstrained = 0;
61 +
62    Constraint *temp_con;
63    Constraint *dummy_plug;
64 <  temp_con = new Constraint[entry_plug->n_SRI];
65 <  n_constrained = 0;
64 >  temp_con = new Constraint[info->n_SRI];
65 >  nConstrained = 0;
66    int constrained = 0;
67    
68    SRI** theArray;
# Line 76 | Line 76 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
76        if(constrained){
77          
78          dummy_plug = theArray[j]->get_constraint();
79 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
80 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
81 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
79 >        temp_con[nConstrained].set_a( dummy_plug->get_a() );
80 >        temp_con[nConstrained].set_b( dummy_plug->get_b() );
81 >        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
82          
83 <        n_constrained++;
83 >        nConstrained++;
84          constrained = 0;
85        }
86      }
# Line 93 | Line 93 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
93        if(constrained){
94          
95          dummy_plug = theArray[j]->get_constraint();
96 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
97 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
98 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
96 >        temp_con[nConstrained].set_a( dummy_plug->get_a() );
97 >        temp_con[nConstrained].set_b( dummy_plug->get_b() );
98 >        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
99          
100 <        n_constrained++;
100 >        nConstrained++;
101          constrained = 0;
102        }
103      }
# Line 110 | Line 110 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
110        if(constrained){
111          
112          dummy_plug = theArray[j]->get_constraint();
113 <        temp_con[n_constrained].set_a( dummy_plug->get_a() );
114 <        temp_con[n_constrained].set_b( dummy_plug->get_b() );
115 <        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
113 >        temp_con[nConstrained].set_a( dummy_plug->get_a() );
114 >        temp_con[nConstrained].set_b( dummy_plug->get_b() );
115 >        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
116          
117 <        n_constrained++;
117 >        nConstrained++;
118          constrained = 0;
119        }
120      }
121    }
122  
123 <  if(n_constrained > 0){
123 >  if(nConstrained > 0){
124      
125 <    is_constrained = 1;
126 <    constrained_i = new int[n_constrained];
127 <    constrained_j = new int[n_constrained];
128 <    constrained_dsqr = new double[n_constrained];
125 >    isConstrained = 1;
126 >
127 >    if(constrainedA != NULL )    delete[] constrainedA;
128 >    if(constrainedB != NULL )    delete[] constrainedB;
129 >    if(constrainedDsqr != NULL ) delete[] constrainedDsqr;
130 >
131 >    constrainedA =    new int[nConstrained];
132 >    constrainedB =    new int[nConstrained];
133 >    constrainedDsqr = new double[nConstrained];
134      
135 <    for( int i = 0; i < n_constrained; i++){
135 >    for( int i = 0; i < nConstrained; i++){
136        
137 <      /* add 1 to the index for the fortran arrays. */
138 <
139 <      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();
137 >      constrainedA[i] = temp_con[i].get_a();
138 >      constrainedB[i] = temp_con[i].get_b();
139 >      constrainedDsqr[i] = temp_con[i].get_dsqr();
140      }
141 +
142 +    
143 +    // save oldAtoms to check for lode balanceing later on.
144 +    
145 +    oldAtoms = nAtoms;
146 +    
147 +    moving = new int[nAtoms];
148 +    moved  = new int[nAtoms];
149 +
150 +    prePos = new double[nAtoms*3];
151    }
152    
153    delete[] temp_con;
154   }
155  
143 – Symplectic::~Symplectic() {
144 –  
145 –  if( n_constrained ){
146 –    delete[] constrained_i;
147 –    delete[] constrained_j;
148 –    delete[] constrained_dsqr;
149 –  }
150 –  
151 – }
156  
157 + void Integrator::integrate( void ){
158  
154 – void Symplectic::integrate( void ){
155 –
156 –  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
157 –
159    int i, j;                         // loop counters
159 –  int nAtoms = entry_plug->n_atoms; // the number of atoms
160    double kE = 0.0;                  // the kinetic energy  
161    double rot_kE;
162    double trans_kE;
# Line 164 | Line 164 | void Symplectic::integrate( void ){
164    double dt2;                       // half the dt
165  
166    double vx, vy, vz;    // the velocities
167 < //  double vx2, vy2, vz2; // the square of the velocities
167 >  double vx2, vy2, vz2; // the square of the velocities
168    double rx, ry, rz;    // the postitions
169    
170    double ji[3];   // the body frame angular momentum
# Line 172 | Line 172 | void Symplectic::integrate( void ){
172    double Tb[3];   // torque in the body frame
173    double angle;   // the angle through which to rotate the rotation matrix
174    double A[3][3]; // the rotation matrix
175 +  double press[9];
176  
177    int time;
178  
179 <  double dt          = entry_plug->dt;
180 <  double runTime     = entry_plug->run_time;
181 <  double sampleTime  = entry_plug->sampleTime;
182 <  double statusTime  = entry_plug->statusTime;
183 <  double thermalTime = entry_plug->thermalTime;
179 >  double dt          = info->dt;
180 >  double runTime     = info->run_time;
181 >  double sampleTime  = info->sampleTime;
182 >  double statusTime  = info->statusTime;
183 >  double thermalTime = info->thermalTime;
184  
185    int n_loops  = (int)( runTime / dt );
186    int sample_n = (int)( sampleTime / dt );
# Line 187 | Line 188 | void Symplectic::integrate( void ){
188    int vel_n    = (int)( thermalTime / dt );
189  
190    int calcPot, calcStress;
191 +  int isError;
192  
193 <  Thermo *tStats;
194 <  StatWriter*  e_out;
195 <  DumpWriter*  dump_out;
193 >  tStats   = new Thermo( info );
194 >  e_out    = new StatWriter( info );
195 >  dump_out = new DumpWriter( info );
196  
197 <  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;
197 >  Atom** atoms = info->atoms;
198    DirectionalAtom* dAtom;
199    dt2 = 0.5 * dt;
200  
201 <  // initialize the forces the before the first step
201 >  // initialize the forces before the first step
202  
203    myFF->doForces(1,1);
204    
205 <  if( entry_plug->setTemp ){
205 >  if( info->setTemp ){
206      
207      tStats->velocitize();
208    }
# Line 219 | Line 212 | void Symplectic::integrate( void ){
212    
213    calcPot = 0;
214  
215 <  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
223 <    calcStress = 1;
224 <  } else {
225 <    calcStress = 0;
226 <  }
215 >  for( tl=0; tl<nLoops; tl++){
216  
217 <  if( n_constrained ){
218 <
219 <    double *Rx = new double[nAtoms];
231 <    double *Ry = new double[nAtoms];
232 <    double *Rz = new double[nAtoms];
217 >    integrateStep( calcPot, calcStress );
218 >      
219 >    time = tl + 1;
220      
221 <    double *Vx = new double[nAtoms];
222 <    double *Vy = new double[nAtoms];
223 <    double *Vz = new double[nAtoms];
224 <    
225 <    double *Fx = new double[nAtoms];
226 <    double *Fy = new double[nAtoms];
227 <    double *Fz = new double[nAtoms];
228 <    
221 >    if( info->setTemp ){
222 >      if( !(time % vel_n) ) tStats->velocitize();
223 >    }
224 >    if( !(time % sample_n) ) dump_out->writeDump( time * dt );
225 >    if( !((time+1) % status_n) ) {
226 >      calcPot = 1;
227 >      calcStress = 1;
228 >    }
229 >    if( !(time % status_n) ){
230 >      e_out->writeStat( time * dt );
231 >      calcPot = 0;
232 >      if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
233 >      else calcStress = 0;
234 >    }    
235  
236 <    for( tl=0; tl < n_loops; tl++ ){
236 >  
237 >  }
238  
239 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
246 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
247 <      
248 <      for( j=0; j<nAtoms; j++ ){
239 >  dump_out->writeFinal();
240  
241 <        Rx[j] = atoms[j]->getX();
242 <        Ry[j] = atoms[j]->getY();
243 <        Rz[j] = atoms[j]->getZ();
241 >  delete dump_out;
242 >  delete e_out;
243 > }
244  
254 –        Vx[j] = atoms[j]->get_vx();
255 –        Vy[j] = atoms[j]->get_vy();
256 –        Vz[j] = atoms[j]->get_vz();
245  
246 <        Fx[j] = atoms[j]->getFx();
247 <        Fy[j] = atoms[j]->getFy();
248 <        Fz[j] = atoms[j]->getFz();
246 > void Integrator::moveA( void ){
247 >  
248 >  int i,j,k;
249 >  int atomIndex, aMatIndex;
250 >  DirectionalAtom* dAtom;
251 >  double Tb[3];
252 >  double ji[3];
253  
254 <      }
255 <        
256 <      v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
257 <                      Fx, Fy, Fz,
258 <                      n_constrained, constrained_dsqr,
259 <                      constrained_i, constrained_j,
260 <                      entry_plug->box_x,
269 <                      entry_plug->box_y,
270 <                      entry_plug->box_z );
271 <      
272 <      for( j=0; j<nAtoms; j++ ){
254 >  for( i=0; i<nAtoms; i++ ){
255 >    atomIndex = i * 3;
256 >    aMatIndex = i * 9;
257 >    
258 >    // velocity half step
259 >    for( j=atomIndex; j<(atomIndex+3); j++ )
260 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
261  
262 <        atoms[j]->setX(Rx[j]);
263 <        atoms[j]->setY(Ry[j]);
264 <        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 <      }
262 >    // position whole step    
263 >    for( j=atomIndex; j<(atomIndex+3); j++ )
264 >      pos[j] += dt * vel[j];
265  
266 +  
267 +    if( atoms[i]->isDirectional() ){
268  
269 < //       for( i=0; i<nAtoms; i++ ){
285 < // //   if( atoms[i]->isDirectional() ){
286 <                  
287 < // //     dAtom = (DirectionalAtom *)atoms[i];
269 >      dAtom = (DirectionalAtom *)atoms[i];
270            
271 < // //     // 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 < //       }
271 >      // get and convert the torque to body frame
272        
273 <      // calculate the forces
273 >      Tb[0] = dAtom->getTx();
274 >      Tb[1] = dAtom->getTy();
275 >      Tb[2] = dAtom->getTz();
276        
277 <      myFF->doForces(calcPot, calcStress);
277 >      dAtom->lab2Body( Tb );
278        
279 <      // 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 );
279 >      // get the angular momentum, and propagate a half step
280        
281 <      for( j=0; j<nAtoms; j++ ){
282 <
283 <        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 <      }
281 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
282 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
283 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
284        
285 < //       for( i=0; i< nAtoms; i++ ){
285 >      // use the angular velocities to propagate the rotation matrix a
286 >      // full time step
287 >      
288 >      // rotate about the x-axis      
289 >      angle = dt2 * ji[0] / dAtom->getIxx();
290 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
291 >      
292 >      // rotate about the y-axis
293 >      angle = dt2 * ji[1] / dAtom->getIyy();
294 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
295 >      
296 >      // rotate about the z-axis
297 >      angle = dt * ji[2] / dAtom->getIzz();
298 >      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
299 >      
300 >      // rotate about the y-axis
301 >      angle = dt2 * ji[1] / dAtom->getIyy();
302 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
303 >      
304 >       // rotate about the x-axis
305 >      angle = dt2 * ji[0] / dAtom->getIxx();
306 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
307 >      
308 >      dAtom->setJx( ji[0] );
309 >      dAtom->setJy( ji[1] );
310 >      dAtom->setJz( ji[2] );
311 >    }
312 >    
313 >  }
314 > }
315  
387 – //      if( atoms[i]->isDirectional() ){
316  
317 < //        dAtom = (DirectionalAtom *)atoms[i];
318 <          
319 < //        // get and convert the torque to body frame
320 <          
321 < //        Tb[0] = dAtom->getTx();
322 < //        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 <    
317 > void Integrator::moveB( void ){
318 >  int i,j,k;
319 >  int atomIndex;
320 >  DirectionalAtom* dAtom;
321 >  double Tb[3];
322 >  double ji[3];
323  
324 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
325 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
324 >  for( i=0; i<nAtoms; i++ ){
325 >    atomIndex = i * 3;
326  
327 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
328 <        myES->NoseHooverAndersonNPT( dt,
329 <                                     tStats->getKinetic(),
419 <                                     tStats->getPressure());
327 >    // velocity half step
328 >    for( j=atomIndex; j<(atomIndex+3); j++ )
329 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
330  
331 <      time = tl + 1;
331 >    if( atoms[i]->isDirectional() ){
332        
333 <      if( entry_plug->setTemp ){
334 <        if( !(time % vel_n) ) tStats->velocitize();
335 <      }
336 <      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
337 <      if( !((time+1) % status_n) ) {
338 <        calcPot = 1;
339 <        // bitwise masking in case we need it for NPT
340 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
341 <      }
342 <      if( !(time % status_n) ){
343 <        e_out->writeStat( time * dt );
344 <        calcPot = 0;
345 <        // bitwise masking in case we need it for NPT
346 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
347 <      }
333 >      dAtom = (DirectionalAtom *)atoms[i];
334 >      
335 >      // get and convert the torque to body frame
336 >      
337 >      Tb[0] = dAtom->getTx();
338 >      Tb[1] = dAtom->getTy();
339 >      Tb[2] = dAtom->getTz();
340 >      
341 >      dAtom->lab2Body( Tb );
342 >      
343 >      // get the angular momentum, and complete the angular momentum
344 >      // half step
345 >      
346 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
347 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
348 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
349 >      
350 >      jx2 = ji[0] * ji[0];
351 >      jy2 = ji[1] * ji[1];
352 >      jz2 = ji[2] * ji[2];
353 >      
354 >      dAtom->setJx( ji[0] );
355 >      dAtom->setJy( ji[1] );
356 >      dAtom->setJz( ji[2] );
357      }
358    }
440 –  else{
359  
360 <    for( tl=0; tl<n_loops; tl++ ){
443 <      
444 <      kE = 0.0;
445 <      rot_kE= 0.0;
446 <      trans_kE = 0.0;
360 > }
361  
362 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
363 <        myES->NoseHooverNVT( dt , 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() ){
362 > void Integrator::preMove( void ){
363 >  int i;
364  
365 < //        dAtom = (DirectionalAtom *)atoms[i];
366 <          
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 <      }
365 >  if( nConstrained ){
366 >    if( oldAtoms != nAtoms ){
367        
368 <      // calculate the forces
368 >      // save oldAtoms to check for lode balanceing later on.
369        
370 <      myFF->doForces(calcPot,calcStress);
370 >      oldAtoms = nAtoms;
371        
372 <      for( i=0; i< nAtoms; i++ ){
373 <        
374 <        // complete the velocity half step
375 <        
376 <        vx = atoms[i]->get_vx() +
377 <          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
378 <        vy = atoms[i]->get_vy() +
379 <          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
380 <        vz = atoms[i]->get_vz() +
381 <          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
382 <        
383 <        atoms[i]->set_vx( vx );
384 <        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() ){
372 >      delete[] moving;
373 >      delete[] moved;
374 >      delete[] oldPos;
375 >      
376 >      moving = new int[nAtoms];
377 >      moved  = new int[nAtoms];
378 >      
379 >      oldPos = new double[nAtoms*3];
380 >    }
381 >  
382 >    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
383 >  }
384 > }  
385  
386 < //        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 , tStats->getKinetic() );
386 > void Integrator::constrainA(){
387  
388 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
389 <        myES->NoseHooverAndersonNPT( dt,
390 <                                     tStats->getKinetic(),
391 <                                     tStats->getPressure());
388 >  int i,j,k;
389 >  int done;
390 >  double pxab, pyab, pzab;
391 >  double rxab, ryab, rzab;
392 >  int a, b;
393 >  double rma, rmb;
394 >  double dx, dy, dz;
395 >  double rabsq, pabsq, rpabsq;
396 >  double diffsq;
397 >  double gab;
398 >  int iteration;
399 >
400 >
401    
402 <      time = tl + 1;
403 <      
404 <      if( entry_plug->setTemp ){
405 <        if( !(time % vel_n) ) tStats->velocitize();
402 >  for( i=0; i<nAtoms; i++){
403 >    
404 >    moving[i] = 0;
405 >    moved[i]  = 1;
406 >  }
407 >  
408 >  
409 >  iteration = 0;
410 >  done = 0;
411 >  while( !done && (iteration < maxIteration )){
412 >
413 >    done = 1;
414 >    for(i=0; i<nConstrained; i++){
415 >
416 >      a = constrainedA[i];
417 >      b = constrainedB[i];
418 >    
419 >      if( moved[a] || moved[b] ){
420 >        
421 >        pxab = pos[3*a+0] - pos[3*b+0];
422 >        pyab = pos[3*a+1] - pos[3*b+1];
423 >        pzab = pos[3*a+2] - pos[3*b+2];
424 >
425 >        //periodic boundary condition
426 >        pxab = pxab - info->box_x * copysign(1, pxab)
427 >          * int(pxab / info->box_x + 0.5);
428 >        pyab = pyab - info->box_y * copysign(1, pyab)
429 >          * int(pyab / info->box_y + 0.5);
430 >        pzab = pzab - info->box_z * copysign(1, pzab)
431 >          * int(pzab / info->box_z + 0.5);
432 >      
433 >        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
434 >        rabsq = constraintedDsqr[i];
435 >        diffsq = pabsq - rabsq;
436 >
437 >        // the original rattle code from alan tidesley
438 >        if (fabs(diffsq) > tol*rabsq*2) {
439 >          rxab = oldPos[3*a+0] - oldPos[3*b+0];
440 >          ryab = oldPos[3*a+1] - oldPos[3*b+1];
441 >          rzab = oldPos[3*a+2] - oldPos[3*b+2];
442 >
443 >          rxab = rxab - info->box_x * copysign(1, rxab)
444 >            * int(rxab / info->box_x + 0.5);
445 >          ryab = ryab - info->box_y * copysign(1, ryab)
446 >            * int(ryab / info->box_y + 0.5);
447 >          rzab = rzab - info->box_z * copysign(1, rzab)
448 >            * int(rzab / info->box_z + 0.5);
449 >
450 >          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
451 >          rpabsq = rpab * rpab;
452 >
453 >
454 >          if (rpabsq < (rabsq * -diffsq)){
455 > #ifdef IS_MPI
456 >            a = atoms[a]->getGlobalIndex();
457 >            b = atoms[b]->getGlobalIndex();
458 > #endif //is_mpi
459 >            sprintf( painCave.errMsg,
460 >                     "Constraint failure in constrainA at atom %d and %d\n.",
461 >                     a, b );
462 >            painCave.isFatal = 1;
463 >            simError();
464 >          }
465 >
466 >          rma = 1.0 / atoms[a]->getMass();
467 >          rmb = 1.0 / atoms[b]->getMass();
468 >          
469 >          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
470 >          dx = rxab * gab;
471 >          dy = ryab * gab;
472 >          dz = rzab * gab;
473 >
474 >          pos[3*a+0] += rma * dx;
475 >          pos[3*a+1] += rma * dy;
476 >          pos[3*a+2] += rma * dz;
477 >
478 >          pos[3*b+0] -= rmb * dx;
479 >          pos[3*b+1] -= rmb * dy;
480 >          pos[3*b+2] -= rmb * dz;
481 >
482 >          dx = dx / dt;
483 >          dy = dy / dt;
484 >          dz = dz / dt;
485 >
486 >          vel[3*a+0] += rma * dx;
487 >          vel[3*a+1] += rma * dy;
488 >          vel[3*a+2] += rma * dz;
489 >
490 >          vel[3*b+0] -= rmb * dx;
491 >          vel[3*b+1] -= rmb * dy;
492 >          vel[3*b+2] -= rmb * dz;
493 >
494 >          moving[a] = 1;
495 >          moving[b] = 1;
496 >          done = 0;
497 >        }
498        }
499 <      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
500 <      if( !((time+1) % status_n) ) {
501 <        calcPot = 1;
502 <        // bitwise masking in case we need it for NPT
503 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
499 >    }
500 >    
501 >    for(i=0; i<nAtoms; i++){
502 >      
503 >      moved[i] = moving[i];
504 >      moving[i] = 0;
505 >    }
506 >
507 >    iteration++;
508 >  }
509 >
510 >  if( !done ){
511 >
512 >    sprintf( painCae.errMsg,
513 >             "Constraint failure in constrainA, too many iterations: %d\n",
514 >             iterations );
515 >    painCave.isFatal = 1;
516 >    simError();
517 >  }
518 >
519 > }
520 >
521 > void Integrator::constrainB( void ){
522 >  
523 >  int i,j,k;
524 >  int done;
525 >  double vxab, vyab, vzab;
526 >  double rxab, ryab, rzab;
527 >  int a, b;
528 >  double rma, rmb;
529 >  double dx, dy, dz;
530 >  double rabsq, pabsq, rvab;
531 >  double diffsq;
532 >  double gab;
533 >  int iteration;
534 >
535 >  for(i=0; i<nAtom; i++){
536 >    moving[i] = 0;
537 >    moved[i] = 1;
538 >  }
539 >
540 >  done = 0;
541 >  while( !done && (iteration < maxIteration ) ){
542 >
543 >    for(i=0; i<nConstrained; i++){
544 >      
545 >      a = constrainedA[i];
546 >      b = constrainedB[i];
547 >
548 >      if( moved[a] || moved[b] ){
549 >        
550 >        vxab = vel[3*a+0] - vel[3*b+0];
551 >        vyab = vel[3*a+1] - vel[3*b+1];
552 >        vzab = vel[3*a+2] - vel[3*b+2];
553 >
554 >        rxab = pos[3*a+0] - pos[3*b+0];q
555 >        ryab = pos[3*a+1] - pos[3*b+1];
556 >        rzab = pos[3*a+2] - pos[3*b+2];
557 >        
558 >        rxab = rxab - info->box_x * copysign(1, rxab)
559 >          * int(rxab / info->box_x + 0.5);
560 >        ryab = ryab - info->box_y * copysign(1, ryab)
561 >          * int(ryab / info->box_y + 0.5);
562 >        rzab = rzab - info->box_z * copysign(1, rzab)
563 >          * int(rzab / info->box_z + 0.5);
564 >
565 >        rma = 1.0 / atoms[a]->getMass();
566 >        rmb = 1.0 / atoms[b]->getMass();
567 >
568 >        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
569 >          
570 >        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
571 >
572 >        if (fabs(gab) > tol) {
573 >          
574 >          dx = rxab * gab;
575 >          dy = ryab * gab;
576 >          dz = rzab * gab;
577 >          
578 >          vel[3*a+0] += rma * dx;
579 >          vel[3*a+1] += rma * dy;
580 >          vel[3*a+2] += rma * dz;
581 >
582 >          vel[3*b+0] -= rmb * dx;
583 >          vel[3*b+1] -= rmb * dy;
584 >          vel[3*b+2] -= rmb * dz;
585 >          
586 >          moving[a] = 1;
587 >          moving[b] = 1;
588 >          done = 0;
589 >        }
590        }
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 –      }      
591      }
592 +
593 +    for(i=0; i<nAtoms; i++){
594 +      moved[i] = moving[i];
595 +      moving[i] = 0;
596 +    }
597 +    
598 +    iteration++;
599    }
600  
601 <  dump_out->writeFinal();
601 >  if( !done ){
602  
603 <  delete dump_out;
604 <  delete e_out;
603 >  
604 >    sprintf( painCae.errMsg,
605 >             "Constraint failure in constrainB, too many iterations: %d\n",
606 >             iterations );
607 >    painCave.isFatal = 1;
608 >    simError();
609 >  }
610 >
611   }
612  
613 < void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
613 >
614 >
615 >
616 >
617 >
618 >
619 > void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
620                           double A[3][3] ){
621  
622    int i,j,k;

Diff Legend

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