--- branches/development/src/primitives/DirectionalAtom.cpp 2012/05/18 21:44:02 1710 +++ branches/development/src/primitives/DirectionalAtom.cpp 2012/09/13 14:10:11 1798 @@ -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) { @@ -117,7 +147,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; @@ -129,7 +160,7 @@ namespace OpenMD { phi = myEuler[0]; theta = myEuler[1]; - psi = myEuler[2]; + // psi = myEuler[2]; cphi = cos(phi); sphi = sin(phi);