--- branches/development/src/nonbonded/SHAPES.cpp	2010/10/03 22:18:59	1505
+++ trunk/src/nonbonded/SHAPES.cpp	2014/11/01 14:12:16	2033
@@ -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 <stdio.h>
@@ -85,22 +86,23 @@ namespace OpenMD {
     // add it to the map:
     AtomTypeProperties atp = atomType->getATP();    
 
-    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, sAtomType) );
-    
-  } else if (atomType->isLennardJones()) {
-      d1 = getLJSigma(atomType) / sqrt(2.0);
-      e1 = getLJEpsilon(atomType);
+    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()) {
+      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"
@@ -261,43 +263,40 @@ namespace OpenMD {
       // sin(theta) is near 0.0:
 
       RealType sti2 = 1.0 - cti*cti;
+      RealType proji;
+      Vector3d dcpidr, dcpidu, dspidr, dspidu;
       if (fabs(sti2) < 1.0e-12) {
-        RealType proji = sqrt(r * 1.0e-12);
-        Vector3d dcpidx(1.0 / proji,
-                        0.0,
-                        
-        dcpidx = 1.0_dp / proji
-          dcpidy = 0.0_dp
-          dcpidux = xi / proji
-          dcpiduy = 0.0_dp
-          dspidx = 0.0_dp
-          dspidy = 1.0_dp / proji
-          dspidux = 0.0_dp
-          dspiduy = yi / proji
-       else
-          proji = sqrt(xi2 + yi2)
-          proji3 = proji*proji*proji
-          dcpidx = 1.0_dp / proji - xi2 / proji3
-          dcpidy = - xi * yi / proji3
-          dcpidux = xi / proji - (xi2 * xi) / proji3
-          dcpiduy = - (xi * yi2) / proji3
-          dspidx = - xi * yi / proji3
-          dspidy = 1.0_dp / proji - yi2 / proji3
-          dspidux = - (yi * xi2) / proji3
-          dspiduy = yi / proji - (yi2 * yi) / proji3
-       endif
+        proji = sqrt(r * 1.0e-12);
+        dcpidr = Vector3d( 1.0 / proji, 0.0,         0.0);         
+        dcpidu = Vector3d( xi / proji,  0.0,         0.0);
+        dspidr = Vector3d( 0.0,         1.0 / proji, 0.0);
+        dspidu = Vector3d (0.0,         yi / proji,  0.0);
+      } else {
+        proji = sqrt(xi2 + yi2);
+        RealType proji3 = proji*proji*proji;
+        dcpidr = Vector3d(1.0_dp / proji - xi2 / proji3,
+                          - xi * yi / proji3,
+                          0.0);
+        dcpidu = Vector3d( xi / proji - (xi2 * xi) / proji3,
+                           - (xi * yi2) / proji3,
+                           0.0);
+        dspidr = Vector3d(- xi * yi / proji3,
+                          1.0_dp / proji - yi2 / proji3,
+                          0.0);
+        dspidu = Vector3d(- (yi * xi2) / proji3,
+                          yi / proji - (yi2 * yi) / proji3,
+                          0.0);
+      }
 
-       cpi = xi / proji
-       dcpidz = 0.0_dp
-       dcpiduz = 0.0_dp
+      cpi = xi / proji;
+      dcpidr.z() = 0.0;
+      dcpidu.z() = 0.0;
+      
+      spi = yi / proji;
+      dspidr.z() = 0.0;
+      dspidu.z() = 0.0;
 
-       spi = yi / proji
-       dspidz = 0.0_dp
-       dspiduz = 0.0_dp
 
-
-
-
     RealType sigma0 = mixer.sigma0;
     RealType dw     = mixer.dw;
     RealType eps0   = mixer.eps0;