--- trunk/OOPSE/libmdtools/Molecule.cpp	2003/03/27 21:07:14	428
+++ trunk/OOPSE/libmdtools/Molecule.cpp	2004/08/23 15:11:36	1452
@@ -1,7 +1,8 @@
-#include <cstdlib>
+#include <stdlib.h>
 
 
 #include "Molecule.hpp"
+#include "simError.h"
 
 
 
@@ -11,13 +12,12 @@ Molecule::Molecule( void ){
   myBonds = NULL;
   myBends = NULL;
   myTorsions = NULL;
-
 }
 
-
-
 Molecule::~Molecule( void ){
   int i;
+  CutoffGroup* cg;
+  vector<CutoffGroup*>::iterator iter;
   
   if( myAtoms != NULL ){
     for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
@@ -39,35 +39,79 @@ Molecule::~Molecule( void ){
     delete[] myTorsions;
   }
 
-  if( myExcludes != NULL ){
-    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
-    delete[] myExcludes;
-  }
+  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
+    delete cg;
+  myCutoffGroups.clear();
+  
 }
 
 
 void Molecule::initialize( molInit &theInit ){
 
+  CutoffGroup* curCutoffGroup;
+  vector<CutoffGroup*>::iterator iterCutoff;
+  Atom* cutoffAtom;
+  vector<Atom*>::iterator iterAtom;
+  int atomIndex;
+  GenericData* gdata;
+  ConsRbData* rbData;
+  RigidBody* oldRb;
+  
   nAtoms = theInit.nAtoms;
   nMembers = nAtoms;
   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;
+
+  myIntegrableObjects = theInit.myIntegrableObjects;
+
+  for (int i = 0; i < myRigidBodies.size(); i++){
+      myRigidBodies[i]->calcRefCoords();
+    //just a quick hack
     
+    gdata = myRigidBodies[i]->getProperty("OldState");
+    if(gdata != NULL){
+      rbData = dynamic_cast<ConsRbData*>(gdata);
+      if(rbData ==NULL)
+        cerr << "dynamic_cast to ConsRbData Error in Molecule::initialize()" << endl;
+      else{
+        oldRb = rbData->getData();
+        oldRb->calcRefCoords();
+      }
+    }//end if(gata != NULL)
+    
+  }//end for(int i = 0; i < myRigidBodies.size(); i++)
+
+  myCutoffGroups = theInit.myCutoffGroups;
+  nCutoffGroups = myCutoffGroups.size();
+
+  myConstraintPairs = theInit.myConstraintPairs;
+  
 }
 
 void Molecule::calcForces( void ){
   
   int i;
+  double com[3];
 
+  for(i=0; i<myRigidBodies.size(); i++) {
+    myRigidBodies[i]->updateAtoms();
+  }
+
+  //calculate the center of mass of the molecule
+  //getCOM(com);  
+  //for(int i = 0; i < nAtoms; i ++)
+  //  myAtoms[i]->setRc(com);  
+  
+
   for(i=0; i<nBonds; i++){
     myBonds[i]->calc_forces();
   }
@@ -79,6 +123,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
+
 }
 
 
@@ -86,6 +133,10 @@ double Molecule::getPotential( void ){
   
   int i;
   double myPot = 0.0;
+
+  for(i=0; i<myRigidBodies.size(); i++) {
+    myRigidBodies[i]->updateAtoms();
+  }
    
   for(i=0; i<nBonds; i++){
     myPot += myBonds[i]->get_potential();
@@ -101,3 +152,130 @@ double Molecule::getPotential( void ){
 
   return myPot;
 }
+
+void Molecule::printMe( void ){
+  
+  int i;
+
+  for(i=0; i<nBonds; i++){
+    myBonds[i]->printMe();
+  }
+
+  for(i=0; i<nBends; i++){
+    myBends[i]->printMe();
+  }
+
+  for(i=0; i<nTorsions; i++){
+    myTorsions[i]->printMe();
+  }
+
+}
+
+void Molecule::moveCOM(double delta[3]){
+  double aPos[3];
+  int i, j;
+
+  for(i=0; i<myIntegrableObjects.size(); i++) {
+    if(myIntegrableObjects[i] != NULL ) {
+      
+      myIntegrableObjects[i]->getPos( aPos );
+      
+      for (j=0; j< 3; j++) 
+        aPos[j] += delta[j];
+
+      myIntegrableObjects[i]->setPos( aPos );
+    }
+  }
+
+  for(i=0; i<myRigidBodies.size(); i++) {
+
+      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 < myRigidBodies.size(); i++) {
+    myRigidBodies[i]->calcForcesAndTorques();   
+  }
+}
+
+void Molecule::getCOM( double COM[3] ) {
+
+  double mass, mtot;
+  double aPos[3];
+  int i, j;
+
+  for (j=0; j<3; j++) 
+    COM[j] = 0.0;
+
+  mtot   = 0.0;
+
+  for (i=0; i < myIntegrableObjects.size(); i++) {
+    if (myIntegrableObjects[i] != NULL) {
+
+      mass = myIntegrableObjects[i]->getMass();
+      mtot   += mass;
+      
+      myIntegrableObjects[i]->getPos( aPos );
+
+      for( j = 0; j < 3; j++) 
+        COM[j] += aPos[j] * mass;
+
+    }
+  }
+
+  for (j = 0; j < 3; j++) 
+    COM[j] /= mtot; 
+}
+
+double Molecule::getCOMvel( double COMvel[3] ) {
+
+  double mass, mtot;
+  double aVel[3];
+  int i, j;
+
+
+  for (j=0; j<3; j++) 
+    COMvel[j] = 0.0;
+
+  mtot   = 0.0;
+
+  for (i=0; i < myIntegrableObjects.size(); i++) {
+    if (myIntegrableObjects[i] != NULL) {
+
+      mass = myIntegrableObjects[i]->getMass();
+      mtot   += mass;
+
+      myIntegrableObjects[i]->getVel(aVel);
+
+      for (j=0; j<3; j++) 
+        COMvel[j] += aVel[j]*mass;
+
+    }
+  }
+
+  for (j=0; j<3; j++) 
+    COMvel[j] /= mtot;
+ 
+  return mtot;
+
+}
+
+double Molecule::getTotalMass()
+{
+
+  double totalMass;
+  
+  totalMass = 0;
+  for(int i =0; i < myIntegrableObjects.size(); i++){
+    totalMass += myIntegrableObjects[i]->getMass();
+  }
+
+  return totalMass;
+}