--- branches/development/src/UseTheForce/ForceField.cpp	2010/07/20 15:43:00	1473
+++ branches/development/src/UseTheForce/ForceField.cpp	2012/05/18 21:44:02	1710
@@ -36,7 +36,8 @@
  * [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).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
 /**
@@ -51,9 +52,8 @@
 #include "UseTheForce/ForceField.hpp"
 #include "utils/simError.h"
 #include "utils/Tuple.hpp"
-#include "UseTheForce/DarkSide/atype_interface.h"
-#include "UseTheForce/DarkSide/fForceOptions_interface.h"
-#include "UseTheForce/DarkSide/switcheroo_interface.h"
+#include "types/LennardJonesAdapter.hpp"
+
 namespace OpenMD {
 
   ForceField::ForceField() { 
@@ -69,18 +69,29 @@ namespace OpenMD {
     }
   }
 
-
-  ForceField::~ForceField() {
-    deleteAtypes();
-    deleteSwitch();   
-  }
-
+  /**
+   * getAtomType by string
+   *
+   * finds the requested atom type in this force field using the string
+   * name of the atom type.
+   */
   AtomType* ForceField::getAtomType(const std::string &at) {
     std::vector<std::string> keys;
     keys.push_back(at);
     return atomTypeCont_.find(keys);
   }
 
+  /**
+   * getAtomType by ident
+   *
+   * finds the requested atom type in this force field using the
+   * integer ident instead of the string name of the atom type.
+   */
+  AtomType* ForceField::getAtomType(int ident) {   
+    std::string at = atypeIdentToName.find(ident)->second;
+    return getAtomType(at);
+  }
+
   BondType* ForceField::getBondType(const std::string &at1, 
 				    const std::string &at2) {
     std::vector<std::string> keys;
@@ -440,6 +451,7 @@ namespace OpenMD {
   }
   
   NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
+    
     std::vector<std::string> keys;
     keys.push_back(at1);
     keys.push_back(at2);    
@@ -449,9 +461,63 @@ namespace OpenMD {
     if (nbiType) {
       return nbiType;
     } else {
-      //if no exact match found, try wild card match
-      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
-    }    
+      AtomType* atype1;
+      AtomType* atype2;
+      std::vector<std::string> at1key;
+      at1key.push_back(at1);
+      atype1 = atomTypeCont_.find(at1key);
+      
+      std::vector<std::string> at2key;
+      at2key.push_back(at2);
+      atype2 = atomTypeCont_.find(at2key);
+      
+      // query atom types for their chains of responsibility
+      std::vector<AtomType*> at1Chain = atype1->allYourBase();
+      std::vector<AtomType*> at2Chain = atype2->allYourBase();
+      
+      std::vector<AtomType*>::iterator i;
+      std::vector<AtomType*>::iterator j;
+      
+      int ii = 0;
+      int jj = 0;
+      int nbiTypeScore;
+      
+      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
+      
+      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
+        jj = 0;
+        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
+          
+          nbiTypeScore = ii + jj;
+          
+          std::vector<std::string> myKeys;
+          myKeys.push_back((*i)->getName());
+          myKeys.push_back((*j)->getName());
+          
+          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
+          if (nbiType) {
+            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
+          }
+          jj++;
+        }
+        ii++;
+      }
+      
+      
+      if (foundNBI.size() > 0) {
+        // sort the foundNBI by the score:
+        std::sort(foundNBI.begin(), foundNBI.end());
+        
+        int bestScore = foundNBI[0].first;
+        std::vector<std::string> theKeys = foundNBI[0].second;
+        
+        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
+        return bestType;
+      } else {
+        //if no exact match found, try wild card match
+        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
+      }
+    }
   }
   
   BondType* ForceField::getExactBondType(const std::string &at1, 
@@ -507,12 +573,14 @@ namespace OpenMD {
   bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
     std::vector<std::string> keys;
     keys.push_back(at);
+    atypeIdentToName[atomType->getIdent()] = at;
     return atomTypeCont_.add(keys, atomType);
   }
 
   bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
     std::vector<std::string> keys;
     keys.push_back(at);
+    atypeIdentToName[atomType->getIdent()] = at;
     return atomTypeCont_.replace(keys, atomType);
   }
 
@@ -569,34 +637,11 @@ namespace OpenMD {
   }
   
   RealType ForceField::getRcutFromAtomType(AtomType* at) {
-    /**@todo */
-    GenericData* data;
     RealType rcut = 0.0;
     
-    if (at->isLennardJones()) {
-      data = at->getPropertyByName("LennardJones");
-      if (data != NULL) {
-	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
-	
-	if (ljData != NULL) {
-	  LJParam ljParam = ljData->getData();
-	  
-	  //by default use 2.5*sigma as cutoff radius
-	  rcut = 2.5 * ljParam.sigma;
-          
-	} else {
-	  sprintf( painCave.errMsg,
-		   "Can not cast GenericData to LJParam\n");
-	  painCave.severity = OPENMD_ERROR;
-	  painCave.isFatal = 1;
-	  simError();          
-	}            
-      } else {
-	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
-	painCave.severity = OPENMD_ERROR;
-	painCave.isFatal = 1;
-	simError();          
-      }
+    LennardJonesAdapter lja = LennardJonesAdapter(at);
+    if (lja.isLennardJones()) {
+      rcut = 2.5 * lja.getSigma();
     }
     return rcut;    
   }
@@ -631,9 +676,4 @@ namespace OpenMD {
     return ffStream;
   }
 
-  void ForceField::setFortranForceOptions(){
-    ForceOptions theseFortranOptions;
-    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
-    setfForceOptions(&theseFortranOptions);
-  }
 } //end namespace OpenMD