--- trunk/src/brains/SimInfo.cpp	2008/09/10 19:51:45	1290
+++ branches/development/src/brains/SimInfo.cpp	2012/09/13 14:10:11	1798
@@ -6,19 +6,10 @@
  * redistribute this software in source and binary code form, provided
  * that the following conditions are met:
  *
- * 1. Acknowledgement of the program authors must be made in any
- *    publication of scientific results based in part on use of the
- *    program.  An acceptable form of acknowledgement is citation of
- *    the article in which the program was described (Matthew
- *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
- *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
- *    Parallel Simulation Engine for Molecular Dynamics,"
- *    J. Comput. Chem. 26, pp. 252-271 (2005))
- *
- * 2. Redistributions of source code must retain the above copyright
+ * 1. Redistributions of source code must retain the above copyright
  *    notice, this list of conditions and the following disclaimer.
  *
- * 3. Redistributions in binary form must reproduce the above copyright
+ * 2. Redistributions in binary form must reproduce the above copyright
  *    notice, this list of conditions and the following disclaimer in the
  *    documentation and/or other materials provided with the
  *    distribution.
@@ -37,6 +28,16 @@
  * arising out of the use of or inability to use software, even if the
  * University of Notre Dame has been advised of the possibility of
  * such damages.
+ *
+ * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
+ * research, please cite the appropriate papers when you publish your
+ * work.  Good starting points are:
+ *                                                                      
+ * [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).          
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
 /**
@@ -54,115 +55,98 @@
 #include "math/Vector3.hpp"
 #include "primitives/Molecule.hpp"
 #include "primitives/StuntDouble.hpp"
-#include "UseTheForce/fCutoffPolicy.h"
-#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
-#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
-#include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
-#include "UseTheForce/doForces_interface.h"
-#include "UseTheForce/DarkSide/neighborLists_interface.h"
-#include "UseTheForce/DarkSide/electrostatic_interface.h"
-#include "UseTheForce/DarkSide/switcheroo_interface.h"
 #include "utils/MemoryUtils.hpp"
 #include "utils/simError.h"
 #include "selection/SelectionManager.hpp"
 #include "io/ForceFieldOptions.hpp"
-#include "UseTheForce/ForceField.hpp"
-
-
+#include "brains/ForceField.hpp"
+#include "nonbonded/SwitchingFunction.hpp"
 #ifdef IS_MPI
-#include "UseTheForce/mpiComponentPlan.h"
-#include "UseTheForce/DarkSide/simParallel_interface.h"
-#endif 
+#include <mpi.h>
+#endif
 
-namespace oopse {
-  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
-    std::map<int, std::set<int> >::iterator i = container.find(index);
-    std::set<int> result;
-    if (i != container.end()) {
-        result = i->second;
-    }
-
-    return result;
-  }
+using namespace std;
+namespace OpenMD {
   
   SimInfo::SimInfo(ForceField* ff, Globals* simParams) : 
     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),
+    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0),
     nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0), 
     nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), 
-    nConstraints_(0), sman_(NULL), fortranInitialized_(false), 
-    calcBoxDipole_(false), useAtomicVirial_(true) {
-
-
-      MoleculeStamp* molStamp;
-      int nMolWithSameStamp;
-      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 nRigidAtoms = 0;
-
-      std::vector<Component*> components = simParams->getComponents();
+    nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), 
+    calcBoxDipole_(false), useAtomicVirial_(true) {    
+    
+    MoleculeStamp* molStamp;
+    int nMolWithSameStamp;
+    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 nRigidAtoms = 0;
+    
+    vector<Component*> components = simParams->getComponents();
+    
+    for (vector<Component*>::iterator i = components.begin(); 
+         i !=components.end(); ++i) {
+      molStamp = (*i)->getMoleculeStamp();
+      nMolWithSameStamp = (*i)->getNMol();
       
-      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
-        molStamp = (*i)->getMoleculeStamp();
-        nMolWithSameStamp = (*i)->getNMol();
-        
-        addMoleculeStamp(molStamp, nMolWithSameStamp);
-
-        //calculate atoms in molecules
-        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;   
-
-        //calculate atoms in cutoff groups
-        int nAtomsInGroups = 0;
-        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
-        
-        for (int j=0; j < nCutoffGroupsInStamp; j++) {
-	  cgStamp = molStamp->getCutoffGroupStamp(j);
-	  nAtomsInGroups += cgStamp->getNMembers();
-        }
-
-        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
-
-        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
-
-        //calculate atoms in rigid bodies
-        int nAtomsInRigidBodies = 0;
-        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
-        
-        for (int j=0; j < nRigidBodiesInStamp; j++) {
-	  rbStamp = molStamp->getRigidBodyStamp(j);
-	  nAtomsInRigidBodies += rbStamp->getNMembers();
-        }
-
-        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
-        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
-        
+      addMoleculeStamp(molStamp, nMolWithSameStamp);
+      
+      //calculate atoms in molecules
+      nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;   
+      
+      //calculate atoms in cutoff groups
+      int nAtomsInGroups = 0;
+      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
+      
+      for (int j=0; j < nCutoffGroupsInStamp; j++) {
+        cgStamp = molStamp->getCutoffGroupStamp(j);
+        nAtomsInGroups += cgStamp->getNMembers();
       }
-
-      //every free atom (atom does not belong to cutoff groups) is a cutoff 
-      //group therefore the total number of cutoff groups in the system is 
-      //equal to the total number of atoms minus number of atoms belong to 
-      //cutoff group defined in meta-data 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 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 
-	                                        + nGlobalRigidBodies_;
-  
-      nGlobalMols_ = molStampIds_.size();
-      molToProcMap_.resize(nGlobalMols_);
+      
+      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
+      
+      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
+      
+      //calculate atoms in rigid bodies
+      int nAtomsInRigidBodies = 0;
+      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
+      
+      for (int j=0; j < nRigidBodiesInStamp; j++) {
+        rbStamp = molStamp->getRigidBodyStamp(j);
+        nAtomsInRigidBodies += rbStamp->getNMembers();
+      }
+      
+      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
+      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
+      
     }
+    
+    //every free atom (atom does not belong to cutoff groups) is a cutoff 
+    //group therefore the total number of cutoff groups in the system is 
+    //equal to the total number of atoms minus number of atoms belong to 
+    //cutoff group defined in meta-data 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 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 
+      + nGlobalRigidBodies_;
+    
+    nGlobalMols_ = molStampIds_.size();
+    molToProcMap_.resize(nGlobalMols_);
+  }
+  
   SimInfo::~SimInfo() {
-    std::map<int, Molecule*>::iterator i;
+    map<int, Molecule*>::iterator i;
     for (i = molecules_.begin(); i != molecules_.end(); ++i) {
       delete i->second;
     }
@@ -173,25 +157,15 @@ namespace oopse {
     delete forceField_;
   }
 
-  int SimInfo::getNGlobalConstraints() {
-    int nGlobalConstraints;
-#ifdef IS_MPI
-    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
-                  MPI_COMM_WORLD);    
-#else
-    nGlobalConstraints =  nConstraints_;
-#endif
-    return nGlobalConstraints;
-  }
 
   bool SimInfo::addMolecule(Molecule* mol) {
     MoleculeIterator i;
-
+    
     i = molecules_.find(mol->getGlobalIndex());
     if (i == molecules_.end() ) {
-
-      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
-        
+      
+      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
+      
       nAtoms_ += mol->getNAtoms();
       nBonds_ += mol->getNBonds();
       nBends_ += mol->getNBends();
@@ -201,15 +175,15 @@ namespace oopse {
       nIntegrableObjects_ += mol->getNIntegrableObjects();
       nCutoffGroups_ += mol->getNCutoffGroups();
       nConstraints_ += mol->getNConstraintPairs();
-
+      
       addInteractionPairs(mol);
-  
+      
       return true;
     } else {
       return false;
     }
   }
-
+  
   bool SimInfo::removeMolecule(Molecule* mol) {
     MoleculeIterator i;
     i = molecules_.find(mol->getGlobalIndex());
@@ -237,8 +211,6 @@ namespace oopse {
     } else {
       return false;
     }
-
-
   }    
 
         
@@ -254,38 +226,54 @@ namespace oopse {
 
 
   void SimInfo::calcNdf() {
-    int ndf_local;
+    int ndf_local, nfq_local;
     MoleculeIterator i;
-    std::vector<StuntDouble*>::iterator j;
+    vector<StuntDouble*>::iterator j;
+    vector<Atom*>::iterator k;
+
     Molecule* mol;
-    StuntDouble* integrableObject;
+    StuntDouble* sd;
+    Atom* atom;
 
     ndf_local = 0;
+    nfq_local = 0;
     
     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; 
-	   integrableObject = mol->nextIntegrableObject(j)) {
 
+      for (sd = mol->beginIntegrableObject(j); sd != NULL; 
+	   sd = mol->nextIntegrableObject(j)) {
+
 	ndf_local += 3;
 
-	if (integrableObject->isDirectional()) {
-	  if (integrableObject->isLinear()) {
+	if (sd->isDirectional()) {
+	  if (sd->isLinear()) {
 	    ndf_local += 2;
 	  } else {
 	    ndf_local += 3;
 	  }
 	}
-            
       }
+
+      for (atom = mol->beginFluctuatingCharge(k); atom != NULL;
+           atom = mol->nextFluctuatingCharge(k)) {
+        if (atom->isFluctuatingCharge()) {
+          nfq_local++;
+        }
+      }
     }
     
+    ndfLocal_ = ndf_local;
+
     // n_constraints is local, so subtract them on each processor
     ndf_local -= nConstraints_;
 
 #ifdef IS_MPI
-    MPI_Allreduce(&ndf_local,&ndf_,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;
 #endif
 
     // nZconstraints_ is global, as are the 3 COM translations for the 
@@ -296,32 +284,52 @@ namespace oopse {
 
   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
     return fdf_;
   }
+  
+  unsigned int SimInfo::getNLocalCutoffGroups(){
+    int nLocalCutoffAtoms = 0;
+    Molecule* mol;
+    MoleculeIterator mi;
+    CutoffGroup* cg;
+    Molecule::CutoffGroupIterator ci;
     
+    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
+      
+      for (cg = mol->beginCutoffGroup(ci); cg != NULL; 
+           cg = mol->nextCutoffGroup(ci)) {
+        nLocalCutoffAtoms += cg->getNumAtom();
+        
+      }        
+    }
+    
+    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
+  }
+    
   void SimInfo::calcNdfRaw() {
     int ndfRaw_local;
 
     MoleculeIterator i;
-    std::vector<StuntDouble*>::iterator j;
+    vector<StuntDouble*>::iterator j;
     Molecule* mol;
-    StuntDouble* integrableObject;
+    StuntDouble* sd;
 
     // Raw degrees of freedom that we have to set
     ndfRaw_local = 0;
     
     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
-	   integrableObject = mol->nextIntegrableObject(j)) {
 
+      for (sd = mol->beginIntegrableObject(j); sd != NULL;
+	   sd = mol->nextIntegrableObject(j)) {
+
 	ndfRaw_local += 3;
 
-	if (integrableObject->isDirectional()) {
-	  if (integrableObject->isLinear()) {
+	if (sd->isDirectional()) {
+	  if (sd->isLinear()) {
 	    ndfRaw_local += 2;
 	  } else {
 	    ndfRaw_local += 3;
@@ -332,7 +340,7 @@ namespace oopse {
     }
     
 #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
@@ -345,7 +353,8 @@ namespace oopse {
 
 
 #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
@@ -356,10 +365,10 @@ namespace oopse {
 
   void SimInfo::addInteractionPairs(Molecule* mol) {
     ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
-    std::vector<Bond*>::iterator bondIter;
-    std::vector<Bend*>::iterator bendIter;
-    std::vector<Torsion*>::iterator torsionIter;
-    std::vector<Inversion*>::iterator inversionIter;
+    vector<Bond*>::iterator bondIter;
+    vector<Bend*>::iterator bendIter;
+    vector<Torsion*>::iterator torsionIter;
+    vector<Inversion*>::iterator inversionIter;
     Bond* bond;
     Bend* bend;
     Torsion* torsion;
@@ -377,30 +386,29 @@ namespace oopse {
     // always be excluded.  These are done at the bottom of this
     // function.
 
-    std::map<int, std::set<int> > atomGroups;
+    map<int, set<int> > atomGroups;
     Molecule::RigidBodyIterator rbIter;
     RigidBody* rb;
     Molecule::IntegrableObjectIterator ii;
-    StuntDouble* integrableObject;
+    StuntDouble* sd;
     
-    for (integrableObject = mol->beginIntegrableObject(ii); 
-         integrableObject != NULL;
-         integrableObject = mol->nextIntegrableObject(ii)) {
+    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
+         sd = mol->nextIntegrableObject(ii)) {
       
-      if (integrableObject->isRigidBody()) {
-        rb = static_cast<RigidBody*>(integrableObject);
-        std::vector<Atom*> atoms = rb->getAtoms();
-        std::set<int> rigidAtoms;
+      if (sd->isRigidBody()) {
+        rb = static_cast<RigidBody*>(sd);
+        vector<Atom*> atoms = rb->getAtoms();
+        set<int> rigidAtoms;
         for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
           rigidAtoms.insert(atoms[i]->getGlobalIndex());
         }
         for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
-          atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
+          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
         }      
       } else {
-        std::set<int> oneAtomSet;
-        oneAtomSet.insert(integrableObject->getGlobalIndex());
-        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
+        set<int> oneAtomSet;
+        oneAtomSet.insert(sd->getGlobalIndex());
+        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
       }
     }  
            
@@ -503,7 +511,7 @@ namespace oopse {
 
     for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
          rb = mol->nextRigidBody(rbIter)) {
-      std::vector<Atom*> atoms = rb->getAtoms();
+      vector<Atom*> atoms = rb->getAtoms();
       for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
 	for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
 	  a = atoms[i]->getGlobalIndex();
@@ -517,10 +525,10 @@ namespace oopse {
 
   void SimInfo::removeInteractionPairs(Molecule* mol) {
     ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
-    std::vector<Bond*>::iterator bondIter;
-    std::vector<Bend*>::iterator bendIter;
-    std::vector<Torsion*>::iterator torsionIter;
-    std::vector<Inversion*>::iterator inversionIter;
+    vector<Bond*>::iterator bondIter;
+    vector<Bend*>::iterator bendIter;
+    vector<Torsion*>::iterator torsionIter;
+    vector<Inversion*>::iterator inversionIter;
     Bond* bond;
     Bend* bend;
     Torsion* torsion;
@@ -530,30 +538,29 @@ namespace oopse {
     int c;
     int d;
 
-    std::map<int, std::set<int> > atomGroups;
+    map<int, set<int> > atomGroups;
     Molecule::RigidBodyIterator rbIter;
     RigidBody* rb;
     Molecule::IntegrableObjectIterator ii;
-    StuntDouble* integrableObject;
+    StuntDouble* sd;
     
-    for (integrableObject = mol->beginIntegrableObject(ii); 
-         integrableObject != NULL;
-         integrableObject = mol->nextIntegrableObject(ii)) {
+    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
+         sd = mol->nextIntegrableObject(ii)) {
       
-      if (integrableObject->isRigidBody()) {
-        rb = static_cast<RigidBody*>(integrableObject);
-        std::vector<Atom*> atoms = rb->getAtoms();
-        std::set<int> rigidAtoms;
+      if (sd->isRigidBody()) {
+        rb = static_cast<RigidBody*>(sd);
+        vector<Atom*> atoms = rb->getAtoms();
+        set<int> rigidAtoms;
         for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
           rigidAtoms.insert(atoms[i]->getGlobalIndex());
         }
         for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
-          atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
+          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
         }      
       } else {
-        std::set<int> oneAtomSet;
-        oneAtomSet.insert(integrableObject->getGlobalIndex());
-        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
+        set<int> oneAtomSet;
+        oneAtomSet.insert(sd->getGlobalIndex());
+        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
       }
     }  
 
@@ -656,7 +663,7 @@ namespace oopse {
 
     for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
          rb = mol->nextRigidBody(rbIter)) {
-      std::vector<Atom*> atoms = rb->getAtoms();
+      vector<Atom*> atoms = rb->getAtoms();
       for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
 	for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
 	  a = atoms[i]->getGlobalIndex();
@@ -679,226 +686,202 @@ namespace oopse {
     molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
   }
 
-  void SimInfo::update() {
 
-    setupSimType();
-
-#ifdef IS_MPI
-    setupFortranParallel();
-#endif
-
-    setupFortranSim();
-
-    //setup fortran force field
-    /** @deprecate */    
-    int isError = 0;
-    
-    setupCutoff();
-    
-    setupElectrostaticSummationMethod( isError );
-    setupSwitchingFunction();
-    setupAccumulateBoxDipole();
-
-    if(isError){
-      sprintf( painCave.errMsg,
-	       "ForceField error: There was an error initializing the forceField in fortran.\n" );
-      painCave.isFatal = 1;
-      simError();
-    }
-
+  /**
+   * update
+   *
+   *  Performs the global checks and variable settings after the
+   *  objects have been created.
+   * 
+   */
+  void SimInfo::update() {   
+    setupSimVariables();
     calcNdf();
     calcNdfRaw();
     calcNdfTrans();
-
-    fortranInitialized_ = true;
   }
