--- trunk/OOPSE/libmdtools/Molecule.cpp	2003/07/15 14:28:54	607
+++ trunk/OOPSE/libmdtools/Molecule.cpp	2004/04/12 20:32:20	1097
@@ -1,4 +1,4 @@
-#include <cstdlib>
+#include <stdlib.h>
 
 
 #include "Molecule.hpp"
@@ -12,6 +12,7 @@ Molecule::Molecule( void ){
   myBonds = NULL;
   myBends = NULL;
   myTorsions = NULL;
+  myRigidBodies = NULL;
 
 }
 
@@ -40,10 +41,12 @@ Molecule::~Molecule( void ){
     delete[] myTorsions;
   }
 
-  if( myExcludes != NULL ){
-    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
-    delete[] myExcludes;
+  if( myRigidBodies != NULL ){
+    for(i=0; i<nRigidBodies; i++) if(myRigidBodies[i] != NULL ) 
+      delete myRigidBodies[i];
+    delete[] myRigidBodies;
   }
+  
 }
 
 
@@ -54,14 +57,14 @@ void Molecule::initialize( molInit &theInit ){
   nBonds = theInit.nBonds;
   nBends = theInit.nBends;
   nTorsions = theInit.nTorsions;
-  nExcludes = theInit.nExcludes;
+  nRigidBodies = theInit.nRigidBodies;
   nOriented = theInit.nOriented;
 
   myAtoms = theInit.myAtoms;
   myBonds = theInit.myBonds;
   myBends = theInit.myBends;
   myTorsions = theInit.myTorsions;
-  myExcludes = theInit.myExcludes;
+  myRigidBodies = theInit.myRigidBodies;
     
 }
 
@@ -69,6 +72,10 @@ void Molecule::calcForces( void ){
   
   int i;
 
+  for(i=0; i<nRigidBodies; i++) {
+    myRigidBodies[i]->updateAtoms();
+  }
+
   for(i=0; i<nBonds; i++){
     myBonds[i]->calc_forces();
   }
@@ -80,6 +87,9 @@ void Molecule::calcForces( void ){
   for(i=0; i<nTorsions; i++){
     myTorsions[i]->calc_forces();
   }
+
+  // Rigid Body forces and torques are done after the fortran force loop
+
 }
 
 
@@ -118,34 +128,57 @@ void Molecule::printMe( void ){
   for(i=0; i<nTorsions; i++){
     myTorsions[i]->printMe();
   }
+
 }
 
 void Molecule::moveCOM(double delta[3]){
-  double x, y, z;
-  int i;
+  double aPos[3];
+  int i, j;
 
   for(i=0; i<nAtoms; i++) {
     if(myAtoms[i] != NULL ) {
+      
+      myAtoms[i]->getPos( aPos );
+      
+      for (j=0; j< 3; j++) 
+        aPos[j] += delta[j];
 
-      x = myAtoms[i]->getX() + delta[0];
-      y = myAtoms[i]->getY() + delta[1];
-      z = myAtoms[i]->getZ() + delta[2];
+      myAtoms[i]->setPos( aPos );
+    }
+  }
 
-      myAtoms[i]->setX(x);
-      myAtoms[i]->setY(y);
-      myAtoms[i]->setZ(z);
+  for(i=0; i<nRigidBodies; i++) {
+
+    if (myRigidBodies[i] != NULL) {
+      
+      myRigidBodies[i]->getPos( aPos );
+
+      for (j=0; j< 3; j++) 
+        aPos[j] += delta[j];
+      
+      myRigidBodies[i]->setPos( aPos );
     }
+  }  
+}
+
+void Molecule::atoms2rigidBodies( void ) {
+  int i;
+  for (i = 0; i < nRigidBodies; i++) {
+    if (myRigidBodies[i] != NULL) {
+      myRigidBodies[i]->calcForcesAndTorques();   
+    }
   }
 }
 
 void Molecule::getCOM( double COM[3] ) {
 
   double mass, mtot;
-  int i;
+  double aPos[3];
+  int i, j;
 
-  COM[0] = 0.0;
-  COM[1] = 0.0;
-  COM[2] = 0.0;
+  for (j=0; j<3; j++) 
+    COM[j] = 0.0;
+
   mtot   = 0.0;
 
   for (i=0; i < nAtoms; i++) {
@@ -153,17 +186,17 @@ void Molecule::getCOM( double COM[3] ) {
 
       mass = myAtoms[i]->getMass();
       mtot   += mass;
-      COM[0] += myAtoms[i]->getX() * mass;
-      COM[1] += myAtoms[i]->getY() * mass;
-      COM[2] += myAtoms[i]->getZ() * mass;
+      
+      myAtoms[i]->getPos( aPos );
 
+      for( j = 0; j < 3; j++) 
+        COM[j] += aPos[j] * mass;
+
     }
   }
 
-  COM[0] /= mtot;
-  COM[1] /= mtot;
-  COM[2] /= mtot;
- 
+  for (j = 0; j < 3; j++) 
+    COM[j] /= mtot; 
 }
 
 double Molecule::getCOMvel( double COMvel[3] ) {
@@ -198,3 +231,19 @@ double Molecule::getCOMvel( double COMvel[3] ) {
   return mtot;
 
 }
+
+double Molecule::getTotalMass()
+{
+  int natoms;
+  Atom** atoms;
+  double totalMass;
+  
+  natoms = getNAtoms();
+  atoms = getMyAtoms();
+  totalMass = 0;
+  for(int i =0; i < natoms; i++){
+    totalMass += atoms[i]->getMass();
+  }
+
+  return totalMass;
+}