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 559 by mmeineke, Thu Jun 19 22:02:44 2003 UTC vs.
Revision 561 by mmeineke, Fri Jun 20 20:29:36 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 >    delete[] oldPos;
53    }
54    
55   }
# Line 146 | Line 147 | void Integrator::checkConstraints( void ){
147      moving = new int[nAtoms];
148      moved  = new int[nAtoms];
149  
150 <    prePos = new double[nAtoms*3];
150 >    oldPos = new double[nAtoms*3];
151    }
152    
153    delete[] temp_con;
# Line 156 | Line 157 | void Integrator::integrate( void ){
157   void Integrator::integrate( void ){
158  
159    int i, j;                         // loop counters
159 –  double kE = 0.0;                  // the kinetic energy  
160 –  double rot_kE;
161 –  double trans_kE;
162 –  int tl;                        // the time loop conter
163 –  double dt2;                       // half the dt
160  
165 –  double vx, vy, vz;    // the velocities
166 –  double vx2, vy2, vz2; // the square of the velocities
167 –  double rx, ry, rz;    // the postitions
168 –  
169 –  double ji[3];   // the body frame angular momentum
170 –  double jx2, jy2, jz2; // the square of the angular momentums
171 –  double Tb[3];   // torque in the body frame
172 –  double angle;   // the angle through which to rotate the rotation matrix
173 –  double A[3][3]; // the rotation matrix
174 –  double press[9];
175 –
176 –  double dt          = info->dt;
161    double runTime     = info->run_time;
162    double sampleTime  = info->sampleTime;
163    double statusTime  = info->statusTime;
# Line 187 | Line 171 | void Integrator::integrate( void ){
171    int calcPot, calcStress;
172    int isError;
173  
174 +
175 +
176    tStats   = new Thermo( info );
177 <  e_out    = new StatWriter( info );
178 <  dump_out = new DumpWriter( info );
177 >  statOut  = new StatWriter( info );
178 >  dumpOut  = new DumpWriter( info );
179  
180 <  Atom** atoms = info->atoms;
180 >  atoms = info->atoms;
181    DirectionalAtom* dAtom;
182 +
183 +  dt = info->dt;
184    dt2 = 0.5 * dt;
185  
186    // initialize the forces before the first step
# Line 204 | Line 192 | void Integrator::integrate( void ){
192      tStats->velocitize();
193    }
194    
195 <  dump_out->writeDump( 0.0 );
196 <  e_out->writeStat( 0.0 );
195 >  dumpOut->writeDump( 0.0 );
196 >  statOut->writeStat( 0.0 );
197    
198    calcPot     = 0;
199    calcStress  = 0;
# Line 223 | Line 211 | void Integrator::integrate( void ){
211    MPIcheckPoint();
212   #endif // is_mpi
213  
214 +
215 +  pos  = Atom::getPosArray();
216 +  vel  = Atom::getVelArray();
217 +  frc  = Atom::getFrcArray();
218 +  trq  = Atom::getTrqArray();
219 +  Amat = Atom::getAmatArray();
220 +
221    while( currTime < runTime ){
222  
223      if( (currTime+dt) >= currStatus ){
224        calcPot = 1;
225        calcStress = 1;
226      }
227 <    
227 >
228      integrateStep( calcPot, calcStress );
229        
230      currTime += dt;
# Line 242 | Line 237 | void Integrator::integrate( void ){
237      }
238  
239      if( currTime >= currSample ){
240 <      dump_out->writeDump( currTime );
240 >      dumpOut->writeDump( currTime );
241        currSample += sampleTime;
242      }
243  
244      if( currTime >= currStatus ){
245 <      e_out->writeStat( time * dt );
245 >      statOut->writeStat( currTime );
246        calcPot = 0;
247        calcStress = 0;
248        currStatus += statusTime;
# Line 261 | Line 256 | void Integrator::integrate( void ){
256  
257    }
258  
259 <  dump_out->writeFinal();
259 >  dumpOut->writeFinal();
260  
261 <  delete dump_out;
262 <  delete e_out;
261 >  delete dumpOut;
262 >  delete statOut;
263   }
264  
265   void Integrator::integrateStep( int calcPot, int calcStress ){
266  
267 +
268 +      
269    // Position full step, and velocity half step
270  
271 <  //preMove();
271 >  preMove();
272    moveA();
273    if( nConstrained ) constrainA();
274  
# Line 294 | Line 291 | void Integrator::moveA( void ){
291    DirectionalAtom* dAtom;
292    double Tb[3];
293    double ji[3];
294 +  double angle;
295  
296    for( i=0; i<nAtoms; i++ ){
297      atomIndex = i * 3;
# Line 331 | Line 329 | void Integrator::moveA( void ){
329        
330        // rotate about the x-axis      
331        angle = dt2 * ji[0] / dAtom->getIxx();
332 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
332 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
333        
334        // rotate about the y-axis
335        angle = dt2 * ji[1] / dAtom->getIyy();
336 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
336 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
337        
338        // rotate about the z-axis
339        angle = dt * ji[2] / dAtom->getIzz();
340 <      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
340 >      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
341        
342        // rotate about the y-axis
343        angle = dt2 * ji[1] / dAtom->getIyy();
344 <      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
344 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
345        
346         // rotate about the x-axis
347        angle = dt2 * ji[0] / dAtom->getIxx();
348 <      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
348 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
349        
350        dAtom->setJx( ji[0] );
351        dAtom->setJy( ji[1] );
# Line 391 | Line 389 | void Integrator::moveB( void ){
389        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
390        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
391        
394 –      jx2 = ji[0] * ji[0];
395 –      jy2 = ji[1] * ji[1];
396 –      jz2 = ji[2] * ji[2];
397 –      
392        dAtom->setJx( ji[0] );
393        dAtom->setJy( ji[1] );
394        dAtom->setJz( ji[2] );
# Line 407 | Line 401 | void Integrator::preMove( void ){
401    int i;
402  
403    if( nConstrained ){
404 <    if( oldAtoms != nAtoms ){
404 >
405 > //    if( oldAtoms != nAtoms ){
406        
407 <      // save oldAtoms to check for lode balanceing later on.
407 > //       // save oldAtoms to check for lode balanceing later on.
408        
409 <      oldAtoms = nAtoms;
409 > //       oldAtoms = nAtoms;
410        
411 <      delete[] moving;
412 <      delete[] moved;
413 <      delete[] oldPos;
411 > //       delete[] moving;
412 > //       delete[] moved;
413 > //       delete[] oldPos;
414        
415 <      moving = new int[nAtoms];
416 <      moved  = new int[nAtoms];
415 > //       moving = new int[nAtoms];
416 > //       moved  = new int[nAtoms];
417        
418 <      oldPos = new double[nAtoms*3];
419 <    }
418 > //       oldPos = new double[nAtoms*3];
419 > //     }
420    
421      for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
422    }
# Line 436 | Line 431 | void Integrator::constrainA(){
431    int a, b;
432    double rma, rmb;
433    double dx, dy, dz;
434 +  double rpab;
435    double rabsq, pabsq, rpabsq;
436    double diffsq;
437    double gab;
# Line 468 | Line 464 | void Integrator::constrainA(){
464  
465          //periodic boundary condition
466          pxab = pxab - info->box_x * copysign(1, pxab)
467 <          * int(pxab / info->box_x + 0.5);
467 >          * int( fabs(pxab) / info->box_x + 0.5);
468          pyab = pyab - info->box_y * copysign(1, pyab)
469 <          * int(pyab / info->box_y + 0.5);
469 >          * int( fabs(pyab) / info->box_y + 0.5);
470          pzab = pzab - info->box_z * copysign(1, pzab)
471 <          * int(pzab / info->box_z + 0.5);
471 >          * int( fabs(pzab) / info->box_z + 0.5);
472        
473          pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
474 <        rabsq = constraintedDsqr[i];
474 >        rabsq = constrainedDsqr[i];
475          diffsq = pabsq - rabsq;
476  
477          // the original rattle code from alan tidesley
# Line 485 | Line 481 | void Integrator::constrainA(){
481            rzab = oldPos[3*a+2] - oldPos[3*b+2];
482  
483            rxab = rxab - info->box_x * copysign(1, rxab)
484 <            * int(rxab / info->box_x + 0.5);
484 >            * int( fabs(rxab) / info->box_x + 0.5);
485            ryab = ryab - info->box_y * copysign(1, ryab)
486 <            * int(ryab / info->box_y + 0.5);
486 >            * int( fabs(ryab) / info->box_y + 0.5);
487            rzab = rzab - info->box_z * copysign(1, rzab)
488 <            * int(rzab / info->box_z + 0.5);
488 >            * int( fabs(rzab) / info->box_z + 0.5);
489  
490            rpab = rxab * pxab + ryab * pyab + rzab * pzab;
491            rpabsq = rpab * rpab;
# Line 553 | Line 549 | void Integrator::constrainA(){
549  
550    if( !done ){
551  
552 <    sprintf( painCae.errMsg,
552 >    sprintf( painCave.errMsg,
553               "Constraint failure in constrainA, too many iterations: %d\n",
554 <             iterations );
554 >             iteration );
555      painCave.isFatal = 1;
556      simError();
557    }
# Line 576 | Line 572 | void Integrator::constrainB( void ){
572    double gab;
573    int iteration;
574  
575 <  for(i=0; i<nAtom; i++){
575 >  for(i=0; i<nAtoms; i++){
576      moving[i] = 0;
577      moved[i] = 1;
578    }
579  
580    done = 0;
581 +  iteration = 0;
582    while( !done && (iteration < maxIteration ) ){
583  
584      for(i=0; i<nConstrained; i++){
# Line 595 | Line 592 | void Integrator::constrainB( void ){
592          vyab = vel[3*a+1] - vel[3*b+1];
593          vzab = vel[3*a+2] - vel[3*b+2];
594  
595 <        rxab = pos[3*a+0] - pos[3*b+0];q
595 >        rxab = pos[3*a+0] - pos[3*b+0];
596          ryab = pos[3*a+1] - pos[3*b+1];
597          rzab = pos[3*a+2] - pos[3*b+2];
598          
599          rxab = rxab - info->box_x * copysign(1, rxab)
600 <          * int(rxab / info->box_x + 0.5);
600 >          * int( fabs(rxab) / info->box_x + 0.5);
601          ryab = ryab - info->box_y * copysign(1, ryab)
602 <          * int(ryab / info->box_y + 0.5);
602 >          * int( fabs(ryab) / info->box_y + 0.5);
603          rzab = rzab - info->box_z * copysign(1, rzab)
604 <          * int(rzab / info->box_z + 0.5);
604 >          * int( fabs(rzab) / info->box_z + 0.5);
605  
606          rma = 1.0 / atoms[a]->getMass();
607          rmb = 1.0 / atoms[b]->getMass();
608  
609          rvab = rxab * vxab + ryab * vyab + rzab * vzab;
610            
611 <        gab = -rvab / ( ( rma + rmb ) * constraintsDsqr[i] );
611 >        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
612  
613          if (fabs(gab) > tol) {
614            
# Line 645 | Line 642 | void Integrator::constrainB( void ){
642    if( !done ){
643  
644    
645 <    sprintf( painCae.errMsg,
645 >    sprintf( painCave.errMsg,
646               "Constraint failure in constrainB, too many iterations: %d\n",
647 <             iterations );
647 >             iteration );
648      painCave.isFatal = 1;
649      simError();
650    }
# Line 661 | Line 658 | void Integrator::rotate( int axes1, int axes2, double
658  
659  
660   void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
661 <                         double A[3][3] ){
661 >                         double A[9] ){
662  
663    int i,j,k;
664    double sinAngle;
# Line 677 | Line 674 | void Integrator::rotate( int axes1, int axes2, double
674  
675    for(i=0; i<3; i++){
676      for(j=0; j<3; j++){
677 <      tempA[j][i] = A[i][j];
677 >      tempA[j][i] = A[3*i + j];
678      }
679    }
680  
# Line 728 | Line 725 | void Integrator::rotate( int axes1, int axes2, double
725    //            A[][] = A[][] * transpose(rot[][])
726  
727  
728 <  // NOte for as yet unknown reason, we are setting the performing the
728 >  // NOte for as yet unknown reason, we are performing the
729    // calculation as:
730    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
731  
732    for(i=0; i<3; i++){
733      for(j=0; j<3; j++){
734 <      A[j][i] = 0.0;
734 >      A[3*j + i] = 0.0;
735        for(k=0; k<3; k++){
736 <        A[j][i] += tempA[i][k] * rot[j][k];
736 >        A[3*j + i] += tempA[i][k] * rot[j][k];
737        }
738      }
739    }

Diff Legend

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