--- trunk/OOPSE/libmdtools/Integrator.cpp 2004/06/01 17:15:43 1212 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/06/04 03:15:31 1234 @@ -1,7 +1,7 @@ #include #include #include - +#include "Rattle.hpp" #ifdef IS_MPI #include "mpiSimulation.hpp" #include @@ -32,7 +32,17 @@ template Integrator::Integrator(SimInfo 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; @@ -45,10 +55,13 @@ template Integrator::Integrator(SimInfo nConstrained = 0; checkConstraints(); - +*/ } template Integrator::~Integrator(){ + if (rattle != NULL) + delete rattle; +/* if (nConstrained){ delete[] constrainedA; delete[] constrainedB; @@ -57,8 +70,10 @@ template Integrator::~Integrator(){ delete[] moved; delete[] oldPos; } +*/ } +/* template void Integrator::checkConstraints(void){ isConstrained = 0; @@ -151,8 +166,8 @@ template void Integrator::checkConstrai delete[] temp_con; } +*/ - template void Integrator::integrate(void){ double runTime = info->run_time; @@ -192,13 +207,12 @@ template void Integrator::integrate(voi // 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(); @@ -255,8 +269,6 @@ template void Integrator::integrate(voi if (info->getTime() >= currStatus){ statOut->writeStat(info->getTime()); - if (info->useSolidThermInt || info->useLiquidThermInt) - statOut->writeRaw(info->getTime()); calcPot = 0; calcStress = 0; currStatus += statusTime; @@ -296,7 +308,8 @@ template void Integrator::integrateStep startProfile(pro3); #endif //profile - preMove(); + //save old state (position, velocity etc) + rattle->doPreConstraint(); #ifdef PROFILE endProfile(pro3); @@ -392,9 +405,7 @@ template void Integrator::moveA(void){ } } - if (nConstrained){ - constrainA(); - } + rattle->doRattleA(); } @@ -435,11 +446,10 @@ template void Integrator::moveB(void){ } } - if (nConstrained){ - constrainB(); - } + rattle->doRattleB(); } +/* template void Integrator::preMove(void){ int i, j; double pos[3]; @@ -698,7 +708,7 @@ template void Integrator::constrainB(vo simError(); } } - +*/ template void Integrator::rotationPropagation ( StuntDouble* sd, double ji[3] ){