diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md
index 706329fc3..85a6f6da8 100644
--- a/feos-dft/CHANGELOG.md
+++ b/feos-dft/CHANGELOG.md
@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Changed
+- The 3D DFT functionalities (3D pores, solvation free energy, free-energy-averaged potentials) are hidden behind the new `3d_dft` feature. For the Python package, the feature is always enabled. [#51](https://github.com/feos-org/feos/pull/51)
 
 ## [0.3.1] - 2022-08-26
 ### Changed
diff --git a/feos-dft/Cargo.toml b/feos-dft/Cargo.toml
index a806b26ca..0a4a6da8d 100644
--- a/feos-dft/Cargo.toml
+++ b/feos-dft/Cargo.toml
@@ -14,23 +14,24 @@ exclude = ["/.github/*", "*.ipynb"]
 
 [package.metadata.docs.rs]
 rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
+features = [ "3d_dft" ]
 
 [dependencies]
 quantity = { version = "0.5", features = ["linalg"] }
 feos-core = { version = "0.3", path = "../feos-core" }
 num-dual = "0.5"
-ndarray = { version = "0.15", features = ["serde", "rayon"] }
-ndarray-stats = "0.5"
+ndarray = "0.15"
 rustdct = "0.7"
 rustfft = "6.0"
 ang = "0.6"
 num-traits = "0.2"
-libc = "0.2"
-gauss-quad = "0.1"
+libm = "0.2"
+gauss-quad = { version = "0.1", optional = true }
 petgraph = "0.6"
 numpy = { version = "0.16", optional = true }
 pyo3 = { version = "0.16", optional = true }
 
 [features]
 default = []
-python = ["pyo3", "numpy", "feos-core/python"]
+3d_dft = ["gauss-quad", "ndarray/rayon"]
+python = ["pyo3", "numpy", "feos-core/python", "3d_dft"]
diff --git a/feos-dft/src/adsorption/external_potential.rs b/feos-dft/src/adsorption/external_potential.rs
index 9caee2d09..2f01bd973 100644
--- a/feos-dft/src/adsorption/external_potential.rs
+++ b/feos-dft/src/adsorption/external_potential.rs
@@ -1,11 +1,14 @@
+#[cfg(feature = "3d_dft")]
 use crate::adsorption::fea_potential::calculate_fea_potential;
 use crate::functional::HelmholtzEnergyFunctional;
+#[cfg(feature = "3d_dft")]
 use crate::geometry::Geometry;
 use feos_core::EosUnit;
-use libc::c_double;
+use libm::tgamma;
 use ndarray::{Array1, Array2, Axis as Axis_nd};
+#[cfg(feature = "3d_dft")]
 use quantity::{QuantityArray2, QuantityScalar};
-use std::f64::consts::PI;
+use std::{f64::consts::PI, marker::PhantomData};
 
 const DELTA_STEELE: f64 = 3.35;
 
@@ -49,6 +52,7 @@ pub enum ExternalPotential<U> {
         rho_s: f64,
     },
     /// Free-energy averaged potential:
+    #[cfg(feature = "3d_dft")]
     FreeEnergyAveraged {
         coordinates: QuantityArray2<U>,
         sigma_ss: Array1<f64>,
@@ -61,6 +65,10 @@ pub enum ExternalPotential<U> {
 
     /// Custom potential
     Custom(Array2<f64>),
+
+    /// Needed to keep `FreeEnergyAveraged` optional
+    #[doc(hidden)]
+    Phantom(PhantomData<U>),
 }
 
 /// Parameters of the fluid required to evaluate the external potential.
@@ -69,6 +77,7 @@ pub trait FluidParameters: HelmholtzEnergyFunctional {
     fn sigma_ff(&self) -> &Array1<f64>;
 }
 
+#[allow(unused_variables)]
 impl<U: EosUnit> ExternalPotential<U> {
     // Evaluate the external potential in cartesian coordinates for a given grid and fluid parameters.
     pub fn calculate_cartesian_potential<P: FluidParameters>(
@@ -183,6 +192,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                                 - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -211,7 +221,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -348,6 +358,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * sigma_sf[i].powi(3)
                             * *rho_s)
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -376,7 +387,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -537,6 +548,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
                                 - sum_n(4, r_grid, sigma_sf[i], pore_size)))
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -565,7 +577,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -578,8 +590,8 @@ fn phi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
 
     (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(m3n2)) * 4.0 * PI.sqrt() / n2m3
         * sigma_r.powf(n2m3)
-        * gamma(n as f64 - 0.5)
-        / gamma(n as f64)
+        * tgamma(n as f64 - 0.5)
+        / tgamma(n as f64)
         * taylor_2f1_phi(r_r, n)
 }
 
@@ -587,8 +599,8 @@ fn psi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
     (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(2.0 - 2.0 * n as f64))
         * 4.0
         * PI.sqrt()
-        * gamma(n as f64 - 0.5)
-        / gamma(n as f64)
+        * tgamma(n as f64 - 0.5)
+        / tgamma(n as f64)
         * sigma_r.powf(2.0 * n as f64 - 2.0)
         * taylor_2f1_psi(r_r, n)
 }
@@ -637,11 +649,3 @@ fn taylor_2f1_psi(x: &Array1<f64>, n: i32) -> Array1<f64> {
         _ => unreachable!(),
     }
 }
