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 563 by mmeineke, Mon Jun 23 21:24:03 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 | Integrator::Integrator( SimInfo* theInfo, ForceFields*
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 +      cerr << "constraint " << constrainedA[i] << " <-> " << constrainedB[i]
142 +           << " => " << constrainedDsqr[i] << "\n";
143      }
144  
145      
# Line 147 | Line 150 | void Integrator::checkConstraints( void ){
150      moving = new int[nAtoms];
151      moved  = new int[nAtoms];
152  
153 <    prePos = new double[nAtoms*3];
153 >    oldPos = new double[nAtoms*3];
154    }
155    
156    delete[] temp_con;
# Line 157 | Line 160 | void Integrator::integrate( void ){
160   void Integrator::integrate( void ){
161  
162    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
163  
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;
164    double runTime     = info->run_time;
165    double sampleTime  = info->sampleTime;
166    double statusTime  = info->statusTime;
# Line 188 | Line 174 | void Integrator::integrate( void ){
174    int calcPot, calcStress;
175    int isError;
176  
177 +
178 +
179    tStats   = new Thermo( info );
180 <  e_out    = new StatWriter( info );
181 <  dump_out = new DumpWriter( info );
180 >  statOut  = new StatWriter( info );
181 >  dumpOut  = new DumpWriter( info );
182  
183 <  Atom** atoms = info->atoms;
183 >  atoms = info->atoms;
184    DirectionalAtom* dAtom;
185 +
186 +  dt = info->dt;
187    dt2 = 0.5 * dt;
188  
189    // initialize the forces before the first step
# Line 205 | Line 195 | void Integrator::integrate( void ){
195      tStats->velocitize();
196    }
197    
198 <  dump_out->writeDump( 0.0 );
199 <  e_out->writeStat( 0.0 );
198 >  dumpOut->writeDump( 0.0 );
199 >  statOut->writeStat( 0.0 );
200    
201    calcPot     = 0;
202    calcStress  = 0;
# Line 215 | Line 205 | void Integrator::integrate( void ){
205    currStatus  = statusTime;
206    currTime    = 0.0;;
207  
208 +
209 +  readyCheck();
210 +
211 + #ifdef IS_MPI
212 +  strcpy( checkPointMsg,
213 +          "The integrator is ready to go." );
214 +  MPIcheckPoint();
215 + #endif // is_mpi
216 +
217 +
218 +  pos  = Atom::getPosArray();
219 +  vel  = Atom::getVelArray();
220 +  frc  = Atom::getFrcArray();
221 +  trq  = Atom::getTrqArray();
222 +  Amat = Atom::getAmatArray();
223 +
224    while( currTime < runTime ){
225  
226      if( (currTime+dt) >= currStatus ){
227        calcPot = 1;
228        calcStress = 1;
229      }
230 <    
230 >
231      integrateStep( calcPot, calcStress );
232        
233      currTime += dt;
# Line 234 | Line 240 | void Integrator::integrate( void ){
240      }
241  
242      if( currTime >= currSample ){
243 <      dump_out->writeDump( currTime );
243 >      dumpOut->writeDump( currTime );
244        currSample += sampleTime;
245      }
246  
247      if( currTime >= currStatus ){
248 <      e_out->writeStat( time * dt );
248 >      statOut->writeStat( currTime );
249        calcPot = 0;
250        calcStress = 0;
251        currStatus += statusTime;
252      }
247 –  }
253  
254 <  dump_out->writeFinal();
254 > #ifdef IS_MPI
255 >    strcpy( checkPointMsg,
256 >            "successfully took a time step." );
257 >    MPIcheckPoint();
258 > #endif // is_mpi
259  
260 <  delete dump_out;
261 <  delete e_out;
260 >  }
261 >
262 >  dumpOut->writeFinal();
263 >
264 >  delete dumpOut;
265 >  delete statOut;
266   }
267  
268   void Integrator::integrateStep( int calcPot, int calcStress ){
269  
270 +
271 +      
272    // Position full step, and velocity half step
273  
274 <  //preMove();
274 >  preMove();
275    moveA();
276    if( nConstrained ) constrainA();
277  
# Line 279 | Line 294 | void Integrator::moveA( void ){
294    DirectionalAtom* dAtom;
295    double Tb[3];
296    double ji[3];
297 +  double angle;
298  
299    for( i=0; i<nAtoms; i++ ){
300      atomIndex = i * 3;
# Line 316 | Line 332 | void Integrator::moveA( void ){
332        
333        // rotate about the x-axis      
334        angle = dt2 * ji[0] / dAtom->getIxx();
335 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
335 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
336        
337        // rotate about the y-axis
338        angle = dt2 * ji[1] / dAtom->getIyy();
339 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
339 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
340        
341        // rotate about the z-axis
342        angle = dt * ji[2] / dAtom->getIzz();
343 <      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
343 >      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
344        
345        // rotate about the y-axis
346        angle = dt2 * ji[1] / dAtom->getIyy();
347 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
347 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
348        
349         // rotate about the x-axis
350        angle = dt2 * ji[0] / dAtom->getIxx();
351 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
351 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
352        
353        dAtom->setJx( ji[0] );
354        dAtom->setJy( ji[1] );
# Line 376 | Line 392 | void Integrator::moveB( void ){
392        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
393        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
394        
379 –      jx2 = ji[0] * ji[0];
380 –      jy2 = ji[1] * ji[1];
381 –      jz2 = ji[2] * ji[2];
382 –      
395        dAtom->setJx( ji[0] );
396        dAtom->setJy( ji[1] );
397        dAtom->setJz( ji[2] );
# Line 392 | Line 404 | void Integrator::preMove( void ){
404    int i;
405  
406    if( nConstrained ){
407 <    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 <  
407 >
408      for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
409    }
410   }  
# Line 418 | Line 415 | void Integrator::constrainA(){
415    int done;
416    double pxab, pyab, pzab;
417    double rxab, ryab, rzab;
418 <  int a, b;
418 >  int a, b, ax, ay, az, bx, by, bz;
419    double rma, rmb;
420    double dx, dy, dz;
421 +  double rpab;
422    double rabsq, pabsq, rpabsq;
423    double diffsq;
424    double gab;
# Line 444 | Line 442 | void Integrator::constrainA(){
442  
443        a = constrainedA[i];
444        b = constrainedB[i];
445 <    
445 >      
446 >      ax = (a*3) + 0;
447 >      ay = (a*3) + 1;
448 >      az = (a*3) + 2;
449 >
450 >      bx = (b*3) + 0;
451 >      by = (b*3) + 1;
452 >      bz = (b*3) + 2;
453 >
454 >
455        if( moved[a] || moved[b] ){
456          
457 <        pxab = pos[3*a+0] - pos[3*b+0];
458 <        pyab = pos[3*a+1] - pos[3*b+1];
459 <        pzab = pos[3*a+2] - pos[3*b+2];
457 >        pxab = pos[ax] - pos[bx];
458 >        pyab = pos[ay] - pos[by];
459 >        pzab = pos[az] - pos[bz];
460  
461          //periodic boundary condition
462          pxab = pxab - info->box_x * copysign(1, pxab)
463 <          * int(pxab / info->box_x + 0.5);
464 <        pyab = pyab - info->box_y * copysign(1, pyab)
465 <          * int(pyab / info->box_y + 0.5);
463 >          * (int)( fabs(pxab / info->box_x) + 0.5);
464 >        pyab = pyab - info->box_y * copysign(1, pyab)
465 >          * (int)( fabs(pyab / info->box_y) + 0.5);
466          pzab = pzab - info->box_z * copysign(1, pzab)
467 <          * int(pzab / info->box_z + 0.5);
467 >          * (int)( fabs(pzab / info->box_z) + 0.5);
468        
469          pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
470 <        rabsq = constraintedDsqr[i];
470 >        rabsq = constrainedDsqr[i];
471          diffsq = pabsq - rabsq;
472  
473          // the original rattle code from alan tidesley
474 <        if (fabs(diffsq) > tol*rabsq*2) {
475 <          rxab = oldPos[3*a+0] - oldPos[3*b+0];
476 <          ryab = oldPos[3*a+1] - oldPos[3*b+1];
477 <          rzab = oldPos[3*a+2] - oldPos[3*b+2];
474 >        if (fabs(diffsq) > (tol*rabsq*2)) {
475 >          rxab = oldPos[ax] - oldPos[bx];
476 >          ryab = oldPos[ay] - oldPos[by];
477 >          rzab = oldPos[az] - oldPos[bz];
478  
479            rxab = rxab - info->box_x * copysign(1, rxab)
480 <            * int(rxab / info->box_x + 0.5);
480 >            * (int)( fabs(rxab / info->box_x) + 0.5);
481            ryab = ryab - info->box_y * copysign(1, ryab)
482 <            * int(ryab / info->box_y + 0.5);
482 >            * (int)( fabs(ryab / info->box_y) + 0.5);
483            rzab = rzab - info->box_z * copysign(1, rzab)
484 <            * int(rzab / info->box_z + 0.5);
484 >            * (int)( fabs(rzab / info->box_z) + 0.5);
485  
486            rpab = rxab * pxab + ryab * pyab + rzab * pzab;
487            rpabsq = rpab * rpab;
488  
489  
490            if (rpabsq < (rabsq * -diffsq)){
491 +
492 +            cerr << "rpabsq = " << rpabsq << ", rabsq = " << rabsq
493 +                 << ", -diffsq = " << -diffsq << "\n";
494 +
495   #ifdef IS_MPI
496              a = atoms[a]->getGlobalIndex();
497              b = atoms[b]->getGlobalIndex();
498   #endif //is_mpi
499              sprintf( painCave.errMsg,
500 <                     "Constraint failure in constrainA at atom %d and %d\n.",
500 >                     "Constraint failure in constrainA at atom %d and %d.\n",
501                       a, b );
502              painCave.isFatal = 1;
503              simError();
# Line 500 | Line 511 | void Integrator::constrainA(){
511            dy = ryab * gab;
512            dz = rzab * gab;
513  
514 <          pos[3*a+0] += rma * dx;
515 <          pos[3*a+1] += rma * dy;
516 <          pos[3*a+2] += rma * dz;
514 >          pos[ax] += rma * dx;
515 >          pos[ay] += rma * dy;
516 >          pos[az] += rma * dz;
517  
518 <          pos[3*b+0] -= rmb * dx;
519 <          pos[3*b+1] -= rmb * dy;
520 <          pos[3*b+2] -= rmb * dz;
518 >          pos[bx] -= rmb * dx;
519 >          pos[by] -= rmb * dy;
520 >          pos[bz] -= rmb * dz;
521  
522            dx = dx / dt;
523            dy = dy / dt;
524            dz = dz / dt;
525  
526 <          vel[3*a+0] += rma * dx;
527 <          vel[3*a+1] += rma * dy;
528 <          vel[3*a+2] += rma * dz;
526 >          vel[ax] += rma * dx;
527 >          vel[ay] += rma * dy;
528 >          vel[az] += rma * dz;
529  
530 <          vel[3*b+0] -= rmb * dx;
531 <          vel[3*b+1] -= rmb * dy;
532 <          vel[3*b+2] -= rmb * dz;
530 >          vel[bx] -= rmb * dx;
531 >          vel[by] -= rmb * dy;
532 >          vel[bz] -= rmb * dz;
533  
534            moving[a] = 1;
535            moving[b] = 1;
# Line 534 | Line 545 | void Integrator::constrainA(){
545      }
546  
547      iteration++;
548 +    cerr << "iterainA = " << iteration << "\n";
549    }
550  
551    if( !done ){
552  
553 <    sprintf( painCae.errMsg,
553 >    sprintf( painCave.errMsg,
554               "Constraint failure in constrainA, too many iterations: %d\n",
555 <             iterations );
555 >             iteration );
556      painCave.isFatal = 1;
557      simError();
558    }
# Line 553 | Line 565 | void Integrator::constrainB( void ){
565    int done;
566    double vxab, vyab, vzab;
567    double rxab, ryab, rzab;
568 <  int a, b;
568 >  int a, b, ax, ay, az, bx, by, bz;
569    double rma, rmb;
570    double dx, dy, dz;
571    double rabsq, pabsq, rvab;
# Line 561 | Line 573 | void Integrator::constrainB( void ){
573    double gab;
574    int iteration;
575  
576 <  for(i=0; i<nAtom; i++){
576 >  for(i=0; i<nAtoms; i++){
577      moving[i] = 0;
578      moved[i] = 1;
579    }
580  
581    done = 0;
582 +  iteration = 0;
583    while( !done && (iteration < maxIteration ) ){
584  
585      for(i=0; i<nConstrained; i++){
# Line 574 | Line 587 | void Integrator::constrainB( void ){
587        a = constrainedA[i];
588        b = constrainedB[i];
589  
590 +      ax = 3*a +0;
591 +      ay = 3*a +1;
592 +      az = 3*a +2;
593 +
594 +      bx = 3*b +0;
595 +      by = 3*b +1;
596 +      bz = 3*b +2;
597 +
598        if( moved[a] || moved[b] ){
599          
600 <        vxab = vel[3*a+0] - vel[3*b+0];
601 <        vyab = vel[3*a+1] - vel[3*b+1];
602 <        vzab = vel[3*a+2] - vel[3*b+2];
600 >        vxab = vel[ax] - vel[bx];
601 >        vyab = vel[ay] - vel[by];
602 >        vzab = vel[az] - vel[bz];
603  
604 <        rxab = pos[3*a+0] - pos[3*b+0];q
605 <        ryab = pos[3*a+1] - pos[3*b+1];
606 <        rzab = pos[3*a+2] - pos[3*b+2];
604 >        rxab = pos[ax] - pos[bx];
605 >        ryab = pos[ay] - pos[by];
606 >        rzab = pos[az] - pos[bz];
607          
608          rxab = rxab - info->box_x * copysign(1, rxab)
609 <          * int(rxab / info->box_x + 0.5);
609 >          * (int)( fabs(rxab / info->box_x) + 0.5);
610          ryab = ryab - info->box_y * copysign(1, ryab)
611 <          * int(ryab / info->box_y + 0.5);
611 >          * (int)( fabs(ryab / info->box_y) + 0.5);
612          rzab = rzab - info->box_z * copysign(1, rzab)
613 <          * int(rzab / info->box_z + 0.5);
613 >          * (int)( fabs(rzab / info->box_z) + 0.5);
614  
615          rma = 1.0 / atoms[a]->getMass();
616          rmb = 1.0 / atoms[b]->getMass();
617  
618          rvab = rxab * vxab + ryab * vyab + rzab * vzab;
619            
620 <        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
620 >        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
621  
622          if (fabs(gab) > tol) {
623            
# Line 604 | Line 625 | void Integrator::constrainB( void ){
625            dy = ryab * gab;
626            dz = rzab * gab;
627            
628 <          vel[3*a+0] += rma * dx;
629 <          vel[3*a+1] += rma * dy;
630 <          vel[3*a+2] += rma * dz;
628 >          vel[ax] += rma * dx;
629 >          vel[ay] += rma * dy;
630 >          vel[az] += rma * dz;
631  
632 <          vel[3*b+0] -= rmb * dx;
633 <          vel[3*b+1] -= rmb * dy;
634 <          vel[3*b+2] -= rmb * dz;
632 >          vel[bx] -= rmb * dx;
633 >          vel[by] -= rmb * dy;
634 >          vel[bz] -= rmb * dz;
635            
636            moving[a] = 1;
637            moving[b] = 1;
# Line 630 | Line 651 | void Integrator::constrainB( void ){
651    if( !done ){
652  
653    
654 <    sprintf( painCae.errMsg,
654 >    sprintf( painCave.errMsg,
655               "Constraint failure in constrainB, too many iterations: %d\n",
656 <             iterations );
656 >             iteration );
657      painCave.isFatal = 1;
658      simError();
659    }
# Line 646 | Line 667 | void Integrator::rotate( int axes1, int axes2, double
667  
668  
669   void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
670 <                         double A[3][3] ){
670 >                         double A[9] ){
671  
672    int i,j,k;
673    double sinAngle;
# Line 662 | Line 683 | void Integrator::rotate( int axes1, int axes2, double
683  
684    for(i=0; i<3; i++){
685      for(j=0; j<3; j++){
686 <      tempA[j][i] = A[i][j];
686 >      tempA[j][i] = A[3*i + j];
687      }
688    }
689  
# Line 713 | Line 734 | void Integrator::rotate( int axes1, int axes2, double
734    //            A[][] = A[][] * transpose(rot[][])
735  
736  
737 <  // NOte for as yet unknown reason, we are setting the performing the
737 >  // NOte for as yet unknown reason, we are performing the
738    // calculation as:
739    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
740  
741    for(i=0; i<3; i++){
742      for(j=0; j<3; j++){
743 <      A[j][i] = 0.0;
743 >      A[3*j + i] = 0.0;
744        for(k=0; k<3; k++){
745 <        A[j][i] += tempA[i][k] * rot[j][k];
745 >        A[3*j + i] += tempA[i][k] * rot[j][k];
746        }
747      }
748    }

Diff Legend

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