--- trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/12 20:32:20 1097 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/06/21 18:52:21 1284 @@ -1,7 +1,8 @@ #include #include #include - +#include "Rattle.hpp" +#include "Roll.hpp" #ifdef IS_MPI #include "mpiSimulation.hpp" #include @@ -33,6 +34,16 @@ template Integrator::Integrator(SimInfo nAtoms = info->n_atoms; integrableObjects = info->integrableObjects; + consFramework = new RattleFramework(info); + + if(consFramework == NULL){ + sprintf(painCave.errMsg, + "Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); + painCave.isFatal = 1; + simError(); + } + +/* // check for constraints constrainedA = NULL; @@ -45,9 +56,13 @@ template Integrator::Integrator(SimInfo nConstrained = 0; checkConstraints(); +*/ } template Integrator::~Integrator(){ + if (consFramework != NULL) + delete consFramework; +/* if (nConstrained){ delete[] constrainedA; delete[] constrainedB; @@ -56,8 +71,10 @@ template Integrator::~Integrator(){ delete[] moved; delete[] oldPos; } +*/ } +/* template void Integrator::checkConstraints(void){ isConstrained = 0; @@ -92,7 +109,7 @@ template void Integrator::checkConstrai if (constrained){ dummy_plug = theArray[j]->get_constraint(); temp_con[nConstrained].set_a(dummy_plug->get_a()); - temp_con[nConstrained].set_b(dummy_plug->get_b()); + temp_con[nConstrained].set_b(Dummy_plug->get_b()); temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr()); nConstrained++; @@ -150,8 +167,8 @@ template void Integrator::checkConstrai delete[] temp_con; } +*/ - template void Integrator::integrate(void){ double runTime = info->run_time; @@ -160,7 +177,7 @@ template void Integrator::integrate(voi double thermalTime = info->thermalTime; double resetTime = info->resetTime; - + double difference; double currSample; double currThermal; double currStatus; @@ -178,21 +195,26 @@ template void Integrator::integrate(voi 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); - //temp test - tStats->getPotential(); + //execute constraint algorithm to make sure at the very beginning the system is constrained + //consFramework->doPreConstraint(); + //consFramework->doConstrainA(); + //calcForce(1, 1); + //consFramework->doConstrainB(); - if (nConstrained){ - preMove(); - constrainA(); - calcForce(1, 1); - constrainB(); - } - if (info->setTemp){ thermalize(); } @@ -213,8 +235,9 @@ template void Integrator::integrate(voi 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; } @@ -269,6 +292,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(); + + delete dumpOut; delete statOut; } @@ -281,7 +309,8 @@ template void Integrator::integrateStep startProfile(pro3); #endif //profile - preMove(); + //save old state (position, velocity etc) + consFramework->doPreConstraint(); #ifdef PROFILE endProfile(pro3); @@ -303,9 +332,7 @@ template void Integrator::integrateStep MPIcheckPoint(); #endif // is_mpi - // calc forces - calcForce(calcPot, calcStress); #ifdef IS_MPI @@ -340,6 +367,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); @@ -378,9 +406,7 @@ template void Integrator::moveA(void){ } } - if (nConstrained){ - constrainA(); - } + consFramework->doConstrainA(); } @@ -421,11 +447,10 @@ template void Integrator::moveB(void){ } } - if (nConstrained){ - constrainB(); - } + consFramework->doConstrainB(); } +/* template void Integrator::preMove(void){ int i, j; double pos[3]; @@ -684,12 +709,13 @@ template void Integrator::constrainB(vo simError(); } } - +*/ 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 @@ -697,26 +723,43 @@ template void Integrator::rotationPropa sd->getA(A); sd->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 ); + 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 ); - // rotate about the z-axis - angle = dt * ji[2] / I[2][2]; - this->rotate( 0, 1, angle, ji, A); + angle = dt * ji[k] / I[k][k]; + this->rotate( i, j, 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 ); - + 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 ); }