diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H
index 80ee376..94c6fba 100644
--- a/Source/Solver/DerivativeAlgorithm.H
+++ b/Source/Solver/DerivativeAlgorithm.H
@@ -6,9 +6,9 @@ using namespace FerroX;
   * Perform first derivative dphi/dz */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DphiDz (
-    amrex::Array4<amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& F,
     amrex::Real const z_hi, amrex::Real const z_lo, 
-    int const i, int const j, int const k,  amrex::GpuArray<amrex::Real, 3> dx, 
+    int const i, int const j, int const k,  amrex::GpuArray<amrex::Real, 3> const dx, 
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>const& prob_lo,
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>const& prob_hi) {
 
@@ -25,8 +25,8 @@ using namespace FerroX;
   * Perform first derivative dF/dx */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DFDx (
-    amrex::Array4<amrex::Real> const& F,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+    amrex::Array4<const amrex::Real> const& F,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx) {
     return (F(i+1,j,k) - F(i-1,j,k))/(2.*dx[0]);
  }
 
@@ -34,8 +34,8 @@ using namespace FerroX;
   * Perform first derivative dF/dy */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DFDy (
-    amrex::Array4<amrex::Real> const& F,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
+    amrex::Array4<const amrex::Real> const& F,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx) {
     return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]);
  }
 
