--- branches/development/src/nonbonded/Electrostatic.cpp	2013/01/08 15:29:03	1822
+++ branches/development/src/nonbonded/Electrostatic.cpp	2013/02/20 15:39:39	1850
@@ -35,7 +35,7 @@
  *                                                                      
  * [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).          
+ * [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).
  */
@@ -106,8 +106,8 @@ namespace OpenMD {
     angstromToM_ = 1.0e-10;
     debyeToCm_ = 3.33564095198e-30;
     
-    // number of points for electrostatic splines
-    np_ = 1000;
+    // Default number of points for electrostatic splines
+    np_ = 140;
     
     // variables to handle different summation methods for long-range 
     // electrostatics:
@@ -246,7 +246,7 @@ namespace OpenMD {
       b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
       b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
       b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
-      selfMult_ = b0c  +  2.0 * a2 * invArootPi;
+      selfMult_ = b0c + a2 * invArootPi;
     } else {
       a2 = 0.0;
       b0c = 1.0 / r;
@@ -279,6 +279,11 @@ namespace OpenMD {
     // working variables for Taylor expansion:
     RealType rmRc, rmRc2, rmRc3, rmRc4;
 
+    // Approximate using splines using a maximum of 0.1 Angstroms
+    // between points.
+    int nptest = int((cutoffRadius_ + 2.0) / 0.1);
+    np_ = (np_ > nptest) ? np_ : nptest;
+  
     // Add a 2 angstrom safety window to deal with cutoffGroups that
     // have charged atoms longer than the cutoffRadius away from each
     // other.  Splining is almost certainly the best choice here.
@@ -685,7 +690,7 @@ namespace OpenMD {
     Vector3d D_a, D_b;  // Dipoles (space-fixed)
     Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
 
-    RealType ri, ri2, ri3, ri4;                  // Distance utility scalars
+    RealType ri;                                 // Distance utility scalar
     RealType rdDa, rdDb;                         // Dipole utility scalars
     Vector3d rxDa, rxDb;                         // Dipole utility vectors
     RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
@@ -730,7 +735,6 @@ namespace OpenMD {
 
     ri = 1.0 /  *(idat.rij);
     Vector3d rhat =  *(idat.d)  * ri;
-    ri2 = ri * ri;
       
     // logicals
 
@@ -804,15 +808,18 @@ namespace OpenMD {
       
       if (idat.excluded) {
         *(idat.skippedCharge2) += C_a;
+      } else {
+        // only do the field if we're not excluded:
+        Eb -= C_a *  pre11_ * v02 * rhat;
       }
-      Eb -= C_a *  pre11_ * v02 * rhat;
     }
     
     if (a_is_Dipole) {
       D_a = *(idat.dipole1);
       rdDa = dot(rhat, D_a);
       rxDa = cross(rhat, D_a);
-      Eb -=  pre12_ * (v13 * rdDa * rhat + v12 * D_a);
+      if (!idat.excluded) 
+        Eb -=  pre12_ * (v13 * rdDa * rhat + v12 * D_a);
     }
     
     if (a_is_Quadrupole) {
@@ -822,7 +829,8 @@ namespace OpenMD {
       rQa = rhat * Q_a;
       rdQar = dot(rhat, Qar);
       rxQar = cross(rhat, Qar);
-      Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24);
+      if (!idat.excluded) 
+        Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24);
     }
     
     if (b_is_Charge) {
@@ -833,15 +841,18 @@ namespace OpenMD {
       
       if (idat.excluded) {
         *(idat.skippedCharge1) += C_b;
+      } else {
+        // only do the field if we're not excluded:
+        Ea += C_b *  pre11_ * v02 * rhat;
       }
-      Ea += C_b *  pre11_ * v02 * rhat;
     }
     
     if (b_is_Dipole) {
       D_b = *(idat.dipole2);
       rdDb = dot(rhat, D_b);
       rxDb = cross(rhat, D_b);
-      Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b);
+      if (!idat.excluded) 
+        Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b);
     }
     
     if (b_is_Quadrupole) {
@@ -851,7 +862,8 @@ namespace OpenMD {
       rQb = rhat * Q_b;
       rdQbr = dot(rhat, Qbr);
       rxQbr = cross(rhat, Qbr);
-      Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24);
+      if (!idat.excluded) 
+        Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24);
     }
     
     if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
@@ -1088,7 +1100,7 @@ namespace OpenMD {
         *(idat.t2) += *(idat.sw) * indirect_Tb;
     }
     return;
-  }  
+  }
     
   void Electrostatic::calcSelfCorrection(SelfData &sdat) {
 
@@ -1131,8 +1143,8 @@ namespace OpenMD {
       
     case esm_SHIFTED_FORCE:
     case esm_SHIFTED_POTENTIAL:
-      if (i_is_Charge) {        
-        self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
+      if (i_is_Charge) {
+        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
         (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
       }
       break;