--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/09/15 16:52:02 763 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/06/02 14:56:18 1221 @@ -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,10 +29,10 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - info->the_integrator = this; nAtoms = info->n_atoms; - + integrableObjects = info->integrableObjects; + // check for constraints constrainedA = NULL; @@ -41,6 +45,7 @@ template Integrator::Integrator(SimInfo nConstrained = 0; checkConstraints(); + } template Integrator::~Integrator(){ @@ -65,7 +70,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(); @@ -111,6 +117,7 @@ template void Integrator::checkConstrai } } + if (nConstrained > 0){ isConstrained = 1; @@ -132,7 +139,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; @@ -147,7 +154,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; @@ -155,29 +161,45 @@ template void Integrator::integrate(voi double thermalTime = info->thermalTime; double resetTime = info->resetTime; - + double difference; double currSample; 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(); + + // remove center of mass drift velocity (in case we passed in a configuration + // that was drifting + tStats->removeCOMdrift(); + + // initialize the retraints if necessary + if (info->useSolidThermInt && !info->useLiquidThermInt) { + myFF->initRestraints(); + } // initialize the forces before the first step calcForce(1, 1); + if (nConstrained){ + preMove(); + constrainA(); + calcForce(1, 1); + constrainB(); + } + if (info->setTemp){ thermalize(); } @@ -192,21 +214,31 @@ 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."); MPIcheckPoint(); #endif // is_mpi - while (info->getTime() < runTime){ - if ((info->getTime() + dt) >= currStatus){ + while (info->getTime() < runTime && !stopIntegrator()){ + difference = info->getTime() + dt - currStatus; + if (difference > 0 || fabs(difference) < 1e-4 ){ 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){ @@ -222,11 +254,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){ @@ -234,6 +266,10 @@ template void Integrator::integrate(voi currReset += resetTime; } } + +#ifdef PROFILE + endProfile( pro2 ); +#endif //profile #ifdef IS_MPI strcpy(checkPointMsg, "successfully took a time step."); @@ -241,7 +277,10 @@ template void Integrator::integrate(voi #endif // is_mpi } - dumpOut->writeFinal(info->getTime()); + // dump out a file containing the omega values for the final configuration + if (info->useSolidThermInt && !info->useLiquidThermInt) + myFF->dumpzAngle(); + delete dumpOut; delete statOut; @@ -250,13 +289,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(); +#ifdef PROFILE + endProfile(pro3); + + startProfile(pro4); +#endif // profile + moveA(); - if (nConstrained){ - constrainA(); - } +#ifdef PROFILE + endProfile(pro4); + + startProfile(pro5); +#endif//profile #ifdef IS_MPI @@ -264,9 +316,7 @@ template void Integrator::integrateStep MPIcheckPoint(); #endif // is_mpi - // calc forces - calcForce(calcPot, calcStress); #ifdef IS_MPI @@ -274,14 +324,19 @@ template void Integrator::integrateStep MPIcheckPoint(); #endif // is_mpi +#ifdef PROFILE + endProfile( pro5 ); + startProfile( pro6 ); +#endif //profile + // finish the velocity half step moveB(); - if (nConstrained){ - constrainB(); - } +#ifdef PROFILE + endProfile(pro6); +#endif // profile #ifdef IS_MPI strcpy(checkPointMsg, "Succesful moveB\n"); @@ -291,21 +346,20 @@ 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 A[3][3], I[3][3]; - double angle; double vel[3], pos[3], frc[3]; double mass; + double omega; + + 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; @@ -313,95 +367,75 @@ 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; - // use the angular velocities to propagate the rotation matrix a - // full time step + this->rotationPropagation( integrableObjects[i], ji ); - 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); - - // 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->setJ(ji); - dAtom->setA(A); + integrableObjects[i]->setJ(ji); } } + + if (nConstrained){ + constrainA(); + } } 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); } } + + if (nConstrained){ + constrainB(); + } } template void Integrator::preMove(void){ @@ -420,7 +454,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]; @@ -560,10 +594,11 @@ template void Integrator::constrainA(){ painCave.isFatal = 1; simError(); } + } template void Integrator::constrainB(void){ - int i, j, k; + int i, j; int done; double posA[3], posB[3]; double velA[3], velB[3]; @@ -572,8 +607,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; @@ -663,6 +697,59 @@ template void Integrator::rotate(int ax } } +template void Integrator::rotationPropagation +( 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 + + sd->getA(A); + sd->getI(I); + + 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 ); + + angle = dt * ji[k] / I[k][k]; + this->rotate( i, j, angle, ji, A); + + angle = dt2 * ji[j] / I[j][j]; + this->rotate( k, i, angle, ji, 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]; + sd->addZangle(angle); + 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, double angle, double ji[3], double A[3][3]){ @@ -728,7 +815,7 @@ template void Integrator::rotate(int ax } } - // rotate the Rotation matrix acording to: + // rotate the Rotation matrix acording to: // A[][] = A[][] * transpose(rot[][]) @@ -756,4 +843,10 @@ template double Integrator::getConserve template double Integrator::getConservedQuantity(void){ return tStats->getTotalE(); -} \ No newline at end of file +} +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(); +}