--- branches/development/src/nonbonded/SHAPES.cpp	2010/09/15 19:32:10	1501
+++ branches/development/src/nonbonded/SHAPES.cpp	2011/11/22 20:38:56	1665
@@ -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).
  */
 
 #include <stdio.h>
@@ -49,23 +50,14 @@ namespace OpenMD {
 
 using namespace std;
 namespace OpenMD {
-  
-  bool SHAPES::initialized_ = false;
-  int SHAPES::lMax_ = 64;
-  int SHAPES::mMax_ = 64;
-  ForceField* SHAPES::forceField_ = NULL;
-  map<int, AtomType*> SHAPES::ShapesMap;
-  map<pair<AtomType*, AtomType*>, SHAPESInteractionData> SHAPES::MixingMap;
-  
-  SHAPES* SHAPES::_instance = NULL;
-
-  SHAPES* SHAPES::Instance() {
-    if (!_instance) {
-      _instance = new SHAPES();
-    }
-    return _instance;
+ 
+  SHAPES::SHAPES() {
+    initialized_ = false;
+    lMax_ = 64;
+    mMax_ = 64;
+    forceField_ = NULL;
   }
-
+  
   void SHAPES::initialize() {    
     
     ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
@@ -75,122 +67,98 @@ namespace OpenMD {
 
     // SHAPES handles all of the SHAPES-SHAPES interactions as well as
     // SHAPES-LJ cross interactions:
-
+    
     for (at = atomTypes->beginType(i); at != NULL;
          at = atomTypes->nextType(i)) {
       
-      if (at->isShape() || at->isLennardJones())
-        addType(at);
+      if (at->isShape()) 
+        addShape(dynamic_cast<ShapeAtomType*>(at));
+
+      if (at->isLennardJones())
+        addLJ(at);
+      
     }
-   
+    
     initialized_ = true;
   }
-      
-  void SHAPES::addType(AtomType* atomType){
+  
+  void SHAPES::addShape(ShapeAtomType* atomType){
     // add it to the map:
     AtomTypeProperties atp = atomType->getATP();    
 
-    pair<map<int,AtomType*>::iterator,bool> ret;    
-    ret = ShapesMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
-    if (ret.second == false) {
-      sprintf( painCave.errMsg,
-               "SHAPES already had a previous entry with ident %d\n",
-               atp.ident);
-      painCave.severity = OPENMD_INFO;
-      painCave.isFatal = 0;
-      simError();         
-    }   
-    
-    if (atomType->isShape()) {
-      ShapeAtomType* sAtomType = dynamic_cast<ShapeAtomType*>(atomType);
-      if (sAtomType == NULL) {
-        sprintf(painCave.errMsg,
-                "SHAPES:: Can't cast to ShapeAtomType");
-        painCave.severity = OPENMD_ERROR;
-        painCave.isFatal = 1;
-        simError();
-      } 
-      ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, sAtomType) );
-
+    if (atomType->isShape() ) {
+      pair<map<int,ShapeAtomType*>::iterator, bool> ret;
+      ret = ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, atomType));
+      if (ret.second == false) {
+        sprintf( painCave.errMsg,
+                 "SHAPES already had a previous entry with ident %d\n",
+                 atp.ident);
+        painCave.severity = OPENMD_INFO;
+        painCave.isFatal = 0;
+        simError();         
+      }   
+      
+      ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, static_cast<ShapeAtomType*>(atomType)) );
+      
     } else if (atomType->isLennardJones()) {
-      d1 = LJ::Instance()->getSigma(atomType) / sqrt(2.0);
-      e1 = LJ::Instance()->getEpsilon(atomType);
+      RealType d1 = getLJSigma(atomType) / sqrt(2.0);
+      RealType e1 = getLJEpsilon(atomType);
     } else {
       sprintf( painCave.errMsg,
                "SHAPES::addType was passed an atomType (%s) that does not\n"
-               "\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
+               "\tappear to be a SHAPES or Lennard-Jones atom.\n",
                atomType->getName().c_str());
       painCave.severity = OPENMD_ERROR;
       painCave.isFatal = 1;
       simError();
-    }
-      
+    }   
+  }
+  
 
