--- trunk/src/primitives/RigidBody.cpp 2005/10/03 15:54:23 642 +++ branches/development/src/primitives/RigidBody.cpp 2013/01/30 14:43:08 1844 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,118 +28,130 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include #include #include "primitives/RigidBody.hpp" #include "utils/simError.h" #include "utils/NumericConstant.hpp" -namespace oopse { - - RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){ - +namespace OpenMD { + + RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), + inertiaTensor_(0.0){ } - + void RigidBody::setPrevA(const RotMat3x3d& a) { ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; - //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; - - for (int i =0 ; i < atoms_.size(); ++i){ + + for (unsigned int i = 0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setPrevA(a * refOrients_[i]); + atoms_[i]->setPrevA(refOrients_[i].transpose() * a); } } - + } - - + + void RigidBody::setA(const RotMat3x3d& a) { ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; - //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; - for (int i =0 ; i < atoms_.size(); ++i){ + for (unsigned int i = 0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setA(a * refOrients_[i]); + atoms_[i]->setA(refOrients_[i].transpose() * a); } } } - + void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) { ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a; - //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_; - - for (int i =0 ; i < atoms_.size(); ++i){ + + for (unsigned int i = 0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { - atoms_[i]->setA(a * refOrients_[i], snapshotNo); + atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo); } } - + } - + Mat3x3d RigidBody::getI() { return inertiaTensor_; } - - std::vector RigidBody::getGrad() { - std::vector grad(6, 0.0); + + std::vector RigidBody::getGrad() { + std::vector grad(6, 0.0); Vector3d force; Vector3d torque; Vector3d myEuler; - double phi, theta, psi; - double cphi, sphi, ctheta, stheta; + RealType phi, theta; + // RealType psi; + RealType cphi, sphi, ctheta, stheta; Vector3d ephi; Vector3d etheta; Vector3d epsi; - + force = getFrc(); torque =getTrq(); myEuler = getA().toEulerAngles(); - + phi = myEuler[0]; theta = myEuler[1]; - psi = myEuler[2]; - + // psi = myEuler[2]; + cphi = cos(phi); sphi = sin(phi); ctheta = cos(theta); stheta = sin(theta); - + // get unit vectors along the phi, theta and psi rotation axes - + ephi[0] = 0.0; ephi[1] = 0.0; ephi[2] = 1.0; - + + //etheta[0] = -sphi; + //etheta[1] = cphi; + //etheta[2] = 0.0; + etheta[0] = cphi; etheta[1] = sphi; - etheta[2] = 0.0; - + etheta[2] = 0.0; + epsi[0] = stheta * cphi; epsi[1] = stheta * sphi; epsi[2] = ctheta; - + //gradient is equal to -force for (int j = 0 ; j<3; j++) grad[j] = -force[j]; - + for (int j = 0; j < 3; j++ ) { - + grad[3] += torque[j]*ephi[j]; grad[4] += torque[j]*etheta[j]; grad[5] += torque[j]*epsi[j]; - + } return grad; } - + void RigidBody::accept(BaseVisitor* v) { v->visit(this); } /**@todo need modification */ void RigidBody::calcRefCoords() { - double mtmp; + RealType mtmp; Vector3d refCOM(0.0); mass_ = 0.0; for (std::size_t i = 0; i < atoms_.size(); ++i) { @@ -157,7 +160,7 @@ namespace oopse { refCOM += refCoords_[i]*mtmp; } refCOM /= mass_; - + // Next, move the origin of the reference coordinate system to the COM: for (std::size_t i = 0; i < atoms_.size(); ++i) { refCoords_[i] -= refCOM; @@ -169,20 +172,20 @@ namespace oopse { Mat3x3d IAtom(0.0); mtmp = atoms_[i]->getMass(); IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; - double r2 = refCoords_[i].lengthSquare(); + RealType r2 = refCoords_[i].lengthSquare(); IAtom(0, 0) += mtmp * r2; IAtom(1, 1) += mtmp * r2; IAtom(2, 2) += mtmp * r2; - + Itmp += IAtom; + //project the inertial moment of directional atoms into this rigid body if (atoms_[i]->isDirectional()) { - IAtom += atoms_[i]->getI(); - Itmp += refOrients_[i].transpose() * IAtom * refOrients_[i]; - } else { - Itmp += IAtom; - } + Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; + } } + // std::cout << Itmp << std::endl; + //diagonalize Vector3d evals; Mat3x3d::diagonalize(Itmp, evals, sU_); @@ -194,7 +197,7 @@ namespace oopse { int nLinearAxis = 0; for (int i = 0; i < 3; i++) { - if (fabs(evals[i]) < oopse::epsilon) { + if (fabs(evals[i]) < OpenMD::epsilon) { linear_ = true; linearAxis_ = i; ++ nLinearAxis; @@ -204,7 +207,7 @@ namespace oopse { if (nLinearAxis > 1) { sprintf( painCave.errMsg, "RigidBody error.\n" - "\tOOPSE found more than one axis in this rigid body with a vanishing \n" + "\tOpenMD found more than one axis in this rigid body with a vanishing \n" "\tmoment of inertia. This can happen in one of three ways:\n" "\t 1) Only one atom was specified, or \n" "\t 2) All atoms were specified at the same location, or\n" @@ -224,9 +227,17 @@ namespace oopse { Vector3d rpos; Vector3d frc(0.0); Vector3d trq(0.0); + Vector3d ef(0.0); Vector3d pos = this->getPos(); - for (int i = 0; i < atoms_.size(); i++) { + AtomType* atype; + int eCount = 0; + + int sl = ((snapshotMan_->getCurrentSnapshot())->*storage_).getStorageLayout(); + + for (unsigned int i = 0; i < atoms_.size(); i++) { + atype = atoms_[i]->getAtomType(); + afrc = atoms_[i]->getFrc(); apos = atoms_[i]->getPos(); rpos = apos - pos; @@ -244,12 +255,83 @@ namespace oopse { atrq = atoms_[i]->getTrq(); trq += atrq; } + + if ((sl & DataStorage::dslElectricField) && (atype->isElectrostatic())) { + ef += atoms_[i]->getElectricField(); + eCount++; + } + } + addFrc(frc); + addTrq(trq); + + if (sl & DataStorage::dslElectricField) { + ef /= eCount; + setElectricField(ef); + } + + } + + Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() { + Vector3d afrc; + Vector3d atrq; + Vector3d apos; + Vector3d rpos; + Vector3d dfrc; + Vector3d frc(0.0); + Vector3d trq(0.0); + Vector3d ef(0.0); + AtomType* atype; + int eCount = 0; + + Vector3d pos = this->getPos(); + Mat3x3d tau_(0.0); + + int sl = ((snapshotMan_->getCurrentSnapshot())->*storage_).getStorageLayout(); + + for (unsigned int i = 0; i < atoms_.size(); i++) { + + afrc = atoms_[i]->getFrc(); + apos = atoms_[i]->getPos(); + rpos = apos - pos; + frc += afrc; + + trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1]; + trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2]; + trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0]; + + // If the atom has a torque associated with it, then we also need to + // migrate the torques onto the center of mass: + + if (atoms_[i]->isDirectional()) { + atrq = atoms_[i]->getTrq(); + trq += atrq; + } + if ((sl & DataStorage::dslElectricField) && (atype->isElectrostatic())) { + ef += atoms_[i]->getElectricField(); + eCount++; + } + + tau_(0,0) -= rpos[0]*afrc[0]; + tau_(0,1) -= rpos[0]*afrc[1]; + tau_(0,2) -= rpos[0]*afrc[2]; + tau_(1,0) -= rpos[1]*afrc[0]; + tau_(1,1) -= rpos[1]*afrc[1]; + tau_(1,2) -= rpos[1]*afrc[2]; + tau_(2,0) -= rpos[2]*afrc[0]; + tau_(2,1) -= rpos[2]*afrc[1]; + tau_(2,2) -= rpos[2]*afrc[2]; + } - - setFrc(frc); - setTrq(trq); - + addFrc(frc); + addTrq(trq); + + if (sl & DataStorage::dslElectricField) { + ef /= eCount; + setElectricField(ef); + } + + return tau_; } void RigidBody::updateAtoms() { @@ -271,7 +353,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(refOrients_[i] * a); + dAtom->setA(refOrients_[i].transpose() * a); } } @@ -298,7 +380,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(refOrients_[i] * a, frame); + dAtom->setA(refOrients_[i].transpose() * a, frame); } } @@ -328,7 +410,7 @@ namespace oopse { Vector3d velRot; - for (int i =0 ; i < refCoords_.size(); ++i) { + for (unsigned int i = 0 ; i < refCoords_.size(); ++i) { atoms_[i]->setVel(rbVel + mat * refCoords_[i]); } @@ -357,7 +439,7 @@ namespace oopse { Vector3d velRot; - for (int i =0 ; i < refCoords_.size(); ++i) { + for (unsigned int i = 0 ; i < refCoords_.size(); ++i) { atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame); } @@ -483,7 +565,7 @@ namespace oopse { "RigidBody error.\n" "\tAtom %s does not have a position specified.\n" "\tThis means RigidBody cannot set up reference coordinates.\n", - ats->getType() ); + ats->getType().c_str() ); painCave.isFatal = 1; simError(); } @@ -503,7 +585,7 @@ namespace oopse { "RigidBody error.\n" "\tAtom %s does not have an orientation specified.\n" "\tThis means RigidBody cannot set up reference orientations.\n", - ats->getType() ); + ats->getType().c_str() ); painCave.isFatal = 1; simError(); }