--- branches/development/src/nonbonded/Electrostatic.cpp	2011/08/18 20:18:19	1613
+++ branches/development/src/nonbonded/Electrostatic.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>
@@ -64,6 +65,7 @@ namespace OpenMD {
     Globals* simParams_ = info_->getSimParams();
 
     summationMap_["HARD"]               = esm_HARD;
+    summationMap_["NONE"]               = esm_HARD;
     summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
     summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL; 
     summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;     
@@ -116,7 +118,7 @@ namespace OpenMD {
         sprintf( painCave.errMsg,
                  "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
                  "\t(Input file specified %s .)\n"
-                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
+                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
                  "\t\"shifted_potential\", \"shifted_force\", or \n"
                  "\t\"reaction_field\".\n", myMethod.c_str() );
         painCave.isFatal = 1;
@@ -583,9 +585,13 @@ 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));
+
           c1 = erfcVal * riji;
           c2 = (-derfcVal + c1) * riji;
         } else {
@@ -672,9 +678,11 @@ namespace OpenMD {
 
           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));
             c1 = erfcVal * ri;
             c2 = (-derfcVal + c1) * ri;
             c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
@@ -709,9 +717,11 @@ namespace OpenMD {
           
         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));
           c1 = erfcVal * riji;
           c2 = (-derfcVal + c1) * riji;
           c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
@@ -798,9 +808,11 @@ namespace OpenMD {
             
           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));
             c1 = erfcVal * ri;
             c2 = (-derfcVal + c1) * ri;
             c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
@@ -877,9 +889,11 @@ namespace OpenMD {
           }
           if (screeningMethod_ == DAMPED) {
             // assemble 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));
             c1 = erfcVal * ri;
             c2 = (-derfcVal + c1) * ri;
             c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
@@ -928,9 +942,11 @@ namespace OpenMD {
 
         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));
           c1 = erfcVal * riji;
           c2 = (-derfcVal + c1) * riji;
           c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
@@ -995,6 +1011,7 @@ namespace OpenMD {
 
       // only accumulate the forces and torques resulting from the
       // indirect reaction field terms.
+
       *(idat.vpair) += indirect_vpair;
       (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
       *(idat.f1) += indirect_dVdr;