--- branches/development/src/primitives/DirectionalAtom.cpp	2012/05/18 21:44:02	1710
+++ branches/development/src/primitives/DirectionalAtom.cpp	2013/02/20 15:39:39	1850
@@ -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) {
@@ -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);