--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/29 00:19:10 837 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/20 16:56:40 1127 @@ -7,6 +7,10 @@ #include #endif //is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif // profile + #include "Integrator.hpp" #include "simError.h" @@ -27,6 +31,7 @@ template Integrator::Integrator(SimInfo } 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; @@ -172,17 +179,21 @@ template void Integrator::integrate(voi readyCheck(); + // remove center of mass drift velocity (in case we passed in a configuration + // that was drifting + tStats->removeCOMdrift(); + // initialize the forces before the first step calcForce(1, 1); - + if (nConstrained){ preMove(); constrainA(); calcForce(1, 1); constrainB(); } - + if (info->setTemp){ thermalize(); } @@ -198,20 +209,29 @@ template void Integrator::integrate(voi statOut->writeStat(info->getTime()); - #ifdef IS_MPI strcpy(checkPointMsg, "The integrator is ready to go."); MPIcheckPoint(); #endif // is_mpi - while (info->getTime() < runTime){ + while (info->getTime() < runTime && !stopIntegrator()){ if ((info->getTime() + dt) >= currStatus){ calcPot = 1; calcStress = 1; } +#ifdef PROFILE + startProfile( pro1 ); +#endif + integrateStep(calcPot, calcStress); +#ifdef PROFILE + endProfile( pro1 ); + + startProfile( pro2 ); +#endif // profile + info->incrTime(dt); if (info->setTemp){ @@ -239,6 +259,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 +270,6 @@ template void Integrator::integrate(voi #endif // is_mpi } - - // write the last frame - dumpOut->writeDump(info->getTime()); - delete dumpOut; delete statOut; } @@ -257,11 +277,26 @@ 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 @@ -279,13 +314,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,19 +336,19 @@ 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 < integrableObjects.size() ; i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getPos(pos); + integrableObjects[i]->getFrc(frc); + + mass = integrableObjects[i]->getMass(); - for (i = 0; i < nAtoms; i++){ - atoms[i]->getVel(vel); - atoms[i]->getPos(pos); - atoms[i]->getFrc(frc); - - mass = atoms[i]->getMass(); - for (j = 0; j < 3; j++){ // velocity half step vel[j] += (dt2 * frc[j] / mass) * eConvert; @@ -314,27 +356,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 +387,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 - 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,38 +687,55 @@ 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]; + int i, j, k; // use the angular velocities to propagate the rotation matrix a // full time step - - dAtom->getA(A); - dAtom->getI(I); - // rotate about the x-axis - angle = dt2 * ji[0] / I[0][0]; - this->rotate( 1, 2, angle, ji, A ); + sd->getA(A); + sd->getI(I); - // 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); + if (sd->isLinear()) { + i = sd->linearAxis(); + j = (i+1)%3; + k = (i+2)%3; + + angle = dt2 * ji[j] / I[j][j]; + this->rotate( k, i, angle, ji, A ); - // rotate about the y-axis - angle = dt2 * ji[1] / I[1][1]; - this->rotate( 2, 0, angle, ji, A ); + angle = dt * ji[k] / I[k][k]; + this->rotate( i, j, angle, ji, A); - // rotate about the x-axis - angle = dt2 * ji[0] / I[0][0]; - this->rotate( 1, 2, angle, ji, A ); + angle = dt2 * ji[j] / I[j][j]; + this->rotate( k, i, angle, ji, A ); - dAtom->setA( A ); + } else { + // rotate about the x-axis + angle = dt2 * ji[0] / I[0][0]; + 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 ); + + } + sd->setA( A ); } template void Integrator::rotate(int axes1, int axes2,