--- branches/development/src/nonbonded/InteractionManager.cpp	2011/06/21 06:34:35	1586
+++ branches/development/src/nonbonded/InteractionManager.cpp	2013/05/15 15:09:35	1874
@@ -35,8 +35,9 @@
  *                                                                      
  * [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]  Vardeman & Gezelter, in progress (2009).                        
+ * [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).
  */
  
 #include "nonbonded/InteractionManager.hpp"
@@ -51,14 +52,29 @@ namespace OpenMD {
     gb_ = new GB();
     sticky_ = new Sticky();
     morse_ = new Morse();
+    repulsivePower_ = new RepulsivePower();
     eam_ = new EAM();
     sc_ = new SC();
     electrostatic_ = new Electrostatic();
     maw_ = new MAW();
   }
 
+  InteractionManager::~InteractionManager() {
+    delete lj_;
+    delete gb_;
+    delete sticky_;
+    delete morse_;
+    delete repulsivePower_;
+    delete eam_;
+    delete sc_;
+    delete electrostatic_;
+    delete maw_;
+  }
+
   void InteractionManager::initialize() {
-    
+
+    if (initialized_) return; 
+
     ForceField* forceField_ = info_->getForceField();
     
     lj_->setForceField(forceField_);
@@ -70,26 +86,25 @@ namespace OpenMD {
     electrostatic_->setSimInfo(info_);
     electrostatic_->setForceField(forceField_);
     maw_->setForceField(forceField_);
+    repulsivePower_->setForceField(forceField_);
 
     ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
     ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
     AtomType* atype1;
     AtomType* atype2;
     pair<AtomType*, AtomType*> key;
-    pair<set<NonBondedInteraction*>::iterator, bool> ret;
     
     for (atype1 = atomTypes->beginType(i1); atype1 != NULL; 
          atype1 = atomTypes->nextType(i1)) {
       
       // add it to the map:
-      AtomTypeProperties atp = atype1->getATP();    
       
       pair<map<int,AtomType*>::iterator,bool> ret;    
-      ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
+      ret = typeMap_.insert( pair<int, AtomType*>(atype1->getIdent(), atype1) );
       if (ret.second == false) {
         sprintf( painCave.errMsg,
                  "InteractionManager already had a previous entry with ident %d\n",
-                 atp.ident);
+                 atype1->getIdent());
         painCave.severity = OPENMD_INFO;
         painCave.isFatal = 0;
         simError();                 
@@ -104,11 +119,7 @@ namespace OpenMD {
 
       for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {        
         atype2 = (*it2).second;
-        
-        bool vdwExplicit = false;
-        bool metExplicit = false;
-        bool hbExplicit = false;
-                       
+                               
         key = make_pair(atype1, atype2);
         
         if (atype1->isLennardJones() && atype2->isLennardJones()) {
@@ -143,6 +154,10 @@ namespace OpenMD {
         
         if (nbiType != NULL) {
 
+          bool vdwExplicit = false;
+          bool metExplicit = false;
+          // bool hbExplicit = false;
+
           if (nbiType->isLennardJones()) {
             // We found an explicit Lennard-Jones interaction.  
             // override all other vdw entries for this pair of atom types:
@@ -176,9 +191,33 @@ namespace OpenMD {
               if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
             }
             interactions_[key].insert(morse_);
+            vdwExplicit = true;
+          }
+
+          if (nbiType->isRepulsivePower()) {
+            if (vdwExplicit) {
+              sprintf( painCave.errMsg,
+                       "InteractionManager::initialize found more than one "
+                       "explicit \n"
+                       "\tvan der Waals interaction for atom types %s - %s\n",
+                       atype1->getName().c_str(), atype2->getName().c_str());
+              painCave.severity = OPENMD_ERROR;
+              painCave.isFatal = 1;
+              simError();
+            }
+            // We found an explicit RepulsivePower interaction.  
+            // override all other vdw entries for this pair of atom types:
+            set<NonBondedInteraction*>::iterator it;
+            for (it = interactions_[key].begin(); 
+                 it != interactions_[key].end(); ++it) {
+              InteractionFamily ifam = (*it)->getFamily();
+              if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
+            }
+            interactions_[key].insert(repulsivePower_);
             vdwExplicit = true;
           }
           
+          
           if (nbiType->isEAM()) {
             // We found an explicit EAM interaction.  
             // override all other metallic entries for this pair of atom types:
@@ -271,18 +310,18 @@ namespace OpenMD {
   }
 
   void InteractionManager::setCutoffRadius(RealType rcut) {
+    
     electrostatic_->setCutoffRadius(rcut);
     eam_->setCutoffRadius(rcut);
   }
 
-  void InteractionManager::setSwitchingRadius(RealType rswitch) {
-    electrostatic_->setSwitchingRadius(rswitch);
-  }
-  
   void InteractionManager::doPrePair(InteractionData idat){
     
     if (!initialized_) initialize();
 	 
+    // excluded interaction, so just return
+    if (idat.excluded) return;
+
     set<NonBondedInteraction*>::iterator it;
 
     for (it = interactions_[ idat.atypes ].begin(); 
@@ -314,26 +353,17 @@ namespace OpenMD {
   void InteractionManager::doPair(InteractionData idat){
     
     if (!initialized_) initialize();
-   
+
     set<NonBondedInteraction*>::iterator it;
 
     for (it = interactions_[ idat.atypes ].begin(); 
-         it != interactions_[ idat.atypes ].end(); ++it)
-      (*it)->calcForce(idat);
-    
-    return;    
-  }
+         it != interactions_[ idat.atypes ].end(); ++it) {
 
-  void InteractionManager::doSkipCorrection(InteractionData idat){
+      // electrostatics still has to worry about indirect
+      // contributions from excluded pairs of atoms:
 
-    if (!initialized_) initialize();  
-
-    set<NonBondedInteraction*>::iterator it;
-
-    for (it = interactions_[ idat.atypes ].begin(); 
-         it != interactions_[ idat.atypes ].end(); ++it){
-      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
-        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat);
+      if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) {
+        (*it)->calcForce(idat);
       }
     }
     
@@ -358,7 +388,7 @@ namespace OpenMD {
 
   RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
     if (!initialized_) initialize();
-    
+
     AtomType* atype = typeMap_[*atid];
 
     pair<AtomType*, AtomType*> key = make_pair(atype, atype);