diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs
index 4b1a2b7..b03f49a 100644
--- a/src/adsorption/mod.rs
+++ b/src/adsorption/mod.rs
@@ -4,7 +4,7 @@ use super::solver::DFTSolver;
 use feos_core::{
     Contributions, EosError, EosResult, EosUnit, EquationOfState, SolverOptions, StateBuilder,
 };
-use ndarray::{arr1, Array1, Dimension, Ix1, Ix3};
+use ndarray::{arr1, Array1, Dimension, Ix1, Ix3, RemoveAxis};
 use quantity::{QuantityArray1, QuantityArray2, QuantityScalar};
 use std::rc::Rc;
 
@@ -55,12 +55,17 @@ where
         })
     }
 
-    fn equilibrium<D: Dimension, F: HelmholtzEnergyFunctional + FluidParameters>(
+    fn equilibrium<
+        D: Dimension + RemoveAxis + 'static,
+        F: HelmholtzEnergyFunctional + FluidParameters,
+    >(
         &self,
         equilibrium: &Adsorption<U, D, F>,
     ) -> EosResult<(QuantityArray1<U>, QuantityArray1<U>)>
     where
         D::Larger: Dimension<Smaller = D>,
+        D::Smaller: Dimension<Larger = D>,
+        <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
     {
         let p_eq = equilibrium.pressure().get(0);
         match self {
@@ -111,10 +116,16 @@ pub type Adsorption1D<U, F> = Adsorption<U, Ix1, F>;
 /// Container structure for adsorption isotherms in 3D pores.
 pub type Adsorption3D<U, F> = Adsorption<U, Ix3, F>;
 
-impl<U: EosUnit, D: Dimension, F: HelmholtzEnergyFunctional + FluidParameters> Adsorption<U, D, F>
+impl<
+        U: EosUnit,
+        D: Dimension + RemoveAxis + 'static,
+        F: HelmholtzEnergyFunctional + FluidParameters,
+    > Adsorption<U, D, F>
 where
     QuantityScalar<U>: std::fmt::Display,
     D::Larger: Dimension<Smaller = D>,
+    D::Smaller: Dimension<Larger = D>,
+    <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
     fn new<S: PoreSpecification<U, D>>(
         functional: &Rc<DFT<F>>,
diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs
index cd1a849..e4b441c 100644
--- a/src/adsorption/pore.rs
+++ b/src/adsorption/pore.rs
@@ -8,7 +8,7 @@ use crate::solver::DFTSolver;
 use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder};
 use ndarray::prelude::*;
 use ndarray::Axis as Axis_nd;
-use ndarray::Zip;
+use ndarray::{RemoveAxis, Zip};
 use ndarray_stats::QuantileExt;
 use quantity::{QuantityArray, QuantityArray2, QuantityArray4, QuantityScalar};
 use std::rc::Rc;
@@ -133,22 +133,19 @@ impl<U: Copy, D: Dimension, F> Clone for PoreProfile<U, D, F> {
     }
 }
 
-impl<U: EosUnit, D: Dimension, F: HelmholtzEnergyFunctional> PoreProfile<U, D, F>
+impl<U: EosUnit, D: Dimension + RemoveAxis + 'static, F: HelmholtzEnergyFunctional>
+    PoreProfile<U, D, F>
 where
     D::Larger: Dimension<Smaller = D>,
+    D::Smaller: Dimension<Larger = D>,
+    <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
     pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> EosResult<()> {
         // Solve the profile
         self.profile.solve(solver, debug)?;
 
         // calculate grand potential density
-        let omega = self
-            .profile
-            .integrate(&self.profile.dft.grand_potential_density(
-                self.profile.temperature,
-                &self.profile.density,
-                &self.profile.convolver,
-            )?);
+        let omega = self.profile.grand_potential()?;
         self.grand_potential = Some(omega);
 
         // calculate interfacial tension
diff --git a/src/functional.rs b/src/functional.rs
index 47114cb..8fe1818 100644
--- a/src/functional.rs
+++ b/src/functional.rs
@@ -213,7 +213,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         let rho = density.to_reduced(U::reference_density())?;
         let (mut f, dfdrho) = self.functional_derivative(t, &rho, convolver)?;
 
-        // calculate the grand potential density
+        // Calculate the grand potential density
         for ((rho, dfdrho), &m) in rho
             .outer_iter()
             .zip(dfdrho.outer_iter())
diff --git a/src/interface/mod.rs b/src/interface/mod.rs
index b535433..0c98417 100644
--- a/src/interface/mod.rs
+++ b/src/interface/mod.rs
@@ -40,11 +40,8 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
 
         // postprocess
         self.surface_tension = Some(self.profile.integrate(
-            &(self.profile.dft.grand_potential_density(
-                self.profile.temperature,
-                &self.profile.density,
-                &self.profile.convolver,
-            )? + self.vle.vapor().pressure(Contributions::Total)),
+            &(self.profile.grand_potential_density()?
+                + self.vle.vapor().pressure(Contributions::Total)),
         ));
         let delta_rho = self.vle.liquid().density - self.vle.vapor().density;
         self.equimolar_radius = Some(
diff --git a/src/profile.rs b/src/profile.rs
index e83ef86..1d0b625 100644
--- a/src/profile.rs
+++ b/src/profile.rs
@@ -564,6 +564,15 @@ where
         Ok(self.integrate(&self.entropy_density(contributions)?))
     }
 
+    pub fn grand_potential_density(&self) -> EosResult<QuantityArray<U, D>> {
+        self.dft
+            .grand_potential_density(self.temperature, &self.density, &self.convolver)
+    }
+
+    pub fn grand_potential(&self) -> EosResult<QuantityScalar<U>> {
+        Ok(self.integrate(&self.grand_potential_density()?))
+    }
+
     pub fn internal_energy(&self, contributions: Contributions) -> EosResult<QuantityScalar<U>> {
         // initialize convolver
         let t = self.temperature.to_reduced(U::reference_temperature())?;
diff --git a/src/python/profile.rs b/src/python/profile.rs
index a73049b..5829ba1 100644
--- a/src/python/profile.rs
+++ b/src/python/profile.rs
@@ -171,6 +171,13 @@ macro_rules! impl_profile {
                     self.0.profile.internal_energy(contributions)?,
                 ))
             }
+
+            #[getter]
+            fn get_grand_potential_density(&self) -> PyResult<$si_arr> {
+                Ok($si_arr::from(
+                    self.0.profile.grand_potential_density()?,
+                ))
+            }
         }
     };
 }
