--- branches/development/src/nonbonded/Electrostatic.cpp	2012/06/07 12:53:46	1750
+++ branches/development/src/nonbonded/Electrostatic.cpp	2012/06/21 19:26:46	1760
@@ -578,12 +578,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;
@@ -1049,6 +1049,7 @@ namespace OpenMD {
       // indirect reaction field terms.
 
       *(idat.vpair) += indirect_vpair;
+      (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] += epot;
       (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
       *(idat.f1) += indirect_dVdr;
       
@@ -1077,6 +1078,7 @@ namespace OpenMD {
       chg1 += *(sdat.flucQ);
       // dVdFQ is really a force, so this is negative the derivative
       *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
+      cerr << "dVdFQ harmonic part = " << *(sdat.dVdFQ) << "\n";
     }
 
     if (summationMethod_ == esm_REACTION_FIELD) {