-    // Now, iterate over all known types and add to the mixing map:
+  LJParam SHAPES::getLJParam(AtomType* atomType) {
     
-    map<int, AtomType*>::iterator it;
-    for( it = ShapesMap.begin(); it != SHAPESMap.end(); ++it) {
-      
-      AtomType* atype2 = (*it).second;
-      
-      RealType d2, l2, e2, er2, dw2;
-      
-      if (atype2->isGayBerne()) {
-        GayBerneParam gb2 = getGayBerneParam(atype2);
-        d2 = gb2.SHAPES_d;
-        l2 = gb2.SHAPES_l;
-        e2 = gb2.SHAPES_eps;
-        er2 = gb2.SHAPES_eps_ratio;
-        dw2 = gb2.SHAPES_dw;
-      } else if (atype2->isLennardJones()) {
-        d2 = LJ::Instance()->getSigma(atype2) / sqrt(2.0);
-        e2 = LJ::Instance()->getEpsilon(atype2);
-        l2 = d2;
-        er2 = 1.0;
-        dw2 = 1.0;
-      } 
-                       
-      SHAPESInteractionData mixer;        
-      
-      //  Cleaver paper uses sqrt of squares to get sigma0 for
-      //  mixed interactions.
-            
-      mixer.sigma0 = sqrt(d1*d1 + d2*d2);
-      mixer.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
-      mixer.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
-      mixer.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
-        ((l2*l2 + d1*d1) * (l1*l1 + d2*d2));
- 
-      // assumed LB mixing rules for now:
-
-      mixer.dw = 0.5 * (dw1 + dw2);
-      mixer.eps0 = sqrt(e1 * e2);
-      
-      RealType er = sqrt(er1 * er2);
-      RealType ermu = pow(er,(1.0 / mu_));
-      RealType xp = (1.0 - ermu) / (1.0 + ermu);
-      RealType ap2 = 1.0 / (1.0 + ermu);
-      
-      mixer.xp2 = xp * xp;
-      mixer.xpap2 = xp * ap2;
-      mixer.xpapi2 = xp / ap2;
-
-      // only add this pairing if at least one of the atoms is a Gay-Berne atom
-
-      if (atomType->isGayBerne() || atype2->isGayBerne()) {
-
-        pair<AtomType*, AtomType*> key1, key2;
-        key1 = make_pair(atomType, atype2);
-        key2 = make_pair(atype2, atomType);
-        
-        MixingMap[key1] = mixer;
-        if (key2 != key1) {
-          MixingMap[key2] = mixer;
-        }
-      }
-    }      
+    // Do sanity checking on the AtomType we were passed before
+    // building any data structures:
+    if (!atomType->isLennardJones()) {
+      sprintf( painCave.errMsg,
+               "SHAPES::getLJParam was passed an atomType (%s) that does not\n"
+               "\tappear to be a Lennard-Jones atom.\n",
+               atomType->getName().c_str());
+      painCave.severity = OPENMD_ERROR;
+      painCave.isFatal = 1;
+      simError();
+    }
+    
+    GenericData* data = atomType->getPropertyByName("LennardJones");
+    if (data == NULL) {
+      sprintf( painCave.errMsg, "SHAPES::getLJParam could not find Lennard-Jones\n"
+               "\tparameters for atomType %s.\n", atomType->getName().c_str());
+      painCave.severity = OPENMD_ERROR;
+      painCave.isFatal = 1;
+      simError(); 
+    }
+    
+    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
+    if (ljData == NULL) {
+      sprintf( painCave.errMsg,
+               "SHAPES::getLJParam could not convert GenericData to LJParam for\n"
+               "\tatom type %s\n", atomType->getName().c_str());
+      painCave.severity = OPENMD_ERROR;
+      painCave.isFatal = 1;
+      simError();          
+    }
+    
+    return ljData->getData();
   }
   
+  RealType SHAPES::getLJEpsilon(AtomType* atomType) {    
+    LJParam ljParam = getLJParam(atomType);
+    return ljParam.epsilon;
+  }
+  RealType SHAPES::getLJSigma(AtomType* atomType) {    
+    LJParam ljParam = getLJParam(atomType);
+    return ljParam.sigma;
+  }
 
   RealType SHAPES::getGayBerneCut(int atid) { 
     if (!initialized_) initialize();