-
-  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
+  
+  /**
+   * getSimulatedAtomTypes
+   *
+   * Returns an STL set of AtomType* that are actually present in this
+   * simulation.  Must query all processors to assemble this information.
+   * 
+   */
+  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
     SimInfo::MoleculeIterator mi;
     Molecule* mol;
     Molecule::AtomIterator ai;
     Atom* atom;
-    std::set<AtomType*> atomTypes;
-
+    set<AtomType*> atomTypes;
+    
     for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
-
-      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
+      for(atom = mol->beginAtom(ai); atom != NULL;
+          atom = mol->nextAtom(ai)) {
 	atomTypes.insert(atom->getAtomType());
-      }
-        
-    }
-
-    return atomTypes;        
-  }
-
-  void SimInfo::setupSimType() {
-    std::set<AtomType*>::iterator i;
-    std::set<AtomType*> atomTypes;
-    atomTypes = getUniqueAtomTypes();
+      }      
+    }    
     
-    int useLennardJones = 0;
-    int useElectrostatic = 0;
-    int useEAM = 0;
-    int useSC = 0;
-    int useCharge = 0;
-    int useDirectional = 0;
-    int useDipole = 0;
-    int useGayBerne = 0;
-    int useSticky = 0;
-    int useStickyPower = 0;
-    int useShape = 0; 
-    int useFLARB = 0; //it is not in AtomType yet
-    int useDirectionalAtom = 0;    
-    int useElectrostatics = 0;
-    //usePBC and useRF are from simParams
-    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
-    int useRF;
-    int useSF;
-    int useSP;
-    int useBoxDipole;
+#ifdef IS_MPI
 