-
-extern "C" {
-    fn tgamma(x: c_double) -> c_double;
-}
-
-fn gamma(x: f64) -> f64 {
-    unsafe { tgamma(x) }
-}
diff --git a/feos-dft/src/adsorption/fea_potential.rs b/feos-dft/src/adsorption/fea_potential.rs
index e1ef4a0c9..39a026446 100644
--- a/feos-dft/src/adsorption/fea_potential.rs
+++ b/feos-dft/src/adsorption/fea_potential.rs
@@ -1,9 +1,9 @@
+use super::pore3d::{calculate_distance2, evaluate_lj_potential};
 use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL};
 use crate::Geometry;
 use feos_core::EosUnit;
 use gauss_quad::GaussLegendre;
 use ndarray::{Array1, Array2, Zip};
-use ndarray_stats::SummaryStatisticsExt;
 use quantity::{QuantityArray2, QuantityScalar};
 use std::f64::consts::PI;
 use std::usize;
@@ -114,7 +114,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
                 let distance2 = calculate_distance2(point, &coordinates, system_size);
                 let potential_sum: f64 = (0..sigma_sf.len())
                     .map(|alpha| {
-                        mi * evaluate(
+                        mi * evaluate_lj_potential(
                             distance2[alpha],
                             sigma_sf[alpha],
                             epsilon_k_sf[alpha],
@@ -125,44 +125,8 @@ pub fn calculate_fea_potential<U: EosUnit>(
                 potential_2d[[i1, i2]] = (-potential_sum.min(MAX_POTENTIAL)).exp();
             }
         }
-        *f = potential_2d.weighted_sum(&weights).unwrap();
+        *f = (potential_2d * &weights).sum();
     });
 
     -temperature * potential.map(|p| (p / weights_sum).ln())
 }
-
-// -temperature * potential.map(|p| (p / weights_sum).ln())
-
-/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
-fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
-    let sigma_r = sigma.powi(2) / distance2;
-
-    let potential: f64 = if distance2 > cutoff_radius2 {
-        0.0
-    } else if distance2 == 0.0 {
-        f64::INFINITY
-    } else {
-        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
-    };
-
-    potential
-}
-
-/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
-fn calculate_distance2(
-    point: [f64; 3],
-    coordinates: &Array2<f64>,
-    system_size: [f64; 3],
-) -> Array1<f64> {
-    Array1::from_shape_fn(coordinates.ncols(), |i| {
-        let mut rx = coordinates[[0, i]] - point[0];
-        let mut ry = coordinates[[1, i]] - point[1];
-        let mut rz = coordinates[[2, i]] - point[2];
-
-        rx = rx - system_size[0] * (rx / system_size[0]).round();
-        ry = ry - system_size[1] * (ry / system_size[1]).round();
-        rz = rz - system_size[2] * (rz / system_size[2]).round();
-
-        rx.powi(2) + ry.powi(2) + rz.powi(2)
-    })
-}
diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs
index 9f8bc953c..e03d005bd 100644
--- a/feos-dft/src/adsorption/mod.rs
+++ b/feos-dft/src/adsorption/mod.rs
@@ -11,10 +11,16 @@ use std::iter;
 use std::rc::Rc;
 
 mod external_potential;
+#[cfg(feature = "3d_dft")]
 mod fea_potential;
 mod pore;
 pub use external_potential::{ExternalPotential, FluidParameters};
-pub use pore::{Pore1D, Pore3D, PoreProfile, PoreProfile1D, PoreProfile3D, PoreSpecification};
+pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
+
+#[cfg(feature = "3d_dft")]
+mod pore3d;
+#[cfg(feature = "3d_dft")]
+pub use pore3d::{Pore3D, PoreProfile3D};
 
 const MAX_ITER_ADSORPTION_EQUILIBRIUM: usize = 50;
 const TOL_ADSORPTION_EQUILIBRIUM: f64 = 1e-8;
diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs
index 8b131b62c..4bbc20a48 100644
--- a/feos-dft/src/adsorption/pore.rs
+++ b/feos-dft/src/adsorption/pore.rs
@@ -3,14 +3,13 @@ use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT};
 use crate::functional_contribution::FunctionalContribution;
 use crate::geometry::{Axis, Geometry, Grid};
-use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
+use crate::profile::{DFTProfile, MAX_POTENTIAL};
 use crate::solver::DFTSolver;
 use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder};
 use ndarray::prelude::*;
 use ndarray::Axis as Axis_nd;
-use ndarray::{RemoveAxis, Zip};
-use ndarray_stats::QuantileExt;
-use quantity::{QuantityArray, QuantityArray2, QuantityArray4, QuantityScalar};
+use ndarray::RemoveAxis;
+use quantity::{QuantityArray, QuantityArray2, QuantityScalar};
 use std::rc::Rc;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
@@ -43,39 +42,6 @@ impl<U: EosUnit> Pore1D<U> {
     }
 }
 
