--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp	2004/12/02 16:04:19	1832
+++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp	2004/12/03 21:30:30	1843
@@ -42,32 +42,24 @@ SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
 namespace oopse {
 
 SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, 
-                                ForceField* ff, Globals* globals) : 
-                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0), 
-                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), 
-                                nCutoffGroups_(0), nConstraints_(0), nZconstraint_(0), sman_(NULL),
-                                fortranInitialized_(false) {
+                                ForceField* ff, Globals* simParams) : 
+                                forceField_(ff), simParams_(simParams), 
+                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
+                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), 
+                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
+                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
+                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
+                                sman_(NULL), fortranInitialized_(false) {
 
+            
     std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
     MoleculeStamp* molStamp;
     int nMolWithSameStamp;
-    int nCutoffAtoms; // number of atoms belong to cutoff groups
-    int nGroups;          //total cutoff groups defined in meta-data file
-    CutoffGroupStamp* cgStamp;
-    int nAtomsInGroups;
-    int nCutoffGroupsInStamp;
-    
+    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
+    int nGroups = 0;          //total cutoff groups defined in meta-data file
+    CutoffGroupStamp* cgStamp;    
     RigidBodyStamp* rbStamp;
-    int nAtomsInRigidBodies;
-    int nRigidBodiesInStamp;
-    int nRigidAtoms;
-    int nRigidBodies; 
-        
-    nGlobalAtoms_ =  0;
-    
-    nGroups = 0;
-    nCutoffAtoms = 0;
-    nRigidBodies = 0;
+    int nRigidAtoms = 0;
     
     for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
         molStamp = i->first;
@@ -80,8 +72,8 @@ SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
 
 
         //calculate atoms in cutoff groups
-        nAtomsInGroups = 0;
-        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
+        int nAtomsInGroups = 0;
+        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
         
         for (int j=0; j < nCutoffGroupsInStamp; j++) {
             cgStamp = molStamp->getCutoffGroup(j);
@@ -92,15 +84,15 @@ SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
         nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
 
         //calculate atoms in rigid bodies
-        nAtomsInRigidBodies = 0;
-        nRigidBodiesInStamp = molStamp->getNCutoffGroups();
+        int nAtomsInRigidBodies = 0;
+        int nRigidBodiesInStamp = molStamp->getNCutoffGroups();
         
         for (int j=0; j < nRigidBodiesInStamp; j++) {
             rbStamp = molStamp->getRigidBody(j);
-            nRigidBodiesInStamp += rbStamp->getNMembers();
+            nAtomsInRigidBodies += rbStamp->getNMembers();
         }
 
-        nRigidBodies += nRigidBodiesInStamp * nMolWithSameStamp;
+        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
         nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
         
     }
@@ -111,14 +103,11 @@ SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
     //file plus the number of cutoff groups defined in meta-data file
     nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
 
-    //every free atom (atom does not belong to rigid bodies) is a rigid body
-    //therefore the total number of cutoff groups in the system is equal to 
+    //every free atom (atom does not belong to rigid bodies) is an integrable object
+    //therefore the total number of  integrable objects in the system is equal to 
     //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
     //file plus the number of  rigid bodies defined in meta-data file
-    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nRigidBodies;
-
-    //initialize globalGroupMembership_, every element of this array will be 0
-    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
+    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
 
     nGlobalMols_ = molStampIds_.size();
 
@@ -134,7 +123,7 @@ SimInfo::~SimInfo() {
     MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
     
     delete sman_;
-    delete globals_;
+    delete simParams_;
     delete forceField_;
 
 }
@@ -144,7 +133,7 @@ bool SimInfo::addMolecule(Molecule* mol) {
     MoleculeIterator i;
 
     i = molecules_.find(mol->getGlobalIndex());
-    if (i != molecules_.end() ) {
+    if (i == molecules_.end() ) {
 
         molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
         
@@ -462,9 +451,9 @@ void SimInfo::setupSimType() {
     int useFLARB = 0; //it is not in AtomType yet
     int useDirectionalAtom = 0;    
     int useElectrostatics = 0;
-    //usePBC and useRF are from globals
-    bool usePBC = globals_->getPBC();
-    bool useRF = globals_->getUseRF();
+    //usePBC and useRF are from simParams
+    bool usePBC = simParams_->getPBC();
+    bool useRF = simParams_->getUseRF();
 
     //loop over all of the atom types
     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
@@ -543,8 +532,8 @@ void SimInfo::setupSimType() {
 
     if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
 
-        if (globals_->haveDielectric()) {
-            fInfo_.dielect = globals_->getDielectric();
+        if (simParams_->haveDielectric()) {
+            fInfo_.dielect = simParams_->getDielectric();
         } else {
             sprintf(painCave.errMsg,
                     "SimSetup Error: No Dielectric constant was set.\n"
@@ -732,7 +721,7 @@ void SimInfo::setupCutoff() {
     
     if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
         
-        if (!globals_->haveRcut()){
+        if (!simParams_->haveRcut()){
             sprintf(painCave.errMsg,
                 "SimCreator Warning: No value was set for the cutoffRadius.\n"
                 "\tOOPSE will use a default value of 15.0 angstroms"
@@ -741,10 +730,10 @@ void SimInfo::setupCutoff() {
             simError();
             rcut_ = 15.0;
         } else{
-            rcut_ = globals_->getRcut();
+            rcut_ = simParams_->getRcut();
         }
 
-        if (!globals_->haveRsw()){
+        if (!simParams_->haveRsw()){
             sprintf(painCave.errMsg,
                 "SimCreator Warning: No value was set for switchingRadius.\n"
                 "\tOOPSE will use a default value of\n"
@@ -753,22 +742,22 @@ void SimInfo::setupCutoff() {
             simError();
             rsw_ = 0.95 * rcut_;
         } else{
-            rsw_ = globals_->getRsw();
+            rsw_ = simParams_->getRsw();
         }
 
     } else {
         // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
         //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
         
-        if (globals_->haveRcut()) {
-            rcut_ = globals_->getRcut();
+        if (simParams_->haveRcut()) {
+            rcut_ = simParams_->getRcut();
         } else {
             //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
             rcut_ = calcMaxCutoffRadius();
         }
 
-        if (globals_->haveRsw()) {
-            rsw_  = globals_->getRsw();
+        if (simParams_->haveRsw()) {
+            rsw_  = simParams_->getRsw();
         } else {
             rsw_ = rcut_;
         }
@@ -828,9 +817,48 @@ std::ostream& operator <<(std::ostream& o, SimInfo& in
     
 }
 
+Vector3d SimInfo::getComVel(){ 
+    SimInfo::MoleculeIterator i;
+    Molecule* mol;
+
+    Vector3d comVel(0.0);
+    double totalMass = 0.0;
+    
+ 
+    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
+        double mass = mol->getMass();
+        totalMass += mass;
+        comVel += mass * mol->getComVel();
+    }  
+
+    comVel /= totalMass;
+
+    return comVel;
+}
+
+Vector3d SimInfo::getCom(){ 
+    SimInfo::MoleculeIterator i;
+    Molecule* mol;
+
+    Vector3d com(0.0);
+    double totalMass = 0.0;
+     
+    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
+        double mass = mol->getMass();
+        totalMass += mass;
+        com += mass * mol->getCom();
+    }  
+
+    com /= totalMass;
+
+    return com;
+
+}        
+
 std::ostream& operator <<(std::ostream& o, SimInfo& info) {
 
     return o;
 }
 
 }//end namespace oopse
+