-    std::string myMethod;
-
-    // set the useRF logical
-    useRF = 0;
-    useSF = 0;
-    useSP = 0;
-
-
-    if (simParams_->haveElectrostaticSummationMethod()) {
-      std::string myMethod = simParams_->getElectrostaticSummationMethod();
-      toUpper(myMethod);
-      if (myMethod == "REACTION_FIELD"){
-        useRF = 1;
-      } else if (myMethod == "SHIFTED_FORCE"){
-	useSF = 1;
-      } else if (myMethod == "SHIFTED_POTENTIAL"){
-	useSP = 1;
-      }
-    }
+    // loop over the found atom types on this processor, and add their
+    // numerical idents to a vector:
     
-    if (simParams_->haveAccumulateBoxDipole()) 
-      if (simParams_->getAccumulateBoxDipole())
-	useBoxDipole = 1;
+    vector<int> foundTypes;
+    set<AtomType*>::iterator i;
+    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) 
+      foundTypes.push_back( (*i)->getIdent() );
 
-    useAtomicVirial_ = simParams_->getUseAtomicVirial();
+    // count_local holds the number of found types on this processor
+    int count_local = foundTypes.size();
 
-    //loop over all of the atom types
-    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
-      useLennardJones |= (*i)->isLennardJones();
-      useElectrostatic |= (*i)->isElectrostatic();
-      useEAM |= (*i)->isEAM();
-      useSC |= (*i)->isSC();
-      useCharge |= (*i)->isCharge();
-      useDirectional |= (*i)->isDirectional();
-      useDipole |= (*i)->isDipole();
-      useGayBerne |= (*i)->isGayBerne();
-      useSticky |= (*i)->isSticky();
-      useStickyPower |= (*i)->isStickyPower();
-      useShape |= (*i)->isShape(); 
-    }
+    int nproc = MPI::COMM_WORLD.Get_size();
 
