--- branches/development/src/nonbonded/Electrostatic.cpp	2012/06/05 17:49:36	1734
+++ branches/development/src/nonbonded/Electrostatic.cpp	2012/07/06 22:01:58	1767
@@ -53,8 +53,8 @@
 #include "io/Globals.hpp"
 #include "nonbonded/SlaterIntegrals.hpp"
 #include "utils/PhysicalConstants.hpp"
+#include "math/erfc.hpp"
 
-
 namespace OpenMD {
   
   Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
@@ -193,8 +193,9 @@ namespace OpenMD {
         
         // throw warning
         sprintf( painCave.errMsg,
-                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
-                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", 
+                 "Electrostatic::initialize: dampingAlpha was not specified in the\n"
+                 "\tinput file.  A default value of %f (1/ang) will be used for the\n"
+                 "\tcutoff of %f (ang).\n", 
                  dampingAlpha_, cutoffRadius_);
         painCave.severity = OPENMD_INFO;
         painCave.isFatal = 0;
@@ -360,12 +361,13 @@ namespace OpenMD {
         vector<RealType> rvals;
         vector<RealType> J1vals;
         vector<RealType> J2vals;
-        for (int i = 0; i < np_; i++) {
+        // don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals.
+        for (int i = 1; i < np_+1; i++) {
           rval = RealType(i) * dr;
           rvals.push_back(rval);
-          J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
+          J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
           // may not be necessary if Slater coulomb integral is symmetric
-          J2vals.push_back(eaData2.hardness *  sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
+          J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
         }
 
         CubicSpline* J1 = new CubicSpline();
@@ -577,12 +579,12 @@ namespace OpenMD {
       if (j_is_Charge) {
         if (screeningMethod_ == DAMPED) {
           // assemble the damping variables
-          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
-          //erfcVal = res.first;
-          //derfcVal = res.second;
+          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
+          erfcVal = res.first;
+          derfcVal = res.second;
 
-          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
-          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
+          //erfcVal = erfc(dampingAlpha_ * *(idat.rij));
+          //derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
 
           c1 = erfcVal * riji;
           c2 = (-derfcVal + c1) * riji;
@@ -1048,6 +1050,9 @@ namespace OpenMD {
       // indirect reaction field terms.
 
       *(idat.vpair) += indirect_vpair;
+      
+      (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=   (*(idat.sw) * vterm + 
+                                                        vFluc1 ) * q_i * q_j;
       (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
       *(idat.f1) += indirect_dVdr;
       
@@ -1076,6 +1081,8 @@ namespace OpenMD {
       chg1 += *(sdat.flucQ);
       // dVdFQ is really a force, so this is negative the derivative
       *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
+      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * 
+        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
     }
 
     if (summationMethod_ == esm_REACTION_FIELD) {