--- branches/development/src/nonbonded/Electrostatic.cpp	2012/05/24 01:29:59	1718
+++ branches/development/src/nonbonded/Electrostatic.cpp	2012/05/24 14:17:42	1721
@@ -48,6 +48,7 @@
 #include "utils/simError.h"
 #include "types/NonBondedInteractionType.hpp"
 #include "types/FixedChargeAdapter.hpp"
+#include "types/FluctuatingChargeAdapter.hpp"
 #include "types/MultipoleAdapter.hpp"
 #include "io/Globals.hpp"
 #include "nonbonded/SlaterIntegrals.hpp"
@@ -288,7 +289,7 @@ namespace OpenMD {
 
     if (fca.isFixedCharge()) {
       electrostaticAtomData.is_Charge = true;
-      electrostaticAtomData.charge = fca.getCharge();
+      electrostaticAtomData.fixedCharge = fca.getCharge();
     }
 
     MultipoleAdapter ma = MultipoleAdapter(atomType);
@@ -315,11 +316,13 @@ namespace OpenMD {
     FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
 
     if (fqa.isFluctuatingCharge()) {
-      electrostaticAtomData.is_FluctuatingCharge = true;
-      electrostaticAtomData.electronegativity = fca.getElectronegativity();
-      electrostaticAtomData.hardness = fca.getHardness();
-      electrostaticAtomData.slaterN = fca.getSlaterN();
-      electrostaticAtomData.slaterZeta = fca.getSlaterZeta();
+      electrostaticAtomData.is_Fluctuating = true;
+      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
+      electrostaticAtomData.hardness = fqa.getHardness();
+      electrostaticAtomData.slaterN = fqa.getSlaterN();
+      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
+    } else {
+      electrostaticAtomData.is_Fluctuating = false;
     }
 
     pair<map<int,AtomType*>::iterator,bool> ret;    
@@ -341,19 +344,20 @@ namespace OpenMD {
     map<AtomType*, ElectrostaticAtomData>::iterator it;
     for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
       AtomType* atype2 = (*it).first;
-      
-      if ((*it).is_FluctuatingCharge && electrostaticAtomData.is_FluctuatingCharge) {
+      ElectrostaticAtomData eaData2 = (*it).second;
+      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
         
         RealType a = electrostaticAtomData.slaterZeta;
-        RealType b = (*it).slaterZeta;
+        RealType b = eaData2.slaterZeta;
         int m = electrostaticAtomData.slaterN;
-        int n = (*it).slaterN;
+        int n = eaData2.slaterN;
 
         // Create the spline of the coulombic integral for s-type
         // Slater orbitals.  Add a 2 angstrom safety window to deal
         // with cutoffGroups that have charged atoms longer than the
         // cutoffRadius away from each other.
 
+        RealType rval;
         RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
         vector<RealType> rvals;
         vector<RealType> J1vals;
@@ -362,12 +366,13 @@ 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 ) );
         }
 
-        CubicSpline J1 = new CubicSpline();
+        CubicSpline* J1 = new CubicSpline();
         J1->addPoints(rvals, J1vals);
-        CubicSpline J2 = new CubicSpline();
+        CubicSpline* J2 = new CubicSpline();
         J2->addPoints(rvals, J2vals);
         
         pair<AtomType*, AtomType*> key1, key2;
@@ -446,6 +451,10 @@ namespace OpenMD {
 
     pair<RealType, RealType> res;
     
+    // splines for coulomb integrals
+    CubicSpline* J1;
+    CubicSpline* J2;
+    
     if (!initialized_) initialize();
     
     ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
@@ -462,14 +471,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.charge;
+      q_i = data1.fixedCharge;
+
+      if (i_is_Fluctuating) {
+        q_i += *(idat.flucQ1);
+      }
+      
       if (idat.excluded) {
         *(idat.skippedCharge2) += q_i;
       }
@@ -507,7 +523,11 @@ namespace OpenMD {
     }
 
     if (j_is_Charge) {
-      q_j = data2.charge;
+      q_j = data2.fixedCharge;
+
+      if (i_is_Fluctuating) 
+        q_j += *(idat.flucQ2);
+
       if (idat.excluded) {
         *(idat.skippedCharge1) += q_j;
       }
@@ -545,6 +565,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;
     
@@ -596,7 +621,25 @@ namespace OpenMD {
 
           vterm = preVal * riji * erfcVal;           
           dudr  = -  *(idat.sw)  * preVal * c2;
+          
+        }
 
+        
+        if (i_is_Fluctuating) {
+          if (!idat.excluded)
+            *(idat.dVdFQ1) += *(idat.sw) * vterm / q_i;
+          else {
+            res = J1->getValueAndDerivativeAt( *(idat.rij) );
+            *(idat.dVdFQ1) += pre11_ * res.first * q_j;
+          }
+        }
+        if (j_is_Fluctuating) {
+          if (!idat.excluded)
+            *(idat.dVdFQ2) += *(idat.sw) * vterm / q_j;
+          else {
+            res = J2->getValueAndDerivativeAt( *(idat.rij) );
+            *(idat.dVdFQ2) += pre11_ * res.first * q_i;
+          }
         }
 
         vpair += vterm;
@@ -1020,7 +1063,7 @@ namespace OpenMD {
       }
     } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
       if (i_is_Charge) {        
-        chg1 = data.charge;
+        chg1 = data.fixedCharge;
         if (screeningMethod_ == DAMPED) {
           self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
         } else {