-    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
-      useDirectionalAtom = 1;
-    }
+    // we need arrays to hold the counts and displacement vectors for
+    // all processors
+    vector<int> counts(nproc, 0);
+    vector<int> disps(nproc, 0);
 
-    if (useCharge || useDipole) {
-      useElectrostatics = 1;
+    // fill the counts array
+    MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
+                              1, MPI::INT);
+  
+    // use the processor counts to compute the displacement array
+    disps[0] = 0;    
+    int totalCount = counts[0];
+    for (int iproc = 1; iproc < nproc; iproc++) {
+      disps[iproc] = disps[iproc-1] + counts[iproc-1];
+      totalCount += counts[iproc];
     }
 
-#ifdef IS_MPI    
-    int temp;
+    // we need a (possibly redundant) set of all found types:
+    vector<int> ftGlobal(totalCount);
+    
+    // now spray out the foundTypes to all the other processors:    
+    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, 
+                               &ftGlobal[0], &counts[0], &disps[0], 
+                               MPI::INT);
 
-    temp = usePBC;
-    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    vector<int>::iterator j;
 
-    temp = useDirectionalAtom;
-    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    // foundIdents is a stl set, so inserting an already found ident
+    // will have no effect.
+    set<int> foundIdents;
 
-    temp = useLennardJones;
-    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
+      foundIdents.insert((*j));
+    
+    // now iterate over the foundIdents and get the actual atom types 
+    // that correspond to these:
+    set<int>::iterator it;
+    for (it = foundIdents.begin(); it != foundIdents.end(); ++it) 
+      atomTypes.insert( forceField_->getAtomType((*it)) );
+ 
+#endif
 
-    temp = useElectrostatics;
-    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    return atomTypes;        
+  }
 
-    temp = useCharge;
-    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
-
-    temp = useDipole;
-    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
-
-    temp = useSticky;
-    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
-
-    temp = useStickyPower;
-    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+  void SimInfo::setupSimVariables() {
+    useAtomicVirial_ = simParams_->getUseAtomicVirial();
+    // we only call setAccumulateBoxDipole if the accumulateBoxDipole
+    // parameter is true
+    calcBoxDipole_ = false;
+    if ( simParams_->haveAccumulateBoxDipole() ) 
+      if ( simParams_->getAccumulateBoxDipole() ) {
+	calcBoxDipole_ = true;       
+      }
     
-    temp = useGayBerne;
-    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    set<AtomType*>::iterator i;
+    set<AtomType*> atomTypes;
+    atomTypes = getSimulatedAtomTypes();    
+    bool usesElectrostatic = false;
+    bool usesMetallic = false;
+    bool usesDirectional = false;
+    bool usesFluctuatingCharges =  false;
+    //loop over all of the atom types
+    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
+      usesElectrostatic |= (*i)->isElectrostatic();
+      usesMetallic |= (*i)->isMetal();
+      usesDirectional |= (*i)->isDirectional();
+      usesFluctuatingCharges |= (*i)->isFluctuatingCharge();
+    }
 
-    temp = useEAM;
-    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
-
-    temp = useSC;
-    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
+#ifdef IS_MPI
+    bool temp;
+    temp = usesDirectional;
+    MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, 
+                              MPI::LOR);
+        
+    temp = usesMetallic;
+    MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, 
+                              MPI::LOR);
     
-    temp = useShape;
-    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);   
+    temp = usesElectrostatic;
+    MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, 
+                              MPI::LOR);
 
-    temp = useFLARB;
-    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    temp = usesFluctuatingCharges;
+    MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, 
+                              MPI::LOR);
+#else
 
-    temp = useRF;
-    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    usesDirectionalAtoms_ = usesDirectional;
+    usesMetallicAtoms_ = usesMetallic;
+    usesElectrostaticAtoms_ = usesElectrostatic;
+    usesFluctuatingCharges_ = usesFluctuatingCharges;
 
-    temp = useSF;
-    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);   
+#endif
+    
+    requiresPrepair_ = usesMetallicAtoms_ ? true : false; 
+    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
+    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
+  }
 
-    temp = useSP;
-    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
 
-    temp = useBoxDipole;
-    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); 
+  vector<int> SimInfo::getGlobalAtomIndices() {
+    SimInfo::MoleculeIterator mi;
+    Molecule* mol;
+    Molecule::AtomIterator ai;
+    Atom* atom;
 
-    temp = useAtomicVirial_;
-    MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); 
+    vector<int> GlobalAtomIndices(getNAtoms(), 0);
+    
+    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
+      
+      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
+	GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
+      }
+    }
+    return GlobalAtomIndices;
+  }
 
-#endif
 
