--- branches/development/src/nonbonded/Electrostatic.cpp	2012/05/18 21:44:02	1710
+++ branches/development/src/nonbonded/Electrostatic.cpp	2012/05/24 01:48:29	1720
@@ -48,8 +48,12 @@
 #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"
+#include "utils/PhysicalConstants.hpp"
+
 
 namespace OpenMD {
   
@@ -285,7 +289,7 @@ namespace OpenMD {
 
     if (fca.isFixedCharge()) {
       electrostaticAtomData.is_Charge = true;
-      electrostaticAtomData.charge = fca.getCharge();
+      electrostaticAtomData.fixedCharge = fca.getCharge();
     }
 
     MultipoleAdapter ma = MultipoleAdapter(atomType);
@@ -309,7 +313,16 @@ namespace OpenMD {
       }
     }
     
+    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
 
+    if (fqa.isFluctuatingCharge()) {
+      electrostaticAtomData.is_Fluctuating = true;
+      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
+      electrostaticAtomData.hardness = fqa.getHardness();
+      electrostaticAtomData.slaterN = fqa.getSlaterN();
+      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
+    }
+
     pair<map<int,AtomType*>::iterator,bool> ret;    
     ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
                                                         atomType) );
@@ -322,7 +335,52 @@ namespace OpenMD {
       simError();         
     }
     
-    ElectrostaticMap[atomType] = electrostaticAtomData;    
+    ElectrostaticMap[atomType] = electrostaticAtomData;   
+
+    // Now, iterate over all known types and add to the mixing map:
+    
+    map<AtomType*, ElectrostaticAtomData>::iterator it;
+    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
+      AtomType* atype2 = (*it).first;
+      ElectrostaticAtomData eaData2 = (*it).second;
+      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
+        
+        RealType a = electrostaticAtomData.slaterZeta;
+        RealType b = eaData2.slaterZeta;
+        int m = electrostaticAtomData.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;
+        vector<RealType> J2vals;
+        for (int i = 0; i < np_; i++) {
+          rval = RealType(i) * dr;
+          rvals.push_back(rval);
+          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
+          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
+        }
+
+        CubicSpline* J1 = new CubicSpline();
+        J1->addPoints(rvals, J1vals);
+        CubicSpline* J2 = new CubicSpline();
+        J2->addPoints(rvals, J2vals);
+        
+        pair<AtomType*, AtomType*> key1, key2;
+        key1 = make_pair(atomType, atype2);
+        key2 = make_pair(atype2, atomType);
+        
+        Jij[key1] = J1;
+        Jij[key2] = J2;
+      }
+    }
+ 
     return;
   }
   
@@ -413,7 +471,7 @@ namespace OpenMD {
     bool j_is_Quadrupole = data2.is_Quadrupole;
     
     if (i_is_Charge) {
-      q_i = data1.charge;
+      q_i = data1.fixedCharge;
       if (idat.excluded) {
         *(idat.skippedCharge2) += q_i;
       }
@@ -451,7 +509,7 @@ namespace OpenMD {
     }
 
     if (j_is_Charge) {
-      q_j = data2.charge;
+      q_j = data2.fixedCharge;
       if (idat.excluded) {
         *(idat.skippedCharge1) += q_j;
       }
@@ -964,7 +1022,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 {