--- trunk/src/primitives/DirectionalAtom.cpp 2013/06/13 14:26:09 1878 +++ trunk/src/primitives/DirectionalAtom.cpp 2013/06/16 15:15:42 1879 @@ -35,7 +35,7 @@ * * [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). + * [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). */ @@ -54,9 +54,12 @@ namespace OpenMD { I_ = da.getI(); MultipoleAdapter ma = MultipoleAdapter(dAtomType); - if (ma.isMultipole()) { - electroBodyFrame_ = ma.getElectroBodyFrame(); + if (ma.isDipole()) { + dipole_ = ma.getDipole(); } + if (ma.isQuadrupole()) { + quadrupole_ = ma.getQuadrupole(); + } // Check if one of the diagonal inertia tensor of this directional // atom is zero: @@ -86,26 +89,53 @@ namespace OpenMD { void DirectionalAtom::setPrevA(const RotMat3x3d& a) { ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a; + if (atomType_->isMultipole()) { - ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; + RotMat3x3d atrans = a.transpose(); + + if (atomType_->isDipole()) { + ((snapshotMan_->getPrevSnapshot())->*storage_).dipole[localIndex_] = atrans * dipole_; + } + + if (atomType_->isQuadrupole()) { + ((snapshotMan_->getPrevSnapshot())->*storage_).quadrupole[localIndex_] = atrans * quadrupole_ * a; + } } } void DirectionalAtom::setA(const RotMat3x3d& a) { ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a; - + if (atomType_->isMultipole()) { - ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; + RotMat3x3d atrans = a.transpose(); + + if (atomType_->isDipole()) { + ((snapshotMan_->getCurrentSnapshot())->*storage_).dipole[localIndex_] = atrans * dipole_; + } + + if (atomType_->isQuadrupole()) { + ((snapshotMan_->getCurrentSnapshot())->*storage_).quadrupole[localIndex_] = atrans * quadrupole_ * a; + } } + } void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) { ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a; - + if (atomType_->isMultipole()) { - ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_; + RotMat3x3d atrans = a.transpose(); + + if (atomType_->isDipole()) { + ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).dipole[localIndex_] = atrans * dipole_; + } + + if (atomType_->isQuadrupole()) { + ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).quadrupole[localIndex_] = atrans * quadrupole_ * a; + } } + } void DirectionalAtom::rotateBy(const RotMat3x3d& m) {