-    fInfo_.SIM_uses_PBC = usePBC;    
-    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
-    fInfo_.SIM_uses_LennardJones = useLennardJones;
-    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
-    fInfo_.SIM_uses_Charges = useCharge;
-    fInfo_.SIM_uses_Dipoles = useDipole;
-    fInfo_.SIM_uses_Sticky = useSticky;
-    fInfo_.SIM_uses_StickyPower = useStickyPower;
-    fInfo_.SIM_uses_GayBerne = useGayBerne;
-    fInfo_.SIM_uses_EAM = useEAM;
-    fInfo_.SIM_uses_SC = useSC;
-    fInfo_.SIM_uses_Shapes = useShape;
-    fInfo_.SIM_uses_FLARB = useFLARB;
-    fInfo_.SIM_uses_RF = useRF;
-    fInfo_.SIM_uses_SF = useSF;
-    fInfo_.SIM_uses_SP = useSP;
-    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
-    fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
-  }
+  vector<int> SimInfo::getGlobalGroupIndices() {
+    SimInfo::MoleculeIterator mi;
+    Molecule* mol;
+    Molecule::CutoffGroupIterator ci;
+    CutoffGroup* cg;
 
-  void SimInfo::setupFortranSim() {
-    int isError;
-    int nExclude, nOneTwo, nOneThree, nOneFour;
-    std::vector<int> fortranGlobalGroupMembership;
+    vector<int> GlobalGroupIndices;
     
-    isError = 0;
-
-    //globalGroupMembership_ is filled by SimCreator    
-    for (int i = 0; i < nGlobalAtoms_; i++) {
-      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
+    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
+      
+      //local index of cutoff group is trivial, it only depends on the
+      //order of travesing
+      for (cg = mol->beginCutoffGroup(ci); cg != NULL; 
+           cg = mol->nextCutoffGroup(ci)) {
+	GlobalGroupIndices.push_back(cg->getGlobalIndex());
+      }        
     }
+    return GlobalGroupIndices;
+  }
 
+
+  void SimInfo::prepareTopology() {
+    int nExclude, nOneTwo, nOneThree, nOneFour;
+
     //calculate mass ratio of cutoff group
-    std::vector<RealType> mfact;
     SimInfo::MoleculeIterator mi;
     Molecule* mol;
     Molecule::CutoffGroupIterator ci;
@@ -907,452 +890,78 @@ namespace oopse {
     Atom* atom;
     RealType totalMass;
 
-    //to avoid memory reallocation, reserve enough space for mfact
-    mfact.reserve(getNCutoffGroups());
+    /**
+     * The mass factor is the relative mass of an atom to the total
+     * mass of the cutoff group it belongs to.  By default, all atoms
+     * are their own cutoff groups, and therefore have mass factors of
+     * 1.  We need some special handling for massless atoms, which
+     * will be treated as carrying the entire mass of the cutoff
+     * group.
+     */
+    massFactors_.clear();
+    massFactors_.resize(getNAtoms(), 1.0);
     
     for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
-      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
+      for (cg = mol->beginCutoffGroup(ci); cg != NULL; 
+           cg = mol->nextCutoffGroup(ci)) {
 
 	totalMass = cg->getMass();
 	for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
 	  // Check for massless groups - set mfact to 1 if true
-	  if (totalMass != 0)
-	    mfact.push_back(atom->getMass()/totalMass);
+	  if (totalMass != 0) 
+            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
 	  else
-	    mfact.push_back( 1.0 );
+	    massFactors_[atom->getLocalIndex()] = 1.0;
 	}
       }       
     }
 
-    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
-    std::vector<int> identArray;
+    // Build the identArray_
 
-    //to avoid memory reallocation, reserve enough space identArray
-    identArray.reserve(getNAtoms());
-    
+    identArray_.clear();
+    identArray_.reserve(getNAtoms());    
     for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
       for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
-	identArray.push_back(atom->getIdent());
+	identArray_.push_back(atom->getIdent());
       }
     }    
-
-    //fill molMembershipArray
-    //molMembershipArray is filled by SimCreator    
-    std::vector<int> molMembershipArray(nGlobalAtoms_);
-    for (int i = 0; i < nGlobalAtoms_; i++) {
-      molMembershipArray[i] = globalMolMembership_[i] + 1;
-    }
     
-    //setup fortran simulation
+    //scan topology 
 
     nExclude = excludedInteractions_.getSize();
     nOneTwo = oneTwoInteractions_.getSize();
     nOneThree = oneThreeInteractions_.getSize();
     nOneFour = oneFourInteractions_.getSize();
 
-    std::cerr << "excludedInteractions contains: " << excludedInteractions_.getSize() << " pairs \n";
-    std::cerr << "oneTwoInteractions contains: " << oneTwoInteractions_.getSize() << " pairs \n";
-    std::cerr << "oneThreeInteractions contains: " << oneThreeInteractions_.getSize() << " pairs \n";
-    std::cerr << "oneFourInteractions contains: " << oneFourInteractions_.getSize() << " pairs \n";
-
     int* excludeList = excludedInteractions_.getPairList();
     int* oneTwoList = oneTwoInteractions_.getPairList();
     int* oneThreeList = oneThreeInteractions_.getPairList();
     int* oneFourList = oneFourInteractions_.getPairList();
 
-    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], 
-                   &nExclude, excludeList, 
-                   &nOneTwo, oneTwoList,
-                   &nOneThree, oneThreeList,
-                   &nOneFour, oneFourList,
-                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_, 
-                   &fortranGlobalGroupMembership[0], &isError); 
-    
-    if( isError ){
-      
-      sprintf( painCave.errMsg,
-	       "There was an error setting the simulation information in fortran.\n" );
-      painCave.isFatal = 1;
-      painCave.severity = OOPSE_ERROR;
-      simError();
-    }
-    
-    
-    sprintf( checkPointMsg,
-	     "succesfully sent the simulation information to fortran.\n");
-    
-    errorCheckPoint();
-    
-    // Setup number of neighbors in neighbor list if present
-    if (simParams_->haveNeighborListNeighbors()) {
-      int nlistNeighbors = simParams_->getNeighborListNeighbors();
-      setNeighbors(&nlistNeighbors);
-    }
-   
+    topologyDone_ = true;
+  }
 
+  void SimInfo::addProperty(GenericData* genData) {
+    properties_.addProperty(genData);  
   }
 
+  void SimInfo::removeProperty(const string& propName) {
+    properties_.removeProperty(propName);  
+  }
 
