--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/07/11 22:34:48 594 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/05/01 18:52:38 1144 @@ -1,18 +1,22 @@ #include -#include -#include +#include +#include #ifdef IS_MPI #include "mpiSimulation.hpp" #include #endif //is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif // profile + #include "Integrator.hpp" #include "simError.h" -Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){ - +template Integrator::Integrator(SimInfo* theInfo, + ForceFields* the_ff){ info = theInfo; myFF = the_ff; isFirst = 1; @@ -21,31 +25,30 @@ Integrator::Integrator( SimInfo *theInfo, ForceFields* nMols = info->n_mol; // give a little love back to the SimInfo object - - if( info->the_integrator != NULL ) delete info->the_integrator; - info->the_integrator = this; + if (info->the_integrator != NULL){ + delete info->the_integrator; + } + nAtoms = info->n_atoms; + integrableObjects = info->integrableObjects; - std::cerr << "integ nAtoms = " << nAtoms << "\n"; - // check for constraints - - constrainedA = NULL; - constrainedB = NULL; + + constrainedA = NULL; + constrainedB = NULL; constrainedDsqr = NULL; - moving = NULL; - moved = NULL; - oldPos = NULL; - + moving = NULL; + moved = NULL; + oldPos = NULL; + nConstrained = 0; checkConstraints(); } -Integrator::~Integrator() { - - if( nConstrained ){ +template Integrator::~Integrator(){ + if (nConstrained){ delete[] constrainedA; delete[] constrainedB; delete[] constrainedDsqr; @@ -53,406 +56,399 @@ Integrator::~Integrator() { delete[] moved; delete[] oldPos; } - } -void Integrator::checkConstraints( void ){ - - +template void Integrator::checkConstraints(void){ isConstrained = 0; - Constraint *temp_con; - Constraint *dummy_plug; + Constraint* temp_con; + Constraint* dummy_plug; temp_con = new Constraint[info->n_SRI]; nConstrained = 0; int constrained = 0; - + SRI** theArray; - for(int i = 0; i < nMols; i++){ - - theArray = (SRI**) molecules[i].getMyBonds(); - for(int j=0; jis_constrained(); - std::cerr << "Is the folowing bond constrained \n"; - theArray[j]->printMe(); - - if(constrained){ - - std::cerr << "Yes\n"; + 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_dsqr(dummy_plug->get_dsqr()); - 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_dsqr( dummy_plug->get_dsqr() ); - - nConstrained++; - constrained = 0; - } - else std::cerr << "No.\n"; + nConstrained++; + constrained = 0; + } } - theArray = (SRI**) molecules[i].getMyBends(); - for(int j=0; jis_constrained(); - - 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_dsqr( dummy_plug->get_dsqr() ); - - nConstrained++; - constrained = 0; + + 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_dsqr(dummy_plug->get_dsqr()); + + nConstrained++; + constrained = 0; } } - theArray = (SRI**) molecules[i].getMyTorsions(); - for(int j=0; jis_constrained(); - - 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_dsqr( dummy_plug->get_dsqr() ); - - nConstrained++; - constrained = 0; + + 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_dsqr(dummy_plug->get_dsqr()); + + nConstrained++; + constrained = 0; } } } - if(nConstrained > 0){ - + + if (nConstrained > 0){ isConstrained = 1; - if(constrainedA != NULL ) delete[] constrainedA; - if(constrainedB != NULL ) delete[] constrainedB; - if(constrainedDsqr != NULL ) delete[] constrainedDsqr; + if (constrainedA != NULL) + delete[] constrainedA; + if (constrainedB != NULL) + delete[] constrainedB; + if (constrainedDsqr != NULL) + delete[] constrainedDsqr; - constrainedA = new int[nConstrained]; - constrainedB = new int[nConstrained]; + constrainedA = new int[nConstrained]; + constrainedB = new int[nConstrained]; constrainedDsqr = new double[nConstrained]; - - for( int i = 0; i < nConstrained; i++){ - + + for (int i = 0; i < nConstrained; i++){ constrainedA[i] = temp_con[i].get_a(); constrainedB[i] = temp_con[i].get_b(); constrainedDsqr[i] = temp_con[i].get_dsqr(); - } - - // save oldAtoms to check for lode balanceing later on. - + + // save oldAtoms to check for lode balancing later on. + oldAtoms = nAtoms; - + moving = new int[nAtoms]; - moved = new int[nAtoms]; + moved = new int[nAtoms]; - oldPos = new double[nAtoms*3]; + oldPos = new double[nAtoms * 3]; } - + delete[] temp_con; } -void Integrator::integrate( void ){ +template void Integrator::integrate(void){ - int i, j; // loop counters - - double runTime = info->run_time; - double sampleTime = info->sampleTime; - double statusTime = info->statusTime; + double runTime = info->run_time; + double sampleTime = info->sampleTime; + double statusTime = info->statusTime; double thermalTime = info->thermalTime; + double resetTime = info->resetTime; + double currSample; double currThermal; double currStatus; - double currTime; + double currReset; int calcPot, calcStress; - int isError; + tStats = new Thermo(info); + statOut = new StatWriter(info); + dumpOut = new DumpWriter(info); - - 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(); + + // remove center of mass drift velocity (in case we passed in a configuration + // that was drifting + tStats->removeCOMdrift(); + // initialize the forces before the first step - myFF->doForces(1,1); + calcForce(1, 1); - if( info->setTemp ){ - - tStats->velocitize(); + if (nConstrained){ + preMove(); + constrainA(); + calcForce(1, 1); + constrainB(); } - dumpOut->writeDump( 0.0 ); - statOut->writeStat( 0.0 ); - + if (info->setTemp){ + thermalize(); + } + calcPot = 0; calcStress = 0; - currSample = sampleTime; - currThermal = thermalTime; - currStatus = statusTime; - currTime = 0.0;; + currSample = sampleTime + info->getTime(); + currThermal = thermalTime+ info->getTime(); + currStatus = statusTime + info->getTime(); + currReset = resetTime + info->getTime(); + dumpOut->writeDump(info->getTime()); + statOut->writeStat(info->getTime()); - readyCheck(); #ifdef IS_MPI - strcpy( checkPointMsg, - "The integrator is ready to go." ); + strcpy(checkPointMsg, "The integrator is ready to go."); MPIcheckPoint(); #endif // is_mpi - - pos = Atom::getPosArray(); - vel = Atom::getVelArray(); - frc = Atom::getFrcArray(); - trq = Atom::getTrqArray(); - Amat = Atom::getAmatArray(); - - while( currTime < runTime ){ - - if( (currTime+dt) >= currStatus ){ + while (info->getTime() < runTime && !stopIntegrator()){ + if ((info->getTime() + dt) >= currStatus){ calcPot = 1; calcStress = 1; } - std::cerr << "calcPot = " << calcPot << "; calcStress = " - << calcStress << "\n"; +#ifdef PROFILE + startProfile( pro1 ); +#endif + + integrateStep(calcPot, calcStress); - integrateStep( calcPot, calcStress ); - - currTime += dt; +#ifdef PROFILE + endProfile( pro1 ); - if( info->setTemp ){ - if( currTime >= currThermal ){ - tStats->velocitize(); - currThermal += thermalTime; + startProfile( pro2 ); +#endif // profile + + info->incrTime(dt); + + if (info->setTemp){ + if (info->getTime() >= currThermal){ + thermalize(); + currThermal += thermalTime; } } - if( currTime >= currSample ){ - dumpOut->writeDump( currTime ); + if (info->getTime() >= currSample){ + dumpOut->writeDump(info->getTime()); currSample += sampleTime; } - if( currTime >= currStatus ){ - statOut->writeStat( currTime ); - calcPot = 0; + if (info->getTime() >= currStatus){ + statOut->writeStat(info->getTime()); + calcPot = 0; calcStress = 0; currStatus += statusTime; - } + } + if (info->resetIntegrator){ + if (info->getTime() >= currReset){ + this->resetIntegrator(); + currReset += resetTime; + } + } + +#ifdef PROFILE + endProfile( pro2 ); +#endif //profile + #ifdef IS_MPI - strcpy( checkPointMsg, - "successfully took a time step." ); + strcpy(checkPointMsg, "successfully took a time step."); MPIcheckPoint(); #endif // is_mpi - } - dumpOut->writeFinal(currTime); - delete dumpOut; delete statOut; } -void Integrator::integrateStep( int calcPot, int calcStress ){ - - - +template void Integrator::integrateStep(int calcPot, + int calcStress){ // Position full step, and velocity half step +#ifdef PROFILE + startProfile(pro3); +#endif //profile + preMove(); + +#ifdef PROFILE + endProfile(pro3); + + startProfile(pro4); +#endif // profile + moveA(); - if( nConstrained ) constrainA(); +#ifdef PROFILE + endProfile(pro4); + + startProfile(pro5); +#endif//profile + + +#ifdef IS_MPI + strcpy(checkPointMsg, "Succesful moveA\n"); + MPIcheckPoint(); +#endif // is_mpi + + // calc forces - myFF->doForces(calcPot,calcStress); + calcForce(calcPot, calcStress); +#ifdef IS_MPI + strcpy(checkPointMsg, "Succesful doForces\n"); + MPIcheckPoint(); +#endif // is_mpi + +#ifdef PROFILE + endProfile( pro5 ); + + startProfile( pro6 ); +#endif //profile + // finish the velocity half step - + moveB(); - if( nConstrained ) constrainB(); - + +#ifdef PROFILE + endProfile(pro6); +#endif // profile + +#ifdef IS_MPI + strcpy(checkPointMsg, "Succesful moveB\n"); + MPIcheckPoint(); +#endif // is_mpi } -void Integrator::moveA( void ){ - - int i,j,k; - int atomIndex, aMatIndex; +template void Integrator::moveA(void){ + size_t i, j; DirectionalAtom* dAtom; - double Tb[3]; - double ji[3]; - double angle; - double A[3][3]; + double Tb[3], ji[3]; + double vel[3], pos[3], frc[3]; + double mass; + + for (i = 0; i < integrableObjects.size() ; i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getPos(pos); + integrableObjects[i]->getFrc(frc); + std::cerr << "i =\t" << i << "\t" << frc[0] << "\t" << frc[1]<< "\t" << frc[2] << "\n"; + + mass = integrableObjects[i]->getMass(); - for( i=0; igetMass() ) * eConvert; + integrableObjects[i]->setVel(vel); + integrableObjects[i]->setPos(pos); - std::cerr<< "MoveA vel[" << i << "] = " - << vel[atomIndex] << "\t" - << vel[atomIndex+1]<< "\t" - << vel[atomIndex+2]<< "\n"; + if (integrableObjects[i]->isDirectional()){ - // position whole step - for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j]; - + // get and convert the torque to body frame - std::cerr<< "MoveA pos[" << i << "] = " - << pos[atomIndex] << "\t" - << pos[atomIndex+1]<< "\t" - << pos[atomIndex+2]<< "\n"; + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); - if( atoms[i]->isDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[i]; - - // get and convert the torque to body frame - - Tb[0] = dAtom->getTx(); - Tb[1] = dAtom->getTy(); - Tb[2] = dAtom->getTz(); - - dAtom->lab2Body( Tb ); - // get the angular momentum, and propagate a half step - - ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert; - ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert; - ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert; - - // use the angular velocities to propagate the rotation matrix a - // full time step - - // get the atom's rotation matrix - - A[0][0] = dAtom->getAxx(); - A[0][1] = dAtom->getAxy(); - A[0][2] = dAtom->getAxz(); - - A[1][0] = dAtom->getAyx(); - A[1][1] = dAtom->getAyy(); - A[1][2] = dAtom->getAyz(); - - A[2][0] = dAtom->getAzx(); - A[2][1] = dAtom->getAzy(); - A[2][2] = dAtom->getAzz(); - - // rotate about the x-axis - angle = dt2 * ji[0] / dAtom->getIxx(); - this->rotate( 1, 2, angle, ji, A ); - - // rotate about the y-axis - angle = dt2 * ji[1] / dAtom->getIyy(); - this->rotate( 2, 0, angle, ji, A ); - - // rotate about the z-axis - angle = dt * ji[2] / dAtom->getIzz(); - this->rotate( 0, 1, angle, ji, A ); - - // rotate about the y-axis - angle = dt2 * ji[1] / dAtom->getIyy(); - this->rotate( 2, 0, angle, ji, A ); - - // rotate about the x-axis - angle = dt2 * ji[0] / dAtom->getIxx(); - this->rotate( 1, 2, angle, ji, A ); - - dAtom->setJx( ji[0] ); - dAtom->setJy( ji[1] ); - dAtom->setJz( ji[2] ); + + integrableObjects[i]->getJ(ji); + + for (j = 0; j < 3; j++) + ji[j] += (dt2 * Tb[j]) * eConvert; + + this->rotationPropagation( integrableObjects[i], ji ); + + integrableObjects[i]->setJ(ji); } - } + + if (nConstrained){ + constrainA(); + } } -void Integrator::moveB( void ){ - int i,j,k; - int atomIndex; - DirectionalAtom* dAtom; - double Tb[3]; - double ji[3]; +template void Integrator::moveB(void){ + int i, j; + double Tb[3], ji[3]; + double vel[3], frc[3]; + double mass; - for( i=0; igetVel(vel); + integrableObjects[i]->getFrc(frc); + mass = integrableObjects[i]->getMass(); + // velocity half step - for( j=atomIndex; j<(atomIndex+3); j++ ) - vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert; + for (j = 0; j < 3; j++) + vel[j] += (dt2 * frc[j] / mass) * eConvert; - std::cerr<< "MoveB vel[" << i << "] = " - << vel[atomIndex] << "\t" - << vel[atomIndex+1]<< "\t" - << vel[atomIndex+2]<< "\n"; + integrableObjects[i]->setVel(vel); + if (integrableObjects[i]->isDirectional()){ - if( atoms[i]->isDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[i]; - // get and convert the torque to body frame - - Tb[0] = dAtom->getTx(); - Tb[1] = dAtom->getTy(); - Tb[2] = dAtom->getTz(); - - dAtom->lab2Body( Tb ); - - // get the angular momentum, and complete the angular momentum - // half step - - ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert; - ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert; - ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert; - - dAtom->setJx( ji[0] ); - dAtom->setJy( ji[1] ); - dAtom->setJz( ji[2] ); + + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); + + // get the angular momentum, and propagate a half step + + integrableObjects[i]->getJ(ji); + + for (j = 0; j < 3; j++) + ji[j] += (dt2 * Tb[j]) * eConvert; + + + integrableObjects[i]->setJ(ji); } } + if (nConstrained){ + constrainB(); + } } -void Integrator::preMove( void ){ - int i; +template void Integrator::preMove(void){ + int i, j; + double pos[3]; - if( nConstrained ){ + if (nConstrained){ + for (i = 0; i < nAtoms; i++){ + atoms[i]->getPos(pos); - for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i]; + for (j = 0; j < 3; j++){ + oldPos[3 * i + j] = pos[j]; + } + } } -} +} -void Integrator::constrainA(){ - - int i,j,k; +template void Integrator::constrainA(){ + int i, j; int done; + double posA[3], posB[3]; + double velA[3], velB[3]; double pab[3]; double rab[3]; int a, b, ax, ay, az, bx, by, bz; @@ -464,109 +460,117 @@ void Integrator::constrainA(){ double gab; int iteration; - for( i=0; igetPos(posA); + atoms[b]->getPos(posB); - info->wrapVector( pab ); + for (j = 0; j < 3; j++) + pab[j] = posA[j] - posB[j]; - pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; + //periodic boundary condition - rabsq = constrainedDsqr[i]; - diffsq = rabsq - pabsq; + info->wrapVector(pab); - // the original rattle code from alan tidesley - if (fabs(diffsq) > (tol*rabsq*2)) { - rab[0] = oldPos[ax] - oldPos[bx]; - rab[1] = oldPos[ay] - oldPos[by]; - rab[2] = oldPos[az] - oldPos[bz]; + pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; - info->wrapVector( rab ); + rabsq = constrainedDsqr[i]; + diffsq = rabsq - pabsq; - rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; + // the original rattle code from alan tidesley + if (fabs(diffsq) > (tol * rabsq * 2)){ + rab[0] = oldPos[ax] - oldPos[bx]; + rab[1] = oldPos[ay] - oldPos[by]; + rab[2] = oldPos[az] - oldPos[bz]; - rpabsq = rpab * rpab; + info->wrapVector(rab); + rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; - if (rpabsq < (rabsq * -diffsq)){ + rpabsq = rpab * rpab; + + if (rpabsq < (rabsq * -diffsq)){ #ifdef IS_MPI - a = atoms[a]->getGlobalIndex(); - b = atoms[b]->getGlobalIndex(); + a = atoms[a]->getGlobalIndex(); + b = atoms[b]->getGlobalIndex(); #endif //is_mpi - sprintf( painCave.errMsg, - "Constraint failure in constrainA at atom %d and %d.\n", - a, b ); - painCave.isFatal = 1; - simError(); - } + sprintf(painCave.errMsg, + "Constraint failure in constrainA at atom %d and %d.\n", a, + b); + painCave.isFatal = 1; + simError(); + } - rma = 1.0 / atoms[a]->getMass(); - rmb = 1.0 / atoms[b]->getMass(); + rma = 1.0 / atoms[a]->getMass(); + rmb = 1.0 / atoms[b]->getMass(); - gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab ); + gab = diffsq / (2.0 * (rma + rmb) * rpab); dx = rab[0] * gab; dy = rab[1] * gab; dz = rab[2] * gab; - pos[ax] += rma * dx; - pos[ay] += rma * dy; - pos[az] += rma * dz; + posA[0] += rma * dx; + posA[1] += rma * dy; + posA[2] += rma * dz; - pos[bx] -= rmb * dx; - pos[by] -= rmb * dy; - pos[bz] -= rmb * dz; + atoms[a]->setPos(posA); + posB[0] -= rmb * dx; + posB[1] -= rmb * dy; + posB[2] -= rmb * dz; + + atoms[b]->setPos(posB); + dx = dx / dt; dy = dy / dt; dz = dz / dt; - vel[ax] += rma * dx; - vel[ay] += rma * dy; - vel[az] += rma * dz; + atoms[a]->getVel(velA); - vel[bx] -= rmb * dx; - vel[by] -= rmb * dy; - vel[bz] -= rmb * dz; + velA[0] += rma * dx; + velA[1] += rma * dy; + velA[2] += rma * dz; - moving[a] = 1; - moving[b] = 1; - done = 0; - } + atoms[a]->setVel(velA); + + atoms[b]->getVel(velB); + + velB[0] -= rmb * dx; + velB[1] -= rmb * dy; + velB[2] -= rmb * dz; + + atoms[b]->setVel(velB); + + moving[a] = 1; + moving[b] = 1; + done = 0; + } } } - - for(i=0; i void Integrator::constrainB(void){ + int i, j; int done; + double posA[3], posB[3]; + double velA[3], velB[3]; double vxab, vyab, vzab; double rab[3]; 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; - for(i=0; igetVel(velA); + atoms[b]->getVel(velB); - rab[0] = pos[ax] - pos[bx]; - rab[1] = pos[ay] - pos[by]; - rab[2] = pos[az] - pos[bz]; - - info->wrapVector( rab ); - - rma = 1.0 / atoms[a]->getMass(); - rmb = 1.0 / atoms[b]->getMass(); + vxab = velA[0] - velB[0]; + vyab = velA[1] - velB[1]; + vzab = velA[2] - velB[2]; - rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab; - - gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] ); + atoms[a]->getPos(posA); + atoms[b]->getPos(posB); - if (fabs(gab) > tol) { - - dx = rab[0] * gab; - dy = rab[1] * gab; - dz = rab[2] * gab; - - vel[ax] += rma * dx; - vel[ay] += rma * dy; - vel[az] += rma * dz; + for (j = 0; j < 3; j++) + rab[j] = posA[j] - posB[j]; - vel[bx] -= rmb * dx; - vel[by] -= rmb * dy; - vel[bz] -= rmb * dz; - - moving[a] = 1; - moving[b] = 1; - done = 0; - } + info->wrapVector(rab); + + rma = 1.0 / atoms[a]->getMass(); + rmb = 1.0 / atoms[b]->getMass(); + + rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab; + + gab = -rvab / ((rma + rmb) * constrainedDsqr[i]); + + if (fabs(gab) > tol){ + dx = rab[0] * gab; + dy = rab[1] * gab; + dz = rab[2] * gab; + + velA[0] += rma * dx; + velA[1] += rma * dy; + velA[2] += rma * dz; + + atoms[a]->setVel(velA); + + velB[0] -= rmb * dx; + velB[1] -= rmb * dy; + velB[2] -= rmb * dz; + + atoms[b]->setVel(velB); + + moving[a] = 1; + moving[b] = 1; + done = 0; + } } } - for(i=0; i 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 + sd->getA(A); + sd->getI(I); + 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 ); + angle = dt * ji[k] / I[k][k]; + this->rotate( i, j, angle, ji, A); -void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], - double A[3][3] ){ + angle = dt2 * ji[j] / I[j][j]; + this->rotate( k, i, angle, ji, A ); - int i,j,k; + } 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]; + 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 ); +} + +template void Integrator::rotate(int axes1, int axes2, + double angle, double ji[3], + double A[3][3]){ + int i, j, k; double sinAngle; double cosAngle; double angleSqr; @@ -704,16 +755,17 @@ void Integrator::rotate( int axes1, int axes2, double // initialize the tempA - for(i=0; i<3; i++){ - for(j=0; j<3; j++){ + for (i = 0; i < 3; i++){ + for (j = 0; j < 3; j++){ tempA[j][i] = A[i][j]; } } // initialize the tempJ - for( i=0; i<3; i++) tempJ[i] = ji[i]; - + for (i = 0; i < 3; i++) + tempJ[i] = ji[i]; + // initalize rot as a unit matrix rot[0][0] = 1.0; @@ -723,14 +775,14 @@ void Integrator::rotate( int axes1, int axes2, double rot[1][0] = 0.0; rot[1][1] = 1.0; rot[1][2] = 0.0; - + rot[2][0] = 0.0; rot[2][1] = 0.0; rot[2][2] = 1.0; - + // use a small angle aproximation for sin and cosine - angleSqr = angle * angle; + angleSqr = angle * angle; angleSqrOver4 = angleSqr / 4.0; top = 1.0 - angleSqrOver4; bottom = 1.0 + angleSqrOver4; @@ -743,17 +795,17 @@ void Integrator::rotate( int axes1, int axes2, double rot[axes1][axes2] = sinAngle; rot[axes2][axes1] = -sinAngle; - + // rotate the momentum acoording to: ji[] = rot[][] * ji[] - - for(i=0; i<3; i++){ + + for (i = 0; i < 3; i++){ ji[i] = 0.0; - for(k=0; k<3; k++){ + for (k = 0; k < 3; k++){ ji[i] += rot[i][k] * tempJ[k]; } } - // rotate the Rotation matrix acording to: + // rotate the Rotation matrix acording to: // A[][] = A[][] * transpose(rot[][]) @@ -761,12 +813,30 @@ void Integrator::rotate( int axes1, int axes2, double // calculation as: // transpose(A[][]) = transpose(A[][]) * transpose(rot[][]) - for(i=0; i<3; i++){ - for(j=0; j<3; j++){ + for (i = 0; i < 3; i++){ + for (j = 0; j < 3; j++){ A[j][i] = 0.0; - for(k=0; k<3; k++){ - A[j][i] += tempA[i][k] * rot[j][k]; + for (k = 0; k < 3; k++){ + A[j][i] += tempA[i][k] * rot[j][k]; } } } } + +template void Integrator::calcForce(int calcPot, int calcStress){ + myFF->doForces(calcPot, calcStress); +} + +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(); +}