--- branches/development/src/primitives/RigidBody.cpp	2012/09/13 14:10:11	1798
+++ branches/development/src/primitives/RigidBody.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).
  */
@@ -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();
+    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);
 
+    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_;
   }