-  void SimInfo::setupFortranParallel() {
-#ifdef IS_MPI    
-    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
-    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
-    std::vector<int> localToGlobalCutoffGroupIndex;
-    SimInfo::MoleculeIterator mi;
-    Molecule::AtomIterator ai;
-    Molecule::CutoffGroupIterator ci;
-    Molecule* mol;
-    Atom* atom;
-    CutoffGroup* cg;
-    mpiSimData parallelData;
-    int isError;
-
-    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
-
-      //local index(index in DataStorge) of atom is important
-      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
-	localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
-      }
-
-      //local index of cutoff group is trivial, it only depends on the order of travesing
-      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
-	localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
-      }        
-        
-    }
-
-    //fill up mpiSimData struct
-    parallelData.nMolGlobal = getNGlobalMolecules();
-    parallelData.nMolLocal = getNMolecules();
-    parallelData.nAtomsGlobal = getNGlobalAtoms();
-    parallelData.nAtomsLocal = getNAtoms();
-    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
-    parallelData.nGroupsLocal = getNCutoffGroups();
-    parallelData.myNode = worldRank;
-    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
-
-    //pass mpiSimData struct and index arrays to fortran
-    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
-                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
-                    &localToGlobalCutoffGroupIndex[0], &isError);
-
-    if (isError) {
-      sprintf(painCave.errMsg,
-	      "mpiRefresh errror: fortran didn't like something we gave it.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    sprintf(checkPointMsg, " mpiRefresh successful.\n");
-    errorCheckPoint();
-
-#endif
+  void SimInfo::clearProperties() {
+    properties_.clearProperties(); 
   }
 
-  void SimInfo::setupCutoff() {           
-    
-    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
-
-    // Check the cutoff policy
-    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
-
-    // Set LJ shifting bools to false
-    ljsp_ = false;
-    ljsf_ = false;
-
-    std::string myPolicy;
-    if (forceFieldOptions_.haveCutoffPolicy()){
-      myPolicy = forceFieldOptions_.getCutoffPolicy();
-    }else if (simParams_->haveCutoffPolicy()) {
-      myPolicy = simParams_->getCutoffPolicy();
-    }
-
-    if (!myPolicy.empty()){
-      toUpper(myPolicy);
-      if (myPolicy == "MIX") {
-        cp = MIX_CUTOFF_POLICY;
-      } else {
-        if (myPolicy == "MAX") {
-          cp = MAX_CUTOFF_POLICY;
-        } else {
-          if (myPolicy == "TRADITIONAL") {            
-            cp = TRADITIONAL_CUTOFF_POLICY;
-          } else {
-            // throw error        
-            sprintf( painCave.errMsg,
-                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
-            painCave.isFatal = 1;
-            simError();
-          }     
-        }           
-      }
-    }           
-    notifyFortranCutoffPolicy(&cp);
-
-    // Check the Skin Thickness for neighborlists
-    RealType skin;
-    if (simParams_->haveSkinThickness()) {
-      skin = simParams_->getSkinThickness();
-      notifyFortranSkinThickness(&skin);
-    }            
-        
-    // Check if the cutoff was set explicitly:
-    if (simParams_->haveCutoffRadius()) {
-      rcut_ = simParams_->getCutoffRadius();
-      if (simParams_->haveSwitchingRadius()) {
-	rsw_  = simParams_->getSwitchingRadius();
-      } else {
-	if (fInfo_.SIM_uses_Charges | 
-	    fInfo_.SIM_uses_Dipoles | 
-	    fInfo_.SIM_uses_RF) {
-	  
-	  rsw_ = 0.85 * rcut_;
-	  sprintf(painCave.errMsg,
-		  "SimCreator Warning: No value was set for the switchingRadius.\n"
-		  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
-		  "\tswitchingRadius = %f. for this simulation\n", rsw_);
-        painCave.isFatal = 0;
-	simError();
-	} else {
-	  rsw_ = rcut_;
-	  sprintf(painCave.errMsg,
-		  "SimCreator Warning: No value was set for the switchingRadius.\n"
-		  "\tOOPSE will use the same value as the cutoffRadius.\n"
-		  "\tswitchingRadius = %f. for this simulation\n", rsw_);
-	  painCave.isFatal = 0;
-	  simError();
-	}
-      }
-
-      if (simParams_->haveElectrostaticSummationMethod()) {
-	std::string myMethod = simParams_->getElectrostaticSummationMethod();
-	toUpper(myMethod);
-        
-	if (myMethod == "SHIFTED_POTENTIAL") {
-	  ljsp_ = true;
-	} else if (myMethod == "SHIFTED_FORCE") {
-	  ljsf_ = true;
-	}
-      }
-      notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
+  vector<string> SimInfo::getPropertyNames() {
+    return properties_.getPropertyNames();  
+  }
       
-    } else {
-      
-      // For electrostatic atoms, we'll assume a large safe value:
-      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
-        sprintf(painCave.errMsg,
-                "SimCreator Warning: No value was set for the cutoffRadius.\n"
-                "\tOOPSE will use a default value of 15.0 angstroms"
-                "\tfor the cutoffRadius.\n");
-        painCave.isFatal = 0;
-	simError();
-	rcut_ = 15.0;
-      
-        if (simParams_->haveElectrostaticSummationMethod()) {
-          std::string myMethod = simParams_->getElectrostaticSummationMethod();
-          toUpper(myMethod);
-      
-      // For the time being, we're tethering the LJ shifted behavior to the
-      // electrostaticSummationMethod keyword options
-	  if (myMethod == "SHIFTED_POTENTIAL") {
-	    ljsp_ = true;
-	  } else if (myMethod == "SHIFTED_FORCE") {
-	    ljsf_ = true;
-	  }
-	  if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
-            if (simParams_->haveSwitchingRadius()){
-              sprintf(painCave.errMsg,
-                      "SimInfo Warning: A value was set for the switchingRadius\n"
-                      "\teven though the electrostaticSummationMethod was\n"
-                      "\tset to %s\n", myMethod.c_str());
-              painCave.isFatal = 1;
-              simError();            
-            } 
-          }
-        }
-      
-        if (simParams_->haveSwitchingRadius()){
-          rsw_ = simParams_->getSwitchingRadius();
-        } else {        
-          sprintf(painCave.errMsg,
-                  "SimCreator Warning: No value was set for switchingRadius.\n"
-                  "\tOOPSE will use a default value of\n"
-                  "\t0.85 * cutoffRadius for the switchingRadius\n");
-          painCave.isFatal = 0;
-          simError();
-          rsw_ = 0.85 * rcut_;
-        }
-
-        notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
-
-      } else {
-        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
-        // We'll punt and let fortran figure out the cutoffs later.
-        
-        notifyFortranYouAreOnYourOwn();
-
-      }
-    }
+  vector<GenericData*> SimInfo::getProperties() { 
+    return properties_.getProperties(); 
   }
 
-  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
-     
-    int errorOut;
-    int esm =  NONE;
-    int sm = UNDAMPED;
-    RealType alphaVal;
-    RealType dielectric;
-    
-    errorOut = isError;
-
-    if (simParams_->haveElectrostaticSummationMethod()) {
-      std::string myMethod = simParams_->getElectrostaticSummationMethod();
-      toUpper(myMethod);
-      if (myMethod == "NONE") {
-        esm = NONE;
-      } else {
-        if (myMethod == "SWITCHING_FUNCTION") {
-          esm = SWITCHING_FUNCTION;
-        } else {
-	  if (myMethod == "SHIFTED_POTENTIAL") {
-	    esm = SHIFTED_POTENTIAL;
-	  } else {
-	    if (myMethod == "SHIFTED_FORCE") {            
-	      esm = SHIFTED_FORCE;
-	    } else {
-	      if (myMethod == "REACTION_FIELD") {
-		esm = REACTION_FIELD;
-		dielectric = simParams_->getDielectric();
-		if (!simParams_->haveDielectric()) {
-		  // throw warning
-		  sprintf( painCave.errMsg,
-			   "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
-			   "\tA default value of %f will be used for the dielectric.\n", dielectric);
-		  painCave.isFatal = 0;
-		  simError();
-		}
-	      } else {
-		// throw error        
-		sprintf( painCave.errMsg,
-			 "SimInfo error: Unknown electrostaticSummationMethod.\n"
-                         "\t(Input file specified %s .)\n"
-                         "\telectrostaticSummationMethod must be one of: \"none\",\n"
-                         "\t\"shifted_potential\", \"shifted_force\", or \n"
-                         "\t\"reaction_field\".\n", myMethod.c_str() );
-		painCave.isFatal = 1;
-		simError();
-	      }     
-	    }           
-	  }
-	}
-      }
-    }
-    
-    if (simParams_->haveElectrostaticScreeningMethod()) {
-      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
-      toUpper(myScreen);
-      if (myScreen == "UNDAMPED") {
-	sm = UNDAMPED;
-      } else {
-	if (myScreen == "DAMPED") {
-	  sm = DAMPED;
-	  if (!simParams_->haveDampingAlpha()) {
-	    // first set a cutoff dependent alpha value
-	    // we assume alpha depends linearly with rcut from 0 to 20.5 ang
-	    alphaVal = 0.5125 - rcut_* 0.025;
-	    // for values rcut > 20.5, alpha is zero
-	    if (alphaVal < 0) alphaVal = 0;
-
-	    // throw warning
-	    sprintf( painCave.errMsg,
-		     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
-                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
-	    painCave.isFatal = 0;
-	    simError();
-	  } else {
-	    alphaVal = simParams_->getDampingAlpha();
-	  }
-	  
-	} else {
-	  // throw error        
-	  sprintf( painCave.errMsg,
-		   "SimInfo error: Unknown electrostaticScreeningMethod.\n"
-                   "\t(Input file specified %s .)\n"
-                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
-                   "or \"damped\".\n", myScreen.c_str() );
-	  painCave.isFatal = 1;
-	  simError();
-	}
-      }
-    }
-    
-    // let's pass some summation method variables to fortran
-    setElectrostaticSummationMethod( &esm );
-    setFortranElectrostaticMethod( &esm );
-    setScreeningMethod( &sm );
-    setDampingAlpha( &alphaVal );
-    setReactionFieldDielectric( &dielectric );
-    initFortranFF( &errorOut );
-  }
-
-  void SimInfo::setupSwitchingFunction() {    
-    int ft = CUBIC;
-
-    if (simParams_->haveSwitchingFunctionType()) {
-      std::string funcType = simParams_->getSwitchingFunctionType();
-      toUpper(funcType);
-      if (funcType == "CUBIC") {
-        ft = CUBIC;
-      } else {
-        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
-          ft = FIFTH_ORDER_POLY;
-	} else {
-	  // throw error        
-	  sprintf( painCave.errMsg,
-		   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
-	  painCave.isFatal = 1;
-	  simError();
-        }           
-      }
-    }
-
-    // send switching function notification to switcheroo
-    setFunctionType(&ft);
-
-  }
-
-  void SimInfo::setupAccumulateBoxDipole() {    
-
-    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
-    if ( simParams_->haveAccumulateBoxDipole() ) 
-      if ( simParams_->getAccumulateBoxDipole() ) {
-	setAccumulateBoxDipole();
-	calcBoxDipole_ = true;
-      }
-
-  }
-
-  void SimInfo::addProperty(GenericData* genData) {
-    properties_.addProperty(genData);  
-  }
-
-  void SimInfo::removeProperty(const std::string& propName) {
-    properties_.removeProperty(propName);  
-  }
-
-  void SimInfo::clearProperties() {
-    properties_.clearProperties(); 
-  }
-
-  std::vector<std::string> SimInfo::getPropertyNames() {
-    return properties_.getPropertyNames();  
-  }
-      
-  std::vector<GenericData*> SimInfo::getProperties() { 
-    return properties_.getProperties(); 
-  }
-
-  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
+  GenericData* SimInfo::getPropertyByName(const string& propName) {
     return properties_.getPropertyByName(propName); 
   }
 
@@ -1366,276 +975,65 @@ namespace oopse {
     Molecule* mol;
     RigidBody* rb;
     Atom* atom;
+    CutoffGroup* cg;
     SimInfo::MoleculeIterator mi;
     Molecule::RigidBodyIterator rbIter;
-    Molecule::AtomIterator atomIter;;
+    Molecule::AtomIterator atomIter;
+    Molecule::CutoffGroupIterator cgIter;
  
     for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
         
-      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
+      for (atom = mol->beginAtom(atomIter); atom != NULL; 
+           atom = mol->nextAtom(atomIter)) {
 	atom->setSnapshotManager(sman_);
       }
         
-      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
+      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_);
+      }
     }    
     
   }
 
