--- trunk/OOPSE/libmdtools/DirectionalAtom.cpp 2003/07/14 21:48:43 599 +++ trunk/OOPSE/libmdtools/DirectionalAtom.cpp 2003/08/07 21:47:18 670 @@ -1,16 +1,76 @@ #include #include "Atom.hpp" +#include "simError.h" +void DirectionalAtom::zeroForces() { + if( hasCoords ){ + frc[offsetX] = 0.0; + frc[offsetY] = 0.0; + frc[offsetZ] = 0.0; + + trq[offsetX] = 0.0; + trq[offsetY] = 0.0; + trq[offsetZ] = 0.0; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to zero frc and trq for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } +} +double DirectionalAtom::getMu( void ) { + if( hasCoords ){ + return mu[index]; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to get Mu for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } + return 0; +} + +void DirectionalAtom::setMu( double the_mu ) { + + if( hasCoords ){ + mu[index] = the_mu; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to set Mu for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } +} + void DirectionalAtom::setA( double the_A[3][3] ){ - - Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2]; - Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2]; - Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2]; - this->updateU(); + if( hasCoords ){ + Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2]; + Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2]; + Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2]; + + this->updateU(); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to set Amat for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } void DirectionalAtom::setI( double the_I[3][3] ){ @@ -24,48 +84,80 @@ void DirectionalAtom::setQ( double the_q[4] ){ double q0Sqr, q1Sqr, q2Sqr, q3Sqr; - q0Sqr = the_q[0] * the_q[0]; - q1Sqr = the_q[1] * the_q[1]; - q2Sqr = the_q[2] * the_q[2]; - q3Sqr = the_q[3] * the_q[3]; - + if( hasCoords ){ + q0Sqr = the_q[0] * the_q[0]; + q1Sqr = the_q[1] * the_q[1]; + q2Sqr = the_q[2] * the_q[2]; + q3Sqr = the_q[3] * the_q[3]; + + + Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr; + Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] ); + Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] ); + + Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] ); + Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr; + Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] ); + + Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] ); + Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] ); + Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr; + + this->updateU(); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to set Q for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } - Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr; - Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] ); - Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] ); - - Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] ); - Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr; - Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] ); - - Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] ); - Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] ); - Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr; - - this->updateU(); } void DirectionalAtom::getA( double the_A[3][3] ){ - the_A[0][0] = Amat[Axx]; - the_A[0][1] = Amat[Axy]; - the_A[0][2] = Amat[Axz]; + if( hasCoords ){ + the_A[0][0] = Amat[Axx]; + the_A[0][1] = Amat[Axy]; + the_A[0][2] = Amat[Axz]; + + the_A[1][0] = Amat[Ayx]; + the_A[1][1] = Amat[Ayy]; + the_A[1][2] = Amat[Ayz]; + + the_A[2][0] = Amat[Azx]; + the_A[2][1] = Amat[Azy]; + the_A[2][2] = Amat[Azz]; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to get Amat for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } - the_A[1][0] = Amat[Ayx]; - the_A[1][1] = Amat[Ayy]; - the_A[1][2] = Amat[Ayz]; - - the_A[2][0] = Amat[Azx]; - the_A[2][1] = Amat[Azy]; - the_A[2][2] = Amat[Azz]; } void DirectionalAtom::printAmatIndex( void ){ - std::cerr << "Atom[" << index << "] index =>\n" - << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n" - << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n" - << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n"; + if( hasCoords ){ + std::cerr << "Atom[" << index << "] index =>\n" + << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n" + << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n" + << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n"; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to print Amat indices for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } @@ -83,64 +175,85 @@ void DirectionalAtom::getQ( double q[4] ){ double t, s; double ad1, ad2, ad3; - t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0; - if( t > 0.0 ){ + if( hasCoords ){ - s = 0.5 / sqrt( t ); - q[0] = 0.25 / s; - q[1] = (Amat[Ayz] - Amat[Azy]) * s; - q[2] = (Amat[Azx] - Amat[Axz]) * s; - q[3] = (Amat[Axy] - Amat[Ayx]) * s; - } - else{ - - ad1 = fabs( Amat[Axx] ); - ad2 = fabs( Amat[Ayy] ); - ad3 = fabs( Amat[Azz] ); - - if( ad1 >= ad2 && ad1 >= ad3 ){ + t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0; + if( t > 0.0 ){ - s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] ); - q[0] = (Amat[Ayz] + Amat[Azy]) / s; - q[1] = 0.5 / s; - q[2] = (Amat[Axy] + Amat[Ayx]) / s; - q[3] = (Amat[Axz] + Amat[Azx]) / s; + s = 0.5 / sqrt( t ); + q[0] = 0.25 / s; + q[1] = (Amat[Ayz] - Amat[Azy]) * s; + q[2] = (Amat[Azx] - Amat[Axz]) * s; + q[3] = (Amat[Axy] - Amat[Ayx]) * s; } - else if( ad2 >= ad1 && ad2 >= ad3 ){ - - s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0; - q[0] = (Amat[Axz] + Amat[Azx]) / s; - q[1] = (Amat[Axy] + Amat[Ayx]) / s; - q[2] = 0.5 / s; - q[3] = (Amat[Ayz] + Amat[Azy]) / s; - } else{ - s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0; - q[0] = (Amat[Axy] + Amat[Ayx]) / s; - q[1] = (Amat[Axz] + Amat[Azx]) / s; - q[2] = (Amat[Ayz] + Amat[Azy]) / s; - q[3] = 0.5 / s; + ad1 = fabs( Amat[Axx] ); + ad2 = fabs( Amat[Ayy] ); + ad3 = fabs( Amat[Azz] ); + + if( ad1 >= ad2 && ad1 >= ad3 ){ + + s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] ); + q[0] = (Amat[Ayz] + Amat[Azy]) / s; + q[1] = 0.5 / s; + q[2] = (Amat[Axy] + Amat[Ayx]) / s; + q[3] = (Amat[Axz] + Amat[Azx]) / s; + } + else if( ad2 >= ad1 && ad2 >= ad3 ){ + + s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0; + q[0] = (Amat[Axz] + Amat[Azx]) / s; + q[1] = (Amat[Axy] + Amat[Ayx]) / s; + q[2] = 0.5 / s; + q[3] = (Amat[Ayz] + Amat[Azy]) / s; + } + else{ + + s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0; + q[0] = (Amat[Axy] + Amat[Ayx]) / s; + q[1] = (Amat[Axz] + Amat[Azx]) / s; + q[2] = (Amat[Ayz] + Amat[Azy]) / s; + q[3] = 0.5 / s; + } } } + else{ + + sprintf( painCave.errMsg, + "Attempt to get Q for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } void DirectionalAtom::setEuler( double phi, double theta, double psi ){ - Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); - Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); - Amat[Axz] = sin(theta) * sin(psi); - - Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); - Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); - Amat[Ayz] = sin(theta) * cos(psi); - - Amat[Azx] = sin(phi) * sin(theta); - Amat[Azy] = -cos(phi) * sin(theta); - Amat[Azz] = cos(theta); - - this->updateU(); + if( hasCoords ){ + Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); + Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); + Amat[Axz] = sin(theta) * sin(psi); + + Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); + Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); + Amat[Ayz] = sin(theta) * cos(psi); + + Amat[Azx] = sin(phi) * sin(theta); + Amat[Azy] = -cos(phi) * sin(theta); + Amat[Azz] = cos(theta); + + this->updateU(); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to set Euler angles for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } @@ -148,33 +261,64 @@ void DirectionalAtom::lab2Body( double r[3] ){ double rl[3]; // the lab frame vector - rl[0] = r[0]; - rl[1] = r[1]; - rl[2] = r[2]; + if( hasCoords ){ + rl[0] = r[0]; + rl[1] = r[1]; + rl[2] = r[2]; + + r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]); + r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]); + r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to convert lab2body for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } - r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]); - r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]); - r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]); } void DirectionalAtom::body2Lab( double r[3] ){ double rb[3]; // the body frame vector - rb[0] = r[0]; - rb[1] = r[1]; - rb[2] = r[2]; - - r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]); - r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]); - r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]); + if( hasCoords ){ + rb[0] = r[0]; + rb[1] = r[1]; + rb[2] = r[2]; + + r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]); + r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]); + r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to convert body2lab for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } void DirectionalAtom::updateU( void ){ - ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz); - ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz); - ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz); + if( hasCoords ){ + ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz); + ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz); + ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz); + } + else{ + + sprintf( painCave.errMsg, + "Attempt to updateU for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } void DirectionalAtom::getJ( double theJ[3] ){ @@ -193,16 +337,36 @@ void DirectionalAtom::getTrq( double theT[3] ){ void DirectionalAtom::getTrq( double theT[3] ){ - theT[0] = trq[offsetX]; - theT[1] = trq[offsetY]; - theT[2] = trq[offsetZ]; + if( hasCoords ){ + theT[0] = trq[offsetX]; + theT[1] = trq[offsetY]; + theT[2] = trq[offsetZ]; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to get Trq for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } } void DirectionalAtom::addTrq( double theT[3] ){ - trq[offsetX] += theT[0]; - trq[offsetY] += theT[1]; - trq[offsetZ] += theT[2]; + if( hasCoords ){ + trq[offsetX] += theT[0]; + trq[offsetY] += theT[1]; + trq[offsetZ] += theT[2]; + } + else{ + + sprintf( painCave.errMsg, + "Attempt to add Trq for atom %d before coords set.\n", + index ); + painCave.isFatal = 1; + simError(); + } }