--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/09/23 20:34:31 782 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/01/30 15:01:09 999 @@ -1,12 +1,16 @@ #include -#include -#include +#include +#include #ifdef IS_MPI #include "mpiSimulation.hpp" #include #endif //is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif // profile + #include "Integrator.hpp" #include "simError.h" @@ -25,7 +29,7 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - + nAtoms = info->n_atoms; // check for constraints @@ -131,7 +135,7 @@ template void Integrator::checkConstrai } - // save oldAtoms to check for lode balanceing later on. + // save oldAtoms to check for lode balancing later on. oldAtoms = nAtoms; @@ -146,7 +150,6 @@ template void Integrator::integrate(voi template void Integrator::integrate(void){ - int i, j; // loop counters double runTime = info->run_time; double sampleTime = info->sampleTime; @@ -159,20 +162,20 @@ template void Integrator::integrate(voi double currThermal; double currStatus; double currReset; - + int calcPot, calcStress; - int isError; tStats = new Thermo(info); statOut = new StatWriter(info); dumpOut = new DumpWriter(info); atoms = info->atoms; - DirectionalAtom* dAtom; dt = info->dt; dt2 = 0.5 * dt; + readyCheck(); + // initialize the forces before the first step calcForce(1, 1); @@ -180,7 +183,7 @@ template void Integrator::integrate(voi if (nConstrained){ preMove(); constrainA(); - calcForce(1, 1); + calcForce(1, 1); constrainB(); } @@ -198,7 +201,6 @@ template void Integrator::integrate(voi dumpOut->writeDump(info->getTime()); statOut->writeStat(info->getTime()); - readyCheck(); #ifdef IS_MPI strcpy(checkPointMsg, "The integrator is ready to go."); @@ -211,8 +213,18 @@ template void Integrator::integrate(voi calcStress = 1; } +#ifdef PROFILE + startProfile( pro1 ); +#endif + integrateStep(calcPot, calcStress); + +#ifdef PROFILE + endProfile( pro1 ); + startProfile( pro2 ); +#endif // profile + info->incrTime(dt); if (info->setTemp){ @@ -228,11 +240,11 @@ template void Integrator::integrate(voi } if (info->getTime() >= currStatus){ - statOut->writeStat(info->getTime()); - calcPot = 0; + statOut->writeStat(info->getTime()); + calcPot = 0; calcStress = 0; currStatus += statusTime; - } + } if (info->resetIntegrator){ if (info->getTime() >= currReset){ @@ -240,6 +252,10 @@ template void Integrator::integrate(voi currReset += resetTime; } } + +#ifdef PROFILE + endProfile( pro2 ); +#endif //profile #ifdef IS_MPI strcpy(checkPointMsg, "successfully took a time step."); @@ -247,8 +263,6 @@ template void Integrator::integrate(voi #endif // is_mpi } - dumpOut->writeFinal(info->getTime()); - delete dumpOut; delete statOut; } @@ -256,13 +270,28 @@ template void Integrator::integrateStep template void Integrator::integrateStep(int calcPot, int calcStress){ // Position full step, and velocity half step + +#ifdef PROFILE + startProfile(pro3); +#endif //profile + preMove(); - moveA(); +#ifdef PROFILE + endProfile(pro3); + startProfile(pro4); +#endif // profile + moveA(); +#ifdef PROFILE + endProfile(pro4); + + startProfile(pro5); +#endif//profile + #ifdef IS_MPI strcpy(checkPointMsg, "Succesful moveA\n"); MPIcheckPoint(); @@ -278,12 +307,19 @@ template void Integrator::integrateStep MPIcheckPoint(); #endif // is_mpi +#ifdef PROFILE + endProfile( pro5 ); + startProfile( pro6 ); +#endif //profile + // finish the velocity half step moveB(); - +#ifdef PROFILE + endProfile(pro6); +#endif // profile #ifdef IS_MPI strcpy(checkPointMsg, "Succesful moveB\n"); @@ -365,7 +401,7 @@ template void Integrator::moveB(void){ if (atoms[i]->isDirectional()){ dAtom = (DirectionalAtom *) atoms[i]; - // get and convert the torque to body frame + // get and convert the torque to body frame dAtom->getTrq(Tb); dAtom->lab2Body(Tb); @@ -403,7 +439,7 @@ template void Integrator::constrainA(){ } template void Integrator::constrainA(){ - int i, j, k; + int i, j; int done; double posA[3], posB[3]; double velA[3], velB[3]; @@ -547,7 +583,7 @@ template void Integrator::constrainB(vo } template void Integrator::constrainB(void){ - int i, j, k; + int i, j; int done; double posA[3], posB[3]; double velA[3], velB[3]; @@ -556,8 +592,7 @@ template void Integrator::constrainB(vo int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; - double rabsq, pabsq, rvab; - double diffsq; + double rvab; double gab; int iteration; @@ -658,28 +693,28 @@ template void Integrator::rotationPropa dAtom->getA(A); dAtom->getI(I); - - // rotate about the x-axis + + // rotate about the x-axis angle = dt2 * ji[0] / I[0][0]; - this->rotate( 1, 2, angle, ji, A ); - + this->rotate( 1, 2, angle, ji, A ); + // rotate about the y-axis angle = dt2 * ji[1] / I[1][1]; this->rotate( 2, 0, angle, ji, A ); - + // rotate about the z-axis angle = dt * ji[2] / I[2][2]; this->rotate( 0, 1, angle, ji, A); - + // rotate about the y-axis angle = dt2 * ji[1] / I[1][1]; this->rotate( 2, 0, angle, ji, A ); - + // rotate about the x-axis angle = dt2 * ji[0] / I[0][0]; this->rotate( 1, 2, angle, ji, A ); - - dAtom->setA( A ); + + dAtom->setA( A ); } template void Integrator::rotate(int axes1, int axes2, @@ -747,7 +782,7 @@ template void Integrator::rotate(int ax } } - // rotate the Rotation matrix acording to: + // rotate the Rotation matrix acording to: // A[][] = A[][] * transpose(rot[][]) @@ -776,3 +811,9 @@ template double Integrator::getConserve template double Integrator::getConservedQuantity(void){ return tStats->getTotalE(); } +template string Integrator::getAdditionalParameters(void){ + //By default, return a null string + //The reason we use string instead of char* is that if we use char*, we will + //return a pointer point to local variable which might cause problem + return string(); +}