diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/Deformation/MaterialForces.h
index 3852be91f4626ba4065a579ea2fd234c94261ed4..89b5c89bbdb051ee7e4942885fe93cd5f9aa7738 100644
--- a/ProcessLib/Deformation/MaterialForces.h
+++ b/ProcessLib/Deformation/MaterialForces.h
@@ -37,12 +37,12 @@ template <int DisplacementDim, typename ShapeFunction,
           typename ShapeMatricesType, typename NodalForceVectorType,
           typename NodalDisplacementVectorType, typename GradientVectorType,
           typename GradientMatrixType, typename IPData, typename StressData,
-          typename IntegrationMethod>
+          typename OutputData, typename IntegrationMethod>
 std::vector<double> const& getMaterialForces(
     std::vector<double> const& local_x, std::vector<double>& nodal_values,
     IntegrationMethod const& _integration_method, IPData const& _ip_data,
-    StressData const& stress_data, MeshLib::Element const& element,
-    bool const is_axially_symmetric)
+    StressData const& stress_data, OutputData const& output_data,
+    MeshLib::Element const& element, bool const is_axially_symmetric)
 {
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
@@ -57,7 +57,8 @@ std::vector<double> const& getMaterialForces(
         auto const& N = _ip_data[ip].N;
         auto const& dNdx = _ip_data[ip].dNdx;
 
-        auto const& psi = _ip_data[ip].free_energy_density;
+        auto const& psi =
+            output_data[ip].free_energy_density_data.free_energy_density;
 
         auto const x_coord =
             NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h
index b5e3e124151c883a94234994d52d2badfe317094..f6525c7d18f5b7f8aa09425803a3f8f28f1a9841 100644
--- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveData.h
@@ -9,6 +9,7 @@
 
 #pragma once
 
+#include "FreeEnergyDensity.h"
 #include "Gravity.h"
 #include "ProcessLib/Reflection/ReflectionData.h"
 #include "SolidDensity.h"
@@ -52,12 +53,14 @@ template <int DisplacementDim>
 struct OutputData
 {
     StrainData<DisplacementDim> eps_data;
+    FreeEnergyDensityData free_energy_density_data;
 
     static auto reflect()
     {
         using Self = OutputData<DisplacementDim>;
 
-        return Reflection::reflectWithoutName(&Self::eps_data);
+        return Reflection::reflectWithoutName(&Self::eps_data,
+                                              &Self::free_energy_density_data);
     }
 };
 
diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
index d0eafec3f1eaad5a39d71a5777d94e4e16effe6c..f5071e1c1a33e5a6a5ab17d2747145045912344e 100644
--- a/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
@@ -40,6 +40,8 @@ void ConstitutiveSetting<DisplacementDim>::eval(
     auto& s_mech_data = cd.s_mech_data;
     auto& volumetric_body_force = cd.volumetric_body_force;
 
+    auto& free_energy_density_data = out.free_energy_density_data;
+
     Temperature const T{T_ref};
     SpaceTimeData const x_t{x_position, t, dt};
     MediaData const media_data{medium};
@@ -47,7 +49,7 @@ void ConstitutiveSetting<DisplacementDim>::eval(
     assertEvalArgsUnique(models.s_mech_model);
     models.s_mech_model.eval(x_t, T, eps_data, eps_data_prev, mat_state,
                              prev_state.stress_data, state.stress_data,
-                             s_mech_data);
+                             s_mech_data, free_energy_density_data);
 
     assertEvalArgsUnique(models.rho_S_model);
     models.rho_S_model.eval(x_t, media_data, T, rho_SR);
diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h
new file mode 100644
index 0000000000000000000000000000000000000000..c52ee5ae847037101fc035bc79aa0d41dfd963f5
--- /dev/null
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/FreeEnergyDensity.h
@@ -0,0 +1,26 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include "Base.h"
+
+namespace ProcessLib::SmallDeformation
+{
+struct FreeEnergyDensityData
+{
+    double free_energy_density;
+
+    static auto reflect()
+    {
+        return ProcessLib::Reflection::reflectWithName(
+            "free_energy_density", &FreeEnergyDensityData::free_energy_density);
+    }
+};
+}  // namespace ProcessLib::SmallDeformation
diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp
index 00eb985ae1f045533edf234917f1ea4ebfbbdf83..283e195768148345aae3b054eeb8f1c4b8979a5d 100644
--- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.cpp
@@ -22,7 +22,8 @@ void SolidMechanicsModel<DisplacementDim>::eval(
     MaterialStateData<DisplacementDim>& mat_state,
     PrevState<StressData<DisplacementDim>> const& stress_data_prev,
     StressData<DisplacementDim>& stress_data,
-    SolidMechanicsDataStateless<DisplacementDim>& current_stateless) const
+    SolidMechanicsDataStateless<DisplacementDim>& current_stateless,
+    FreeEnergyDensityData& free_energy_density_data) const
 {
     namespace MPL = MaterialPropertyLib;
 
@@ -59,6 +60,11 @@ void SolidMechanicsModel<DisplacementDim>::eval(
 
     std::tie(stress_data.sigma, mat_state.material_state_variables,
              current_stateless.stiffness_tensor) = std::move(*solution);
+
+    free_energy_density_data.free_energy_density =
+        solid_material_.computeFreeEnergyDensity(
+            x_t.t, x_t.x, x_t.dt, eps_data.eps, stress_data.sigma,
+            *mat_state.material_state_variables);
 }
 
 template struct SolidMechanicsModel<2>;
diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h
index 572ac2391de9982ad4eb4b4f87e6a648b0cdb7fd..95eafc247ab868845ea12ca1052325ae4eaa6997 100644
--- a/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/SolidMechanics.h
@@ -9,6 +9,7 @@
 
 #pragma once
 
+#include "FreeEnergyDensity.h"
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "MaterialState.h"
 #include "StressData.h"
@@ -36,15 +37,15 @@ struct SolidMechanicsModel
     {
     }
 
-    void eval(
-        SpaceTimeData const& x_t,
-        Temperature const& temperature,
-        StrainData<DisplacementDim> const& eps_data,
-        PrevState<StrainData<DisplacementDim>> const& eps_data_prev,
-        MaterialStateData<DisplacementDim>& mat_state,
-        PrevState<StressData<DisplacementDim>> const& stress_data_prev,
-        StressData<DisplacementDim>& stress_data,
-        SolidMechanicsDataStateless<DisplacementDim>& current_stateless) const;
+    void eval(SpaceTimeData const& x_t,
+              Temperature const& temperature,
+              StrainData<DisplacementDim> const& eps_data,
+              PrevState<StrainData<DisplacementDim>> const& eps_data_prev,
+              MaterialStateData<DisplacementDim>& mat_state,
+              PrevState<StressData<DisplacementDim>> const& stress_data_prev,
+              StressData<DisplacementDim>& stress_data,
+              SolidMechanicsDataStateless<DisplacementDim>& current_stateless,
+              FreeEnergyDensityData& free_energy_density_data) const;
 
     auto getInternalVariables() const
     {
diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index 5ac074a06eefdd74b3dd5ec02e4c20a71c64ef01..50626d70315af314ba61224d402c2d3b0b2d15b3 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
@@ -154,12 +154,6 @@ struct SmallDeformationLocalAssemblerInterface
         }
     }
 
-    virtual std::vector<double> const& getIntPtFreeEnergyDensity(
-        const double t,
-        std::vector<GlobalVector*> const& x,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
-        std::vector<double>& cache) const = 0;
-
     // TODO move to NumLib::ExtrapolatableElement
     unsigned getNumberOfIntegrationPoints() const
     {
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index c9c5899d3d3408372c2df829e78a3224eb062b79..888044b7fa358565432045aec664e61c08b80baa 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -40,8 +40,6 @@ template <typename BMatricesType, typename ShapeMatricesType,
           int DisplacementDim>
 struct IntegrationPointData final
 {
-    double free_energy_density = 0;
-
     double integration_weight;
     typename ShapeMatricesType::NodalRowVectorType N;
     typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
@@ -315,15 +313,6 @@ public:
                 this->prev_states_[ip], this->material_states_[ip],
                 this->output_data_[ip]);
 
-            auto const& eps = this->output_data_[ip].eps_data.eps;
-            auto const& sigma = this->current_states_[ip].stress_data.sigma;
-            auto& state = *this->material_states_[ip].material_state_variables;
-
-            // Update free energy density needed for material forces.
-            _ip_data[ip].free_energy_density =
-                this->solid_material_.computeFreeEnergyDensity(
-                    t, x_position, dt, eps, sigma, state);
-
             this->material_states_[ip].pushBackState();
         }
 
@@ -341,9 +330,10 @@ public:
             DisplacementDim, ShapeFunction, ShapeMatricesType,
             typename BMatricesType::NodalForceVectorType,
             NodalDisplacementVectorType, GradientVectorType,
-            GradientMatrixType>(
-            local_x, nodal_values, this->integration_method_, _ip_data,
-            this->current_states_, this->element_, this->is_axially_symmetric_);
+            GradientMatrixType>(local_x, nodal_values,
+                                this->integration_method_, _ip_data,
+                                this->current_states_, this->output_data_,
+                                this->element_, this->is_axially_symmetric_);
     }
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
@@ -355,22 +345,6 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
-    std::vector<double> const& getIntPtFreeEnergyDensity(
-        const double /*t*/,
-        std::vector<GlobalVector*> const& /*x*/,
-        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
-        std::vector<double>& cache) const override
-    {
-        cache.clear();
-        cache.reserve(_ip_data.size());
-
-        transform(
-            cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
-            [](auto const& ip_data) { return ip_data.free_energy_density; });
-
-        return cache;
-    }
-
 private:
     static constexpr auto localDOF(std::vector<double> const& x)
     {
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index e9ae759b7ff377ca92534944f0486c6de6f2fa7e..bed27c80494cb1229d67facbedb942a48c053e2b 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -104,10 +104,6 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
                              std::move(get_ip_values_function)));
     };
 
-    add_secondary_variable("free_energy_density",
-                           1,
-                           &LocalAssemblerInterface::getIntPtFreeEnergyDensity);
-
     ProcessLib::Reflection::addReflectedSecondaryVariables<DisplacementDim>(
         SmallDeformationLocalAssemblerInterface<
             DisplacementDim>::getReflectionDataForOutput(),