-  Vector3d SimInfo::getComVel(){ 
-    SimInfo::MoleculeIterator i;
-    Molecule* mol;
 
-    Vector3d comVel(0.0);
-    RealType totalMass = 0.0;
-    
- 
-    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-      RealType mass = mol->getMass();
-      totalMass += mass;
-      comVel += mass * mol->getComVel();
-    }  
+  ostream& operator <<(ostream& o, SimInfo& info) {
 
-#ifdef IS_MPI
-    RealType tmpMass = totalMass;
-    Vector3d tmpComVel(comVel);    
-    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-#endif
-
-    comVel /= totalMass;
-
-    return comVel;
-  }
-
-  Vector3d SimInfo::getCom(){ 
-    SimInfo::MoleculeIterator i;
-    Molecule* mol;
-
-    Vector3d com(0.0);
-    RealType totalMass = 0.0;
-     
-    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-      RealType mass = mol->getMass();
-      totalMass += mass;
-      com += mass * mol->getCom();
-    }  
-
-#ifdef IS_MPI
-    RealType tmpMass = totalMass;
-    Vector3d tmpCom(com);    
-    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-#endif
-
-    com /= totalMass;
-
-    return com;
-
-  }        
-
-  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
-
     return o;
   }
    
-   
-   /* 
-   Returns center of mass and center of mass velocity in one function call.
-   */
-   
-   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ 
-      SimInfo::MoleculeIterator i;
-      Molecule* mol;
-      
-    
-      RealType totalMass = 0.0;
-    
-
-      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-         RealType mass = mol->getMass();
-         totalMass += mass;
-         com += mass * mol->getCom();
-         comVel += mass * mol->getComVel();           
-      }  
-      
-#ifdef IS_MPI
-      RealType tmpMass = totalMass;
-      Vector3d tmpCom(com);  
-      Vector3d tmpComVel(comVel);
-      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-#endif
-      
-      com /= totalMass;
-      comVel /= totalMass;
-   }        
-   
-   /* 
-   Return intertia tensor for entire system and angular momentum Vector.
-
-
-       [  Ixx -Ixy  -Ixz ]
-  J =| -Iyx  Iyy  -Iyz |
-       [ -Izx -Iyz   Izz ]
-    */
-
-   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
-      
- 
-      RealType xx = 0.0;
-      RealType yy = 0.0;
-      RealType zz = 0.0;
-      RealType xy = 0.0;
-      RealType xz = 0.0;
-      RealType yz = 0.0;
-      Vector3d com(0.0);
-      Vector3d comVel(0.0);
-      
-      getComAll(com, comVel);
-      
-      SimInfo::MoleculeIterator i;
-      Molecule* mol;
-      
-      Vector3d thisq(0.0);
-      Vector3d thisv(0.0);
-
-      RealType thisMass = 0.0;
-     
-      
-      
-   
-      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
-        
-         thisq = mol->getCom()-com;
-         thisv = mol->getComVel()-comVel;
-         thisMass = mol->getMass();
-         // Compute moment of intertia coefficients.
-         xx += thisq[0]*thisq[0]*thisMass;
-         yy += thisq[1]*thisq[1]*thisMass;
-         zz += thisq[2]*thisq[2]*thisMass;
-         
-         // compute products of intertia
-         xy += thisq[0]*thisq[1]*thisMass;
-         xz += thisq[0]*thisq[2]*thisMass;
-         yz += thisq[1]*thisq[2]*thisMass;
-            
-         angularMomentum += cross( thisq, thisv ) * thisMass;
-            
-      }  
-      
-      
-      inertiaTensor(0,0) = yy + zz;
-      inertiaTensor(0,1) = -xy;
-      inertiaTensor(0,2) = -xz;
-      inertiaTensor(1,0) = -xy;
-      inertiaTensor(1,1) = xx + zz;
-      inertiaTensor(1,2) = -yz;
-      inertiaTensor(2,0) = -xz;
-      inertiaTensor(2,1) = -yz;
-      inertiaTensor(2,2) = xx + yy;
-      
-#ifdef IS_MPI
-      Mat3x3d tmpI(inertiaTensor);
-      Vector3d tmpAngMom;
-      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-#endif
-               
-      return;
-   }
-
-   //Returns the angular momentum of the system
-   Vector3d SimInfo::getAngularMomentum(){
-      
-      Vector3d com(0.0);
-      Vector3d comVel(0.0);
-      Vector3d angularMomentum(0.0);
-      
-      getComAll(com,comVel);
-      
-      SimInfo::MoleculeIterator i;
-      Molecule* mol;
-      
-      Vector3d thisr(0.0);
-      Vector3d thisp(0.0);
-      
-      RealType thisMass;
-      
-      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {         
-        thisMass = mol->getMass(); 
-	thisr = mol->getCom()-com;
-	thisp = (mol->getComVel()-comVel)*thisMass;
-         
-	angularMomentum += cross( thisr, thisp );
-         
-      }  
-       
-#ifdef IS_MPI
-      Vector3d tmpAngMom;
-      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
-#endif
-      
-      return angularMomentum;
-   }
-   
+  
   StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
