--- trunk/src/primitives/RigidBody.cpp	2005/10/12 21:00:59	663
+++ trunk/src/primitives/RigidBody.cpp	2007/04/06 21:53:43	1126
@@ -51,11 +51,10 @@ namespace oopse {
 
   void RigidBody::setPrevA(const RotMat3x3d& a) {
     ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
-    //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
 
     for (int i =0 ; i < atoms_.size(); ++i){
       if (atoms_[i]->isDirectional()) {
-	atoms_[i]->setPrevA(a * refOrients_[i]);
+	atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
       }
     }
 
@@ -64,11 +63,10 @@ namespace oopse {
       
   void RigidBody::setA(const RotMat3x3d& a) {
     ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
-    //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
 
     for (int i =0 ; i < atoms_.size(); ++i){
       if (atoms_[i]->isDirectional()) {
-	atoms_[i]->setA(a * refOrients_[i]);
+	atoms_[i]->setA(refOrients_[i].transpose() * a);
       }
     }
   }    
@@ -79,7 +77,7 @@ namespace oopse {
 
     for (int i =0 ; i < atoms_.size(); ++i){
       if (atoms_[i]->isDirectional()) {
-	atoms_[i]->setA(a * refOrients_[i], snapshotNo);
+	atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
       }
     }
 
@@ -89,13 +87,13 @@ namespace oopse {
     return inertiaTensor_;
   }    
 
-  std::vector<double> RigidBody::getGrad() {
-    std::vector<double> grad(6, 0.0);
+  std::vector<RealType> RigidBody::getGrad() {
+    std::vector<RealType> grad(6, 0.0);
     Vector3d force;
     Vector3d torque;
     Vector3d myEuler;
-    double phi, theta, psi;
-    double cphi, sphi, ctheta, stheta;
+    RealType phi, theta, psi;
+    RealType cphi, sphi, ctheta, stheta;
     Vector3d ephi;
     Vector3d etheta;
     Vector3d epsi;
@@ -148,7 +146,7 @@ namespace oopse {
 
   /**@todo need modification */
   void  RigidBody::calcRefCoords() {
-    double mtmp;
+    RealType mtmp;
     Vector3d refCOM(0.0);
     mass_ = 0.0;
     for (std::size_t i = 0; i < atoms_.size(); ++i) {
@@ -169,7 +167,7 @@ namespace oopse {
       Mat3x3d IAtom(0.0);  
       mtmp = atoms_[i]->getMass();
       IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
-      double r2 = refCoords_[i].lengthSquare();
+      RealType r2 = refCoords_[i].lengthSquare();
       IAtom(0, 0) += mtmp * r2;
       IAtom(1, 1) += mtmp * r2;
       IAtom(2, 2) += mtmp * r2;
@@ -181,7 +179,7 @@ namespace oopse {
       } 
     }
 
-    std::cout << Itmp << std::endl;
+    //    std::cout << Itmp << std::endl;
 
     //diagonalize 
     Vector3d evals;
@@ -223,7 +221,7 @@ namespace oopse {
     Vector3d apos;
     Vector3d rpos;
     Vector3d frc(0.0);
-    Vector3d trq(0.0);
+    Vector3d trq(0.0);    
     Vector3d pos = this->getPos();
     for (int i = 0; i < atoms_.size(); i++) {
 
@@ -243,13 +241,56 @@ namespace oopse {
       if (atoms_[i]->isDirectional()) {
 	atrq = atoms_[i]->getTrq();
 	trq += atrq;
-      }
+      }      
+    }         
+    addFrc(frc);
+    addTrq(trq);    
+  }
+
+  Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
+    Vector3d afrc;
+    Vector3d atrq;
+    Vector3d apos;
+    Vector3d rpos;
+    Vector3d frc(0.0);
+    Vector3d trq(0.0);    
+    Vector3d pos = this->getPos();
+    Mat3x3d tau_(0.0);
+
+    for (int i = 0; i < atoms_.size(); i++) {
+
+      afrc = atoms_[i]->getFrc();
+      apos = atoms_[i]->getPos();
+      rpos = apos - pos;
         
-    }
-    
-    setFrc(frc);
-    setTrq(trq);
-    
+      frc += afrc;
+
+      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
+      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
+      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
+
+      // If the atom has a torque associated with it, then we also need to 
+      // migrate the torques onto the center of mass:
+
+      if (atoms_[i]->isDirectional()) {
+	atrq = atoms_[i]->getTrq();
+	trq += atrq;
+      }
+      
+      tau_(0,0) -= rpos[0]*afrc[0];
+      tau_(0,1) -= rpos[0]*afrc[1];
+      tau_(0,2) -= rpos[0]*afrc[2];
+      tau_(1,0) -= rpos[1]*afrc[0];
+      tau_(1,1) -= rpos[1]*afrc[1];
+      tau_(1,2) -= rpos[1]*afrc[2];
+      tau_(2,0) -= rpos[2]*afrc[0];
+      tau_(2,1) -= rpos[2]*afrc[1];
+      tau_(2,2) -= rpos[2]*afrc[2];
+      
+    }         
+    addFrc(frc);
+    addTrq(trq);
+    return tau_;
   }
 
   void  RigidBody::updateAtoms() {
@@ -271,7 +312,7 @@ namespace oopse {
       if (atoms_[i]->isDirectional()) {
           
 	dAtom = (DirectionalAtom *) atoms_[i];
-	dAtom->setA(refOrients_[i] * a);
+	dAtom->setA(refOrients_[i].transpose() * a);
       }
 
     }
@@ -298,7 +339,7 @@ namespace oopse {
       if (atoms_[i]->isDirectional()) {
           
 	dAtom = (DirectionalAtom *) atoms_[i];
-	dAtom->setA(refOrients_[i] * a, frame);
+	dAtom->setA(refOrients_[i].transpose() * a, frame);
       }
 
     }
@@ -483,7 +524,7 @@ namespace oopse {
 	       "RigidBody error.\n"
 	       "\tAtom %s does not have a position specified.\n"
 	       "\tThis means RigidBody cannot set up reference coordinates.\n",
-	       ats->getType() );
+	       ats->getType().c_str() );
       painCave.isFatal = 1;
       simError();
     }
@@ -503,7 +544,7 @@ namespace oopse {
 		 "RigidBody error.\n"
 		 "\tAtom %s does not have an orientation specified.\n"
 		 "\tThis means RigidBody cannot set up reference orientations.\n",
-		 ats->getType() );
+		 ats->getType().c_str() );
 	painCave.isFatal = 1;
 	simError();
       }