-/// Parameters required to specify a 3D pore.
-pub struct Pore3D<U> {
-    system_size: [QuantityScalar<U>; 3],
-    n_grid: [usize; 3],
-    coordinates: QuantityArray2<U>,
-    sigma_ss: Array1<f64>,
-    epsilon_k_ss: Array1<f64>,
-    potential_cutoff: Option<f64>,
-    cutoff_radius: Option<QuantityScalar<U>>,
-}
-
-impl<U> Pore3D<U> {
-    pub fn new(
-        system_size: [QuantityScalar<U>; 3],
-        n_grid: [usize; 3],
-        coordinates: QuantityArray2<U>,
-        sigma_ss: Array1<f64>,
-        epsilon_k_ss: Array1<f64>,
-        potential_cutoff: Option<f64>,
-        cutoff_radius: Option<QuantityScalar<U>>,
-    ) -> Self {
-        Self {
-            system_size,
-            n_grid,
-            coordinates,
-            sigma_ss,
-            epsilon_k_ss,
-            potential_cutoff,
-            cutoff_radius,
-        }
-    }
-}
-
 /// Trait for the generic implementation of adsorption applications.
 pub trait PoreSpecification<U: EosUnit, D: Dimension> {
     /// Initialize a new single pore.
@@ -119,8 +85,6 @@ pub struct PoreProfile<U, D: Dimension, F> {
 
 /// Density profile and properties of a 1D confined system.
 pub type PoreProfile1D<U, F> = PoreProfile<U, Ix1, F>;
-/// Density profile and properties of a 3D confined system.
-pub type PoreProfile3D<U, F> = PoreProfile<U, Ix3, F>;
 
 impl<U: Copy, D: Dimension, F> Clone for PoreProfile<U, D, F> {
     fn clone(&self) -> Self {
@@ -179,7 +143,13 @@ impl<U: EosUnit> PoreSpecification<U, Ix1> for Pore1D<U> {
 
         let axis = match self.geometry {
             Geometry::Cartesian => {
-                let potential_offset = POTENTIAL_OFFSET * bulk.eos.sigma_ff().max().unwrap();
+                let potential_offset = POTENTIAL_OFFSET
+                    * bulk
+                        .eos
+                        .sigma_ff()
+                        .iter()
+                        .max_by(|a, b| a.total_cmp(b))
+                        .unwrap();
                 Axis::new_cartesian(n_grid, 0.5 * self.pore_size, Some(potential_offset))?
             }
             Geometry::Cylindrical => Axis::new_polar(n_grid, self.pore_size)?,
@@ -219,65 +189,6 @@ impl<U: EosUnit> PoreSpecification<U, Ix1> for Pore1D<U> {
     }
 }
 
-impl<U: EosUnit> PoreSpecification<U, Ix3> for Pore3D<U> {
-    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
-        &self,
-        bulk: &State<U, DFT<F>>,
-        density: Option<&QuantityArray4<U>>,
-        external_potential: Option<&Array4<f64>>,
-    ) -> EosResult<PoreProfile3D<U, F>> {
-        let dft: &F = &bulk.eos;
-
-        // generate grid
-        let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None)?;
-        let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None)?;
-        let z = Axis::new_cartesian(self.n_grid[2], self.system_size[2], None)?;
-
-        // move center of geometry of solute to box center
-        let coordinates = Array2::from_shape_fn(self.coordinates.raw_dim(), |(i, j)| {
-            (self.coordinates.get((i, j)))
-                .to_reduced(U::reference_length())
-                .unwrap()
-        });
-
-        // temperature
-        let t = bulk.temperature.to_reduced(U::reference_temperature())?;
-
-        // calculate external potential
-        let external_potential = external_potential.map_or_else(
-            || {
-                external_potential_3d(
-                    dft,
-                    [&x, &y, &z],
-                    self.system_size,
-                    coordinates,
-                    &self.sigma_ss,
-                    &self.epsilon_k_ss,
-                    self.cutoff_radius,
-                    self.potential_cutoff,
-                    t,
-                )
-            },
-            |e| Ok(e.clone()),
-        )?;
-
-        // initialize convolver
-        let grid = Grid::Periodical3(x, y, z);
-        let weight_functions = dft.weight_functions(t);
-        let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
-
-        Ok(PoreProfile {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density)?,
-            grand_potential: None,
-            interfacial_tension: None,
-        })
-    }
-
-    fn dimension(&self) -> i32 {
-        3
-    }
-}
-
 fn external_potential_1d<U: EosUnit, P: FluidParameters>(
     pore_width: QuantityScalar<U>,
     temperature: QuantityScalar<U>,
@@ -334,108 +245,6 @@ fn external_potential_1d<U: EosUnit, P: FluidParameters>(
     Ok(external_potential)
 }
 
-pub fn external_potential_3d<U: EosUnit, F: FluidParameters>(
-    functional: &F,
-    axis: [&Axis; 3],
-    system_size: [QuantityScalar<U>; 3],
-    coordinates: Array2<f64>,
-    sigma_ss: &Array1<f64>,
-    epsilon_ss: &Array1<f64>,
-    cutoff_radius: Option<QuantityScalar<U>>,
-    potential_cutoff: Option<f64>,
-    reduced_temperature: f64,
-) -> EosResult<Array4<f64>> {
-    // allocate external potential
-    let m = functional.m();
-    let mut external_potential = Array4::zeros((
-        m.len(),
-        axis[0].grid.len(),
-        axis[1].grid.len(),
-        axis[2].grid.len(),
-    ));
-
-    let system_size = [
-        system_size[0].to_reduced(U::reference_length())?,
-        system_size[1].to_reduced(U::reference_length())?,
-        system_size[2].to_reduced(U::reference_length())?,
-    ];
-
-    let cutoff_radius = cutoff_radius
-        .unwrap_or(CUTOFF_RADIUS * U::reference_length())
-        .to_reduced(U::reference_length())?;
-
-    // square cut-off radius
-    let cutoff_radius2 = cutoff_radius.powi(2);
-
-    // calculate external potential
-    let sigma_ff = functional.sigma_ff();
-    let epsilon_k_ff = functional.epsilon_k_ff();
-
-    Zip::indexed(&mut external_potential).par_for_each(|(i, ix, iy, iz), u| {
-        let distance2 = calculate_distance2(
-            [&axis[0].grid[ix], &axis[1].grid[iy], &axis[2].grid[iz]],
-            &coordinates,
-            system_size,
-        );
-        let sigma_sf = sigma_ss.mapv(|s| (s + sigma_ff[i]) / 2.0);
-        let epsilon_sf = epsilon_ss.mapv(|e| (e * epsilon_k_ff[i]).sqrt());
-        *u = (0..sigma_ss.len())
-            .map(|alpha| {
-                m[i] * evaluate(
-                    distance2[alpha],
-                    sigma_sf[alpha],
-                    epsilon_sf[alpha],
-                    cutoff_radius2,
-                )
-            })
-            .sum::<f64>()
-            / reduced_temperature
-    });
-
-    let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
-    external_potential.map_inplace(|x| {
-        if *x > potential_cutoff {
-            *x = potential_cutoff
-        }
-    });
-
-    Ok(external_potential)
-}
-
-/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
-fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
-    let sigma_r = sigma.powi(2) / distance2;
-
-    let potential: f64 = if distance2 > cutoff_radius2 {
-        0.0
-    } else if distance2 == 0.0 {
-        f64::INFINITY
-    } else {
-        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
-    };
-
-    potential
-}
-
-/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
-fn calculate_distance2(
-    point: [&f64; 3],
-    coordinates: &Array2<f64>,
-    system_size: [f64; 3],
-) -> Array1<f64> {
-    Array1::from_shape_fn(coordinates.ncols(), |i| {
-        let mut rx = coordinates[[0, i]] - point[0];
-        let mut ry = coordinates[[1, i]] - point[1];
-        let mut rz = coordinates[[2, i]] - point[2];
-
-        rx -= system_size[0] * (rx / system_size[0]).round();
-        ry -= system_size[1] * (ry / system_size[1]).round();
-        rz -= system_size[2] * (rz / system_size[2]).round();
-
-        rx.powi(2) + ry.powi(2) + rz.powi(2)
-    })
-}
-
 const EPSILON_HE: f64 = 10.9;
 const SIGMA_HE: f64 = 2.64;
 
diff --git a/feos-dft/src/adsorption/pore3d.rs b/feos-dft/src/adsorption/pore3d.rs
new file mode 100644
index 000000000..47d823bb7
--- /dev/null
+++ b/feos-dft/src/adsorption/pore3d.rs
@@ -0,0 +1,212 @@
+use super::pore::{PoreProfile, PoreSpecification};
+use crate::adsorption::FluidParameters;
+use crate::convolver::ConvolverFFT;
+use crate::functional::{HelmholtzEnergyFunctional, DFT};
+use crate::geometry::{Axis, Grid};
+use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
+use feos_core::{EosResult, EosUnit, State};
+use ndarray::prelude::*;
+use ndarray::Zip;
+use quantity::{QuantityArray2, QuantityArray4, QuantityScalar};
+
+/// Parameters required to specify a 3D pore.
+pub struct Pore3D<U> {
+    system_size: [QuantityScalar<U>; 3],
+    n_grid: [usize; 3],
+    coordinates: QuantityArray2<U>,
+    sigma_ss: Array1<f64>,
+    epsilon_k_ss: Array1<f64>,
+    potential_cutoff: Option<f64>,
+    cutoff_radius: Option<QuantityScalar<U>>,
+}
+
+impl<U> Pore3D<U> {
+    pub fn new(
+        system_size: [QuantityScalar<U>; 3],
+        n_grid: [usize; 3],
+        coordinates: QuantityArray2<U>,
+        sigma_ss: Array1<f64>,
+        epsilon_k_ss: Array1<f64>,
+        potential_cutoff: Option<f64>,
+        cutoff_radius: Option<QuantityScalar<U>>,
+    ) -> Self {
+        Self {
+            system_size,
+            n_grid,
+            coordinates,
+            sigma_ss,
+            epsilon_k_ss,
+            potential_cutoff,
+            cutoff_radius,
+        }
+    }
+}
+
+/// Density profile and properties of a 3D confined system.
+pub type PoreProfile3D<U, F> = PoreProfile<U, Ix3, F>;
+
+impl<U: EosUnit> PoreSpecification<U, Ix3> for Pore3D<U> {
+    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
+        &self,
+        bulk: &State<U, DFT<F>>,
+        density: Option<&QuantityArray4<U>>,
+        external_potential: Option<&Array4<f64>>,
+    ) -> EosResult<PoreProfile3D<U, F>> {
+        let dft: &F = &bulk.eos;
+
+        // generate grid
+        let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None)?;
+        let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None)?;
+        let z = Axis::new_cartesian(self.n_grid[2], self.system_size[2], None)?;
+
+        // move center of geometry of solute to box center
+        let coordinates = Array2::from_shape_fn(self.coordinates.raw_dim(), |(i, j)| {
+            (self.coordinates.get((i, j)))
+                .to_reduced(U::reference_length())
+                .unwrap()
+        });
+
+        // temperature
+        let t = bulk.temperature.to_reduced(U::reference_temperature())?;
+
+        // calculate external potential
+        let external_potential = external_potential.map_or_else(
+            || {
+                external_potential_3d(
+                    dft,
+                    [&x, &y, &z],
+                    self.system_size,
+                    coordinates,
+                    &self.sigma_ss,
+                    &self.epsilon_k_ss,
+                    self.cutoff_radius,
+                    self.potential_cutoff,
+                    t,
+                )
+            },
+            |e| Ok(e.clone()),
+        )?;
+
+        // initialize convolver
+        let grid = Grid::Periodical3(x, y, z);
+        let weight_functions = dft.weight_functions(t);
+        let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
+
+        Ok(PoreProfile {
+            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density)?,
+            grand_potential: None,
+            interfacial_tension: None,
+        })
+    }
+
+    fn dimension(&self) -> i32 {
+        3
+    }
+}
+
+pub fn external_potential_3d<U: EosUnit, F: FluidParameters>(
+    functional: &F,
+    axis: [&Axis; 3],
+    system_size: [QuantityScalar<U>; 3],
+    coordinates: Array2<f64>,
+    sigma_ss: &Array1<f64>,
+    epsilon_ss: &Array1<f64>,
+    cutoff_radius: Option<QuantityScalar<U>>,
+    potential_cutoff: Option<f64>,
+    reduced_temperature: f64,
+) -> EosResult<Array4<f64>> {
+    // allocate external potential
+    let m = functional.m();
+    let mut external_potential = Array4::zeros((
+        m.len(),
+        axis[0].grid.len(),
+        axis[1].grid.len(),
+        axis[2].grid.len(),
+    ));
+
+    let system_size = [
+        system_size[0].to_reduced(U::reference_length())?,
+        system_size[1].to_reduced(U::reference_length())?,
+        system_size[2].to_reduced(U::reference_length())?,
+    ];
+
+    let cutoff_radius = cutoff_radius
+        .unwrap_or(CUTOFF_RADIUS * U::reference_length())
+        .to_reduced(U::reference_length())?;
+
+    // square cut-off radius
+    let cutoff_radius2 = cutoff_radius.powi(2);
+
+    // calculate external potential
+    let sigma_ff = functional.sigma_ff();
+    let epsilon_k_ff = functional.epsilon_k_ff();
+
+    Zip::indexed(&mut external_potential).par_for_each(|(i, ix, iy, iz), u| {
+        let distance2 = calculate_distance2(
+            [axis[0].grid[ix], axis[1].grid[iy], axis[2].grid[iz]],
+            &coordinates,
+            system_size,
+        );
+        let sigma_sf = sigma_ss.mapv(|s| (s + sigma_ff[i]) / 2.0);
+        let epsilon_sf = epsilon_ss.mapv(|e| (e * epsilon_k_ff[i]).sqrt());
+        *u = (0..sigma_ss.len())
+            .map(|alpha| {
+                m[i] * evaluate_lj_potential(
+                    distance2[alpha],
+                    sigma_sf[alpha],
+                    epsilon_sf[alpha],
+                    cutoff_radius2,
+                )
+            })
+            .sum::<f64>()
+            / reduced_temperature
+    });
+
+    let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
+    external_potential.map_inplace(|x| {
+        if *x > potential_cutoff {
+            *x = potential_cutoff
+        }
+    });
+
+    Ok(external_potential)
+}
+
+/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
+pub(super) fn evaluate_lj_potential(
+    distance2: f64,
+    sigma: f64,
+    epsilon: f64,
+    cutoff_radius2: f64,
+) -> f64 {
+    let sigma_r = sigma.powi(2) / distance2;
+
+    let potential: f64 = if distance2 > cutoff_radius2 {
+        0.0
+    } else if distance2 == 0.0 {
+        f64::INFINITY
+    } else {
+        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
+    };
+
+    potential
+}
+
+/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
+pub(super) fn calculate_distance2(
+    point: [f64; 3],
+    coordinates: &Array2<f64>,
+    system_size: [f64; 3],
+) -> Array1<f64> {
+    Array1::from_shape_fn(coordinates.ncols(), |i| {
+        let mut rx = coordinates[[0, i]] - point[0];
+        let mut ry = coordinates[[1, i]] - point[1];
+        let mut rz = coordinates[[2, i]] - point[2];
+
+        rx -= system_size[0] * (rx / system_size[0]).round();
+        ry -= system_size[1] * (ry / system_size[1]).round();
+        rz -= system_size[2] * (rz / system_size[2]).round();
+
+        rx.powi(2) + ry.powi(2) + rz.powi(2)
+    })
+}
diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs
index 1d0b62542..870b2b142 100644
--- a/feos-dft/src/profile.rs
+++ b/feos-dft/src/profile.rs
@@ -16,6 +16,7 @@ use std::ops::MulAssign;
 use std::rc::Rc;
 
 pub(crate) const MAX_POTENTIAL: f64 = 50.0;
