--- trunk/src/brains/ForceManager.cpp	2013/08/05 16:13:46	1923
+++ trunk/src/brains/ForceManager.cpp	2014/10/31 18:57:19	2032
@@ -57,7 +57,7 @@
 #include "primitives/Torsion.hpp"
 #include "primitives/Inversion.hpp"
 #include "nonbonded/NonBondedInteraction.hpp"
-#include "perturbations/ElectricField.hpp"
+#include "perturbations/UniformField.hpp"
 #include "parallel/ForceMatrixDecomposition.hpp"
 
 #include <cstdio>
@@ -405,9 +405,9 @@ namespace OpenMD {
       interactionMan_->setSimInfo(info_);
       interactionMan_->initialize();
 
-      // We want to delay the cutoffs until after the interaction
-      // manager has set up the atom-atom interactions so that we can
-      // query them for suggested cutoff values
+      //! We want to delay the cutoffs until after the interaction
+      //! manager has set up the atom-atom interactions so that we can
+      //! query them for suggested cutoff values
       setupCutoffs();
 
       info_->prepareTopology();      
@@ -417,19 +417,20 @@ namespace OpenMD {
       if (doHeatFlux_) doParticlePot_ = true;
 
       doElectricField_ = info_->getSimParams()->getOutputElectricField();
+      doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
    
     }
 
     ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
     
-    // Force fields can set options on how to scale van der Waals and
-    // electrostatic interactions for atoms connected via bonds, bends
-    // and torsions in this case the topological distance between
-    // atoms is:
-    // 0 = topologically unconnected
-    // 1 = bonded together 
-    // 2 = connected via a bend
-    // 3 = connected via a torsion
+    //! Force fields can set options on how to scale van der Waals and
+    //! electrostatic interactions for atoms connected via bonds, bends
+    //! and torsions in this case the topological distance between
+    //! atoms is:
+    //! 0 = topologically unconnected
+    //! 1 = bonded together 
+    //! 2 = connected via a bend
+    //! 3 = connected via a torsion
     
     vdwScale_.reserve(4);
     fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
@@ -447,8 +448,8 @@ namespace OpenMD {
     electrostaticScale_[2] = fopts.getelectrostatic13scale();
     electrostaticScale_[3] = fopts.getelectrostatic14scale();    
     
-    if (info_->getSimParams()->haveElectricField()) {
-      ElectricField* eField = new ElectricField(info_);
+    if (info_->getSimParams()->haveUniformField()) {
+      UniformField* eField = new UniformField(info_);
       perturbations_.push_back(eField);
     }
 
@@ -637,14 +638,15 @@ namespace OpenMD {
     // Collect from all nodes.  This should eventually be moved into a
     // SystemDecomposition, but this is a better place than in
     // Thermo to do the collection.
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, 
-                              MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, 
-                              MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, 
-                              MPI::REALTYPE, MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, 
-                              MPI::REALTYPE, MPI::SUM);
+
+    MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, 
+                  MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, 
+                  MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
 #endif
 
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
@@ -666,6 +668,7 @@ namespace OpenMD {
     DataStorage* config = &(curSnapshot->atomData);
     DataStorage* cgConfig = &(curSnapshot->cgData);
 
+
     //calculate the center of mass of cutoff group
 
     SimInfo::MoleculeIterator mi;
@@ -709,11 +712,14 @@ namespace OpenMD {
     RealType dVdFQ1(0.0);
     RealType dVdFQ2(0.0);
     potVec longRangePotential(0.0);
-    potVec reciprocalPotential(0.0);
+    RealType reciprocalPotential(0.0);
     potVec workPot(0.0);
     potVec exPot(0.0);
     Vector3d eField1(0.0);
     Vector3d eField2(0.0);
+    RealType sPot1(0.0);
+    RealType sPot2(0.0);
+                   
     vector<int>::iterator ia, jb;
 
     int loopStart, loopEnd;
@@ -729,13 +735,16 @@ namespace OpenMD {
     idat.dVdFQ1 = &dVdFQ1;
     idat.dVdFQ2 = &dVdFQ2;
     idat.eField1 = &eField1;
-    idat.eField2 = &eField2;   
+    idat.eField2 = &eField2; 
+    idat.sPot1 = &sPot1;
+    idat.sPot2 = &sPot2;
     idat.f1 = &f1;
     idat.sw = &sw;
     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
     idat.doParticlePot = doParticlePot_;
     idat.doElectricField = doElectricField_;
+    idat.doSitePotential = doSitePotential_;
     sdat.doParticlePot = doParticlePot_;
     
     loopEnd = PAIR_LOOP;
@@ -770,12 +779,13 @@ namespace OpenMD {
         rgrpsq = d_grp.lengthSquare();
 
         if (rgrpsq < rCutSq) {
-
           if (iLoop == PAIR_LOOP) {
             vij = 0.0;
             fij.zero();
             eField1.zero();
             eField2.zero();
+            sPot1 = 0.0;
+            sPot2 = 0.0;
           }
           
           in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, 
@@ -827,7 +837,7 @@ namespace OpenMD {
                
                 r = sqrt( *(idat.r2) );
                 idat.rij = &r;
-               
+
                 if (iLoop == PREPAIR_LOOP) {
                   interactionMan_->doPrePair(idat);
                 } else {
@@ -928,6 +938,8 @@ namespace OpenMD {
     fDecomp_->collectData();
     if (cutoffMethod_ == EWALD_FULL) {
       interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
+
+      curSnapshot->setReciprocalPotential(reciprocalPotential);
     }
         
     if (info_->requiresSelfCorrection()) {
@@ -941,7 +953,7 @@ namespace OpenMD {
     fDecomp_->collectSelfData();
 
     longRangePotential = *(fDecomp_->getEmbeddingPotential()) + 
-      *(fDecomp_->getPairwisePotential()) + reciprocalPotential;
+      *(fDecomp_->getPairwisePotential());
 
     curSnapshot->setLongRangePotential(longRangePotential);
     
@@ -975,8 +987,8 @@ namespace OpenMD {
     }
     
 #ifdef IS_MPI
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, 
-                              MPI::REALTYPE, MPI::SUM);
+    MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
 #endif
     curSnapshot->setStressTensor(stressTensor);