--- trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/19 22:13:01 1125 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/05/22 18:16:18 1187 @@ -32,7 +32,7 @@ template Integrator::Integrator(SimInfo nAtoms = info->n_atoms; integrableObjects = info->integrableObjects; - + // check for constraints constrainedA = NULL; @@ -45,6 +45,7 @@ template Integrator::Integrator(SimInfo nConstrained = 0; checkConstraints(); + } template Integrator::~Integrator(){ @@ -160,7 +161,7 @@ template void Integrator::integrate(voi double thermalTime = info->thermalTime; double resetTime = info->resetTime; - + double difference; double currSample; double currThermal; double currStatus; @@ -179,6 +180,15 @@ 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->useThermInt) { + myFF->initRestraints(); + } + // initialize the forces before the first step calcForce(1, 1); @@ -211,7 +221,8 @@ template void Integrator::integrate(voi #endif // is_mpi while (info->getTime() < runTime && !stopIntegrator()){ - if ((info->getTime() + dt) >= currStatus){ + difference = info->getTime() + dt - currStatus; + if (difference > 0 || fabs(difference) < 1e-4 ){ calcPot = 1; calcStress = 1; } @@ -244,6 +255,7 @@ template void Integrator::integrate(voi if (info->getTime() >= currStatus){ statOut->writeStat(info->getTime()); + statOut->writeRaw(info->getTime()); calcPot = 0; calcStress = 0; currStatus += statusTime; @@ -266,6 +278,11 @@ template void Integrator::integrate(voi #endif // is_mpi } + // dump out a file containing the omega values for the final configuration + if (info->useThermInt) + myFF->dumpzAngle(); + + delete dumpOut; delete statOut; } @@ -299,10 +316,8 @@ template void Integrator::integrateStep strcpy(checkPointMsg, "Succesful moveA\n"); MPIcheckPoint(); #endif // is_mpi - // calc forces - calcForce(calcPot, calcStress); #ifdef IS_MPI @@ -337,6 +352,7 @@ template void Integrator::moveA(void){ 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); @@ -720,6 +736,7 @@ template void Integrator::rotationPropa // 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