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

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 558 by mmeineke, Thu Jun 19 19:21:23 2003 UTC vs.
Revision 572 by mmeineke, Wed Jul 2 21:26:55 2003 UTC

# Line 1 | Line 1
1   #include <iostream>
2   #include <cstdlib>
3 + #include <cmath>
4  
5   #ifdef IS_MPI
6   #include "mpiSimulation.hpp"
# Line 10 | Line 11
11   #include "simError.h"
12  
13  
14 < Integrator::Integrator( SimInfo* theInfo, ForceFields* the_ff ){
14 > Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){
15    
16    info = theInfo;
17    myFF = the_ff;
# Line 33 | Line 34 | Integrator::Integrator( SimInfo* theInfo, ForceFields*
34    constrainedDsqr = NULL;
35    moving          = NULL;
36    moved           = NULL;
37 <  prePos          = NULL;
37 >  oldPos          = NULL;
38    
39    nConstrained = 0;
40  
# Line 48 | Line 49 | Integrator::~Integrator() {
49      delete[] constrainedDsqr;
50      delete[] moving;
51      delete[] moved;
52 <    delete[] prePos;
52 < k
52 >    delete[] oldPos;
53    }
54    
55   }
# Line 137 | Line 137 | void Integrator::checkConstraints( void ){
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      
# Line 147 | Line 148 | void Integrator::checkConstraints( void ){
148      moving = new int[nAtoms];
149      moved  = new int[nAtoms];
150  
151 <    prePos = new double[nAtoms*3];
151 >    oldPos = new double[nAtoms*3];
152    }
153    
154    delete[] temp_con;
# Line 157 | Line 158 | void Integrator::integrate( void ){
158   void Integrator::integrate( void ){
159  
160    int i, j;                         // loop counters
160 –  double kE = 0.0;                  // the kinetic energy  
161 –  double rot_kE;
162 –  double trans_kE;
163 –  int tl;                        // the time loop conter
164 –  double dt2;                       // half the dt
161  
166 –  double vx, vy, vz;    // 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
171 –  double jx2, jy2, jz2; // the square of the angular momentums
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 –  double dt          = info->dt;
162    double runTime     = info->run_time;
163    double sampleTime  = info->sampleTime;
164    double statusTime  = info->statusTime;
# Line 188 | Line 172 | void Integrator::integrate( void ){
172    int calcPot, calcStress;
173    int isError;
174  
175 +
176 +
177    tStats   = new Thermo( info );
178 <  e_out    = new StatWriter( info );
179 <  dump_out = new DumpWriter( info );
178 >  statOut  = new StatWriter( info );
179 >  dumpOut  = new DumpWriter( info );
180  
181 <  Atom** atoms = info->atoms;
181 >  atoms = info->atoms;
182    DirectionalAtom* dAtom;
183 +
184 +  dt = info->dt;
185    dt2 = 0.5 * dt;
186  
187    // initialize the forces before the first step
# Line 205 | Line 193 | void Integrator::integrate( void ){
193      tStats->velocitize();
194    }
195    
196 <  dump_out->writeDump( 0.0 );
197 <  e_out->writeStat( 0.0 );
196 >  dumpOut->writeDump( 0.0 );
197 >  statOut->writeStat( 0.0 );
198    
199    calcPot     = 0;
200    calcStress  = 0;
# Line 215 | Line 203 | void Integrator::integrate( void ){
203    currStatus  = statusTime;
204    currTime    = 0.0;;
205  
206 +
207 +  readyCheck();
208 +
209 + #ifdef IS_MPI
210 +  strcpy( checkPointMsg,
211 +          "The integrator is ready to go." );
212 +  MPIcheckPoint();
213 + #endif // is_mpi
214 +
215 +
216 +  pos  = Atom::getPosArray();
217 +  vel  = Atom::getVelArray();
218 +  frc  = Atom::getFrcArray();
219 +  trq  = Atom::getTrqArray();
220 +  Amat = Atom::getAmatArray();
221 +
222    while( currTime < runTime ){
223  
224      if( (currTime+dt) >= currStatus ){
225        calcPot = 1;
226        calcStress = 1;
227      }
228 <    
228 >
229      integrateStep( calcPot, calcStress );
230        
231      currTime += dt;
# Line 234 | Line 238 | void Integrator::integrate( void ){
238      }
239  
240      if( currTime >= currSample ){
241 <      dump_out->writeDump( currTime );
241 >      dumpOut->writeDump( currTime );
242        currSample += sampleTime;
243      }
244  
245      if( currTime >= currStatus ){
246 <      e_out->writeStat( time * dt );
246 >      statOut->writeStat( currTime );
247        calcPot = 0;
248        calcStress = 0;
249        currStatus += statusTime;
250      }
251 +
252 + #ifdef IS_MPI
253 +    strcpy( checkPointMsg,
254 +            "successfully took a time step." );
255 +    MPIcheckPoint();
256 + #endif // is_mpi
257 +
258    }
259  
260 <  dump_out->writeFinal();
260 >  dumpOut->writeFinal(currTime);
261  
262 <  delete dump_out;
263 <  delete e_out;
262 >  delete dumpOut;
263 >  delete statOut;
264   }
265  
266   void Integrator::integrateStep( int calcPot, int calcStress ){
267  
268 +
269 +      
270    // Position full step, and velocity half step
271  
272 <  //preMove();
272 >  preMove();
273    moveA();
274    if( nConstrained ) constrainA();
275  
# Line 279 | Line 292 | void Integrator::moveA( void ){
292    DirectionalAtom* dAtom;
293    double Tb[3];
294    double ji[3];
295 +  double angle;
296  
297 +
298 +
299    for( i=0; i<nAtoms; i++ ){
300      atomIndex = i * 3;
301      aMatIndex = i * 9;
302 <    
302 >
303      // velocity half step
304      for( j=atomIndex; j<(atomIndex+3); j++ )
305        vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
306  
307      // position whole step    
308 <    for( j=atomIndex; j<(atomIndex+3); j++ )
309 <      pos[j] += dt * vel[j];
294 <
295 <  
308 >    for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
309 >    
310      if( atoms[i]->isDirectional() ){
311  
312        dAtom = (DirectionalAtom *)atoms[i];
# Line 316 | Line 330 | void Integrator::moveA( void ){
330        
331        // rotate about the x-axis      
332        angle = dt2 * ji[0] / dAtom->getIxx();
333 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
333 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
334        
335        // rotate about the y-axis
336        angle = dt2 * ji[1] / dAtom->getIyy();
337 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
337 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
338        
339        // rotate about the z-axis
340        angle = dt * ji[2] / dAtom->getIzz();
341 <      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
341 >      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
342        
343        // rotate about the y-axis
344        angle = dt2 * ji[1] / dAtom->getIyy();
345 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
345 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
346        
347         // rotate about the x-axis
348        angle = dt2 * ji[0] / dAtom->getIxx();
349 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
349 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
350        
351        dAtom->setJx( ji[0] );
352        dAtom->setJy( ji[1] );
# Line 376 | Line 390 | void Integrator::moveB( void ){
390        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
391        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
392        
379 –      jx2 = ji[0] * ji[0];
380 –      jy2 = ji[1] * ji[1];
381 –      jz2 = ji[2] * ji[2];
382 –      
393        dAtom->setJx( ji[0] );
394        dAtom->setJy( ji[1] );
395        dAtom->setJz( ji[2] );
# Line 392 | Line 402 | void Integrator::preMove( void ){
402    int i;
403  
404    if( nConstrained ){
405 <    if( oldAtoms != nAtoms ){
396 <      
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 <  
405 >
406      for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
407    }
408   }  
# Line 416 | Line 411 | void Integrator::constrainA(){
411  
412    int i,j,k;
413    int done;
414 <  double pxab, pyab, pzab;
415 <  double rxab, ryab, rzab;
416 <  int a, b;
414 >  double pab[3];
415 >  double rab[3];
416 >  int a, b, ax, ay, az, bx, by, bz;
417    double rma, rmb;
418    double dx, dy, dz;
419 +  double rpab;
420    double rabsq, pabsq, rpabsq;
421    double diffsq;
422    double gab;
# Line 433 | Line 429 | void Integrator::constrainA(){
429      moving[i] = 0;
430      moved[i]  = 1;
431    }
432 <  
437 <  
432 >
433    iteration = 0;
434    done = 0;
435    while( !done && (iteration < maxIteration )){
# Line 444 | Line 439 | void Integrator::constrainA(){
439  
440        a = constrainedA[i];
441        b = constrainedB[i];
442 <    
442 >      
443 >      ax = (a*3) + 0;
444 >      ay = (a*3) + 1;
445 >      az = (a*3) + 2;
446 >
447 >      bx = (b*3) + 0;
448 >      by = (b*3) + 1;
449 >      bz = (b*3) + 2;
450 >
451        if( moved[a] || moved[b] ){
452          
453 <        pxab = pos[3*a+0] - pos[3*b+0];
454 <        pyab = pos[3*a+1] - pos[3*b+1];
455 <        pzab = pos[3*a+2] - pos[3*b+2];
453 >        pab[0] = pos[ax] - pos[bx];
454 >        pab[1] = pos[ay] - pos[by];
455 >        pab[2] = pos[az] - pos[bz];
456  
457 <        //periodic boundary condition
458 <        pxab = pxab - info->box_x * copysign(1, pxab)
459 <          * 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 <        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
463 <        rabsq = constraintedDsqr[i];
464 <        diffsq = pabsq - rabsq;
457 >        //periodic boundary condition
458 >
459 >        info->wrapVector( pab );
460  
461 +        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
462 +
463 +        rabsq = constrainedDsqr[i];
464 +        diffsq = rabsq - pabsq;
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);
467 >        if (fabs(diffsq) > (tol*rabsq*2)) {
468 >          rab[0] = oldPos[ax] - oldPos[bx];
469 >          rab[1] = oldPos[ay] - oldPos[by];
470 >          rab[2] = oldPos[az] - oldPos[bz];
471  
472 <          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
472 >          info->wrapVector( rab );
473 >
474 >          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
475 >
476            rpabsq = rpab * rpab;
477  
478  
479            if (rpabsq < (rabsq * -diffsq)){
480 +
481   #ifdef IS_MPI
482              a = atoms[a]->getGlobalIndex();
483              b = atoms[b]->getGlobalIndex();
484   #endif //is_mpi
485              sprintf( painCave.errMsg,
486 <                     "Constraint failure in constrainA at atom %d and %d\n.",
486 >                     "Constraint failure in constrainA at atom %d and %d.\n",
487                       a, b );
488              painCave.isFatal = 1;
489              simError();
# Line 494 | Line 491 | void Integrator::constrainA(){
491  
492            rma = 1.0 / atoms[a]->getMass();
493            rmb = 1.0 / atoms[b]->getMass();
494 <          
494 >
495            gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
499 –          dx = rxab * gab;
500 –          dy = ryab * gab;
501 –          dz = rzab * gab;
496  
497 <          pos[3*a+0] += rma * dx;
498 <          pos[3*a+1] += rma * dy;
499 <          pos[3*a+2] += rma * dz;
497 >          dx = rab[0] * gab;
498 >          dy = rab[1] * gab;
499 >          dz = rab[2] * gab;
500  
501 <          pos[3*b+0] -= rmb * dx;
502 <          pos[3*b+1] -= rmb * dy;
503 <          pos[3*b+2] -= rmb * dz;
501 >          pos[ax] += rma * dx;
502 >          pos[ay] += rma * dy;
503 >          pos[az] += rma * dz;
504  
505 +          pos[bx] -= rmb * dx;
506 +          pos[by] -= rmb * dy;
507 +          pos[bz] -= rmb * dz;
508 +
509            dx = dx / dt;
510            dy = dy / dt;
511            dz = dz / dt;
512  
513 <          vel[3*a+0] += rma * dx;
514 <          vel[3*a+1] += rma * dy;
515 <          vel[3*a+2] += rma * dz;
513 >          vel[ax] += rma * dx;
514 >          vel[ay] += rma * dy;
515 >          vel[az] += rma * dz;
516  
517 <          vel[3*b+0] -= rmb * dx;
518 <          vel[3*b+1] -= rmb * dy;
519 <          vel[3*b+2] -= rmb * dz;
517 >          vel[bx] -= rmb * dx;
518 >          vel[by] -= rmb * dy;
519 >          vel[bz] -= rmb * dz;
520  
521            moving[a] = 1;
522            moving[b] = 1;
# Line 538 | Line 536 | void Integrator::constrainA(){
536  
537    if( !done ){
538  
539 <    sprintf( painCae.errMsg,
539 >    sprintf( painCave.errMsg,
540               "Constraint failure in constrainA, too many iterations: %d\n",
541 <             iterations );
541 >             iteration );
542      painCave.isFatal = 1;
543      simError();
544    }
# Line 552 | Line 550 | void Integrator::constrainB( void ){
550    int i,j,k;
551    int done;
552    double vxab, vyab, vzab;
553 <  double rxab, ryab, rzab;
554 <  int a, b;
553 >  double rab[3];
554 >  int a, b, ax, ay, az, bx, by, bz;
555    double rma, rmb;
556    double dx, dy, dz;
557    double rabsq, pabsq, rvab;
# Line 561 | Line 559 | void Integrator::constrainB( void ){
559    double gab;
560    int iteration;
561  
562 <  for(i=0; i<nAtom; i++){
562 >  for(i=0; i<nAtoms; i++){
563      moving[i] = 0;
564      moved[i] = 1;
565    }
566  
567    done = 0;
568 +  iteration = 0;
569    while( !done && (iteration < maxIteration ) ){
570  
571 +    done = 1;
572 +
573      for(i=0; i<nConstrained; i++){
574        
575        a = constrainedA[i];
576        b = constrainedB[i];
577  
578 +      ax = (a*3) + 0;
579 +      ay = (a*3) + 1;
580 +      az = (a*3) + 2;
581 +
582 +      bx = (b*3) + 0;
583 +      by = (b*3) + 1;
584 +      bz = (b*3) + 2;
585 +
586        if( moved[a] || moved[b] ){
587          
588 <        vxab = vel[3*a+0] - vel[3*b+0];
589 <        vyab = vel[3*a+1] - vel[3*b+1];
590 <        vzab = vel[3*a+2] - vel[3*b+2];
588 >        vxab = vel[ax] - vel[bx];
589 >        vyab = vel[ay] - vel[by];
590 >        vzab = vel[az] - vel[bz];
591  
592 <        rxab = pos[3*a+0] - pos[3*b+0];q
593 <        ryab = pos[3*a+1] - pos[3*b+1];
594 <        rzab = pos[3*a+2] - pos[3*b+2];
592 >        rab[0] = pos[ax] - pos[bx];
593 >        rab[1] = pos[ay] - pos[by];
594 >        rab[2] = pos[az] - pos[bz];
595          
596 <        rxab = rxab - info->box_x * copysign(1, rxab)
597 <          * 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 <
596 >        info->wrapVector( rab );
597 >        
598          rma = 1.0 / atoms[a]->getMass();
599          rmb = 1.0 / atoms[b]->getMass();
600  
601 <        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
601 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
602            
603 <        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
603 >        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
604  
605          if (fabs(gab) > tol) {
606            
607 <          dx = rxab * gab;
608 <          dy = ryab * gab;
609 <          dz = rzab * gab;
607 >          dx = rab[0] * gab;
608 >          dy = rab[1] * gab;
609 >          dz = rab[2] * gab;
610            
611 <          vel[3*a+0] += rma * dx;
612 <          vel[3*a+1] += rma * dy;
613 <          vel[3*a+2] += rma * dz;
611 >          vel[ax] += rma * dx;
612 >          vel[ay] += rma * dy;
613 >          vel[az] += rma * dz;
614  
615 <          vel[3*b+0] -= rmb * dx;
616 <          vel[3*b+1] -= rmb * dy;
617 <          vel[3*b+2] -= rmb * dz;
615 >          vel[bx] -= rmb * dx;
616 >          vel[by] -= rmb * dy;
617 >          vel[bz] -= rmb * dz;
618            
619            moving[a] = 1;
620            moving[b] = 1;
# Line 630 | Line 634 | void Integrator::constrainB( void ){
634    if( !done ){
635  
636    
637 <    sprintf( painCae.errMsg,
637 >    sprintf( painCave.errMsg,
638               "Constraint failure in constrainB, too many iterations: %d\n",
639 <             iterations );
639 >             iteration );
640      painCave.isFatal = 1;
641      simError();
642    }
# Line 646 | Line 650 | void Integrator::rotate( int axes1, int axes2, double
650  
651  
652   void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
653 <                         double A[3][3] ){
653 >                         double A[9] ){
654  
655    int i,j,k;
656    double sinAngle;
# Line 662 | Line 666 | void Integrator::rotate( int axes1, int axes2, double
666  
667    for(i=0; i<3; i++){
668      for(j=0; j<3; j++){
669 <      tempA[j][i] = A[i][j];
669 >      tempA[j][i] = A[3*i + j];
670      }
671    }
672  
# Line 713 | Line 717 | void Integrator::rotate( int axes1, int axes2, double
717    //            A[][] = A[][] * transpose(rot[][])
718  
719  
720 <  // NOte for as yet unknown reason, we are setting the performing the
720 >  // NOte for as yet unknown reason, we are performing the
721    // calculation as:
722    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
723  
724    for(i=0; i<3; i++){
725      for(j=0; j<3; j++){
726 <      A[j][i] = 0.0;
726 >      A[3*j + i] = 0.0;
727        for(k=0; k<3; k++){
728 <        A[j][i] += tempA[i][k] * rot[j][k];
728 >        A[3*j + i] += tempA[i][k] * rot[j][k];
729        }
730      }
731    }

Diff Legend

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