diff --git a/Source/FieldSolver/FerroE/Eff_Field_Grad.H b/Source/FieldSolver/FerroE/Eff_Field_Grad.H
new file mode 100644
index 000000000..d1e711614
--- /dev/null
+++ b/Source/FieldSolver/FerroE/Eff_Field_Grad.H
@@ -0,0 +1,106 @@
+/*
+  Derivative algorithms for calculating the gradient energy terms
+*/
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real UpwardDx (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i+1,j,k) - F(i,j,k))/(dx[0]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DownwardDx (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i,j,k) - F(i-1,j,k))/(dx[0]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real UpwardDy (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i,j+1,k) - F(i,j,k))/(dx[1]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DownwardDy (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i,j,k) - F(i,j-1,k))/(dx[1]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real UpwardDz (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i,j,k+1) - F(i,j,k))/(dx[2]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DownwardDz (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+   return (F(i,j,k) - F(i,j,k-1))/(dx[2]);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DoubleDx (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> const& fe) {
+
+   if (fe(i,j,k) == 1 && fe(i+1,j,k) == 0) {//upper end of fe
+
+      return (0.0 - DownwardDx(F, i, j, k, dx))/(dx[0]);
+
+   } else if (fe(i,j,k) == 1 && fe(i-1,j,k) == 0) { //lower end of fe
+
+      return (UpwardDx(F, i, j, k, dx) - 0.0)/(dx[0]);
+
+   } else {//bulk
+
+      return (UpwardDx(F, i, j, k, dx) - DownwardDx(F, i, j, k, dx))/(dx[0]);
+
+   }
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DoubleDy (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> const& fe) {
+   
+   if (fe(i,j,k) == 1 && fe(i,j+1,k) == 0) {//upper end of fe
+
+      return (0.0 - DownwardDy(F, i, j, k, dx))/(dx[1]);
+
+   } else if (fe(i,j,k) == 1 && fe(i,j-1,k) == 0) { //lower end of fe
+
+      return (UpwardDy(F, i, j, k, dx) - 0.0)/(dx[1]);
+
+   } else {//bulk
+
+      return (UpwardDy(F, i, j, k, dx) - DownwardDy(F, i, j, k, dx))/(dx[1]);
+
+   }
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+static amrex::Real DoubleDz (
+   amrex::Array4<amrex::Real> const& F,
+   int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx, amrex::Array4<amrex::Real> const& fe) {
+   
+   if (fe(i,j,k) == 1 && fe(i,j,k+1) == 0) {//upper end of fe
+
+      return (0.0 - DownwardDz(F, i, j, k, dx))/(dx[2]);
+
+   } else if (fe(i,j,k) == 1 && fe(i,j,k-1) == 0) { //lower end of fe
+
+      return (UpwardDz(F, i, j, k, dx) - 0.0)/(dx[2]);
+
+   } else {//bulk
+
+      return (UpwardDz(F, i, j, k, dx) - DownwardDz(F, i, j, k, dx))/(dx[2]);
+
+   }
+}
+
diff --git a/Source/FieldSolver/FerroE/FerroE.H b/Source/FieldSolver/FerroE/FerroE.H
index fcdcb0778..1b043288f 100644
--- a/Source/FieldSolver/FerroE/FerroE.H
+++ b/Source/FieldSolver/FerroE/FerroE.H
@@ -56,6 +56,8 @@ public:
      const amrex::Real alpha_1122 = 1.637e10;
      const amrex::Real alpha_1123 = 1.367e10;
 
+     //Gradient energy coefficients
+     const amrex::Real G_11 = 5.1e-10;
 };
 
 
diff --git a/Source/FieldSolver/FerroE/FerroE.cpp b/Source/FieldSolver/FerroE/FerroE.cpp
index 64df20e9e..3303da6bd 100644
--- a/Source/FieldSolver/FerroE/FerroE.cpp
+++ b/Source/FieldSolver/FerroE/FerroE.cpp
@@ -1,5 +1,6 @@
 #include "FerroE.H"
 #include "Eff_Field_Landau.H"
+#include "Eff_Field_grad.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 #include "Utils/WarpXUtil.H"
 #include "WarpX.H"
@@ -189,7 +190,9 @@ FerroE::EvolveP (amrex::Real dt)
      amrex::Print() << " evolve P \n";
      auto & warpx = WarpX::GetInstance();
      int include_Landau = warpx.include_Landau;
+     int include_grad = warpx.include_grad;
      const int lev = 0;
+     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = warpx.Geom(lev).CellSizeArray();
 
      //Px, Py, and Pz have 2 components each. Px(i,j,k,0) is Px and Px(i,j,k,1) is dPx/dt and so on..
      amrex::MultiFab * Px = warpx.get_pointer_polarization_fp(lev, 0);
@@ -224,10 +227,14 @@ FerroE::EvolveP (amrex::Real dt)
                    Ex_eff += compute_ex_Landau(Px_arr(i,j,k,0), Py_arr(i,j,k,0), Pz_arr(i,j,k,0));
                 }
 
-                  //get dPx/dt using numerical integration
-                  update_v(Px_arr(i,j,k,1), dt, Ex_eff, gamma, mu);
-                  //get Px 
-                  Px_arr(i,j,k,0) += dt*Px_arr(i,j,k,1);
+                if (include_grad == 1){
+                   Ex_eff += G_11*DoubleDx(Px_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Px_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Px_arr, i, j, k, dx, fe_arr);
+                }
+
+                //get dPx/dt using numerical integration
+                update_v(Px_arr(i,j,k,1), dt, Ex_eff, gamma, mu);
+                //get Px 
+                Px_arr(i,j,k,0) += dt*Px_arr(i,j,k,1);
             }
         },
         [=] AMREX_GPU_DEVICE (int i, int j, int k) {
@@ -239,10 +246,14 @@ FerroE::EvolveP (amrex::Real dt)
                    Ey_eff += compute_ey_Landau(Px_arr(i,j,k,0), Py_arr(i,j,k,0), Pz_arr(i,j,k,0));
                 }
                   