@@ -43,9 +43,9 @@ using namespace FerroX;
   * Perform first derivative dP/dx */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DPDx (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
 ) {
     if (mask(i-1,j,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
       
@@ -111,9 +111,9 @@ using namespace FerroX;
   * Perform first derivative dP/dy */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DPDy (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
 ) {
 
     if (mask(i,j-1,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
@@ -180,9 +180,9 @@ using namespace FerroX;
   * Perform first derivative dP/dz */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DPDz (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
 ) {
 
     if (mask(i,j,k-1) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
@@ -255,9 +255,9 @@ using namespace FerroX;
   * Perform double derivative (d^2)P/dx^2 */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DoubleDPDx (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const dx
     ) {
         
     if (mask(i-1,j,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
@@ -325,9 +325,9 @@ using namespace FerroX;
   * Perform double derivative (d^2)P/dy^2 */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DoubleDPDy (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
     ) {
         
     if (mask(i,j-1,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
@@ -394,9 +394,9 @@ using namespace FerroX;
   * Perform double derivative (d^2)P/dz^2 */
  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
  static amrex::Real DoubleDPDz (
-    amrex::Array4<amrex::Real> const& F,
-    amrex::Array4<amrex::Real> const& mask,
-    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx
+    amrex::Array4<const amrex::Real> const& F,
+    amrex::Array4<const amrex::Real> const& mask,
+    int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx
     ) {
         
     if (mask(i,j,k-1) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary
@@ -469,9 +469,9 @@ using namespace FerroX;
 /**
   * Perform double derivative (d^2)P/dxdy */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-static amrex::Real DoubleDPDxDy (amrex::Array4<amrex::Real> const& F,
-                               amrex::Array4<amrex::Real> const& mask,
-                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
+static amrex::Real DoubleDPDxDy (amrex::Array4<const amrex::Real> const& F,
+                               amrex::Array4<const amrex::Real> const& mask,
+                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
 {
       return (DPDy(F, mask, i+1, j, k, dx) - DPDy(F, mask, i-1, j, k, dx)) / 2. /dx[0]; 
 }
@@ -479,9 +479,9 @@ static amrex::Real DoubleDPDxDy (amrex::Array4<amrex::Real> const& F,
 /**
   * Perform double derivative (d^2)P/dxdz */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-static amrex::Real DoubleDPDxDz (amrex::Array4<amrex::Real> const& F,
-                               amrex::Array4<amrex::Real> const& mask,
-                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
+static amrex::Real DoubleDPDxDz (amrex::Array4<const amrex::Real> const& F,
+                               amrex::Array4<const amrex::Real> const& mask,
+                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
 {
       return (DPDz(F, mask, i+1, j, k, dx) - DPDz(F, mask, i-1, j, k, dx)) / 2. /dx[0]; 
 }
@@ -489,9 +489,9 @@ static amrex::Real DoubleDPDxDz (amrex::Array4<amrex::Real> const& F,
 /**
   * Perform double derivative (d^2)P/dydz */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-static amrex::Real DoubleDPDyDz (amrex::Array4<amrex::Real> const& F,
-                               amrex::Array4<amrex::Real> const& mask,
-                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx)
+static amrex::Real DoubleDPDyDz (amrex::Array4<const amrex::Real> const& F,
+                               amrex::Array4<const amrex::Real> const& mask,
+                               int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> const& dx)
 {
       return (DPDz(F, mask, i, j+1, k, dx) - DPDz(F, mask, i, j-1, k, dx)) / 2. /dx[1];  
 }
diff --git a/Source/Solver/Energy_Calculation.H b/Source/Solver/Energy_Calculation.H
new file mode 100644
index 0000000..f41bed8
--- /dev/null
+++ b/Source/Solver/Energy_Calculation.H
@@ -0,0 +1,29 @@
+#ifndef FERROX_TOTALFREEENERGY_H_
+#define FERROX_TOTALFREEENERGY_H_
+
+#include <AMReX_Array.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_REAL.H>
+
+// Computes the Landau free energy over the ferroelectric region (MaterialMask == 0.0)
+Real ComputeLandauEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                          MultiFab& MaterialMask,
+                         const Geometry& geom,
+                         Real alpha,
+                         Real beta,
+                         Real gamma_coeff);
+
+
+Real ComputeElectrostaticEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                                 Array<MultiFab, AMREX_SPACEDIM>& E,
+                                const MultiFab& MaterialMask,
+                                const Geometry& geom);
+
+
+Real ComputeGradientEnergy(const Array<MultiFab, AMREX_SPACEDIM>& P,
+                            const MultiFab& MaterialMask,
+                            const Geometry& geom,
+                            Real g11,
+                            Real g44);
+#endif
\ No newline at end of file
diff --git a/Source/Solver/Energy_Calculation.cpp b/Source/Solver/Energy_Calculation.cpp
new file mode 100644
index 0000000..ce98041
--- /dev/null
+++ b/Source/Solver/Energy_Calculation.cpp
@@ -0,0 +1,121 @@
+#include <AMReX_MultiFab.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_REAL.H>
+#include <cmath>
+#include "DerivativeAlgorithm.H"
+#include "Energy_Calculation.H"
+
+using namespace amrex;
+
+Real ComputeLandauEnergy(Array<MultiFab, AMREX_SPACEDIM>& P,
+                         MultiFab& MaterialMask,
+                         const Geometry& geom,
+                         Real alpha,
+                         Real beta,
+                         Real gamma_coeff)
+{
+    const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+
+    ReduceOps<ReduceOpSum> reduce_op;
+    ReduceData<Real> reduce_data(reduce_op);
+    using ReduceTuple = typename decltype(reduce_data)::Type;
+
+    for (MFIter mfi(P[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        auto const& px = P[0].array(mfi);
+        auto const& py = P[1].array(mfi);
+        auto const& pz = P[2].array(mfi);
+        auto const& mask = MaterialMask.array(mfi);
+
+        reduce_op.eval(bx, reduce_data,
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
+                if (mask(i,j,k) == 0.0) {
+                    Real P2 = px(i,j,k)*px(i,j,k) + py(i,j,k)*py(i,j,k) + pz(i,j,k)*pz(i,j,k);
+                    Real P4 = P2 * P2;
+                    Real P6 = P4 * P2;
+                    Real f_landau = 0.5 * alpha * P2 + 0.25 * beta * P4 + (1.0/6.0) * gamma_coeff * P6;
+                    return {f_landau * dV};
+                } else {
+                    return {0.0};
+                }
+            });
+    }
+
+    Real landau_energy = amrex::get<0>(reduce_data.value());
+    ParallelDescriptor::ReduceRealSum(landau_energy);
+    return landau_energy;
+}
+
+
+Real ComputeGradientEnergy(const Array<MultiFab, AMREX_SPACEDIM>& P,
+                           const MultiFab& MaterialMask,
+                           const Geometry& geom,
+                           Real g11,
+                           Real g44)
+{
+    const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+    const auto dx = geom.CellSizeArray();  // GpuArray<Real, 3>
+
+    ReduceOps<ReduceOpSum> reduce_op;
+    ReduceData<Real> reduce_data(reduce_op);
+    using ReduceTuple = typename decltype(reduce_data)::Type;
+
+    for (MFIter mfi(P[2], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        auto const& pz = P[2].array(mfi);  // Scalar polarization (Pz)
+        auto const& mask = MaterialMask.array(mfi);
+
+        reduce_op.eval(bx, reduce_data,
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
+                if (mask(i,j,k) == 0.0) {
+                    Real dPdx = DPDx(pz, mask, i, j, k, dx);
+                    Real dPdy = DPDy(pz, mask, i, j, k, dx);
+                    Real dPdz = DPDz(pz, mask, i, j, k, dx);
+
+                    Real f_grad = 0.5 * (g44 * (dPdx*dPdx + dPdy*dPdy) + g11 * dPdz*dPdz);
+                    return {f_grad * dV};
+                } else {
+                    return {0.0};
+                }
+            });
+    }
+
+    Real grad_energy = amrex::get<0>(reduce_data.value());
+    ParallelDescriptor::ReduceRealSum(grad_energy);
+    return grad_energy;
+}
+
+
+
+Real ComputeElectrostaticEnergy( Array<MultiFab, AMREX_SPACEDIM>& P,
+                                 Array<MultiFab, AMREX_SPACEDIM>& E,
+                                const MultiFab& MaterialMask,
+                                const Geometry& geom)
+{
+    const Real dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+
+    ReduceOps<ReduceOpSum> reduce_op;
+    ReduceData<Real> reduce_data(reduce_op);
+    using ReduceTuple = typename decltype(reduce_data)::Type;
+
+    for (MFIter mfi(P[2], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        auto const& pz = P[2].array(mfi);  // scalar polarization: P = Pz
+        auto const& ez = E[2].array(mfi);  // scalar electric field: E = Ez
+        auto const& mask = MaterialMask.array(mfi);
+
+        reduce_op.eval(bx, reduce_data,
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) -> ReduceTuple {
+                if (mask(i,j,k) == 0.0) {
+                    Real f_elec = ez(i,j,k) * pz(i,j,k);  // E · P
+                    return {f_elec * dV};
+                } else {
+                    return {0.0};
+                }
+            });
+    }
+
+    Real electrostatic_energy = amrex::get<0>(reduce_data.value());
+    ParallelDescriptor::ReduceRealSum(electrostatic_energy);
+    return electrostatic_energy;
+}
\ No newline at end of file
diff --git a/Source/Solver/Make.package b/Source/Solver/Make.package
index aeadbd5..5bf6fc1 100644
--- a/Source/Solver/Make.package
+++ b/Source/Solver/Make.package
@@ -2,11 +2,13 @@ CEXE_sources += ElectrostaticSolver.cpp
 CEXE_sources += Initialization.cpp
 CEXE_sources += ChargeDensity.cpp
 CEXE_sources += TotalEnergyDensity.cpp
+CEXE_sources += Energy_Calculation.cpp
 
 CEXE_headers += ElectrostaticSolver.H
 CEXE_headers += Initialization.H
 CEXE_headers += ChargeDensity.H
 CEXE_headers += TotalEnergyDensity.H
+CEXE_headers += Energy_Calculation.H
 
 VPATH_LOCATIONS   += $(CODE_HOME)/Source/Solver
 INCLUDE_LOCATIONS += $(CODE_HOME)/Source/Solver
diff --git a/Source/main.cpp b/Source/main.cpp
index 96b0caa..faee83d 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -18,6 +18,7 @@
 #include "Solver/Initialization.H"
 #include "Solver/ChargeDensity.H"
 #include "Solver/TotalEnergyDensity.H"
+#include "Solver/Energy_Calculation.H"
 #include "Input/BoundaryConditions/BoundaryConditions.H"
 #include "Input/GeometryProperties/GeometryProperties.H"
 #include "Utils/SelectWarpXUtils/WarpXUtil.H"
@@ -258,7 +259,18 @@ void main_main (c_FerroX& rFerroX)
 
             // Calculate E from Phi
             ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);
+            Real Landau_Energy = ComputeLandauEnergy(P_old, MaterialMask, geom, alpha, beta, FerroX::gamma);
+            Real Gradient_Energy = ComputeGradientEnergy(P_old, MaterialMask, geom, g11, g44);
+            Real Electrostatic_Energy = ComputeElectrostaticEnergy(P_old, E, MaterialMask, geom);
+            Real Total_Energy = Landau_Energy + Gradient_Energy + Electrostatic_Energy;
 
+            amrex::Print() << "Landau_Energy value is: " << Landau_Energy << '\n';
+            amrex::Print() << "Gradient_Energy value is: " << Gradient_Energy << '\n';
+            amrex::Print() << "Electrostatic_Energy value is: " << Electrostatic_Energy << '\n';
+            amrex::Print() << "Total_Energy value is: " << Total_Energy << '\n';
+
+
+            
             // compute f^n = f(P^n,Phi^n)
             CalculateTDGL_RHS(GL_rhs, P_old, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);
 
@@ -291,7 +303,6 @@ void main_main (c_FerroX& rFerroX)
 
                 //update E using PoissonPhi computed with P_new_pre
                 ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);
-
                 // compute f^{n+1,*} = f(P^{n+1,*},Phi^{n+1,*})
                 CalculateTDGL_RHS(GL_rhs_pre, P_new_pre, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi);