--- trunk/src/brains/SimInfo.cpp	2012/08/22 02:28:28	1782
+++ trunk/src/brains/SimInfo.cpp	2013/08/19 13:12:00	1929
@@ -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).
  */
@@ -91,6 +91,13 @@ 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);
@@ -268,8 +275,9 @@ namespace OpenMD {
     ndf_local -= nConstraints_;
 
 #ifdef IS_MPI
-    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-    MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM);
+    MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1,
+                              MPI::INT, MPI::SUM);
 #else
     ndf_ = ndf_local;
     nGlobalFluctuatingCharges_ = nfq_local;
@@ -283,7 +291,7 @@ namespace OpenMD {
 
   int SimInfo::getFdf() {
 #ifdef IS_MPI
-    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+    MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM);
 #else
     fdf_ = fdf_local;
 #endif
@@ -339,7 +347,7 @@ namespace OpenMD {
     }
     
 #ifdef IS_MPI
-    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+    MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM);
 #else
     ndfRaw_ = ndfRaw_local;
 #endif
@@ -352,7 +360,8 @@ namespace OpenMD {
 
 
 #ifdef IS_MPI
-    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+    MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1, 
+                              MPI::INT, MPI::SUM);
 #else
     ndfTrans_ = ndfTrans_local;
 #endif
@@ -779,7 +788,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
@@ -877,7 +902,6 @@ namespace OpenMD {
 
 
   void SimInfo::prepareTopology() {
-    int nExclude, nOneTwo, nOneThree, nOneFour;
 
     //calculate mass ratio of cutoff group
     SimInfo::MoleculeIterator mi;
@@ -914,28 +938,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 
-
-    nExclude = excludedInteractions_.getSize();
-    nOneTwo = oneTwoInteractions_.getSize();
-    nOneThree = oneThreeInteractions_.getSize();
-    nOneFour = oneFourInteractions_.getSize();
-
-    int* excludeList = excludedInteractions_.getPairList();
-    int* oneTwoList = oneTwoInteractions_.getPairList();
-    int* oneThreeList = oneThreeInteractions_.getPairList();
-    int* oneFourList = oneFourInteractions_.getPairList();
-
+       
     topologyDone_ = true;
   }
 
@@ -1007,7 +1024,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");
@@ -1025,8 +1042,8 @@ namespace OpenMD {
   int SimInfo::getNGlobalConstraints() {
     int nGlobalConstraints;
 #ifdef IS_MPI
-    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
-                  MPI_COMM_WORLD);    
+    MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, 
+                              MPI::INT, MPI::SUM);
 #else
     nGlobalConstraints =  nConstraints_;
 #endif