--- trunk/src/primitives/RigidBody.cpp 2008/01/23 16:38:22 1211 +++ 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,13 +28,23 @@ * 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 { +namespace OpenMD { RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){ @@ -52,7 +53,7 @@ namespace oopse { void RigidBody::setPrevA(const RotMat3x3d& a) { ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; - for (int i =0 ; i < atoms_.size(); ++i){ + for (unsigned int i = 0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { atoms_[i]->setPrevA(refOrients_[i].transpose() * a); } @@ -64,7 +65,7 @@ namespace oopse { void RigidBody::setA(const RotMat3x3d& a) { ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; - for (int i =0 ; i < atoms_.size(); ++i){ + for (unsigned int i = 0 ; i < atoms_.size(); ++i){ if (atoms_[i]->isDirectional()) { atoms_[i]->setA(refOrients_[i].transpose() * a); } @@ -73,10 +74,8 @@ namespace oopse { 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(refOrients_[i].transpose() * a, snapshotNo); } @@ -93,7 +92,8 @@ namespace oopse { Vector3d force; Vector3d torque; Vector3d myEuler; - RealType phi, theta, psi; + RealType phi, theta; + // RealType psi; RealType cphi, sphi, ctheta, stheta; Vector3d ephi; Vector3d etheta; @@ -105,7 +105,7 @@ namespace oopse { phi = myEuler[0]; theta = myEuler[1]; - psi = myEuler[2]; + // psi = myEuler[2]; cphi = cos(phi); sphi = sin(phi); @@ -118,9 +118,13 @@ namespace oopse { 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; @@ -193,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; @@ -203,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" @@ -222,10 +226,18 @@ namespace oopse { Vector3d apos; Vector3d rpos; Vector3d frc(0.0); - Vector3d trq(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; @@ -242,10 +254,21 @@ namespace oopse { if (atoms_[i]->isDirectional()) { 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() { @@ -255,11 +278,17 @@ namespace oopse { Vector3d rpos; Vector3d dfrc; Vector3d frc(0.0); - Vector3d trq(0.0); + Vector3d trq(0.0); + Vector3d ef(0.0); + AtomType* atype; + int eCount = 0; + Vector3d pos = this->getPos(); Mat3x3d tau_(0.0); - for (int i = 0; i < atoms_.size(); i++) { + int sl = ((snapshotMan_->getCurrentSnapshot())->*storage_).getStorageLayout(); + + for (unsigned int i = 0; i < atoms_.size(); i++) { afrc = atoms_[i]->getFrc(); apos = atoms_[i]->getPos(); @@ -277,6 +306,10 @@ namespace oopse { 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]; @@ -292,6 +325,12 @@ namespace oopse { } addFrc(frc); addTrq(trq); + + if (sl & DataStorage::dslElectricField) { + ef /= eCount; + setElectricField(ef); + } + return tau_; } @@ -371,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]); } @@ -400,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); }