--- branches/development/src/nonbonded/Electrostatic.cpp	2012/05/24 01:48:29	1720
+++ branches/development/src/nonbonded/Electrostatic.cpp	2012/05/24 20:59:54	1723
@@ -217,7 +217,6 @@ namespace OpenMD {
         addType(at);
     }
     
-
     cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
     rcuti_ = 1.0 / cutoffRadius_;
     rcuti2_ = rcuti_ * rcuti_;
@@ -284,6 +283,7 @@ namespace OpenMD {
     electrostaticAtomData.is_Dipole = false;
     electrostaticAtomData.is_SplitDipole = false;
     electrostaticAtomData.is_Quadrupole = false;
+    electrostaticAtomData.is_Fluctuating = false;
 
     FixedChargeAdapter fca = FixedChargeAdapter(atomType);
 
@@ -364,6 +364,7 @@ namespace OpenMD {
           rval = RealType(i) * dr;
           rvals.push_back(rval);
           J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
+          // may not be necessary if Slater coulomb integral is symmetric
           J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
         }
 
@@ -446,7 +447,12 @@ namespace OpenMD {
     Vector3d indirect_dVdr(V3Zero);
     Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
 
+    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
     pair<RealType, RealType> res;
+    
+    // splines for coulomb integrals
+    CubicSpline* J1;
+    CubicSpline* J2;
     
     if (!initialized_) initialize();
     
@@ -464,14 +470,21 @@ namespace OpenMD {
     bool i_is_Dipole = data1.is_Dipole;
     bool i_is_SplitDipole = data1.is_SplitDipole;
     bool i_is_Quadrupole = data1.is_Quadrupole;
+    bool i_is_Fluctuating = data1.is_Fluctuating;
 
     bool j_is_Charge = data2.is_Charge;
     bool j_is_Dipole = data2.is_Dipole;
     bool j_is_SplitDipole = data2.is_SplitDipole;
     bool j_is_Quadrupole = data2.is_Quadrupole;
+    bool j_is_Fluctuating = data2.is_Fluctuating;
     
     if (i_is_Charge) {
       q_i = data1.fixedCharge;
+
+      if (i_is_Fluctuating) {
+        q_i += *(idat.flucQ1);
+      }
+      
       if (idat.excluded) {
         *(idat.skippedCharge2) += q_i;
       }
@@ -510,6 +523,10 @@ namespace OpenMD {
 
     if (j_is_Charge) {
       q_j = data2.fixedCharge;
+
+      if (i_is_Fluctuating) 
+        q_j += *(idat.flucQ2);
+
       if (idat.excluded) {
         *(idat.skippedCharge1) += q_j;
       }
@@ -547,6 +564,11 @@ namespace OpenMD {
       duduz_j = V3Zero;
     }
     
+    if (i_is_Fluctuating && j_is_Fluctuating) {
+      J1 = Jij[idat.atypes];
+      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
+    }
+
     epot = 0.0;
     dVdr = V3Zero;
     
@@ -598,12 +620,38 @@ namespace OpenMD {
 
           vterm = preVal * riji * erfcVal;           
           dudr  = -  *(idat.sw)  * preVal * c2;
-
+          
         }
-
+        
         vpair += vterm;
         epot +=  *(idat.sw)  * vterm;
-        dVdr += dudr * rhat;                 
+        dVdr += dudr * rhat;
+
+        if (i_is_Fluctuating) {
+          if (idat.excluded) {
+            // vFluc1 is the difference between the direct coulomb integral
+            // and the normal 1/r-like  interaction between point charges.
+            coulInt = J1->getValueAt( *(idat.rij) );
+            vFluc1 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
+          } else {
+            vFluc1 = 0.0;
+          }
+          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) / q_i;
+        }
+
+        if (j_is_Fluctuating) {
+          if (idat.excluded) {
+            // vFluc2 is the difference between the direct coulomb integral
+            // and the normal 1/r-like  interaction between point charges.
+            coulInt = J2->getValueAt( *(idat.rij) );
+            vFluc2 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
+          } else {
+            vFluc2 = 0.0;
+          }
+          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) / q_j;
+        }
+          
+ 
       }
 
       if (j_is_Dipole) {
@@ -676,6 +724,9 @@ namespace OpenMD {
           duduz_j += -preSw * pot_term * rhat;
 
         }
+        if (i_is_Fluctuating) {
+          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
+        }
       }
 
       if (j_is_Quadrupole) {
@@ -728,6 +779,10 @@ namespace OpenMD {
         dudux_j += preSw * qxx_j * cx_j * rhatdot2;
         duduy_j += preSw * qyy_j * cy_j * rhatdot2;
         duduz_j += preSw * qzz_j * cz_j * rhatdot2;
+        if (i_is_Fluctuating) {
+          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
+        }
+
       }
     }
     
@@ -803,7 +858,12 @@ namespace OpenMD {
           // calculate derivatives for the forces and torques
           dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
           duduz_i += preSw * pot_term * rhat;
+        }
+
+        if (j_is_Fluctuating) {
+          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
         }
+
       }
 
       if (j_is_Dipole) {
@@ -954,6 +1014,11 @@ namespace OpenMD {
         dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
         duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
         duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
+
+        if (j_is_Fluctuating) {
+          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
+        }
+
       }
     }
 
@@ -992,7 +1057,6 @@ namespace OpenMD {
         *(idat.t2) -= cross(uz_j, indirect_duduz_j);
     }
 
-
     return;
   }