+#[cfg(feature = "3d_dft")]
 pub(crate) const CUTOFF_RADIUS: f64 = 14.0;
 
 /// General specifications for the chemical potential in a DFT calculation.
diff --git a/feos-dft/src/solvation/mod.rs b/feos-dft/src/solvation/mod.rs
index bb23fff09..482d097f8 100644
--- a/feos-dft/src/solvation/mod.rs
+++ b/feos-dft/src/solvation/mod.rs
@@ -1,208 +1,8 @@
 //! Solvation free energies and pair correlaion functions.
-use crate::adsorption::FluidParameters;
-use crate::convolver::ConvolverFFT;
-use crate::functional::{HelmholtzEnergyFunctional, DFT};
-use crate::geometry::{Axis, Grid};
-use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
-use crate::solver::DFTSolver;
-use feos_core::{Contributions, EosResult, EosUnit, State};
-use ndarray::prelude::*;
-use ndarray::Zip;
-use quantity::{QuantityArray2, QuantityScalar};
-
 mod pair_correlation;
 pub use pair_correlation::{PairCorrelation, PairPotential};
 
-/// Density profile and properties of a solute in a inhomogeneous bulk fluid.
-pub struct SolvationProfile<U: EosUnit, F: HelmholtzEnergyFunctional> {
-    pub profile: DFTProfile<U, Ix3, F>,
-    pub grand_potential: Option<QuantityScalar<U>>,
-    pub solvation_free_energy: Option<QuantityScalar<U>>,
-}
-
-impl<U: EosUnit, F: HelmholtzEnergyFunctional> Clone for SolvationProfile<U, F> {
-    fn clone(&self) -> Self {
-        Self {
-            profile: self.profile.clone(),
-            grand_potential: self.grand_potential,
-            solvation_free_energy: self.solvation_free_energy,
-        }
-    }
-}
-
-impl<U: EosUnit, F: HelmholtzEnergyFunctional> SolvationProfile<U, F> {
-    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.grand_potential()?;
-        self.grand_potential = Some(omega);
-
-        // calculate solvation free energy
-        self.solvation_free_energy = Some(
-            (omega + self.profile.bulk.pressure(Contributions::Total) * self.profile.volume())
-                / U::reference_moles(),
-        );
-
-        Ok(())
-    }
-
-    pub fn solve(mut self, solver: Option<&DFTSolver>) -> EosResult<Self> {
-        self.solve_inplace(solver, false)?;
-        Ok(self)
-    }
-}
-
-impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> SolvationProfile<U, F> {
-    pub fn new(
-        bulk: &State<U, DFT<F>>,
-        n_grid: [usize; 3],
-        coordinates: QuantityArray2<U>,
-        sigma_ss: Array1<f64>,
-        epsilon_ss: Array1<f64>,
-        system_size: Option<[QuantityScalar<U>; 3]>,
-        cutoff_radius: Option<QuantityScalar<U>>,
-        potential_cutoff: Option<f64>,
-    ) -> EosResult<Self> {
-        let dft: &F = &bulk.eos;
-
-        let system_size = system_size.unwrap_or([40.0 * U::reference_length(); 3]);
-
-        // generate grid
-        let x = Axis::new_cartesian(n_grid[0], system_size[0], None)?;
-        let y = Axis::new_cartesian(n_grid[1], system_size[1], None)?;
-        let z = Axis::new_cartesian(n_grid[2], system_size[2], None)?;
-
-        // move center of geometry of solute to box center
-        let mut coordinates = Array2::from_shape_fn(coordinates.raw_dim(), |(i, j)| {
-            (coordinates.get((i, j)))
-                .to_reduced(U::reference_length())
-                .unwrap()
-        });
-
-        let center = [
-            system_size[0].to_reduced(U::reference_length())? / 2.0,
-            system_size[1].to_reduced(U::reference_length())? / 2.0,
-            system_size[2].to_reduced(U::reference_length())? / 2.0,
-        ];
-
-        let shift: Array2<f64> = Array2::from_shape_fn((3, 1), |(i, _)| {
-            center[i] - coordinates.row(i).sum() / coordinates.ncols() as f64
-        });
-
-        coordinates = coordinates + shift;
-
-        // temperature
-        let t = bulk.temperature.to_reduced(U::reference_temperature())?;
-
-        // calculate external potential
-        let external_potential = external_potential_3d(
-            dft,
-            [&x, &y, &z],
-            coordinates,
-            sigma_ss,
-            epsilon_ss,
-            cutoff_radius,
-            potential_cutoff,
-            t,
-        )?;
-
-        // initialize convolver
-        let grid = Grid::Cartesian3(x, y, z);
-        let weight_functions = dft.weight_functions(t);
-        let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
-
-        Ok(Self {
-            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None)?,
-            grand_potential: None,
-            solvation_free_energy: None,
-        })
-    }
-}
-
-fn external_potential_3d<U: EosUnit, F: FluidParameters>(
-    functional: &F,
-    axis: [&Axis; 3],
-    coordinates: Array2<f64>,
-    sigma_ss: Array1<f64>,
-    epsilon_ss: Array1<f64>,
-    cutoff_radius: Option<QuantityScalar<U>>,
-    potential_cutoff: Option<f64>,
-    reduced_temperature: f64,
-) -> EosResult<Array4<f64>> {
-    // allocate external potential
-    let m = functional.m();
-    let mut external_potential = Array4::zeros((
-        m.len(),
-        axis[0].grid.len(),
-        axis[1].grid.len(),
-        axis[2].grid.len(),
-    ));
-
-    let cutoff_radius = cutoff_radius
-        .unwrap_or(CUTOFF_RADIUS * U::reference_length())
-        .to_reduced(U::reference_length())?;
-
-    // square cut-off radius
-    let cutoff_radius2 = cutoff_radius.powi(2);
-
-    // calculate external potential
-    let sigma_ff = functional.sigma_ff();
-    let epsilon_k_ff = functional.epsilon_k_ff();
-
-    Zip::indexed(&mut external_potential).par_for_each(|(i, ix, iy, iz), u| {
-        let distance2 = calculate_distance2(
-            [&axis[0].grid[ix], &axis[1].grid[iy], &axis[2].grid[iz]],
-            &coordinates,
-        );
-        let sigma_sf = sigma_ss.mapv(|s| (s + sigma_ff[i]) / 2.0);
-        let epsilon_sf = epsilon_ss.mapv(|e| (e * epsilon_k_ff[i]).sqrt());
-        *u = (0..sigma_ss.len())
-            .map(|alpha| {
-                m[i] * evaluate(
-                    distance2[alpha],
-                    sigma_sf[alpha],
-                    epsilon_sf[alpha],
-                    cutoff_radius2,
-                )
-            })
-            .sum::<f64>()
-            / reduced_temperature
-    });
-
-    let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
-    external_potential.map_inplace(|x| {
-        if *x > potential_cutoff {
-            *x = potential_cutoff
-        }
-    });
-
-    Ok(external_potential)
-}
-
-/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
-fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
-    let sigma_r = sigma.powi(2) / distance2;
-
-    let potential: f64 = if distance2 > cutoff_radius2 {
-        0.0
-    } else if distance2 == 0.0 {
-        f64::INFINITY
-    } else {
-        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
-    };
-
-    potential
-}
-
-/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
-fn calculate_distance2(point: [&f64; 3], coordinates: &Array2<f64>) -> Array1<f64> {
-    Array1::from_shape_fn(coordinates.ncols(), |i| {
-        let rx = coordinates[[0, i]] - point[0];
-        let ry = coordinates[[1, i]] - point[1];
-        let rz = coordinates[[2, i]] - point[2];
-
-        rx.powi(2) + ry.powi(2) + rz.powi(2)
-    })
-}
+#[cfg(feature = "3d_dft")]
+mod solvation_profile;
+#[cfg(feature = "3d_dft")]
+pub use solvation_profile::SolvationProfile;
diff --git a/feos-dft/src/solvation/pair_correlation.rs b/feos-dft/src/solvation/pair_correlation.rs
index dc24a30f1..b8a9e8720 100644
--- a/feos-dft/src/solvation/pair_correlation.rs
+++ b/feos-dft/src/solvation/pair_correlation.rs
@@ -1,9 +1,9 @@
 //! Functionalities for the calculation of pair correlation functions.
