--- trunk/src/brains/ForceManager.cpp	2014/04/17 19:07:31	1987
+++ trunk/src/brains/ForceManager.cpp	2014/12/11 16:16:43	2047
@@ -57,7 +57,8 @@
 #include "primitives/Torsion.hpp"
 #include "primitives/Inversion.hpp"
 #include "nonbonded/NonBondedInteraction.hpp"
-#include "perturbations/ElectricField.hpp"
+#include "perturbations/UniformField.hpp"
+#include "perturbations/UniformGradient.hpp"
 #include "parallel/ForceMatrixDecomposition.hpp"
 
 #include <cstdio>
@@ -394,9 +395,6 @@ namespace OpenMD {
     switcher_->setSwitch(rSwitch_, rCut_);
   }
 
-
-
-  
   void ForceManager::initialize() {
 
     if (!info_->isTopologyDone()) {
@@ -405,9 +403,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 +415,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,11 +446,17 @@ 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);
     }
-
+    if (info_->getSimParams()->haveUniformGradientStrength() ||
+        info_->getSimParams()->haveUniformGradientDirection1() ||
+        info_->getSimParams()->haveUniformGradientDirection2() ) {
+      UniformGradient* eGrad = new UniformGradient(info_);
+      perturbations_.push_back(eGrad);
+    }
+    
     usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
     
     fDecomp_->distributeInitialData();
@@ -716,6 +721,9 @@ namespace OpenMD {
     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;
@@ -731,13 +739,17 @@ 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.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE ||
+                         cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
     idat.doParticlePot = doParticlePot_;
     idat.doElectricField = doElectricField_;
+    idat.doSitePotential = doSitePotential_;
     sdat.doParticlePot = doParticlePot_;
     
     loopEnd = PAIR_LOOP;
@@ -758,7 +770,7 @@ namespace OpenMD {
       }
 
       for (vector<pair<int, int> >::iterator it = neighborList_.begin(); 
-             it != neighborList_.end(); ++it) {
+           it != neighborList_.end(); ++it) {
                 
         cg1 = (*it).first;
         cg2 = (*it).second;
@@ -777,6 +789,8 @@ namespace OpenMD {
             fij.zero();
             eField1.zero();
             eField2.zero();
+            sPot1 = 0.0;
+            sPot2 = 0.0;
           }
           
           in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, 
@@ -949,7 +963,7 @@ namespace OpenMD {
     curSnapshot->setLongRangePotential(longRangePotential);
     
     curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
-                                         *(fDecomp_->getExcludedPotential()));
+                                       *(fDecomp_->getExcludedPotential()));
 
   }
 
@@ -985,37 +999,37 @@ namespace OpenMD {
     
     if (info_->getSimParams()->getUseLongRangeCorrections()) {
       /*
-      RealType vol = curSnapshot->getVolume();
-      RealType Elrc(0.0);
-      RealType Wlrc(0.0);
+        RealType vol = curSnapshot->getVolume();
+        RealType Elrc(0.0);
+        RealType Wlrc(0.0);
 
-      set<AtomType*>::iterator i;
-      set<AtomType*>::iterator j;
+        set<AtomType*>::iterator i;
+        set<AtomType*>::iterator j;
     
-      RealType n_i, n_j;
-      RealType rho_i, rho_j;
-      pair<RealType, RealType> LRI;
+        RealType n_i, n_j;
+        RealType rho_i, rho_j;
+        pair<RealType, RealType> LRI;
       
-      for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
+        for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
         n_i = RealType(info_->getGlobalCountOfType(*i));
         rho_i = n_i /  vol;
         for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
-          n_j = RealType(info_->getGlobalCountOfType(*j));
-          rho_j = n_j / vol;
+        n_j = RealType(info_->getGlobalCountOfType(*j));
+        rho_j = n_j / vol;
           
-          LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
+        LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
 
-          Elrc += n_i   * rho_j * LRI.first;
-          Wlrc -= rho_i * rho_j * LRI.second;
+        Elrc += n_i   * rho_j * LRI.first;
+        Wlrc -= rho_i * rho_j * LRI.second;
         }
-      }
-      Elrc *= 2.0 * NumericConstant::PI;
-      Wlrc *= 2.0 * NumericConstant::PI;
+        }
+        Elrc *= 2.0 * NumericConstant::PI;
+        Wlrc *= 2.0 * NumericConstant::PI;
 
-      RealType lrp = curSnapshot->getLongRangePotential();
-      curSnapshot->setLongRangePotential(lrp + Elrc);
-      stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
-      curSnapshot->setStressTensor(stressTensor);
+        RealType lrp = curSnapshot->getLongRangePotential();
+        curSnapshot->setLongRangePotential(lrp + Elrc);
+        stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
+        curSnapshot->setStressTensor(stressTensor);
       */
      
     }