--- trunk/src/brains/SimInfo.cpp	2012/09/10 18:38:44	1796
+++ trunk/src/brains/SimInfo.cpp	2013/12/05 18:19:26	1953
@@ -35,7 +35,7 @@
  *                                                                      
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
- * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
  * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
@@ -47,6 +47,9 @@
  * @version 1.0
  */
 
+#ifdef IS_MPI
+#include <mpi.h>
+#endif
 #include <algorithm>
 #include <set>
 #include <map>
@@ -61,9 +64,6 @@
 #include "io/ForceFieldOptions.hpp"
 #include "brains/ForceField.hpp"
 #include "nonbonded/SwitchingFunction.hpp"
-#ifdef IS_MPI
-#include <mpi.h>
-#endif
 
 using namespace std;
 namespace OpenMD {
@@ -72,10 +72,12 @@ namespace OpenMD {
     forceField_(ff), simParams_(simParams), 
     ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
     nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), 
-    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0),
-    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0), 
-    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), 
-    nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), 
+    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), 
+    nGlobalFluctuatingCharges_(0), nGlobalBonds_(0), nGlobalBends_(0), 
+    nGlobalTorsions_(0), nGlobalInversions_(0), nAtoms_(0), nBonds_(0), 
+    nBends_(0), nTorsions_(0), nInversions_(0), nRigidBodies_(0), 
+    nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), 
+    nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), 
     calcBoxDipole_(false), useAtomicVirial_(true) {    
     
     MoleculeStamp* molStamp;
@@ -91,12 +93,23 @@ namespace OpenMD {
     for (vector<Component*>::iterator i = components.begin(); 
          i !=components.end(); ++i) {
       molStamp = (*i)->getMoleculeStamp();
+      if ( (*i)->haveRegion() ) {        
+        molStamp->setRegion( (*i)->getRegion() );
+      } else {
+        // set the region to a disallowed value:
+        molStamp->setRegion( -1 );
+      }
+
       nMolWithSameStamp = (*i)->getNMol();
       
       addMoleculeStamp(molStamp, nMolWithSameStamp);
       
       //calculate atoms in molecules
-      nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;   
+      nGlobalAtoms_ += molStamp->getNAtoms() * nMolWithSameStamp;
+      nGlobalBonds_ += molStamp->getNBonds() * nMolWithSameStamp;
+      nGlobalBends_ += molStamp->getNBends() * nMolWithSameStamp;
+      nGlobalTorsions_ += molStamp->getNTorsions() * nMolWithSameStamp;
+      nGlobalInversions_ += molStamp->getNInversions() * nMolWithSameStamp;
       
       //calculate atoms in cutoff groups
       int nAtomsInGroups = 0;
@@ -411,6 +424,7 @@ namespace OpenMD {
         atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
       }
     }  
+
            
     for (bond= mol->beginBond(bondIter); bond != NULL; 
          bond = mol->nextBond(bondIter)) {
@@ -781,7 +795,23 @@ namespace OpenMD {
 
     return atomTypes;        
   }
+
+
+  int getGlobalCountOfType(AtomType* atype) {
+    /*
+    set<AtomType*> atypes = getSimulatedAtomTypes();
+    map<AtomType*, int> counts_;
 
+    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
+      for(atom = mol->beginAtom(ai); atom != NULL;
+          atom = mol->nextAtom(ai)) {
+	atom->getAtomType();
+      }      
+    }    
+    */
+    return 0;
+  }
+
   void SimInfo::setupSimVariables() {
     useAtomicVirial_ = simParams_->getUseAtomicVirial();
     // we only call setAccumulateBoxDipole if the accumulateBoxDipole
@@ -915,23 +945,21 @@ namespace OpenMD {
       }       
     }
 
-    // Build the identArray_
+    // Build the identArray_ and regions_
 
     identArray_.clear();
-    identArray_.reserve(getNAtoms());    
-    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
+    identArray_.reserve(getNAtoms());   
+    regions_.clear();
+    regions_.reserve(getNAtoms());
+ 
+    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
+      int reg = mol->getRegion();      
       for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
 	identArray_.push_back(atom->getIdent());
+        regions_.push_back(reg);
       }
     }    
-    
-    //scan topology 
-
-    int* excludeList = excludedInteractions_.getPairList();
-    int* oneTwoList = oneTwoInteractions_.getPairList();
-    int* oneThreeList = oneThreeInteractions_.getPairList();
-    int* oneFourList = oneFourInteractions_.getPairList();
-
+       
     topologyDone_ = true;
   }
 
@@ -966,33 +994,55 @@ namespace OpenMD {
     delete sman_;
     sman_ = sman;
 
-    Molecule* mol;
-    RigidBody* rb;
-    Atom* atom;
-    CutoffGroup* cg;
     SimInfo::MoleculeIterator mi;
+    Molecule::AtomIterator ai;
     Molecule::RigidBodyIterator rbIter;
-    Molecule::AtomIterator atomIter;
     Molecule::CutoffGroupIterator cgIter;
+    Molecule::BondIterator bondIter;
+    Molecule::BendIterator bendIter;
+    Molecule::TorsionIterator torsionIter;
+    Molecule::InversionIterator inversionIter;
  
+    Molecule* mol;
+    Atom* atom;
+    RigidBody* rb;
+    CutoffGroup* cg;
+    Bond* bond;
+    Bend* bend;
+    Torsion* torsion;
+    Inversion* inversion;    
+
     for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
         
-      for (atom = mol->beginAtom(atomIter); atom != NULL; 
-           atom = mol->nextAtom(atomIter)) {
+      for (atom = mol->beginAtom(ai); atom != NULL; 
+           atom = mol->nextAtom(ai)) {
 	atom->setSnapshotManager(sman_);
-      }
-        
+      }        
       for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
            rb = mol->nextRigidBody(rbIter)) {
 	rb->setSnapshotManager(sman_);
       }
-
       for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; 
            cg = mol->nextCutoffGroup(cgIter)) {
 	cg->setSnapshotManager(sman_);
       }
-    }    
-    
+      for (bond = mol->beginBond(bondIter); bond != NULL; 
+           bond = mol->nextBond(bondIter)) {
+        bond->setSnapshotManager(sman_);
+      }
+      for (bend = mol->beginBend(bendIter); bend != NULL; 
+           bend = mol->nextBend(bendIter)) {
+        bend->setSnapshotManager(sman_);
+      }
+      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
+           torsion = mol->nextTorsion(torsionIter)) {
+        torsion->setSnapshotManager(sman_);
+      }
+      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
+           inversion = mol->nextInversion(inversionIter)) {
+        inversion->setSnapshotManager(sman_);
+      }
+    }
   }
 
 
@@ -1003,7 +1053,7 @@ namespace OpenMD {
    
   
   StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
-    if (index >= IOIndexToIntegrableObject.size()) {
+    if (index >= int(IOIndexToIntegrableObject.size())) {
       sprintf(painCave.errMsg,
               "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n"
               "\tindex exceeds number of known objects!\n");