--- trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/20 16:56:40 1127 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/05/27 20:06:38 1207 @@ -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; @@ -183,6 +184,11 @@ template void Integrator::integrate(voi // 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); @@ -215,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; } @@ -248,6 +255,8 @@ template void Integrator::integrate(voi if (info->getTime() >= currStatus){ statOut->writeStat(info->getTime()); + if (info->useThermInt) + statOut->writeRaw(info->getTime()); calcPot = 0; calcStress = 0; currStatus += statusTime; @@ -270,6 +279,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; } @@ -304,9 +318,7 @@ template void Integrator::integrateStep MPIcheckPoint(); #endif // is_mpi - // calc forces - calcForce(calcPot, calcStress); #ifdef IS_MPI @@ -341,6 +353,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); @@ -724,6 +737,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