--- trunk/OOPSE/libmdtools/Integrator.cpp 2004/02/17 19:23:44 1057 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2004/04/12 20:32:20 1097 @@ -31,6 +31,7 @@ template Integrator::Integrator(SimInfo } nAtoms = info->n_atoms; + integrableObjects = info->integrableObjects; // check for constraints @@ -334,18 +335,18 @@ template void Integrator::moveA(void){ template void Integrator::moveA(void){ - int i, j; + size_t i, j; DirectionalAtom* dAtom; double Tb[3], ji[3]; double vel[3], pos[3], frc[3]; double mass; - - for (i = 0; i < nAtoms; i++){ - atoms[i]->getVel(vel); - atoms[i]->getPos(pos); - atoms[i]->getFrc(frc); - - mass = atoms[i]->getMass(); + + for (i = 0; i < integrableObjects.size() ; i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getPos(pos); + integrableObjects[i]->getFrc(frc); + + mass = integrableObjects[i]->getMass(); for (j = 0; j < 3; j++){ // velocity half step @@ -354,27 +355,26 @@ template void Integrator::moveA(void){ pos[j] += dt * vel[j]; } - atoms[i]->setVel(vel); - atoms[i]->setPos(pos); + integrableObjects[i]->setVel(vel); + integrableObjects[i]->setPos(pos); - if (atoms[i]->isDirectional()){ - dAtom = (DirectionalAtom *) atoms[i]; + if (integrableObjects[i]->isDirectional()){ // get and convert the torque to body frame - dAtom->getTrq(Tb); - dAtom->lab2Body(Tb); + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); // get the angular momentum, and propagate a half step - dAtom->getJ(ji); + integrableObjects[i]->getJ(ji); for (j = 0; j < 3; j++) ji[j] += (dt2 * Tb[j]) * eConvert; - this->rotationPropagation( dAtom, ji ); + this->rotationPropagation( integrableObjects[i], ji ); - dAtom->setJ(ji); + integrableObjects[i]->setJ(ji); } } @@ -386,40 +386,38 @@ template void Integrator::moveB(void){ template void Integrator::moveB(void){ int i, j; - DirectionalAtom* dAtom; double Tb[3], ji[3]; double vel[3], frc[3]; double mass; - for (i = 0; i < nAtoms; i++){ - atoms[i]->getVel(vel); - atoms[i]->getFrc(frc); + for (i = 0; i < integrableObjects.size(); i++){ + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getFrc(frc); - mass = atoms[i]->getMass(); + mass = integrableObjects[i]->getMass(); // velocity half step for (j = 0; j < 3; j++) vel[j] += (dt2 * frc[j] / mass) * eConvert; - atoms[i]->setVel(vel); + integrableObjects[i]->setVel(vel); - if (atoms[i]->isDirectional()){ - dAtom = (DirectionalAtom *) atoms[i]; + if (integrableObjects[i]->isDirectional()){ // get and convert the torque to body frame - dAtom->getTrq(Tb); - dAtom->lab2Body(Tb); + integrableObjects[i]->getTrq(Tb); + integrableObjects[i]->lab2Body(Tb); // get the angular momentum, and propagate a half step - dAtom->getJ(ji); + integrableObjects[i]->getJ(ji); for (j = 0; j < 3; j++) ji[j] += (dt2 * Tb[j]) * eConvert; - dAtom->setJ(ji); + integrableObjects[i]->setJ(ji); } } @@ -688,7 +686,7 @@ template void Integrator::rotationPropa } template void Integrator::rotationPropagation -( DirectionalAtom* dAtom, double ji[3] ){ +( StuntDouble* sd, double ji[3] ){ double angle; double A[3][3], I[3][3]; @@ -696,8 +694,8 @@ template void Integrator::rotationPropa // use the angular velocities to propagate the rotation matrix a // full time step - dAtom->getA(A); - dAtom->getI(I); + sd->getA(A); + sd->getI(I); // rotate about the x-axis angle = dt2 * ji[0] / I[0][0]; @@ -719,7 +717,7 @@ template void Integrator::rotationPropa angle = dt2 * ji[0] / I[0][0]; this->rotate( 1, 2, angle, ji, A ); - dAtom->setA( A ); + sd->setA( A ); } template void Integrator::rotate(int axes1, int axes2,