--- trunk/src/nonbonded/Electrostatic.cpp	2013/08/23 15:59:23	1933
+++ trunk/src/nonbonded/Electrostatic.cpp	2014/04/14 18:32:51	1981
@@ -40,6 +40,10 @@
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 
+#ifdef IS_MPI
+#include <mpi.h>
+#endif
+
 #include <stdio.h>
 #include <string.h>
 
@@ -57,9 +61,7 @@
 #include "math/erfc.hpp"
 #include "math/SquareMatrix.hpp"
 #include "primitives/Molecule.hpp"
-#ifdef IS_MPI
-#include <mpi.h>
-#endif
+#include "flucq/FluctuatingChargeForces.hpp"
 
 namespace OpenMD {
   
@@ -69,8 +71,20 @@ namespace OpenMD {
                                   haveDampingAlpha_(false), 
                                   haveDielectric_(false),
                                   haveElectroSplines_(false)
-  {}
+  {
+    flucQ_ = new FluctuatingChargeForces(info_);
+  }
   
+  void Electrostatic::setForceField(ForceField *ff) {
+    forceField_ = ff;
+    flucQ_->setForceField(forceField_);
+  }
+
+  void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
+    simTypes_ = simtypes;
+    flucQ_->setSimulatedAtomTypes(simTypes_);
+  }
+
   void Electrostatic::initialize() { 
     
     Globals* simParams_ = info_->getSimParams();
@@ -1139,9 +1153,11 @@ namespace OpenMD {
         
     if (i_is_Fluctuating) {
       C_a += *(sdat.flucQ);
-      *(sdat.flucQfrc) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
-      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * 
-        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
+
+      flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),  
+                                 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY], 
+                                 *(sdat.flucQfrc) );
+
     }
 
     switch (summationMethod_) {
@@ -1244,7 +1260,6 @@ namespace OpenMD {
     vector<vector<RealType> > els;
     vector<vector<RealType> > ems;
     vector<vector<RealType> > ens;
-
     
     int nMax = info_->getNAtoms();
     
@@ -1267,7 +1282,6 @@ namespace OpenMD {
     Vector3d t( 2.0 * M_PI );
     t.Vdiv(t, box);
 
-    
     SimInfo::MoleculeIterator mi;
     Molecule::AtomIterator ai;
     int i;
@@ -1440,7 +1454,7 @@ namespace OpenMD {
                   Q = atom->getQuadrupole() * mPoleConverter; 
                   Qk = Q * kVec;                  
                   qk = dot(kVec, Qk);
-                  qxk[i] = cross(kVec, Qk);
+                  qxk[i] = -cross(kVec, Qk);
                   qkc[i] = qk * ckr[i];
                   qks[i] = qk * skr[i];
                 }              
@@ -1457,18 +1471,18 @@ namespace OpenMD {
             qkss = std::accumulate(qks.begin(),qks.end(),0.0);
             
 #ifdef IS_MPI
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
-            MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE, 
-                                      MPI::SUM);
+            MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE, 
+                          MPI_SUM, MPI_COMM_WORLD);
 #endif        
             
             // Accumulate potential energy and virial contribution:
@@ -1502,7 +1516,11 @@ namespace OpenMD {
                                             +skr[i]*(ckss+dkcs-qkss));
                
                 atom->addFrc( 4.0 * rvol * qfrc * kVec );
-                
+
+                if (atom->isFluctuatingCharge()) { 
+                  atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
+                }
+                  
                 if (data.is_Dipole) {
                   atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
                 }