--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/28 16:03:37 829 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/12 20:32:20 1097 @@ -7,6 +7,10 @@ #include #endif //is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif // profile + #include "Integrator.hpp" #include "simError.h" @@ -25,8 +29,9 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - + nAtoms = info->n_atoms; + integrableObjects = info->integrableObjects; // check for constraints @@ -64,7 +69,8 @@ template void Integrator::checkConstrai SRI** theArray; for (int i = 0; i < nMols; i++){ - theArray = (SRI * *) molecules[i].getMyBonds(); + + theArray = (SRI * *) molecules[i].getMyBonds(); for (int j = 0; j < molecules[i].getNBonds(); j++){ constrained = theArray[j]->is_constrained(); @@ -110,6 +116,7 @@ template void Integrator::checkConstrai } } + if (nConstrained > 0){ isConstrained = 1; @@ -131,7 +138,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; @@ -158,7 +165,7 @@ template void Integrator::integrate(voi double currThermal; double currStatus; double currReset; - + int calcPot, calcStress; tStats = new Thermo(info); @@ -176,10 +183,13 @@ template void Integrator::integrate(voi calcForce(1, 1); + //temp test + tStats->getPotential(); + if (nConstrained){ preMove(); constrainA(); - calcForce(1, 1); + calcForce(1, 1); constrainB(); } @@ -198,7 +208,6 @@ template void Integrator::integrate(voi statOut->writeStat(info->getTime()); - #ifdef IS_MPI strcpy(checkPointMsg, "The integrator is ready to go."); MPIcheckPoint(); @@ -210,8 +219,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){ @@ -227,11 +246,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){ @@ -239,6 +258,10 @@ template void Integrator::integrate(voi currReset += resetTime; } } + +#ifdef PROFILE + endProfile( pro2 ); +#endif //profile #ifdef IS_MPI strcpy(checkPointMsg, "successfully took a time step."); @@ -246,10 +269,6 @@ template void Integrator::integrate(voi #endif // is_mpi } - - // write the last frame - dumpOut->writeDump(info->getTime()); - delete dumpOut; delete statOut; } @@ -257,13 +276,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(); @@ -279,13 +313,20 @@ 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"); MPIcheckPoint(); @@ -294,18 +335,18 @@ template void Integrator::moveA(void){ template void Integrator::moveA(void){ - int i, j; + size_t i, j; DirectionalAtom* dAtom; double Tb[3], ji[3]; double vel[3], pos[3], frc[3]; double mass; - - for (i = 0; i < nAtoms; i++){ - atoms[i]->getVel(vel); - atoms[i]->getPos(pos); - atoms[i]->getFrc(frc); - - mass = atoms[i]->getMass(); + + for (i = 0; i < integrableObjects.size() ; i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getPos(pos); + integrableObjects[i]->getFrc(frc); + + mass = integrableObjects[i]->getMass(); for (j = 0; j < 3; j++){ // velocity half step @@ -314,27 +355,26 @@ template void Integrator::moveA(void){ pos[j] += dt * vel[j]; } - atoms[i]->setVel(vel); - atoms[i]->setPos(pos); + integrableObjects[i]->setVel(vel); + integrableObjects[i]->setPos(pos); - if (atoms[i]->isDirectional()){ - dAtom = (DirectionalAtom *) atoms[i]; + if (integrableObjects[i]->isDirectional()){ // get and convert the torque to body frame - dAtom->getTrq(Tb); - dAtom->lab2Body(Tb); + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); // get the angular momentum, and propagate a half step - dAtom->getJ(ji); + integrableObjects[i]->getJ(ji); for (j = 0; j < 3; j++) ji[j] += (dt2 * Tb[j]) * eConvert; - this->rotationPropagation( dAtom, ji ); + this->rotationPropagation( integrableObjects[i], ji ); - dAtom->setJ(ji); + integrableObjects[i]->setJ(ji); } } @@ -346,40 +386,38 @@ template void Integrator::moveB(void){ template void Integrator::moveB(void){ int i, j; - DirectionalAtom* dAtom; double Tb[3], ji[3]; double vel[3], frc[3]; double mass; - for (i = 0; i < nAtoms; i++){ - atoms[i]->getVel(vel); - atoms[i]->getFrc(frc); + for (i = 0; i < integrableObjects.size(); i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getFrc(frc); - mass = atoms[i]->getMass(); + mass = integrableObjects[i]->getMass(); // velocity half step for (j = 0; j < 3; j++) vel[j] += (dt2 * frc[j] / mass) * eConvert; - atoms[i]->setVel(vel); + integrableObjects[i]->setVel(vel); - if (atoms[i]->isDirectional()){ - dAtom = (DirectionalAtom *) atoms[i]; + if (integrableObjects[i]->isDirectional()){ - // get and convert the torque to body frame + // get and convert the torque to body frame - dAtom->getTrq(Tb); - dAtom->lab2Body(Tb); + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); // get the angular momentum, and propagate a half step - dAtom->getJ(ji); + integrableObjects[i]->getJ(ji); for (j = 0; j < 3; j++) ji[j] += (dt2 * Tb[j]) * eConvert; - dAtom->setJ(ji); + integrableObjects[i]->setJ(ji); } } @@ -648,7 +686,7 @@ template void Integrator::rotationPropa } template void Integrator::rotationPropagation -( DirectionalAtom* dAtom, double ji[3] ){ +( StuntDouble* sd, double ji[3] ){ double angle; double A[3][3], I[3][3]; @@ -656,30 +694,30 @@ template void Integrator::rotationPropa // use the angular velocities to propagate the rotation matrix a // full time step - dAtom->getA(A); - dAtom->getI(I); - - // rotate about the x-axis + sd->getA(A); + sd->getI(I); + + // 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 ); + + sd->setA( A ); } template void Integrator::rotate(int axes1, int axes2, @@ -747,7 +785,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 +814,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(); +}