diff --git a/src/solvation/mod.rs b/src/solvation/mod.rs
index e739c8c..bb23fff 100644
--- a/src/solvation/mod.rs
+++ b/src/solvation/mod.rs
@@ -36,13 +36,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> SolvationProfile<U, F> {
         self.profile.solve(solver, debug)?;
 
         // calculate grand potential density
-        let omega = self
-            .profile
-            .integrate(&self.profile.dft.grand_potential_density(
-                self.profile.temperature,
-                &self.profile.density,
-                &self.profile.convolver,
-            )?);
+        let omega = self.profile.grand_potential()?;
         self.grand_potential = Some(omega);
 
         // calculate solvation free energy
diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs
index 2a895c7..a9b5794 100644
--- a/src/solvation/pair_correlation.rs
+++ b/src/solvation/pair_correlation.rs
@@ -76,11 +76,8 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<U
 
         // calculate self solvation free energy
         self.self_solvation_free_energy = Some(self.profile.integrate(
-            &(self.profile.dft.grand_potential_density(
-                self.profile.temperature,
-                &self.profile.density,
-                &self.profile.convolver,
-            )? + self.profile.bulk.pressure(Contributions::Total)),
+            &(self.profile.grand_potential_density()?
+                + self.profile.bulk.pressure(Contributions::Total)),
         ));
 
         // calculate structure factor