--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/09/04 21:48:35 746 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/10/29 20:41:39 843 @@ -1,6 +1,6 @@ #include -#include -#include +#include +#include #ifdef IS_MPI #include "mpiSimulation.hpp" @@ -25,7 +25,6 @@ template Integrator::Integrator(SimInfo if (info->the_integrator != NULL){ delete info->the_integrator; } - info->the_integrator = this; nAtoms = info->n_atoms; @@ -147,7 +146,6 @@ template void Integrator::integrate(voi template void Integrator::integrate(void){ - int i, j; // loop counters double runTime = info->run_time; double sampleTime = info->sampleTime; @@ -160,23 +158,30 @@ template void Integrator::integrate(voi double currThermal; double currStatus; double currReset; - + int calcPot, calcStress; - int isError; tStats = new Thermo(info); statOut = new StatWriter(info); dumpOut = new DumpWriter(info); atoms = info->atoms; - DirectionalAtom* dAtom; dt = info->dt; dt2 = 0.5 * dt; + readyCheck(); + // initialize the forces before the first step calcForce(1, 1); + + if (nConstrained){ + preMove(); + constrainA(); + calcForce(1, 1); + constrainB(); + } if (info->setTemp){ thermalize(); @@ -192,7 +197,6 @@ template void Integrator::integrate(voi dumpOut->writeDump(info->getTime()); statOut->writeStat(info->getTime()); - readyCheck(); #ifdef IS_MPI strcpy(checkPointMsg, "The integrator is ready to go."); @@ -222,11 +226,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){ @@ -241,8 +245,10 @@ template void Integrator::integrate(voi #endif // is_mpi } - dumpOut->writeFinal(info->getTime()); + // write the last frame + dumpOut->writeDump(info->getTime()); + delete dumpOut; delete statOut; } @@ -254,9 +260,7 @@ template void Integrator::integrateStep moveA(); - if (nConstrained){ - constrainA(); - } + #ifdef IS_MPI @@ -279,10 +283,8 @@ template void Integrator::integrateStep moveB(); - if (nConstrained){ - constrainB(); - } + #ifdef IS_MPI strcpy(checkPointMsg, "Succesful moveB\n"); MPIcheckPoint(); @@ -294,8 +296,6 @@ template void Integrator::moveA(void){ int i, j; DirectionalAtom* dAtom; double Tb[3], ji[3]; - double A[3][3], I[3][3]; - double angle; double vel[3], pos[3], frc[3]; double mass; @@ -331,37 +331,15 @@ template void Integrator::moveA(void){ for (j = 0; j < 3; j++) ji[j] += (dt2 * Tb[j]) * eConvert; - // use the angular velocities to propagate the rotation matrix a - // full time step + this->rotationPropagation( dAtom, ji ); - dAtom->getA(A); - dAtom->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); - - // 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->setJ(ji); - dAtom->setA(A); } } + + if (nConstrained){ + constrainA(); + } } @@ -387,7 +365,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); @@ -403,6 +381,10 @@ template void Integrator::moveB(void){ dAtom->setJ(ji); } } + + if (nConstrained){ + constrainB(); + } } template void Integrator::preMove(void){ @@ -421,7 +403,7 @@ template void Integrator::constrainA(){ } template void Integrator::constrainA(){ - int i, j, k; + int i, j; int done; double posA[3], posB[3]; double velA[3], velB[3]; @@ -561,10 +543,11 @@ template void Integrator::constrainA(){ painCave.isFatal = 1; simError(); } + } template void Integrator::constrainB(void){ - int i, j, k; + int i, j; int done; double posA[3], posB[3]; double velA[3], velB[3]; @@ -573,8 +556,7 @@ template void Integrator::constrainB(vo int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; - double rabsq, pabsq, rvab; - double diffsq; + double rvab; double gab; int iteration; @@ -664,6 +646,41 @@ template void Integrator::rotate(int ax } } +template void Integrator::rotationPropagation +( DirectionalAtom* dAtom, double ji[3] ){ + + double angle; + double A[3][3], I[3][3]; + + // use the angular velocities to propagate the rotation matrix a + // full time step + + dAtom->getA(A); + dAtom->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 ); + + // 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 ); +} + template void Integrator::rotate(int axes1, int axes2, double angle, double ji[3], double A[3][3]){ @@ -729,7 +746,7 @@ template void Integrator::rotate(int ax } } - // rotate the Rotation matrix acording to: + // rotate the Rotation matrix acording to: // A[][] = A[][] * transpose(rot[][]) @@ -754,3 +771,13 @@ template void Integrator::thermalize(){ template void Integrator::thermalize(){ tStats->velocitize(); } + +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(); +}