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 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC vs.
Revision 556 by mmeineke, Wed Jun 18 22:20:12 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 <  molecules = entry_plug->molecules;
20 <  nMols = entry_plug->n_mol;
19 >  molecules = info->molecules;
20 >  nMols = info->n_mol;
21  
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  
51 –  mass = new double[entry_plug->n_atoms];
52 –  for(int i = 0; i < entry_plug->n_atoms; i++){
53 –    mass[i] = entry_plug->atoms[i]->getMass();
54 –  }  
55 –
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 74 | 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 91 | 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 108 | 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;
133 <      constrained_j[i] = temp_con[i].get_b() + 1;
134 <      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  
141 – Symplectic::~Symplectic() {
142 –  
143 –  if( n_constrained ){
144 –    delete[] constrained_i;
145 –    delete[] constrained_j;
146 –    delete[] constrained_dsqr;
147 –  }
148 –  
149 – }
156  
157 + void Integrator::integrate( void ){
158  
152 – void Symplectic::integrate( void ){
153 –
154 –  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
155 –
159    int i, j;                         // loop counters
157 –  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 172 | Line 174 | void Symplectic::integrate( void ){
174    double A[3][3]; // the rotation matrix
175    double press[9];
176  
177 <  int time;
177 >  double dt          = info->dt;
178 >  double runTime     = info->run_time;
179 >  double sampleTime  = info->sampleTime;
180 >  double statusTime  = info->statusTime;
181 >  double thermalTime = info->thermalTime;
182  
183 <  double dt          = entry_plug->dt;
184 <  double runTime     = entry_plug->run_time;
185 <  double sampleTime  = entry_plug->sampleTime;
186 <  double statusTime  = entry_plug->statusTime;
181 <  double thermalTime = entry_plug->thermalTime;
183 >  double currSample;
184 >  double currThermal;
185 >  double currStatus;
186 >  double currTime;
187  
183 –  int n_loops  = (int)( runTime / dt );
184 –  int sample_n = (int)( sampleTime / dt );
185 –  int status_n = (int)( statusTime / dt );
186 –  int vel_n    = (int)( thermalTime / dt );
187 –
188    int calcPot, calcStress;
189 +  int isError;
190  
191 <  Thermo *tStats;
192 <  StatWriter*  e_out;
193 <  DumpWriter*  dump_out;
191 >  tStats   = new Thermo( info );
192 >  e_out    = new StatWriter( info );
193 >  dump_out = new DumpWriter( info );
194  
195 <  tStats   = new Thermo( entry_plug );
195 <  e_out    = new StatWriter( entry_plug );
196 <  dump_out = new DumpWriter( entry_plug );
197 <
198 <  Atom** atoms = entry_plug->atoms;
195 >  Atom** atoms = info->atoms;
196    DirectionalAtom* dAtom;
197    dt2 = 0.5 * dt;
198  
199 <  // initialize the forces the before the first step
199 >  // initialize the forces before the first step
200  
201    myFF->doForces(1,1);
202    
203 <  if( entry_plug->setTemp ){
203 >  if( info->setTemp ){
204      
205      tStats->velocitize();
206    }
# Line 211 | Line 208 | void Symplectic::integrate( void ){
208    dump_out->writeDump( 0.0 );
209    e_out->writeStat( 0.0 );
210    
211 <  calcPot = 0;
211 >  calcPot     = 0;
212 >  calcStress  = 0;
213 >  currSample  = sampleTime;
214 >  currThermal = thermalTime;
215 >  currStatus  = statusTime;
216 >  currTime    = 0.0;;
217  
218 <  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
219 <    calcStress = 1;
220 <  } else {
221 <    calcStress = 0;
218 >  while( currTime < runTime ){
219 >
220 >    if( (currTime+dt) >= currStatus ){
221 >      calcPot = 1;
222 >      calcStress = 1;
223 >    }
224 >    
225 >    integrateStep( calcPot, calcStress );
226 >      
227 >    currTime += dt;
228 >
229 >    if( info->setTemp ){
230 >      if( currTime >= currThermal ){
231 >        tStats->velocitize();
232 >        currThermal += thermalTime;
233 >      }
234 >    }
235 >
236 >    if( currTime >= currSample ){
237 >      dump_out->writeDump( currTime );
238 >      currSample += sampleTime;
239 >    }
240 >
241 >    if( currTime >= currStatus ){
242 >      e_out->writeStat( time * dt );
243 >      calcPot = 0;
244 >      calcStress = 0;
245 >      currStatus += statusTime;
246 >    }
247    }
248  
249 <  if( n_constrained ){
249 >  dump_out->writeFinal();
250  
251 <    double *Rx = new double[nAtoms];
252 <    double *Ry = new double[nAtoms];
253 <    double *Rz = new double[nAtoms];
254 <    
255 <    double *Vx = new double[nAtoms];
256 <    double *Vy = new double[nAtoms];
257 <    double *Vz = new double[nAtoms];
258 <    
259 <    double *Fx = new double[nAtoms];
260 <    double *Fy = new double[nAtoms];
261 <    double *Fz = new double[nAtoms];
235 <    
251 >  delete dump_out;
252 >  delete e_out;
253 > }
254 >
255 > void Integrator::integrateStep( int calcPot, int calcStress ){
256 >
257 >  // Position full step, and velocity half step
258 >
259 >  preMove();
260 >  moveA();
261 >  if( nConstrained ) constrainA();
262  
263 <    for( tl=0; tl < n_loops; tl++ ){
263 >  // calc forces
264  
265 <      
240 <      for( j=0; j<nAtoms; j++ ){
265 >  myFF->doForces(calcPot,calcStress);
266  
267 <        Rx[j] = atoms[j]->getX();
268 <        Ry[j] = atoms[j]->getY();
269 <        Rz[j] = atoms[j]->getZ();
267 >  // finish the velocity  half step
268 >  
269 >  moveB();
270 >  if( nConstrained ) constrainB();
271 >  
272 > }
273  
246 –        Vx[j] = atoms[j]->get_vx();
247 –        Vy[j] = atoms[j]->get_vy();
248 –        Vz[j] = atoms[j]->get_vz();
274  
275 <        Fx[j] = atoms[j]->getFx();
276 <        Fy[j] = atoms[j]->getFy();
277 <        Fz[j] = atoms[j]->getFz();
275 > void Integrator::moveA( void ){
276 >  
277 >  int i,j,k;
278 >  int atomIndex, aMatIndex;
279 >  DirectionalAtom* dAtom;
280 >  double Tb[3];
281 >  double ji[3];
282  
283 <      }
284 <        
285 <      v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
286 <                      Fx, Fy, Fz,
287 <                      n_constrained, constrained_dsqr,
288 <                      constrained_i, constrained_j,
289 <                      entry_plug->box_x,
261 <                      entry_plug->box_y,
262 <                      entry_plug->box_z );
263 <      
264 <      for( j=0; j<nAtoms; j++ ){
283 >  for( i=0; i<nAtoms; i++ ){
284 >    atomIndex = i * 3;
285 >    aMatIndex = i * 9;
286 >    
287 >    // velocity half step
288 >    for( j=atomIndex; j<(atomIndex+3); j++ )
289 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
290  
291 <        atoms[j]->setX(Rx[j]);
292 <        atoms[j]->setY(Ry[j]);
293 <        atoms[j]->setZ(Rz[j]);
269 <        
270 <        atoms[j]->set_vx(Vx[j]);
271 <        atoms[j]->set_vy(Vy[j]);
272 <        atoms[j]->set_vz(Vz[j]);
273 <      }
291 >    // position whole step    
292 >    for( j=atomIndex; j<(atomIndex+3); j++ )
293 >      pos[j] += dt * vel[j];
294  
295 +  
296 +    if( atoms[i]->isDirectional() ){
297  
298 <      for( i=0; i<nAtoms; i++ ){
277 <        if( atoms[i]->isDirectional() ){
278 <                  
279 <          dAtom = (DirectionalAtom *)atoms[i];
298 >      dAtom = (DirectionalAtom *)atoms[i];
299            
300 <          // get and convert the torque to body frame
282 <          
283 <          Tb[0] = dAtom->getTx();
284 <          Tb[1] = dAtom->getTy();
285 <          Tb[2] = dAtom->getTz();
286 <          
287 <          dAtom->lab2Body( Tb );
288 <          
289 <          // get the angular momentum, and propagate a half step
290 <          
291 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
292 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
293 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
294 <          
295 <          // get the atom's rotation matrix
296 <          
297 <          A[0][0] = dAtom->getAxx();
298 <          A[0][1] = dAtom->getAxy();
299 <          A[0][2] = dAtom->getAxz();
300 <          
301 <          A[1][0] = dAtom->getAyx();
302 <          A[1][1] = dAtom->getAyy();
303 <          A[1][2] = dAtom->getAyz();
304 <          
305 <          A[2][0] = dAtom->getAzx();
306 <          A[2][1] = dAtom->getAzy();
307 <          A[2][2] = dAtom->getAzz();
308 <          
309 <          
310 <          // use the angular velocities to propagate the rotation matrix a
311 <          // full time step
312 <          
313 <          
314 <          angle = dt2 * ji[0] / dAtom->getIxx();
315 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
316 <          
317 <          angle = dt2 * ji[1] / dAtom->getIyy();
318 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
319 <          
320 <          angle = dt * ji[2] / dAtom->getIzz();
321 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
322 <          
323 <          angle = dt2 * ji[1] / dAtom->getIyy();
324 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
325 <          
326 <          angle = dt2 * ji[0] / dAtom->getIxx();
327 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
328 <          
329 <          
330 <          dAtom->setA( A );
331 <          dAtom->setJx( ji[0] );
332 <          dAtom->setJy( ji[1] );
333 <          dAtom->setJz( ji[2] );
334 <        }
335 <      }
300 >      // get and convert the torque to body frame
301        
302 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
303 <        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
304 <
340 <      // calculate the forces
302 >      Tb[0] = dAtom->getTx();
303 >      Tb[1] = dAtom->getTy();
304 >      Tb[2] = dAtom->getTz();
305        
306 <      myFF->doForces(calcPot, calcStress);
306 >      dAtom->lab2Body( Tb );
307        
308 <      // move b
308 >      // get the angular momentum, and propagate a half step
309 >      
310 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
311 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
312 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
313 >      
314 >      // use the angular velocities to propagate the rotation matrix a
315 >      // full time step
316 >      
317 >      // rotate about the x-axis      
318 >      angle = dt2 * ji[0] / dAtom->getIxx();
319 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
320 >      
321 >      // rotate about the y-axis
322 >      angle = dt2 * ji[1] / dAtom->getIyy();
323 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
324 >      
325 >      // rotate about the z-axis
326 >      angle = dt * ji[2] / dAtom->getIzz();
327 >      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
328 >      
329 >      // rotate about the y-axis
330 >      angle = dt2 * ji[1] / dAtom->getIyy();
331 >      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
332 >      
333 >       // rotate about the x-axis
334 >      angle = dt2 * ji[0] / dAtom->getIxx();
335 >      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
336 >      
337 >      dAtom->setJx( ji[0] );
338 >      dAtom->setJy( ji[1] );
339 >      dAtom->setJz( ji[2] );
340 >    }
341 >    
342 >  }
343 > }
344  
346 –      for( j=0; j<nAtoms; j++ ){
345  
346 <        Rx[j] = atoms[j]->getX();
347 <        Ry[j] = atoms[j]->getY();
348 <        Rz[j] = atoms[j]->getZ();
346 > void Integrator::moveB( void ){
347 >  int i,j,k;
348 >  int atomIndex;
349 >  DirectionalAtom* dAtom;
350 >  double Tb[3];
351 >  double ji[3];
352  
353 <        Vx[j] = atoms[j]->get_vx();
354 <        Vy[j] = atoms[j]->get_vy();
354 <        Vz[j] = atoms[j]->get_vz();
353 >  for( i=0; i<nAtoms; i++ ){
354 >    atomIndex = i * 3;
355  
356 <        Fx[j] = atoms[j]->getFx();
357 <        Fy[j] = atoms[j]->getFy();
358 <        Fz[j] = atoms[j]->getFz();
359 <      }
360 <        
361 <      v_constrain_b_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
362 <                      Fx, Fy, Fz,
363 <                      kE, n_constrained, constrained_dsqr,
364 <                      constrained_i, constrained_j,
365 <                      entry_plug->box_x,
366 <                      entry_plug->box_y,
367 <                      entry_plug->box_z );
368 <      
369 <      for( j=0; j<nAtoms; j++ ){
356 >    // velocity half step
357 >    for( j=atomIndex; j<(atomIndex+3); j++ )
358 >      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
359  
360 <        atoms[j]->setX(Rx[j]);
361 <        atoms[j]->setY(Ry[j]);
362 <        atoms[j]->setZ(Rz[j]);
360 >    if( atoms[i]->isDirectional() ){
361 >      
362 >      dAtom = (DirectionalAtom *)atoms[i];
363 >      
364 >      // get and convert the torque to body frame
365 >      
366 >      Tb[0] = dAtom->getTx();
367 >      Tb[1] = dAtom->getTy();
368 >      Tb[2] = dAtom->getTz();
369 >      
370 >      dAtom->lab2Body( Tb );
371 >      
372 >      // get the angular momentum, and complete the angular momentum
373 >      // half step
374 >      
375 >      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
376 >      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
377 >      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
378 >      
379 >      jx2 = ji[0] * ji[0];
380 >      jy2 = ji[1] * ji[1];
381 >      jz2 = ji[2] * ji[2];
382 >      
383 >      dAtom->setJx( ji[0] );
384 >      dAtom->setJy( ji[1] );
385 >      dAtom->setJz( ji[2] );
386 >    }
387 >  }
388  
389 <        atoms[j]->set_vx(Vx[j]);
390 <        atoms[j]->set_vy(Vy[j]);
391 <        atoms[j]->set_vz(Vz[j]);
392 <      }
389 > }
390 >
391 > void Integrator::preMove( void ){
392 >  int i;
393 >
394 >  if( nConstrained ){
395 >    if( oldAtoms != nAtoms ){
396        
397 <      for( i=0; i< nAtoms; i++ ){
397 >      // save oldAtoms to check for lode balanceing later on.
398 >      
399 >      oldAtoms = nAtoms;
400 >      
401 >      delete[] moving;
402 >      delete[] moved;
403 >      delete[] oldPos;
404 >      
405 >      moving = new int[nAtoms];
406 >      moved  = new int[nAtoms];
407 >      
408 >      oldPos = new double[nAtoms*3];
409 >    }
410 >  
411 >    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
412 >  }
413 > }  
414  
415 <        if( atoms[i]->isDirectional() ){
415 > void Integrator::constrainA(){
416  
417 <          dAtom = (DirectionalAtom *)atoms[i];
418 <          
419 <          // get and convert the torque to body frame
420 <          
421 <          Tb[0] = dAtom->getTx();
422 <          Tb[1] = dAtom->getTy();
423 <          Tb[2] = dAtom->getTz();
424 <          
425 <          dAtom->lab2Body( Tb );
426 <          
427 <          // get the angular momentum, and complete the angular momentum
395 <          // half step
396 <          
397 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
398 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
399 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
400 <          
401 <          dAtom->setJx( ji[0] );
402 <          dAtom->setJy( ji[1] );
403 <          dAtom->setJz( ji[2] );
404 <        }
405 <      }
406 <    
417 >  int i,j,k;
418 >  int done;
419 >  double pxab, pyab, pzab;
420 >  double rxab, ryab, rzab;
421 >  int a, b;
422 >  double rma, rmb;
423 >  double dx, dy, dz;
424 >  double rabsq, pabsq, rpabsq;
425 >  double diffsq;
426 >  double gab;
427 >  int iteration;
428  
408 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
409 –        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
429  
430 <      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
431 <        tStats->getPressureTensor(press);
432 <        myES->NoseHooverAndersonNPT( dt,
433 <                                     tStats->getKinetic(),
434 <                                     press);
435 <      }
430 >  
431 >  for( i=0; i<nAtoms; i++){
432 >    
433 >    moving[i] = 0;
434 >    moved[i]  = 1;
435 >  }
436 >  
437 >  
438 >  iteration = 0;
439 >  done = 0;
440 >  while( !done && (iteration < maxIteration )){
441  
442 <      time = tl + 1;
442 >    done = 1;
443 >    for(i=0; i<nConstrained; i++){
444 >
445 >      a = constrainedA[i];
446 >      b = constrainedB[i];
447 >    
448 >      if( moved[a] || moved[b] ){
449 >        
450 >        pxab = pos[3*a+0] - pos[3*b+0];
451 >        pyab = pos[3*a+1] - pos[3*b+1];
452 >        pzab = pos[3*a+2] - pos[3*b+2];
453 >
454 >        //periodic boundary condition
455 >        pxab = pxab - info->box_x * copysign(1, pxab)
456 >          * int(pxab / info->box_x + 0.5);
457 >        pyab = pyab - info->box_y * copysign(1, pyab)
458 >          * int(pyab / info->box_y + 0.5);
459 >        pzab = pzab - info->box_z * copysign(1, pzab)
460 >          * int(pzab / info->box_z + 0.5);
461        
462 <      if( entry_plug->setTemp ){
463 <        if( !(time % vel_n) ) tStats->velocitize();
462 >        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
463 >        rabsq = constraintedDsqr[i];
464 >        diffsq = pabsq - rabsq;
465 >
466 >        // the original rattle code from alan tidesley
467 >        if (fabs(diffsq) > tol*rabsq*2) {
468 >          rxab = oldPos[3*a+0] - oldPos[3*b+0];
469 >          ryab = oldPos[3*a+1] - oldPos[3*b+1];
470 >          rzab = oldPos[3*a+2] - oldPos[3*b+2];
471 >
472 >          rxab = rxab - info->box_x * copysign(1, rxab)
473 >            * int(rxab / info->box_x + 0.5);
474 >          ryab = ryab - info->box_y * copysign(1, ryab)
475 >            * int(ryab / info->box_y + 0.5);
476 >          rzab = rzab - info->box_z * copysign(1, rzab)
477 >            * int(rzab / info->box_z + 0.5);
478 >
479 >          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
480 >          rpabsq = rpab * rpab;
481 >
482 >
483 >          if (rpabsq < (rabsq * -diffsq)){
484 > #ifdef IS_MPI
485 >            a = atoms[a]->getGlobalIndex();
486 >            b = atoms[b]->getGlobalIndex();
487 > #endif //is_mpi
488 >            sprintf( painCave.errMsg,
489 >                     "Constraint failure in constrainA at atom %d and %d\n.",
490 >                     a, b );
491 >            painCave.isFatal = 1;
492 >            simError();
493 >          }
494 >
495 >          rma = 1.0 / atoms[a]->getMass();
496 >          rmb = 1.0 / atoms[b]->getMass();
497 >          
498 >          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
499 >          dx = rxab * gab;
500 >          dy = ryab * gab;
501 >          dz = rzab * gab;
502 >
503 >          pos[3*a+0] += rma * dx;
504 >          pos[3*a+1] += rma * dy;
505 >          pos[3*a+2] += rma * dz;
506 >
507 >          pos[3*b+0] -= rmb * dx;
508 >          pos[3*b+1] -= rmb * dy;
509 >          pos[3*b+2] -= rmb * dz;
510 >
511 >          dx = dx / dt;
512 >          dy = dy / dt;
513 >          dz = dz / dt;
514 >
515 >          vel[3*a+0] += rma * dx;
516 >          vel[3*a+1] += rma * dy;
517 >          vel[3*a+2] += rma * dz;
518 >
519 >          vel[3*b+0] -= rmb * dx;
520 >          vel[3*b+1] -= rmb * dy;
521 >          vel[3*b+2] -= rmb * dz;
522 >
523 >          moving[a] = 1;
524 >          moving[b] = 1;
525 >          done = 0;
526 >        }
527        }
423 –      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
424 –      if( !((time+1) % status_n) ) {
425 –        calcPot = 1;
426 –        calcStress = 1;
427 –      }
428 –      if( !(time % status_n) ){
429 –        e_out->writeStat( time * dt );
430 –        calcPot = 0;
431 –        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
432 –        else calcStress = 0;
433 –      }
528      }
529 +    
530 +    for(i=0; i<nAtoms; i++){
531 +      
532 +      moved[i] = moving[i];
533 +      moving[i] = 0;
534 +    }
535 +
536 +    iteration++;
537    }
436 –  else{
538  
539 <    for( tl=0; tl<n_loops; tl++ ){
539 >  if( !done ){
540 >
541 >    sprintf( painCae.errMsg,
542 >             "Constraint failure in constrainA, too many iterations: %d\n",
543 >             iterations );
544 >    painCave.isFatal = 1;
545 >    simError();
546 >  }
547 >
548 > }
549 >
550 > void Integrator::constrainB( void ){
551 >  
552 >  int i,j,k;
553 >  int done;
554 >  double vxab, vyab, vzab;
555 >  double rxab, ryab, rzab;
556 >  int a, b;
557 >  double rma, rmb;
558 >  double dx, dy, dz;
559 >  double rabsq, pabsq, rvab;
560 >  double diffsq;
561 >  double gab;
562 >  int iteration;
563 >
564 >  for(i=0; i<nAtom; i++){
565 >    moving[i] = 0;
566 >    moved[i] = 1;
567 >  }
568 >
569 >  done = 0;
570 >  while( !done && (iteration < maxIteration ) ){
571 >
572 >    for(i=0; i<nConstrained; i++){
573        
574 <      kE = 0.0;
575 <      rot_kE= 0.0;
442 <      trans_kE = 0.0;
574 >      a = constrainedA[i];
575 >      b = constrainedB[i];
576  
577 <      for( i=0; i<nAtoms; i++ ){
577 >      if( moved[a] || moved[b] ){
578          
579 <        // velocity half step
579 >        vxab = vel[3*a+0] - vel[3*b+0];
580 >        vyab = vel[3*a+1] - vel[3*b+1];
581 >        vzab = vel[3*a+2] - vel[3*b+2];
582 >
583 >        rxab = pos[3*a+0] - pos[3*b+0];q
584 >        ryab = pos[3*a+1] - pos[3*b+1];
585 >        rzab = pos[3*a+2] - pos[3*b+2];
586          
587 <        vx = atoms[i]->get_vx() +
588 <          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
589 <        vy = atoms[i]->get_vy() +
590 <          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
591 <        vz = atoms[i]->get_vz() +
592 <          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
454 <        
455 <        // position whole step
456 <        
457 <        rx = atoms[i]->getX() + dt * vx;
458 <        ry = atoms[i]->getY() + dt * vy;
459 <        rz = atoms[i]->getZ() + dt * vz;
460 <        
461 <        atoms[i]->setX( rx );
462 <        atoms[i]->setY( ry );
463 <        atoms[i]->setZ( rz );
464 <        
465 <        atoms[i]->set_vx( vx );
466 <        atoms[i]->set_vy( vy );
467 <        atoms[i]->set_vz( vz );
468 <        
469 <        if( atoms[i]->isDirectional() ){
587 >        rxab = rxab - info->box_x * copysign(1, rxab)
588 >          * int(rxab / info->box_x + 0.5);
589 >        ryab = ryab - info->box_y * copysign(1, ryab)
590 >          * int(ryab / info->box_y + 0.5);
591 >        rzab = rzab - info->box_z * copysign(1, rzab)
592 >          * int(rzab / info->box_z + 0.5);
593  
594 <          dAtom = (DirectionalAtom *)atoms[i];
594 >        rma = 1.0 / atoms[a]->getMass();
595 >        rmb = 1.0 / atoms[b]->getMass();
596 >
597 >        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
598            
599 <          // get and convert the torque to body frame
599 >        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
600 >
601 >        if (fabs(gab) > tol) {
602            
603 <          Tb[0] = dAtom->getTx();
604 <          Tb[1] = dAtom->getTy();
605 <          Tb[2] = dAtom->getTz();
603 >          dx = rxab * gab;
604 >          dy = ryab * gab;
605 >          dz = rzab * gab;
606            
607 <          dAtom->lab2Body( Tb );
608 <          
609 <          // get the angular momentum, and propagate a half step
482 <          
483 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
484 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
485 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
486 <          
487 <          // get the atom's rotation matrix
488 <          
489 <          A[0][0] = dAtom->getAxx();
490 <          A[0][1] = dAtom->getAxy();
491 <          A[0][2] = dAtom->getAxz();
492 <          
493 <          A[1][0] = dAtom->getAyx();
494 <          A[1][1] = dAtom->getAyy();
495 <          A[1][2] = dAtom->getAyz();
496 <          
497 <          A[2][0] = dAtom->getAzx();
498 <          A[2][1] = dAtom->getAzy();
499 <          A[2][2] = dAtom->getAzz();
500 <          
501 <          
502 <          // use the angular velocities to propagate the rotation matrix a
503 <          // full time step
504 <          
505 <          
506 <          angle = dt2 * ji[0] / dAtom->getIxx();
507 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
508 <          
509 <          angle = dt2 * ji[1] / dAtom->getIyy();
510 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
511 <          
512 <          angle = dt * ji[2] / dAtom->getIzz();
513 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
514 <          
515 <          angle = dt2 * ji[1] / dAtom->getIyy();
516 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
517 <          
518 <          angle = dt2 * ji[0] / dAtom->getIxx();
519 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
520 <          
521 <          
522 <          dAtom->setA( A );
523 <          dAtom->setJx( ji[0] );
524 <          dAtom->setJy( ji[1] );
525 <          dAtom->setJz( ji[2] );
526 <        }
527 <      }
607 >          vel[3*a+0] += rma * dx;
608 >          vel[3*a+1] += rma * dy;
609 >          vel[3*a+2] += rma * dz;
610  
611 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
612 <        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
613 <      
532 <      
533 <      // calculate the forces
534 <      
535 <      myFF->doForces(calcPot,calcStress);
536 <      
537 <      for( i=0; i< nAtoms; i++ ){
538 <        
539 <        // complete the velocity half step
540 <        
541 <        vx = atoms[i]->get_vx() +
542 <          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
543 <        vy = atoms[i]->get_vy() +
544 <          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
545 <        vz = atoms[i]->get_vz() +
546 <          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
547 <        
548 <        atoms[i]->set_vx( vx );
549 <        atoms[i]->set_vy( vy );
550 <        atoms[i]->set_vz( vz );
551 <        
552 <        vx2 = vx * vx;
553 <        vy2 = vy * vy;
554 <        vz2 = vz * vz;
555 <        
556 <        if( atoms[i]->isDirectional() ){
557 <
558 <          dAtom = (DirectionalAtom *)atoms[i];
611 >          vel[3*b+0] -= rmb * dx;
612 >          vel[3*b+1] -= rmb * dy;
613 >          vel[3*b+2] -= rmb * dz;
614            
615 <          // get and convert the torque to body frame
616 <          
617 <          Tb[0] = dAtom->getTx();
563 <          Tb[1] = dAtom->getTy();
564 <          Tb[2] = dAtom->getTz();
565 <          
566 <          dAtom->lab2Body( Tb );
567 <          
568 <          // get the angular momentum, and complete the angular momentum
569 <          // half step
570 <          
571 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
572 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
573 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
574 <          
575 <          jx2 = ji[0] * ji[0];
576 <          jy2 = ji[1] * ji[1];
577 <          jz2 = ji[2] * ji[2];
578 <          
579 <          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
580 <            + (jz2 / dAtom->getIzz());
581 <          
582 <          dAtom->setJx( ji[0] );
583 <          dAtom->setJy( ji[1] );
584 <          dAtom->setJz( ji[2] );
615 >          moving[a] = 1;
616 >          moving[b] = 1;
617 >          done = 0;
618          }
586 –
619        }
620 <    
589 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
590 <        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
620 >    }
621  
622 <      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
623 <        tStats->getPressureTensor(press);
624 <        myES->NoseHooverAndersonNPT( dt,
595 <                                     tStats->getKinetic(),
596 <                                     press);
597 <      }
598 <  
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 <        calcStress = 1;
608 <      }
609 <      if( !(time % status_n) ){
610 <        e_out->writeStat( time * dt );
611 <        calcPot = 0;
612 <        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
613 <        else calcStress = 0;
614 <      }      
622 >    for(i=0; i<nAtoms; i++){
623 >      moved[i] = moving[i];
624 >      moving[i] = 0;
625      }
626 +    
627 +    iteration++;
628    }
629  
630 <  dump_out->writeFinal();
630 >  if( !done ){
631  
632 <  delete dump_out;
633 <  delete e_out;
632 >  
633 >    sprintf( painCae.errMsg,
634 >             "Constraint failure in constrainB, too many iterations: %d\n",
635 >             iterations );
636 >    painCave.isFatal = 1;
637 >    simError();
638 >  }
639 >
640   }
641  
642 < void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
642 >
643 >
644 >
645 >
646 >
647 >
648 > void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
649                           double A[3][3] ){
650  
651    int i,j,k;

Diff Legend

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