--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/16 19:16:24 804 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/29 17:55:28 841 @@ -1,6 +1,6 @@ #include -#include -#include +#include +#include #ifdef IS_MPI #include "mpiSimulation.hpp" @@ -25,7 +25,7 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - + nAtoms = info->n_atoms; // check for constraints @@ -158,7 +158,7 @@ template void Integrator::integrate(voi double currThermal; double currStatus; double currReset; - + int calcPot, calcStress; tStats = new Thermo(info); @@ -174,15 +174,20 @@ template void Integrator::integrate(voi // initialize the forces before the first step + std::cerr << "Before initial Force calc\n"; + calcForce(1, 1); if (nConstrained){ preMove(); constrainA(); - calcForce(1, 1); + calcForce(1, 1); constrainB(); + std::cerr << "premove done\n"; } - + + + if (info->setTemp){ thermalize(); } @@ -227,11 +232,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 +244,8 @@ template void Integrator::integrate(voi currReset += resetTime; } } + + std::cerr << "done with time = " << info->getTime() << "\n"; #ifdef IS_MPI strcpy(checkPointMsg, "successfully took a time step."); @@ -366,7 +373,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); @@ -658,28 +665,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 +754,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 +783,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(); +}