-    return IOIndexToIntegrableObject.at(index);
+    if (index >= IOIndexToIntegrableObject.size()) {
+      sprintf(painCave.errMsg,
+              "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n"
+              "\tindex exceeds number of known objects!\n");
+      painCave.isFatal = 1;
+      simError();
+      return NULL;
+    } else
+      return IOIndexToIntegrableObject.at(index);
   }
   
-  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
+  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
     IOIndexToIntegrableObject= v;
   }
 
-  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes 
-     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
-     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to 
-     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
-  */
-  void SimInfo::getGyrationalVolume(RealType &volume){
-    Mat3x3d intTensor;
-    RealType det;
-    Vector3d dummyAngMom; 
-    RealType sysconstants;
-    RealType geomCnst;
-
-    geomCnst = 3.0/2.0;
-    /* Get the inertial tensor and angular momentum for free*/
-    getInertiaTensor(intTensor,dummyAngMom);
-    
-    det = intTensor.determinant();
-    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
-    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
-    return;
+  int SimInfo::getNGlobalConstraints() {
+    int nGlobalConstraints;
+#ifdef IS_MPI
+    MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, 
+                              MPI::INT, MPI::SUM);
+#else
+    nGlobalConstraints =  nConstraints_;
+#endif
+    return nGlobalConstraints;
   }
 
-  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
-    Mat3x3d intTensor;
-    Vector3d dummyAngMom; 
-    RealType sysconstants;
-    RealType geomCnst;
+}//end namespace OpenMD
 
-    geomCnst = 3.0/2.0;
-    /* Get the inertial tensor and angular momentum for free*/
-    getInertiaTensor(intTensor,dummyAngMom);
-    
-    detI = intTensor.determinant();
-    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
-    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
-    return;
-  }
-/*
-   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
-      assert( v.size() == nAtoms_ + nRigidBodies_);
-      sdByGlobalIndex_ = v;
-    }
-
-    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
-      //assert(index < nAtoms_ + nRigidBodies_);
-      return sdByGlobalIndex_.at(index);
-    }   
-*/   
-}//end namespace oopse
-