--- trunk/src/primitives/RigidBody.cpp 2010/03/30 15:05:38 1424 +++ trunk/src/primitives/RigidBody.cpp 2015/03/07 21:41:51 2071 @@ -35,8 +35,9 @@ * * [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] Vardeman & Gezelter, in progress (2009). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include #include @@ -52,7 +53,7 @@ namespace OpenMD { 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 OpenMD { 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 OpenMD { 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 OpenMD { 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 OpenMD { phi = myEuler[0]; theta = myEuler[1]; - psi = myEuler[2]; + // psi = myEuler[2]; cphi = cos(phi); sphi = sin(phi); @@ -226,10 +226,18 @@ namespace OpenMD { 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; @@ -246,10 +254,21 @@ namespace OpenMD { 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() { @@ -259,12 +278,20 @@ namespace OpenMD { 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++) { + atype = atoms_[i]->getAtomType(); + afrc = atoms_[i]->getFrc(); apos = atoms_[i]->getPos(); rpos = apos - pos; @@ -282,6 +309,11 @@ namespace OpenMD { 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]; @@ -296,6 +328,12 @@ namespace OpenMD { } addFrc(frc); addTrq(trq); + + if (sl & DataStorage::dslElectricField) { + ef /= eCount; + setElectricField(ef); + } + return tau_; } @@ -317,7 +355,7 @@ namespace OpenMD { if (atoms_[i]->isDirectional()) { - dAtom = (DirectionalAtom *) atoms_[i]; + dAtom = dynamic_cast(atoms_[i]); dAtom->setA(refOrients_[i].transpose() * a); } @@ -344,7 +382,7 @@ namespace OpenMD { if (atoms_[i]->isDirectional()) { - dAtom = (DirectionalAtom *) atoms_[i]; + dAtom = dynamic_cast(atoms_[i]); dAtom->setA(refOrients_[i].transpose() * a, frame); } @@ -353,7 +391,7 @@ namespace OpenMD { } void RigidBody::updateAtomVel() { - Mat3x3d skewMat;; + Mat3x3d skewMat; Vector3d ji = getJ(); Mat3x3d I = getI(); @@ -375,7 +413,7 @@ namespace OpenMD { 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]); } @@ -404,7 +442,7 @@ namespace OpenMD { 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); } @@ -451,25 +489,26 @@ namespace OpenMD { Vector3d ji = getJ(); Mat3x3d I = getI(); - skewMat(0, 0) =0; - skewMat(0, 1) = ji[2] /I(2, 2); - skewMat(0, 2) = -ji[1] /I(1, 1); + skewMat(0, 0) = 0; + skewMat(0, 1) = ji[2] / I(2, 2); + skewMat(0, 2) = -ji[1] / I(1, 1); - skewMat(1, 0) = -ji[2] /I(2, 2); - skewMat(1, 1) = 0; - skewMat(1, 2) = ji[0]/I(0, 0); + skewMat(1, 0) = -ji[2] / I(2, 2); + skewMat(1, 1) = 0; + skewMat(1, 2) = ji[0]/ I(0, 0); - skewMat(2, 0) =ji[1] /I(1, 1); - skewMat(2, 1) = -ji[0]/I(0, 0); - skewMat(2, 2) = 0; + skewMat(2, 0) = ji[1] / I(1, 1); + skewMat(2, 1) = -ji[0] / I(0, 0); + skewMat(2, 2) = 0; velRot = (getA() * skewMat).transpose() * ref; - vel =getVel() + velRot; + vel = getVel() + velRot; return true; } else { - std::cerr << index << " is an invalid index, current rigid body contains " + std::cerr << index + << " is an invalid index, current rigid body contains " << atoms_.size() << "atoms" << std::endl; return false; }