-                  //get dPy/dt using numerical integration
-                  update_v(Py_arr(i,j,k,1), dt, Ey_eff, gamma, mu);
-                  //get Py 
-                  Py_arr(i,j,k,0) += dt*Py_arr(i,j,k,1);
+                if (include_grad == 1){
+                   Ey_eff += G_11*DoubleDx(Py_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Py_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Py_arr, i, j, k, dx, fe_arr);
+                }
+
+                //get dPy/dt using numerical integration
+                update_v(Py_arr(i,j,k,1), dt, Ey_eff, gamma, mu);
+                //get Py 
+                Py_arr(i,j,k,0) += dt*Py_arr(i,j,k,1);
             }
         },
         [=] AMREX_GPU_DEVICE (int i, int j, int k) {
@@ -254,10 +265,14 @@ FerroE::EvolveP (amrex::Real dt)
                    Ez_eff += compute_ez_Landau(Px_arr(i,j,k,0), Py_arr(i,j,k,0), Pz_arr(i,j,k,0));
                 }
                   
-                  //get dPz/dt using numerical integration
-                  update_v(Pz_arr(i,j,k,1), dt, Ez_eff, gamma, mu);
-                  //get Pz 
-                  Pz_arr(i,j,k,0) += dt*Pz_arr(i,j,k,1);
+                if (include_grad == 1){
+                   Ez_eff += G_11*DoubleDx(Pz_arr, i, j, k, dx, fe_arr) + G_11*DoubleDy(Pz_arr, i, j, k, dx, fe_arr) + G_11*DoubleDz(Pz_arr, i, j, k, dx, fe_arr);
+                }
+
+                //get dPz/dt using numerical integration
+                update_v(Pz_arr(i,j,k,1), dt, Ez_eff, gamma, mu);
+                //get Pz 
+                Pz_arr(i,j,k,0) += dt*Pz_arr(i,j,k,1);
             }
         }
     );
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 2cfaca06b..e2b029243 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -367,6 +367,7 @@ public:
 #endif
 #ifdef WARPX_FERROE
     int include_Landau = 1;
+    int include_grad = 0;
 #endif
     //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid
     //! using finite centering of order given by #current_centering_nox, #current_centering_noy,