--- branches/development/src/primitives/RigidBody.cpp	2010/07/09 23:08:25	1465
+++ branches/development/src/primitives/RigidBody.cpp	2013/01/30 14:43:08	1844
@@ -36,7 +36,8 @@
  * [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).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 #include <algorithm>
 #include <math.h>
@@ -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,11 +278,17 @@ 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++) {
       
       afrc = atoms_[i]->getFrc();
       apos = atoms_[i]->getPos();
@@ -282,6 +307,10 @@ 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 +325,12 @@ namespace OpenMD {
     }
     addFrc(frc);
     addTrq(trq);
+
+    if (sl & DataStorage::dslElectricField) {
+      ef /= eCount;
+      setElectricField(ef);
+    }
+
     return tau_;
   }
 
@@ -375,7 +410,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 +439,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);
     }