--- trunk/src/nonbonded/Electrostatic.cpp	2013/08/21 16:52:56	1932
+++ trunk/src/nonbonded/Electrostatic.cpp	2014/02/26 14:14:50	1969
@@ -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();
@@ -987,7 +1001,6 @@ namespace OpenMD {
         F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
         Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
         Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
-
         // Even if we excluded this pair from direct interactions, we
         // still have the reaction-field-mediated dipole-dipole
         // interaction:
@@ -1047,7 +1060,7 @@ namespace OpenMD {
         trQaQb = QaQb.trace();
         rQaQb = rhat * QaQb;
         QaQbr = QaQb * rhat;
-        QaxQb = cross(Q_a, Q_b);
+        QaxQb = mCross(Q_a, Q_b);
         rQaQbr = dot(rQa, Qbr);
         rQaxQbr = cross(rQa, Qbr);
         
@@ -1078,7 +1091,6 @@ namespace OpenMD {
         //             + 4.0 * cross(rhat, QbQar) 
 
         Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
-
       }
     }
 
@@ -1141,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_) {
@@ -1441,8 +1455,8 @@ namespace OpenMD {
                 if (data.is_Quadrupole) {
                   Q = atom->getQuadrupole() * mPoleConverter; 
                   Qk = Q * kVec;                  
-                  qk = dot(Qk, kVec);
-                  qxk[i] = cross(Qk, kVec);
+                  qk = dot(kVec, Qk);
+                  qxk[i] = -cross(kVec, Qk);
                   qkc[i] = qk * ckr[i];
                   qks[i] = qk * skr[i];
                 }              
@@ -1459,18 +1473,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:
@@ -1501,7 +1515,7 @@ namespace OpenMD {
                 RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
                                          -ckr[i]*(ckss+dkcs-qkss));
                 RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
-                                             +skr[i]*(ckss+dkcs-qkss));
+                                            +skr[i]*(ckss+dkcs-qkss));
                
                 atom->addFrc( 4.0 * rvol * qfrc * kVec );