diff --git a/CHANGELOG.md b/CHANGELOG.md
index a795c53..4d69d72 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `Pore1D::initialize` and `Pore3D::initialize` now accept initial values for the density profiles as optional arguments. [#24](https://github.com/feos-org/feos-dft/pull/24)
 - Internally restructured the `DFT` structure to avoid redundant data. [#24](https://github.com/feos-org/feos-dft/pull/24)
 - Removed the `m` function in `FluidParameters`, it is instead inferred from `HelmholtzEnergyFunctional` which is now a supertrait of `FluidParameters`. [#24](https://github.com/feos-org/feos-dft/pull/24)
+- Added optional field `cutoff_radius` to `ExternalPotential::FreeEnergyAveraged`. [#25](https://github.com/feos-org/feos-dft/pull/25)
 
 ### Packaging
 - Updated `pyo3` and `numpy` dependencies to 0.16.
diff --git a/src/adsorption/external_potential.rs b/src/adsorption/external_potential.rs
index d23a9ec..47a14d7 100644
--- a/src/adsorption/external_potential.rs
+++ b/src/adsorption/external_potential.rs
@@ -49,6 +49,7 @@ pub enum ExternalPotential<U> {
         pore_center: [f64; 3],
         system_size: [QuantityScalar<U>; 3],
         n_grid: [usize; 2],
+        cutoff_radius: Option<f64>,
     },
 
     /// Custom potential
@@ -168,6 +169,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -185,6 +187,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         n_grid,
                         temperature,
                         Geometry::Cartesian,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),
@@ -314,6 +317,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -331,6 +335,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         n_grid,
                         temperature,
                         Geometry::Cylindrical,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),
@@ -475,6 +480,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -492,6 +498,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         n_grid,
                         temperature,
                         Geometry::Spherical,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),
diff --git a/src/adsorption/fea_potential.rs b/src/adsorption/fea_potential.rs
index da81444..e1ef4a0 100644
--- a/src/adsorption/fea_potential.rs
+++ b/src/adsorption/fea_potential.rs
@@ -20,12 +20,13 @@ pub fn calculate_fea_potential<U: EosUnit>(
     n_grid: &[usize; 2],
     temperature: f64,
     geometry: Geometry,
+    cutoff_radius: Option<f64>,
 ) -> Array1<f64> {
     // allocate external potential
     let mut potential: Array1<f64> = Array1::zeros(grid.len());
 
     // calculate squared cutoff radius
-    let cutoff_radius2 = CUTOFF_RADIUS.powi(2);
+    let cutoff_radius2 = cutoff_radius.unwrap_or(CUTOFF_RADIUS).powi(2);
 
     // dimensionless solid coordinates
     let coordinates = Array2::from_shape_fn(coordinates.raw_dim(), |(i, j)| {
diff --git a/src/python/adsorption/external_potential.rs b/src/python/adsorption/external_potential.rs
index 6d595b7..9b23e91 100644
--- a/src/python/adsorption/external_potential.rs
+++ b/src/python/adsorption/external_potential.rs
@@ -185,13 +185,15 @@ impl PyExternalPotential {
     ///     The size of the unit cell.
     /// n_grid : [int; 2]
     ///     The number of grid points in each direction.
+    /// cutoff_radius : float, optional
+    ///     The cutoff used in the calculation of fluid/wall interactions.
     /// Returns
     /// -------
     /// ExternalPotential
     ///
     #[staticmethod]
     #[pyo3(
-        text_signature = "(coordinates, sigma_ss, epsilon_k_ss, pore_center, system_size, n_grid)"
+        text_signature = "(coordinates, sigma_ss, epsilon_k_ss, pore_center, system_size, n_grid, cutoff_radius=None)"
     )]
     pub fn FreeEnergyAveraged(
         coordinates: &PySIArray2,
@@ -200,6 +202,7 @@ impl PyExternalPotential {
         pore_center: [f64; 3],
         system_size: [PySINumber; 3],
         n_grid: [usize; 2],
+        cutoff_radius: Option<f64>,
     ) -> Self {
         Self(ExternalPotential::FreeEnergyAveraged {
             coordinates: coordinates.clone().into(),
@@ -212,6 +215,7 @@ impl PyExternalPotential {
                 system_size[2].into(),
             ],
             n_grid,
+            cutoff_radius,
         })
     }
 }