--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/28 16:03:37 829 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/06/04 03:15:31 1234 @@ -1,12 +1,16 @@ #include #include #include - +#include "Rattle.hpp" #ifdef IS_MPI #include "mpiSimulation.hpp" #include #endif //is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif // profile + #include "Integrator.hpp" #include "simError.h" @@ -25,9 +29,20 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - + nAtoms = info->n_atoms; + integrableObjects = info->integrableObjects; + rattle = new RattleFramework(info); + + if(rattle == NULL){ + sprintf(painCave.errMsg, + "Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); + painCave.isFatal = 1; + simError(); + } + +/* // check for constraints constrainedA = NULL; @@ -40,9 +55,13 @@ template Integrator::Integrator(SimInfo nConstrained = 0; checkConstraints(); +*/ } template Integrator::~Integrator(){ + if (rattle != NULL) + delete rattle; +/* if (nConstrained){ delete[] constrainedA; delete[] constrainedB; @@ -51,8 +70,10 @@ template Integrator::~Integrator(){ delete[] moved; delete[] oldPos; } +*/ } +/* template void Integrator::checkConstraints(void){ isConstrained = 0; @@ -64,7 +85,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 +132,7 @@ template void Integrator::checkConstrai } } + if (nConstrained > 0){ isConstrained = 1; @@ -131,7 +154,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; @@ -143,8 +166,8 @@ template void Integrator::checkConstrai delete[] temp_con; } +*/ - template void Integrator::integrate(void){ double runTime = info->run_time; @@ -153,12 +176,12 @@ 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; tStats = new Thermo(info); @@ -172,16 +195,24 @@ 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 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(); - } + //execute constraint algorithm to make sure at the very beginning the system is constrained + rattle->doPreConstraint(); + rattle->doRattleA(); + calcForce(1, 1); + rattle->doRattleB(); if (info->setTemp){ thermalize(); @@ -198,19 +229,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){ - 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); @@ -227,11 +268,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 +280,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 +291,11 @@ template void Integrator::integrate(voi #endif // is_mpi } + // dump out a file containing the omega values for the final configuration + if (info->useSolidThermInt && !info->useLiquidThermInt) + myFF->dumpzAngle(); + - // write the last frame - dumpOut->writeDump(info->getTime()); - delete dumpOut; delete statOut; } @@ -257,21 +303,35 @@ template void Integrator::integrateStep template void Integrator::integrateStep(int calcPot, int calcStress){ // Position full step, and velocity half step - preMove(); - moveA(); +#ifdef PROFILE + startProfile(pro3); +#endif //profile + //save old state (position, velocity etc) + rattle->doPreConstraint(); +#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(); #endif // is_mpi - // calc forces - calcForce(calcPot, calcStress); #ifdef IS_MPI @@ -279,13 +339,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 +361,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 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; @@ -314,80 +382,74 @@ 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); } } - if (nConstrained){ - constrainA(); - } + rattle->doRattleA(); } 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]; - - // get and convert the torque to body frame + if (integrableObjects[i]->isDirectional()){ - dAtom->getTrq(Tb); - dAtom->lab2Body(Tb); + // get and convert the torque to body frame + 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(); - } + rattle->doRattleB(); } +/* template void Integrator::preMove(void){ int i, j; double pos[3]; @@ -646,40 +708,58 @@ template void Integrator::constrainB(vo simError(); } } - +*/ 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 ); - - // 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->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, @@ -747,7 +827,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 +856,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(); +}