--- branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp	2004/12/02 03:12:25	1824
+++ branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp	2004/12/13 22:30:27	1883
@@ -39,35 +39,32 @@ namespace oopse {
 #include "utils/MemoryUtils.hpp"
 #include "utils/simError.h"
 
+#ifdef IS_MPI
+#include "UseTheForce/mpiComponentPlan.h"
+#include "UseTheForce/DarkSide/simParallel_interface.h"
+#endif 
+
 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 +77,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 +89,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,15 +108,12 @@ 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;
+    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
 
-    //initialize globalGroupMembership_, every element of this array will be 0
-    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
-
     nGlobalMols_ = molStampIds_.size();
 
 #ifdef IS_MPI    
@@ -134,7 +128,7 @@ SimInfo::~SimInfo() {
     MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
     
     delete sman_;
-    delete globals_;
+    delete simParams_;
     delete forceField_;
 
 }
@@ -144,9 +138,9 @@ bool SimInfo::addMolecule(Molecule* mol) {
     MoleculeIterator i;
 
     i = molecules_.find(mol->getGlobalIndex());
-    if (i != molecules_.end() ) {
+    if (i == molecules_.end() ) {
 
-        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
+        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
         
         nAtoms_ += mol->getNAtoms();
         nBonds_ += mol->getNBonds();
@@ -157,6 +151,8 @@ bool SimInfo::addMolecule(Molecule* mol) {
         nCutoffGroups_ += mol->getNCutoffGroups();
         nConstraints_ += mol->getNConstraints();
 
+        addExcludePairs(mol);
+        
         return true;
     } else {
         return false;
@@ -180,6 +176,7 @@ bool SimInfo::removeMolecule(Molecule* mol) {
         nCutoffGroups_ -= mol->getNCutoffGroups();
         nConstraints_ -= mol->getNConstraints();
 
+        removeExcludePairs(mol);
         molecules_.erase(mol->getGlobalIndex());
 
         delete mol;
@@ -462,9 +459,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 +540,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"
@@ -613,7 +610,7 @@ void SimInfo::setupFortranSim() {
     //molMembershipArray is filled by SimCreator    
     std::vector<int> molMembershipArray(nGlobalAtoms_);
     for (int i = 0; i < nGlobalAtoms_; i++) {
-        molMembershipArray.push_back(globalMolMembership_[i] + 1);
+        molMembershipArray[i] = globalMolMembership_[i] + 1;
     }
     
     //setup fortran simulation
@@ -681,11 +678,11 @@ void SimInfo::setupFortranParallel() {
     parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
     parallelData.nGroupsLocal = getNCutoffGroups();
     parallelData.myNode = worldRank;
-    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
+    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
 
     //pass mpiSimData struct and index arrays to fortran
-    setFsimParallel(parallelData, &(parallelData->nAtomsLocal),
-                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
+    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
+                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
                     &localToGlobalCutoffGroupIndex[0], &isError);
 
     if (isError) {
@@ -732,7 +729,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 +738,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 +750,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 +825,48 @@ std::ostream& operator <<(ostream& o, SimInfo& info) {
     
 }
 
-std::ostream& operator <<(ostream& o, SimInfo& info) {
+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
+