-use super::{Axis, DFTProfile, Grid};
 use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, DFT};
 use crate::profile::MAX_POTENTIAL;
 use crate::solver::DFTSolver;
+use crate::{Axis, DFTProfile, Grid};
 use feos_core::{Contributions, EosResult, EosUnit, State};
 use ndarray::prelude::*;
 use quantity::QuantityScalar;
diff --git a/feos-dft/src/solvation/solvation_profile.rs b/feos-dft/src/solvation/solvation_profile.rs
new file mode 100644
index 000000000..fd3e5d1ab
--- /dev/null
+++ b/feos-dft/src/solvation/solvation_profile.rs
@@ -0,0 +1,204 @@
+use crate::adsorption::FluidParameters;
+use crate::convolver::ConvolverFFT;
+use crate::functional::{HelmholtzEnergyFunctional, DFT};
+use crate::geometry::{Axis, Grid};
+use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
+use crate::solver::DFTSolver;
+use feos_core::{Contributions, EosResult, EosUnit, State};
+use ndarray::prelude::*;
+use ndarray::Zip;
+use quantity::{QuantityArray2, QuantityScalar};
+
+/// Density profile and properties of a solute in a inhomogeneous bulk fluid.
+pub struct SolvationProfile<U: EosUnit, F: HelmholtzEnergyFunctional> {
+    pub profile: DFTProfile<U, Ix3, F>,
+    pub grand_potential: Option<QuantityScalar<U>>,
+    pub solvation_free_energy: Option<QuantityScalar<U>>,
+}
+
+impl<U: EosUnit, F: HelmholtzEnergyFunctional> Clone for SolvationProfile<U, F> {
+    fn clone(&self) -> Self {
+        Self {
+            profile: self.profile.clone(),
+            grand_potential: self.grand_potential,
+            solvation_free_energy: self.solvation_free_energy,
+        }
+    }
+}
+
+impl<U: EosUnit, F: HelmholtzEnergyFunctional> SolvationProfile<U, F> {
+    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.grand_potential()?;
+        self.grand_potential = Some(omega);
+
+        // calculate solvation free energy
+        self.solvation_free_energy = Some(
+            (omega + self.profile.bulk.pressure(Contributions::Total) * self.profile.volume())
+                / U::reference_moles(),
+        );
+
+        Ok(())
+    }
+
+    pub fn solve(mut self, solver: Option<&DFTSolver>) -> EosResult<Self> {
+        self.solve_inplace(solver, false)?;
+        Ok(self)
+    }
+}
+
+impl<U: EosUnit, F: HelmholtzEnergyFunctional + FluidParameters> SolvationProfile<U, F> {
+    pub fn new(
+        bulk: &State<U, DFT<F>>,
+        n_grid: [usize; 3],
+        coordinates: QuantityArray2<U>,
+        sigma_ss: Array1<f64>,
+        epsilon_ss: Array1<f64>,
+        system_size: Option<[QuantityScalar<U>; 3]>,
+        cutoff_radius: Option<QuantityScalar<U>>,
+        potential_cutoff: Option<f64>,
+    ) -> EosResult<Self> {
+        let dft: &F = &bulk.eos;
+
+        let system_size = system_size.unwrap_or([40.0 * U::reference_length(); 3]);
+
+        // generate grid
+        let x = Axis::new_cartesian(n_grid[0], system_size[0], None)?;
+        let y = Axis::new_cartesian(n_grid[1], system_size[1], None)?;
+        let z = Axis::new_cartesian(n_grid[2], system_size[2], None)?;
+
+        // move center of geometry of solute to box center
+        let mut coordinates = Array2::from_shape_fn(coordinates.raw_dim(), |(i, j)| {
+            (coordinates.get((i, j)))
+                .to_reduced(U::reference_length())
+                .unwrap()
+        });
+
+        let center = [
+            system_size[0].to_reduced(U::reference_length())? / 2.0,
+            system_size[1].to_reduced(U::reference_length())? / 2.0,
+            system_size[2].to_reduced(U::reference_length())? / 2.0,
+        ];
+
+        let shift: Array2<f64> = Array2::from_shape_fn((3, 1), |(i, _)| {
+            center[i] - coordinates.row(i).sum() / coordinates.ncols() as f64
+        });
+
+        coordinates = coordinates + shift;
+
+        // temperature
+        let t = bulk.temperature.to_reduced(U::reference_temperature())?;
+
+        // calculate external potential
+        let external_potential = external_potential_3d(
+            dft,
+            [&x, &y, &z],
+            coordinates,
+            sigma_ss,
+            epsilon_ss,
+            cutoff_radius,
+            potential_cutoff,
+            t,
+        )?;
+
+        // initialize convolver
+        let grid = Grid::Cartesian3(x, y, z);
+        let weight_functions = dft.weight_functions(t);
+        let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
+
+        Ok(Self {
+            profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None)?,
+            grand_potential: None,
+            solvation_free_energy: None,
+        })
+    }
+}
+
+fn external_potential_3d<U: EosUnit, F: FluidParameters>(
+    functional: &F,
+    axis: [&Axis; 3],
+    coordinates: Array2<f64>,
+    sigma_ss: Array1<f64>,
+    epsilon_ss: Array1<f64>,
+    cutoff_radius: Option<QuantityScalar<U>>,
+    potential_cutoff: Option<f64>,
+    reduced_temperature: f64,
+) -> EosResult<Array4<f64>> {
+    // allocate external potential
+    let m = functional.m();
+    let mut external_potential = Array4::zeros((
+        m.len(),
+        axis[0].grid.len(),
+        axis[1].grid.len(),
+        axis[2].grid.len(),
+    ));
+
+    let cutoff_radius = cutoff_radius
+        .unwrap_or(CUTOFF_RADIUS * U::reference_length())
+        .to_reduced(U::reference_length())?;
+
+    // square cut-off radius
+    let cutoff_radius2 = cutoff_radius.powi(2);
+
+    // calculate external potential
+    let sigma_ff = functional.sigma_ff();
+    let epsilon_k_ff = functional.epsilon_k_ff();
+
+    Zip::indexed(&mut external_potential).par_for_each(|(i, ix, iy, iz), u| {
+        let distance2 = calculate_distance2(
+            [&axis[0].grid[ix], &axis[1].grid[iy], &axis[2].grid[iz]],
+            &coordinates,
+        );
+        let sigma_sf = sigma_ss.mapv(|s| (s + sigma_ff[i]) / 2.0);
+        let epsilon_sf = epsilon_ss.mapv(|e| (e * epsilon_k_ff[i]).sqrt());
+        *u = (0..sigma_ss.len())
+            .map(|alpha| {
+                m[i] * evaluate(
+                    distance2[alpha],
+                    sigma_sf[alpha],
+                    epsilon_sf[alpha],
+                    cutoff_radius2,
+                )
+            })
+            .sum::<f64>()
+            / reduced_temperature
+    });
+
+    let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
+    external_potential.map_inplace(|x| {
+        if *x > potential_cutoff {
+            *x = potential_cutoff
+        }
+    });
+
+    Ok(external_potential)
+}
+
+/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
+fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
+    let sigma_r = sigma.powi(2) / distance2;
+
+    let potential: f64 = if distance2 > cutoff_radius2 {
+        0.0
+    } else if distance2 == 0.0 {
+        f64::INFINITY
+    } else {
+        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
+    };
+
+    potential
+}
+
+/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
+fn calculate_distance2(point: [&f64; 3], coordinates: &Array2<f64>) -> Array1<f64> {
+    Array1::from_shape_fn(coordinates.ncols(), |i| {
+        let rx = coordinates[[0, i]] - point[0];
+        let ry = coordinates[[1, i]] - point[1];
+        let rz = coordinates[[2, i]] - point[2];
+
+        rx.powi(2) + ry.powi(2) + rz.powi(2)
+    })
+}