diff --git a/Cargo.toml b/Cargo.toml
index 814347e11..26b35a219 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -22,7 +22,7 @@ keywords = [
 categories = ["science"]
 
 [workspace.dependencies]
-quantity = "0.12"
+quantity = { git = "https://github.com/itt-ustutt/quantity", branch = "consts_instead_of_typenum" }
 num-dual = "0.13"
 ndarray = "0.17"
 nalgebra = "0.34"
@@ -33,7 +33,6 @@ serde = "1.0"
 serde_json = "1.0"
 indexmap = "2.0"
 itertools = "0.14"
-typenum = "1.16"
 rayon = "1.11"
 petgraph = "0.8"
 rustdct = "0.7"
diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml
index 94717bdb3..f33cd1255 100644
--- a/crates/feos-core/Cargo.toml
+++ b/crates/feos-core/Cargo.toml
@@ -25,7 +25,6 @@ serde = { workspace = true, features = ["derive"] }
 serde_json = { workspace = true, features = ["preserve_order"] }
 indexmap = { workspace = true, features = ["serde"] }
 rayon = { workspace = true, optional = true }
-typenum = { workspace = true }
 itertools = { workspace = true }
 
 [dev-dependencies]
diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs
index 1fe7bb56c..326aff883 100644
--- a/crates/feos-core/src/equation_of_state/residual.rs
+++ b/crates/feos-core/src/equation_of_state/residual.rs
@@ -3,9 +3,10 @@ use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, U1, alloca
 use num_dual::{DualNum, Gradients, partial, partial2, second_derivative, third_derivative};
 use quantity::ad::first_derivative;
 use quantity::*;
-use std::ops::Deref;
+use std::ops::{Deref, Div};
 use std::sync::Arc;
-use typenum::Quot;
+
+type Quot<T1, T2> = <T1 as Div<T2>>::Output;
 
 /// Molar weight of all components.
 ///
diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs
index b883193ca..b51ae0653 100644
--- a/crates/feos-core/src/lib.rs
+++ b/crates/feos-core/src/lib.rs
@@ -3,7 +3,6 @@
 #![warn(clippy::allow_attributes)]
 use quantity::{Quantity, SIUnit};
 use std::ops::{Div, Mul};
-use typenum::Integer;
 
 /// Print messages with level `Verbosity::Iter` or higher.
 #[macro_export]
@@ -132,20 +131,20 @@ const fn powi(x: f64, n: i32) -> f64 {
 /// Conversion between reduced units and SI units.
 pub trait ReferenceSystem {
     type Inner;
-    type T: Integer;
-    type L: Integer;
-    type M: Integer;
-    type I: Integer;
-    type THETA: Integer;
-    type N: Integer;
-    type J: Integer;
-    const FACTOR: f64 = powi(REFERENCE_VALUES[0], Self::T::I32)
-        * powi(REFERENCE_VALUES[1], Self::L::I32)
-        * powi(REFERENCE_VALUES[2], Self::M::I32)
-        * powi(REFERENCE_VALUES[3], Self::I::I32)
-        * powi(REFERENCE_VALUES[4], Self::THETA::I32)
-        * powi(REFERENCE_VALUES[5], Self::N::I32)
-        * powi(REFERENCE_VALUES[6], Self::J::I32);
+    const T: i8;
+    const L: i8;
+    const M: i8;
+    const I: i8;
+    const THETA: i8;
+    const N: i8;
+    const J: i8;
+    const FACTOR: f64 = powi(REFERENCE_VALUES[0], Self::T as i32)
+        * powi(REFERENCE_VALUES[1], Self::L as i32)
+        * powi(REFERENCE_VALUES[2], Self::M as i32)
+        * powi(REFERENCE_VALUES[3], Self::I as i32)
+        * powi(REFERENCE_VALUES[4], Self::THETA as i32)
+        * powi(REFERENCE_VALUES[5], Self::N as i32)
+        * powi(REFERENCE_VALUES[6], Self::J as i32);
 
     fn from_reduced(value: Self::Inner) -> Self
     where
@@ -161,17 +160,25 @@ pub trait ReferenceSystem {
 }
 
 /// Conversion to and from reduced units
-impl<Inner, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
-    ReferenceSystem for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
+impl<
+    Inner,
+    const T: i8,
+    const L: i8,
+    const M: i8,
+    const I: i8,
+    const THETA: i8,
+    const N: i8,
+    const J: i8,
+> ReferenceSystem for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
 {
     type Inner = Inner;
-    type T = T;
-    type L = L;
-    type M = M;
-    type I = I;
-    type THETA = THETA;
-    type N = N;
-    type J = J;
+    const T: i8 = T;
+    const L: i8 = L;
+    const M: i8 = M;
+    const I: i8 = I;
+    const THETA: i8 = THETA;
+    const N: i8 = N;
+    const J: i8 = J;
     fn from_reduced(value: Inner) -> Self
     where
         Inner: Mul<f64, Output = Inner>,
diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs
index dfc600545..da8a12b60 100644
--- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs
+++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs
@@ -12,7 +12,6 @@ use ndarray::Array1;
 use num_dual::linalg::LU;
 use num_dual::{DualNum, DualStruct, Gradients};
 use quantity::{Density, Dimensionless, Moles, Pressure, Quantity, RGAS, SIUnit, Temperature};
-use typenum::{N1, N2, P1, Z0};
 
 const MAX_ITER_INNER: usize = 5;
 const TOL_INNER: f64 = 1e-9;
@@ -94,7 +93,7 @@ impl<D: DualNum<f64> + Copy> TemperatureOrPressure<D> for Temperature<D> {
 // used instead of the explicit unit. Maybe the type is too complicated for the
 // compiler?
 impl<D: DualNum<f64> + Copy> TemperatureOrPressure<D>
-    for Quantity<D, SIUnit<N2, N1, P1, Z0, Z0, Z0, Z0>>
+    for Quantity<D, SIUnit<-2, -1, 1, 0, 0, 0, 0>>
 {
     type Other = Temperature<D>;
     const IDENTIFIER: &'static str = "pressure";
diff --git a/crates/feos-core/src/state/builder.rs b/crates/feos-core/src/state/builder.rs
index 02c61c93d..c418a7c1f 100644
--- a/crates/feos-core/src/state/builder.rs
+++ b/crates/feos-core/src/state/builder.rs
@@ -14,25 +14,24 @@ use quantity::*;
 /// # use quantity::*;
 /// # use nalgebra::dvector;
 /// # use approx::assert_relative_eq;
-/// # use typenum::P3;
 /// # fn main() -> FeosResult<()> {
 /// // Create a state for given T,V,N
 /// let eos = &PengRobinson::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?);
 /// let state = StateBuilder::new(&eos)
 ///                 .temperature(300.0 * KELVIN)
-///                 .volume(12.5 * METER.powi::<P3>())
+///                 .volume(12.5 * METER.powi::<3>())
 ///                 .moles(&(dvector![2.5] * MOL))
 ///                 .build()?;
-/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<P3>());
+/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<3>());
 ///
 /// // For a pure component, the composition does not need to be specified.
 /// let eos = &PengRobinson::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?);
 /// let state = StateBuilder::new(&eos)
 ///                 .temperature(300.0 * KELVIN)
-///                 .volume(12.5 * METER.powi::<P3>())
+///                 .volume(12.5 * METER.powi::<3>())
 ///                 .total_moles(2.5 * MOL)
 ///                 .build()?;
-/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<P3>());
+/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<3>());
 ///
 /// // The state can be constructed without providing any extensive property.
 /// let eos = &PengRobinson::new(
@@ -45,10 +44,10 @@ use quantity::*;
 /// );
 /// let state = StateBuilder::new(&eos)
 ///                 .temperature(300.0 * KELVIN)
-///                 .partial_density(&(dvector![0.2, 0.6] * MOL / METER.powi::<P3>()))
+///                 .partial_density(&(dvector![0.2, 0.6] * MOL / METER.powi::<3>()))
 ///                 .build()?;
 /// assert_relative_eq!(state.molefracs, dvector![0.25, 0.75]);
-/// assert_relative_eq!(state.density, 0.8 * MOL / METER.powi::<P3>());
+/// assert_relative_eq!(state.density, 0.8 * MOL / METER.powi::<3>());
 /// # Ok(())
 /// # }
 /// ```
diff --git a/crates/feos-core/src/state/cache.rs b/crates/feos-core/src/state/cache.rs
index 2ee54a737..b9a1326f7 100644
--- a/crates/feos-core/src/state/cache.rs
+++ b/crates/feos-core/src/state/cache.rs
@@ -1,8 +1,10 @@
 use nalgebra::allocator::Allocator;
 use nalgebra::{DefaultAllocator, Dim, OVector, Scalar};
 use quantity::*;
+use std::ops::Sub;
 use std::sync::OnceLock;
-use typenum::Diff;
+
+type Diff<T1, T2> = <T1 as Sub<T2>>::Output;
 
 #[derive(Clone, Debug)]
 #[expect(clippy::type_complexity)]
diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs
index 4f896d72e..5da7a43c5 100644
--- a/crates/feos-core/src/state/mod.rs
+++ b/crates/feos-core/src/state/mod.rs
@@ -798,12 +798,11 @@ mod critical_point;
 mod tests {
     use super::*;
     use nalgebra::dvector;
-    use typenum::P3;
 
     #[test]
     fn test_validate() {
         let temperature = 298.15 * KELVIN;
-        let density = 3000.0 * MOL / METER.powi::<P3>();
+        let density = 3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![0.03, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_ok());
     }
@@ -811,7 +810,7 @@ mod tests {
     #[test]
     fn test_negative_temperature() {
         let temperature = -298.15 * KELVIN;
-        let density = 3000.0 * MOL / METER.powi::<P3>();
+        let density = 3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![0.03, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_err());
     }
@@ -819,7 +818,7 @@ mod tests {
     #[test]
     fn test_nan_temperature() {
         let temperature = f64::NAN * KELVIN;
-        let density = 3000.0 * MOL / METER.powi::<P3>();
+        let density = 3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![0.03, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_err());
     }
@@ -827,7 +826,7 @@ mod tests {
     #[test]
     fn test_negative_mole_number() {
         let temperature = 298.15 * KELVIN;
-        let density = 3000.0 * MOL / METER.powi::<P3>();
+        let density = 3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![-0.03, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_err());
     }
@@ -835,7 +834,7 @@ mod tests {
     #[test]
     fn test_nan_mole_number() {
         let temperature = 298.15 * KELVIN;
-        let density = 3000.0 * MOL / METER.powi::<P3>();
+        let density = 3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![f64::NAN, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_err());
     }
@@ -843,7 +842,7 @@ mod tests {
     #[test]
     fn test_negative_density() {
         let temperature = 298.15 * KELVIN;
-        let density = -3000.0 * MOL / METER.powi::<P3>();
+        let density = -3000.0 * MOL / METER.powi::<3>();
         let molefracs = dvector![0.01, 0.02, 0.05];
         assert!(validate(temperature, density, &molefracs).is_err());
     }
diff --git a/crates/feos-dft/Cargo.toml b/crates/feos-dft/Cargo.toml
index 9a010add0..2e9e52a7f 100644
--- a/crates/feos-dft/Cargo.toml
+++ b/crates/feos-dft/Cargo.toml
@@ -25,7 +25,6 @@ num-traits = { workspace = true }
 libm = { workspace = true }
 gauss-quad = { workspace = true, optional = true }
 petgraph = { workspace = true }
-typenum = { workspace = true }
 
 feos-core = { workspace = true }
 
diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
index c7602b7d5..faf19e62d 100644
--- a/crates/feos-dft/src/adsorption/pore.rs
+++ b/crates/feos-dft/src/adsorption/pore.rs
@@ -19,12 +19,12 @@ use quantity::{
     Temperature, Volume,
 };
 use rustdct::DctNum;
-use typenum::Diff;
+use std::ops::Sub;
 
 const POTENTIAL_OFFSET: f64 = 2.0;
 const DEFAULT_GRID_POINTS: usize = 2048;
 
-pub type _HenryCoefficient = Diff<_Moles, _Pressure>;
+pub type _HenryCoefficient = <_Moles as Sub<_Pressure>>::Output;
 pub type HenryCoefficient<T> = Quantity<T, _HenryCoefficient>;
 
 /// Parameters required to specify a 1D pore.
diff --git a/crates/feos-dft/src/pdgt.rs b/crates/feos-dft/src/pdgt.rs
index d3a9fa30a..47ad359eb 100644
--- a/crates/feos-dft/src/pdgt.rs
+++ b/crates/feos-dft/src/pdgt.rs
@@ -10,7 +10,9 @@ use quantity::{
     Temperature,
 };
 use std::ops::{Add, AddAssign, Sub};
-use typenum::{Diff, P2, Sum};
+
+type Sum<T1, T2> = <T1 as Add<T2>>::Output;
+type Diff<T1, T2> = <T1 as Sub<T2>>::Output;
 
 type _InfluenceParameter = Diff<Sum<_MolarEnergy, _Area>, _Density>;
 type InfluenceParameter<T> = Quantity<T, _InfluenceParameter>;
@@ -218,7 +220,7 @@ pub trait PdgtFunctionalProperties: HelmholtzEnergyFunctional {
 
             // calculate interfacial width
             let w_temp = integrate_trapezoidal(&rho_r * &*z * z_int, dx);
-            *w = (24.0 * (w_temp - 0.5 * ze.powi::<P2>())).sqrt();
+            *w = (24.0 * (w_temp - 0.5 * ze.powi::<2>())).sqrt();
 
             // shift density profile
             *z -= ze;
diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs
index 7bc9b7ac9..edac0e7b2 100644
--- a/crates/feos-dft/src/profile/mod.rs
+++ b/crates/feos-dft/src/profile/mod.rs
@@ -12,7 +12,6 @@ use num_dual::DualNum;
 use quantity::{_Volume, DEGREES, Density, Length, Moles, Quantity, Temperature, Volume};
 use std::ops::{Add, MulAssign};
 use std::sync::Arc;
-use typenum::Sum;
 
 mod properties;
 
@@ -304,7 +303,7 @@ where
     pub fn integrate<S: Data<Elem = f64>, U>(
         &self,
         profile: &Quantity<ArrayBase<S, D>, U>,
-    ) -> Quantity<f64, Sum<_Volume, U>>
+    ) -> Quantity<f64, <_Volume as Add<U>>::Output>
     where
         _Volume: Add<U>,
     {
@@ -322,7 +321,7 @@ where
     pub fn integrate_comp<S: Data<Elem = f64>, U>(
         &self,
         profile: &Quantity<ArrayBase<S, D::Larger>, U>,
-    ) -> Quantity<DVector<f64>, Sum<_Volume, U>>
+    ) -> Quantity<DVector<f64>, <_Volume as Add<U>>::Output>
     where
         _Volume: Add<U>,
     {
@@ -335,7 +334,7 @@ where
     pub fn integrate_segments<S: Data<Elem = f64>, U>(
         &self,
         profile: &Quantity<ArrayBase<S, D::Larger>, U>,
-    ) -> Quantity<DVector<f64>, Sum<_Volume, U>>
+    ) -> Quantity<DVector<f64>, <_Volume as Add<U>>::Output>
     where
         _Volume: Add<U>,
     {
diff --git a/crates/feos/Cargo.toml b/crates/feos/Cargo.toml
index 5965e4354..28e6e6109 100644
--- a/crates/feos/Cargo.toml
+++ b/crates/feos/Cargo.toml
@@ -26,7 +26,6 @@ serde_json = { workspace = true }
 indexmap = { workspace = true }
 rayon = { workspace = true, optional = true }
 itertools = { workspace = true }
-typenum = { workspace = true }
 
 feos-core = { workspace = true }
 feos-derive = { workspace = true, optional = true }
diff --git a/crates/feos/benches/dft_pore.rs b/crates/feos/benches/dft_pore.rs
index 440c3648e..667c5ff12 100644
--- a/crates/feos/benches/dft_pore.rs
+++ b/crates/feos/benches/dft_pore.rs
@@ -10,7 +10,6 @@ use feos::hard_sphere::{FMTFunctional, FMTVersion};
 use feos::pcsaft::{PcSaftFunctional, PcSaftParameters};
 use nalgebra::dvector;
 use quantity::{ANGSTROM, KELVIN, NAV};
-use typenum::P3;
 
 fn fmt(c: &mut Criterion) {
     let mut group = c.benchmark_group("DFT_pore_fmt");
@@ -23,7 +22,7 @@ fn fmt(c: &mut Criterion) {
         None,
         None,
     );
-    let bulk = State::new_pure(&func, KELVIN, 0.75 / NAV / ANGSTROM.powi::<P3>()).unwrap();
+    let bulk = State::new_pure(&func, KELVIN, 0.75 / NAV / ANGSTROM.powi::<3>()).unwrap();
     group.bench_function("liquid", |b| {
         b.iter(|| pore.initialize(&bulk, None, None).unwrap().solve(None))
     });
diff --git a/crates/feos/benches/dual_numbers.rs b/crates/feos/benches/dual_numbers.rs
index 6ec66a63c..454c6828f 100644
--- a/crates/feos/benches/dual_numbers.rs
+++ b/crates/feos/benches/dual_numbers.rs
@@ -13,7 +13,6 @@ use feos_core::parameter::PureRecord;
 use nalgebra::{DVector, Dyn, dvector};
 use num_dual::{Dual2_64, Dual3_64, Dual64, DualNum, HyperDual64};
 use quantity::*;
-use typenum::P3;
 
 /// Helper function to create a state for given parameters.
 /// - temperature is 80% of critical temperature,
@@ -129,7 +128,7 @@ fn methane_co2_pcsaft(c: &mut Criterion) {
 
     // 230 K, 50 bar, x0 = 0.15
     let temperature = 230.0 * KELVIN;
-    let density = 24.16896 * KILO * MOL / METER.powi::<P3>();
+    let density = 24.16896 * KILO * MOL / METER.powi::<3>();
     let volume = 10.0 * MOL / density;
     let x = dvector![0.15, 0.85];
     let moles = &x * 10.0 * MOL;
diff --git a/crates/feos/benches/state_properties.rs b/crates/feos/benches/state_properties.rs
index 401addf03..916f094ed 100644
--- a/crates/feos/benches/state_properties.rs
+++ b/crates/feos/benches/state_properties.rs
@@ -5,7 +5,6 @@ use feos::core::{Contributions, Residual, State};
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use nalgebra::{DVector, dvector};
 use quantity::*;
-use typenum::P3;
 
 /// Evaluate a property of a state given the EoS, the property to compute,
 /// temperature, volume, moles, and the contributions to consider.
@@ -42,7 +41,7 @@ fn properties_pcsaft(c: &mut Criterion) {
     .unwrap();
     let eos = PcSaft::new(parameters);
     let t = 300.0 * KELVIN;
-    let density = 71.18 * KILO * MOL / METER.powi::<P3>();
+    let density = 71.18 * KILO * MOL / METER.powi::<3>();
     let v = 100.0 * MOL / density;
     let x = dvector![1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0];
     let m = &x * 100.0 * MOL;
@@ -92,7 +91,7 @@ fn properties_pcsaft_polar(c: &mut Criterion) {
     .unwrap();
     let eos = PcSaft::new(parameters);
     let t = 300.0 * KELVIN;
-    let density = 71.18 * KILO * MOL / METER.powi::<P3>();
+    let density = 71.18 * KILO * MOL / METER.powi::<3>();
     let v = 100.0 * MOL / density;
     let x = dvector![1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0];
     let m = &x * 100.0 * MOL;
diff --git a/crates/feos/src/association/mod.rs b/crates/feos/src/association/mod.rs
index a12e12e14..72bd97a38 100644
--- a/crates/feos/src/association/mod.rs
+++ b/crates/feos/src/association/mod.rs
@@ -18,7 +18,7 @@ pub use dft::YuWuAssociationFunctional;
 /// [AssociationStrength::association_strength] multiplies the model-specific
 /// site-site association strength [AssociationStrength::association_strength_ij]
 /// with the contact value of the hard-sphere pair correlation function.
-/// 
+///
 /// For implementations that require a different form,
 /// [AssociationStrength::association_strength] can be overwritten.
 pub trait AssociationStrength: HardSphereProperties {
@@ -463,7 +463,6 @@ mod tests_gc_pcsaft {
     use nalgebra::dvector;
     use num_dual::Dual64;
     use quantity::{METER, MOL, PASCAL, Pressure};
-    use typenum::P3;
 
     #[test]
     fn test_assoc_propanol() {
@@ -471,7 +470,7 @@ mod tests_gc_pcsaft {
         let params = GcPcSaftEosParameters::new(&parameters);
         let contrib = Association::new(50, 1e-10);
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
@@ -499,7 +498,7 @@ mod tests_gc_pcsaft {
         let params = GcPcSaftEosParameters::new(&parameters);
         let contrib = Association::new_cross_association(50, 1e-10);
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
@@ -527,7 +526,7 @@ mod tests_gc_pcsaft {
         let params = GcPcSaftEosParameters::new(&parameters);
         let contrib = Association::new(50, 1e-10);
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = (dvector![1.5, 2.5] * MOL).to_reduced().map(Dual64::from_re);
         let total_moles = moles.sum();
         let molar_volume = volume / total_moles;
diff --git a/crates/feos/src/epcsaft/eos/mod.rs b/crates/feos/src/epcsaft/eos/mod.rs
index 06aa87b63..03a2713d0 100644
--- a/crates/feos/src/epcsaft/eos/mod.rs
+++ b/crates/feos/src/epcsaft/eos/mod.rs
@@ -173,13 +173,12 @@ mod tests {
     use approx::assert_relative_eq;
     use feos_core::*;
     use nalgebra::dvector;
-    use typenum::P3;
 
     #[test]
     fn ideal_gas_pressure() {
         let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
         let t = 200.0 * KELVIN;
-        let v = 1e-3 * METER.powi::<P3>();
+        let v = 1e-3 * METER.powi::<3>();
         let n = dvector![1.0] * MOL;
         let s = State::new_nvt(&&e, t, v, &n).unwrap();
         let p_ig = s.total_moles * RGAS * t / v;
@@ -195,7 +194,7 @@ mod tests {
     fn ideal_gas_heat_capacity_joback() {
         let e = ElectrolytePcSaft::new(propane_parameters()).unwrap();
         let t = 200.0 * KELVIN;
-        let v = 1e-3 * METER.powi::<P3>();
+        let v = 1e-3 * METER.powi::<3>();
         let n = dvector![1.0] * MOL;
         let s = State::new_nvt(&&e, t, v, &n).unwrap();
         let p_ig = s.total_moles * RGAS * t / v;
@@ -282,7 +281,7 @@ mod tests {
         let e2 = ElectrolytePcSaft::new(butane_parameters()).unwrap();
         let e12 = ElectrolytePcSaft::new(propane_butane_parameters()).unwrap();
         let t = 300.0 * KELVIN;
-        let v = 0.02456883872966545 * METER.powi::<P3>();
+        let v = 0.02456883872966545 * METER.powi::<3>();
         let m1 = dvector![2.0] * MOL;
         let m1m = dvector![2.0, 0.0] * MOL;
         let m2m = dvector![0.0, 2.0] * MOL;
diff --git a/crates/feos/src/gc_pcsaft/eos/dispersion.rs b/crates/feos/src/gc_pcsaft/eos/dispersion.rs
index a1a8169bf..3f464507d 100644
--- a/crates/feos/src/gc_pcsaft/eos/dispersion.rs
+++ b/crates/feos/src/gc_pcsaft/eos/dispersion.rs
@@ -129,13 +129,12 @@ mod test {
     use nalgebra::dvector;
     use num_dual::Dual64;
     use quantity::{METER, MOL, PASCAL, Pressure};
-    use typenum::P3;
 
     #[test]
     fn test_dispersion_propane() {
         let parameters = propane();
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
@@ -153,7 +152,7 @@ mod test {
     fn test_dispersion_propanol() {
         let parameters = GcPcSaftEosParameters::new(&propanol());
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
diff --git a/crates/feos/src/gc_pcsaft/eos/hard_chain.rs b/crates/feos/src/gc_pcsaft/eos/hard_chain.rs
index 50b3c4afa..f67dffda4 100644
--- a/crates/feos/src/gc_pcsaft/eos/hard_chain.rs
+++ b/crates/feos/src/gc_pcsaft/eos/hard_chain.rs
@@ -42,13 +42,12 @@ mod test {
     use nalgebra::dvector;
     use num_dual::Dual64;
     use quantity::{METER, MOL, PASCAL, Pressure};
-    use typenum::P3;
 
     #[test]
     fn test_hc_propane() {
         let parameters = propane();
         let temperature = 300.0;
-        let volume = METER.powi::<P3>().to_reduced();
+        let volume = METER.powi::<3>().to_reduced();
         let volume = Dual64::from_re(volume).derivative();
         let moles = (1.5 * MOL).to_reduced();
         let state = StateHD::new(
@@ -70,7 +69,7 @@ mod test {
     fn test_hc_propanol() {
         let parameters = GcPcSaftEosParameters::new(&propanol());
         let temperature = 300.0;
-        let volume = METER.powi::<P3>().to_reduced();
+        let volume = METER.powi::<3>().to_reduced();
         let volume = Dual64::from_re(volume).derivative();
         let moles = (1.5 * MOL).to_reduced();
         let state = StateHD::new(
diff --git a/crates/feos/src/gc_pcsaft/eos/mod.rs b/crates/feos/src/gc_pcsaft/eos/mod.rs
index 4198da197..ae6acf960 100644
--- a/crates/feos/src/gc_pcsaft/eos/mod.rs
+++ b/crates/feos/src/gc_pcsaft/eos/mod.rs
@@ -150,13 +150,12 @@ mod test {
     use nalgebra::dvector;
     use num_dual::Dual64;
     use quantity::{METER, MOL, PASCAL, Pressure};
-    use typenum::P3;
 
     #[test]
     fn hs_propane() {
         let parameters = propane();
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
@@ -174,7 +173,7 @@ mod test {
     fn hs_propanol() {
         let parameters = GcPcSaftEosParameters::new(&propanol());
         let temperature = 300.0;
-        let volume = Dual64::from_re(METER.powi::<P3>().to_reduced()).derivative();
+        let volume = Dual64::from_re(METER.powi::<3>().to_reduced()).derivative();
         let moles = Dual64::from_re((1.5 * MOL).to_reduced());
         let molar_volume = volume / moles;
         let state = StateHD::new(
diff --git a/crates/feos/src/ideal_gas/dippr.rs b/crates/feos/src/ideal_gas/dippr.rs
index 927dd2e07..7553e376c 100644
--- a/crates/feos/src/ideal_gas/dippr.rs
+++ b/crates/feos/src/ideal_gas/dippr.rs
@@ -159,7 +159,6 @@ mod tests {
     use feos_core::{Contributions, EquationOfState, StateBuilder};
     use num_dual::first_derivative;
     use quantity::*;
-    use typenum::P3;
 
     use super::*;
 
@@ -173,7 +172,7 @@ mod tests {
         let dippr = Dippr::new(DipprParameters::new_pure(record.clone())?);
         let eos = EquationOfState::ideal_gas(dippr.clone());
         let temperature = 300.0 * KELVIN;
-        let volume = METER.powi::<P3>();
+        let volume = METER.powi::<3>();
         let state = StateBuilder::new(&&eos)
             .temperature(temperature)
             .volume(volume)
@@ -217,7 +216,7 @@ mod tests {
         let dippr = Dippr::new(DipprParameters::new_pure(record.clone())?);
         let eos = EquationOfState::ideal_gas(dippr.clone());
         let temperature = 300.0 * KELVIN;
-        let volume = METER.powi::<P3>();
+        let volume = METER.powi::<3>();
         let state = StateBuilder::new(&&eos)
             .temperature(temperature)
             .volume(volume)
@@ -263,7 +262,7 @@ mod tests {
         let dippr = Dippr::new(DipprParameters::new_pure(record.clone())?);
         let eos = EquationOfState::ideal_gas(dippr.clone());
         let temperature = 20.0 * KELVIN;
-        let volume = METER.powi::<P3>();
+        let volume = METER.powi::<3>();
         let state = StateBuilder::new(&&eos)
             .temperature(temperature)
             .volume(volume)
diff --git a/crates/feos/src/ideal_gas/joback.rs b/crates/feos/src/ideal_gas/joback.rs
index bef674722..742c6ef2e 100644
--- a/crates/feos/src/ideal_gas/joback.rs
+++ b/crates/feos/src/ideal_gas/joback.rs
@@ -179,7 +179,6 @@ mod tests {
     use nalgebra::dvector;
     use quantity::*;
     use std::collections::HashMap;
-    use typenum::P3;
 
     use super::*;
 
@@ -263,7 +262,7 @@ mod tests {
         let pr = PureRecord::new(Identifier::default(), 1.0, jr);
         let joback = Joback::new(JobackParameters::new_pure(pr)?);
         let eos = EquationOfState::ideal_gas(joback);
-        let state = State::new_pure(&&eos, 1000.0 * KELVIN, 1.0 * MOL / METER.powi::<P3>())?;
+        let state = State::new_pure(&&eos, 1000.0 * KELVIN, 1.0 * MOL / METER.powi::<3>())?;
         assert!(
             ((state.molar_isobaric_heat_capacity(Contributions::IdealGas)
                 / (JOULE / MOL / KELVIN))
@@ -294,7 +293,7 @@ mod tests {
         )?);
         let eos = EquationOfState::ideal_gas(joback.clone());
         let temperature = 300.0 * KELVIN;
-        let volume = METER.powi::<P3>();
+        let volume = METER.powi::<3>();
         let moles = &dvector![1.0, 3.0] * MOL;
         let state = StateBuilder::new(&&eos)
             .temperature(temperature)
diff --git a/crates/feos/src/multiparameter/mod.rs b/crates/feos/src/multiparameter/mod.rs
index f8ab47634..68bbd98b8 100644
--- a/crates/feos/src/multiparameter/mod.rs
+++ b/crates/feos/src/multiparameter/mod.rs
@@ -138,7 +138,6 @@ mod test {
     use nalgebra::{Dyn, SVector, U2, dvector};
     use num_dual::{Dual2Vec, hessian};
     use quantity::{GRAM, KELVIN, KILO, KILOGRAM, METER, MOL, RGAS};
-    use typenum::P3;
 
     use super::*;
 
@@ -262,7 +261,7 @@ mod test {
     #[test]
     fn test_ideal_gas_hack() {
         let t = 647. * KELVIN;
-        let rho = 358. * KILOGRAM / METER.powi::<P3>();
+        let rho = 358. * KILOGRAM / METER.powi::<3>();
         let eos = &water();
         let mw = eos.molar_weight.get(0);
         let moles = dvector![1.8] * MOL;
@@ -270,7 +269,7 @@ mod test {
         let phi_feos = (a_feos / RGAS / moles.sum() / t).into_value();
         println!("A:          {a_feos}");
         println!("phi(feos):  {phi_feos}");
-        let delta = (rho / (eos.rhoc * MOL / METER.powi::<P3>() * mw)).into_value();
+        let delta = (rho / (eos.rhoc * MOL / METER.powi::<3>() * mw)).into_value();
         let tau = (eos.tc * KELVIN / t).into_value();
         let phi = eos.ideal_gas[0]
             .terms
diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs
index 0fca1fa6d..b81d3ac7e 100644
--- a/crates/feos/src/pcsaft/eos/mod.rs
+++ b/crates/feos/src/pcsaft/eos/mod.rs
@@ -9,7 +9,6 @@ use num_dual::{DualNum, partial2};
 use quantity::ad::first_derivative;
 use quantity::*;
 use std::f64::consts::{FRAC_PI_6, PI};
-use typenum::P2;
 
 pub(crate) mod dispersion;
 pub(crate) mod hard_chain;
@@ -241,7 +240,7 @@ impl EntropyScaling for PcSaft {
                 let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
                 5.0 / 16.0 * (mw.get(i) * KB / NAV * temperature / PI).sqrt()
                     / omega22(tr)
-                    / (p.sigma[i] * ANGSTROM).powi::<P2>()
+                    / (p.sigma[i] * ANGSTROM).powi::<2>()
             })
             .collect();
         let mut ce_mix = 0.0 * MILLI * PASCAL * SECOND;
@@ -291,7 +290,7 @@ impl EntropyScaling for PcSaft {
         let res: Vec<_> = (0..self.components())
             .map(|i| {
                 let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value();
-                3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<P2>() / omega11(tr) / (density * NAV)
+                3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<2>() / omega11(tr) / (density * NAV)
                     * (temperature * RGAS / PI / mw.get(i) / p.m[i]).sqrt()
             })
             .collect();
@@ -393,13 +392,12 @@ mod tests {
     use feos_core::*;
     use nalgebra::dvector;
     use quantity::{BAR, KELVIN, METER, PASCAL, RGAS};
-    use typenum::{P2, P3};
 
     #[test]
     fn ideal_gas_pressure() {
         let e = &propane_parameters();
         let t = 200.0 * KELVIN;
-        let v = 1e-3 * METER.powi::<P3>();
+        let v = 1e-3 * METER.powi::<3>();
         let n = dvector![1.0] * MOL;
         let s = State::new_nvt(&e, t, v, &n).unwrap();
         let p_ig = s.total_moles * RGAS * t / v;
@@ -415,7 +413,7 @@ mod tests {
     fn ideal_gas_heat_capacity_joback() {
         let e = &propane_parameters();
         let t = 200.0 * KELVIN;
-        let v = 1e-3 * METER.powi::<P3>();
+        let v = 1e-3 * METER.powi::<3>();
         let n = dvector![1.0] * MOL;
         let s = State::new_nvt(&e, t, v, &n).unwrap();
         let p_ig = s.total_moles * RGAS * t / v;
@@ -504,7 +502,7 @@ mod tests {
         let e2 = &butane_parameters();
         let e12 = &propane_butane_parameters();
         let t = 300.0 * KELVIN;
-        let v = 0.02456883872966545 * METER.powi::<P3>();
+        let v = 0.02456883872966545 * METER.powi::<3>();
         let m1 = dvector![2.0] * MOL;
         let m1m = dvector![2.0, 0.0] * MOL;
         let m2m = dvector![0.0, 2.0] * MOL;
@@ -577,7 +575,7 @@ mod tests {
         let s = State::new_npt(&e, t, p, &n, None)?;
         assert_relative_eq!(
             s.diffusion(),
-            0.01505 * (CENTI * METER).powi::<P2>() / SECOND,
+            0.01505 * (CENTI * METER).powi::<2>() / SECOND,
             epsilon = 1e-5
         );
         assert_relative_eq!(
diff --git a/crates/feos/src/pets/eos/mod.rs b/crates/feos/src/pets/eos/mod.rs
index 0ad191726..f1d1d76cb 100644
--- a/crates/feos/src/pets/eos/mod.rs
+++ b/crates/feos/src/pets/eos/mod.rs
@@ -144,13 +144,12 @@ mod tests {
     use feos_core::{Contributions, PhaseEquilibrium, State, StateHD};
     use nalgebra::dvector;
     use quantity::{BAR, KELVIN, METER, MOL, RGAS};
-    use typenum::P3;
 
     #[test]
     fn ideal_gas_pressure() {
         let e = &Pets::new(argon_parameters());
         let t = 200.0 * KELVIN;
-        let v = 1e-3 * METER.powi::<P3>();
+        let v = 1e-3 * METER.powi::<3>();
         let n = dvector![1.0] * MOL;
         let s = State::new_nvt(&e, t, v, &n).unwrap();
         let p_ig = s.total_moles * RGAS * t / v;
@@ -212,7 +211,7 @@ mod tests {
         let e2 = &Pets::new(krypton_parameters());
         let e12 = &Pets::new(argon_krypton_parameters());
         let t = 300.0 * KELVIN;
-        let v = 0.02456883872966545 * METER.powi::<P3>();
+        let v = 0.02456883872966545 * METER.powi::<3>();
         let m1 = dvector![2.0] * MOL;
         let m1m = dvector![2.0, 0.0] * MOL;
         let m2m = dvector![0.0, 2.0] * MOL;
diff --git a/crates/feos/src/uvtheory/eos/mod.rs b/crates/feos/src/uvtheory/eos/mod.rs
index bc38dcd33..9fa963ae2 100644
--- a/crates/feos/src/uvtheory/eos/mod.rs
+++ b/crates/feos/src/uvtheory/eos/mod.rs
@@ -110,7 +110,6 @@ mod test {
     use feos_core::{FeosResult, State};
     use nalgebra::dvector;
     use quantity::{ANGSTROM, KELVIN, MOL, NAV, RGAS};
-    use typenum::P3;
 
     #[test]
     fn helmholtz_energy_pure_wca() -> FeosResult<()> {
@@ -123,7 +122,7 @@ mod test {
         let reduced_density = 1.0;
         let temperature = reduced_temperature * eps_k * KELVIN;
         let moles = dvector![2.0] * MOL;
-        let volume = (sig * ANGSTROM).powi::<P3>() / reduced_density * NAV * 2.0 * MOL;
+        let volume = (sig * ANGSTROM).powi::<3>() / reduced_density * NAV * 2.0 * MOL;
         let s = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
         let a = (s.residual_molar_helmholtz_energy() / (RGAS * temperature)).into_value();
         assert_relative_eq!(a, 2.972986567516, max_relative = 1e-12); //wca
@@ -147,7 +146,7 @@ mod test {
         let reduced_density = 1.0;
         let temperature = reduced_temperature * eps_k * KELVIN;
         let moles = dvector![2.0] * MOL;
-        let volume = (sig * ANGSTROM).powi::<P3>() / reduced_density * NAV * 2.0 * MOL;
+        let volume = (sig * ANGSTROM).powi::<3>() / reduced_density * NAV * 2.0 * MOL;
         let s = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
 
         let a = (s.residual_molar_helmholtz_energy() / (RGAS * temperature)).into_value();
@@ -173,7 +172,7 @@ mod test {
         let reduced_density = 0.5;
         let temperature = reduced_temperature * eps_k * KELVIN;
         let moles = dvector![2.0] * MOL;
-        let volume = (sig * ANGSTROM).powi::<P3>() / reduced_density * NAV * 2.0 * MOL;
+        let volume = (sig * ANGSTROM).powi::<3>() / reduced_density * NAV * 2.0 * MOL;
         let s = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
         let a = (s.residual_molar_helmholtz_energy() / (RGAS * temperature)).into_value();
         dbg!(a);
@@ -209,7 +208,7 @@ mod test {
         let reduced_density = 1.0;
         let moles = dvector![1.7, 0.3] * MOL;
         let total_moles = moles.sum();
-        let volume = (sig_x * ANGSTROM).powi::<P3>() / reduced_density * NAV * total_moles;
+        let volume = (sig_x * ANGSTROM).powi::<3>() / reduced_density * NAV * total_moles;
 
         // EoS
         let options = UVTheoryOptions {
@@ -242,7 +241,7 @@ mod test {
         let reduced_density = 0.9;
         let moles = dvector![0.4, 0.6] * MOL;
         let total_moles = moles.sum();
-        let volume = (p.sigma[0] * ANGSTROM).powi::<P3>() / reduced_density * NAV * total_moles;
+        let volume = (p.sigma[0] * ANGSTROM).powi::<3>() / reduced_density * NAV * total_moles;
 
         // EoS
         let eos_wca = &UVTheory::new(parameters);
@@ -267,7 +266,7 @@ mod test {
         // state
         let reduced_temperature = 1.5;
         let t_x = reduced_temperature * p.epsilon_k[0] * KELVIN;
-        let sigma_x_3 = (0.4 + 0.6 * 8.0) * ANGSTROM.powi::<P3>();
+        let sigma_x_3 = (0.4 + 0.6 * 8.0) * ANGSTROM.powi::<3>();
         let density = 0.52000000000000002 / sigma_x_3;
         let moles = dvector![0.4, 0.6] * MOL;
         let total_moles = moles.sum();
diff --git a/crates/feos/tests/gc_pcsaft/binary.rs b/crates/feos/tests/gc_pcsaft/binary.rs
index 0419a8987..e9b9ed434 100644
--- a/crates/feos/tests/gc_pcsaft/binary.rs
+++ b/crates/feos/tests/gc_pcsaft/binary.rs
@@ -6,7 +6,6 @@ use feos_core::parameter::IdentifierOption;
 use feos_core::{Contributions, FeosResult, State};
 use nalgebra::dvector;
 use quantity::{KELVIN, METER, MOL};
-use typenum::P3;
 
 #[test]
 fn test_binary() -> FeosResult<()> {
@@ -70,9 +69,9 @@ fn test_polar_term() -> FeosResult<()> {
     let eos1 = &GcPcSaft::new(parameters1);
     let eos2 = &GcPcSaft::new(parameters2);
     let moles = dvector![0.5, 0.5] * MOL;
-    let p1 = State::new_nvt(&eos1, 300.0 * KELVIN, METER.powi::<P3>(), &moles)?
+    let p1 = State::new_nvt(&eos1, 300.0 * KELVIN, METER.powi::<3>(), &moles)?
         .pressure(Contributions::Total);
-    let p2 = State::new_nvt(&eos2, 300.0 * KELVIN, METER.powi::<P3>(), &moles)?
+    let p2 = State::new_nvt(&eos2, 300.0 * KELVIN, METER.powi::<3>(), &moles)?
         .pressure(Contributions::Total);
     println!("{p1} {p2}");
     assert_eq!(p1, p2);
diff --git a/crates/feos/tests/gc_pcsaft/dft.rs b/crates/feos/tests/gc_pcsaft/dft.rs
index 5d245bfb1..a9813c56e 100644
--- a/crates/feos/tests/gc_pcsaft/dft.rs
+++ b/crates/feos/tests/gc_pcsaft/dft.rs
@@ -10,7 +10,6 @@ use feos_dft::{DFTSolver, Geometry};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
-use typenum::P3;
 
 #[test]
 #[allow(non_snake_case)]
@@ -40,7 +39,7 @@ fn test_bulk_implementation() -> Result<(), Box<dyn Error>> {
     let eos = GcPcSaft::new(parameters);
     let func = GcPcSaftFunctional::new(parameters_func);
     let t = 200.0 * KELVIN;
-    let v = 0.002 * METER.powi::<P3>() * NAV / NAV_old;
+    let v = 0.002 * METER.powi::<3>() * NAV / NAV_old;
     let n = dvector![1.5] * MOL;
     let state_eos = State::new_nvt(&&eos, t, v, &n)?;
     let state_func = State::new_nvt(&&func, t, v, &n)?;
@@ -123,7 +122,7 @@ fn test_bulk_association() -> Result<(), Box<dyn Error>> {
     let func = GcPcSaftFunctional::new(func_parameters);
 
     let t = 200.0 * KELVIN;
-    let v = 0.002 * METER.powi::<P3>();
+    let v = 0.002 * METER.powi::<3>();
     let n = dvector![1.5] * MOL;
     let state_eos = State::new_nvt(&&eos, t, v, &n)?;
     let state_func = State::new_nvt(&&func, t, v, &n)?;
@@ -171,13 +170,13 @@ fn test_dft() -> Result<(), Box<dyn Error>> {
 
     assert_relative_eq!(
         vle.vapor().density,
-        12.8820179191167643 * MOL / METER.powi::<P3>() * NAV_old / NAV,
+        12.8820179191167643 * MOL / METER.powi::<3>() * NAV_old / NAV,
         max_relative = 1e-13,
     );
 
     assert_relative_eq!(
         vle.liquid().density,
-        13.2705903446123212 * KILO * MOL / METER.powi::<P3>() * NAV_old / NAV,
+        13.2705903446123212 * KILO * MOL / METER.powi::<3>() * NAV_old / NAV,
         max_relative = 1e-13,
     );
 
diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs
index 097dd607a..6e4ee3156 100644
--- a/crates/feos/tests/pcsaft/critical_point.rs
+++ b/crates/feos/tests/pcsaft/critical_point.rs
@@ -5,7 +5,6 @@ use feos_core::parameter::IdentifierOption;
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
-use typenum::P3;
 
 #[test]
 fn test_critical_point_pure() -> Result<(), Box<dyn Error>> {
@@ -21,7 +20,7 @@ fn test_critical_point_pure() -> Result<(), Box<dyn Error>> {
     assert_relative_eq!(cp.temperature, 375.12441 * KELVIN, max_relative = 1e-8);
     assert_relative_eq!(
         cp.density,
-        4733.00377 * MOL / METER.powi::<P3>(),
+        4733.00377 * MOL / METER.powi::<3>(),
         max_relative = 1e-6
     );
     Ok(())
@@ -42,7 +41,7 @@ fn test_critical_point_mix() -> Result<(), Box<dyn Error>> {
     assert_relative_eq!(cp.temperature, 407.93481 * KELVIN, max_relative = 1e-8);
     assert_relative_eq!(
         cp.density,
-        4265.50745 * MOL / METER.powi::<P3>(),
+        4265.50745 * MOL / METER.powi::<3>(),
         max_relative = 1e-6
     );
     Ok(())
diff --git a/crates/feos/tests/pcsaft/dft.rs b/crates/feos/tests/pcsaft/dft.rs
index b6c6f3880..e8d057e47 100644
--- a/crates/feos/tests/pcsaft/dft.rs
+++ b/crates/feos/tests/pcsaft/dft.rs
@@ -12,7 +12,6 @@ use nalgebra::dvector;
 use ndarray::Axis;
 use quantity::*;
 use std::error::Error;
-use typenum::P3;
 
 fn parameters(comp: &str) -> FeosResult<PcSaftParameters> {
     PcSaftParameters::from_json(
@@ -35,7 +34,7 @@ fn test_bulk_implementations() -> Result<(), Box<dyn Error>> {
     let func_full =
         PcSaftFunctional::new_full(parameters("water_np")?, FMTVersion::KierlikRosinberg);
     let t = 300.0 * KELVIN;
-    let v = 0.002 * METER.powi::<P3>() * NAV / NAV_old;
+    let v = 0.002 * METER.powi::<3>() * NAV / NAV_old;
     let n = dvector![1.5] * MOL;
     let state = State::new_nvt(&&eos, t, v, &n)?;
     let state_pure = State::new_nvt(&&func_pure, t, v, &n)?;
@@ -135,7 +134,7 @@ fn test_dft_propane() -> Result<(), Box<dyn Error>> {
         (&func_full_vec).solve_pdgt(&vle_full_vec, 198, 0, None)?.1
     );
 
-    let vapor_density = 12.2557486248527745 * MOL / METER.powi::<P3>() * NAV_old / NAV;
+    let vapor_density = 12.2557486248527745 * MOL / METER.powi::<3>() * NAV_old / NAV;
     assert_relative_eq!(
         vle_pure.vapor().density,
         vapor_density,
@@ -152,7 +151,7 @@ fn test_dft_propane() -> Result<(), Box<dyn Error>> {
         max_relative = 1e-13,
     );
 
-    let liquid_density = 13.8941749145544549 * KILO * MOL / METER.powi::<P3>() * NAV_old / NAV;
+    let liquid_density = 13.8941749145544549 * KILO * MOL / METER.powi::<3>() * NAV_old / NAV;
     assert_relative_eq!(
         vle_pure.liquid().density,
         liquid_density,
@@ -254,7 +253,7 @@ fn test_dft_water() -> Result<(), Box<dyn Error>> {
         vle_full_vec.liquid().density
     );
 
-    let vapor_density = 75.8045715345905222 * MOL / METER.powi::<P3>() * NAV_old / NAV;
+    let vapor_density = 75.8045715345905222 * MOL / METER.powi::<3>() * NAV_old / NAV;
     assert_relative_eq!(
         vle_pure.vapor().density,
         vapor_density,
@@ -266,7 +265,7 @@ fn test_dft_water() -> Result<(), Box<dyn Error>> {
         max_relative = 1e-13,
     );
 
-    let liquid_density = 47.8480850281608454 * KILO * MOL / METER.powi::<P3>() * NAV_old / NAV;
+    let liquid_density = 47.8480850281608454 * KILO * MOL / METER.powi::<3>() * NAV_old / NAV;
     assert_relative_eq!(
         vle_pure.liquid().density,
         liquid_density,
diff --git a/crates/feos/tests/pcsaft/state_creation_mixture.rs b/crates/feos/tests/pcsaft/state_creation_mixture.rs
index e30d4058c..b7b93c73a 100644
--- a/crates/feos/tests/pcsaft/state_creation_mixture.rs
+++ b/crates/feos/tests/pcsaft/state_creation_mixture.rs
@@ -6,7 +6,6 @@ use feos_core::{Contributions, EquationOfState, FeosResult, StateBuilder};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
-use typenum::P3;
 
 fn propane_butane_parameters() -> FeosResult<(PcSaftParameters, Vec<Joback>)> {
     let saft = PcSaftParameters::from_json(
@@ -61,7 +60,7 @@ fn pressure_entropy_molefracs() -> Result<(), Box<dyn Error>> {
 fn volume_temperature_molefracs() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_butane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
-    let volume = 1.5e-3 * METER.powi::<P3>();
+    let volume = 1.5e-3 * METER.powi::<3>();
     let moles = MOL;
     let x = dvector![0.3, 0.7];
     let state = StateBuilder::new(&&saft)
@@ -79,7 +78,7 @@ fn temperature_partial_density() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_butane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
     let x = dvector![0.3, 0.7];
-    let partial_density = x.clone() * MOL / METER.powi::<P3>();
+    let partial_density = x.clone() * MOL / METER.powi::<3>();
     let density = partial_density.sum();
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
@@ -98,7 +97,7 @@ fn temperature_density_molefracs() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_butane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
     let x = dvector![0.3, 0.7];
-    let density = MOL / METER.powi::<P3>();
+    let density = MOL / METER.powi::<3>();
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
         .density(density)
diff --git a/crates/feos/tests/pcsaft/state_creation_pure.rs b/crates/feos/tests/pcsaft/state_creation_pure.rs
index c8fea1780..d5bc9769e 100644
--- a/crates/feos/tests/pcsaft/state_creation_pure.rs
+++ b/crates/feos/tests/pcsaft/state_creation_pure.rs
@@ -7,7 +7,6 @@ use feos_core::{
 };
 use quantity::*;
 use std::error::Error;
-use typenum::P3;
 
 fn propane_parameters() -> FeosResult<(PcSaftParameters, Vec<Joback>)> {
     let saft = PcSaftParameters::from_json(
@@ -29,7 +28,7 @@ fn propane_parameters() -> FeosResult<(PcSaftParameters, Vec<Joback>)> {
 fn temperature_volume() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
-    let volume = 1.5e-3 * METER.powi::<P3>();
+    let volume = 1.5e-3 * METER.powi::<3>();
     let moles = MOL;
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
@@ -44,7 +43,7 @@ fn temperature_volume() -> Result<(), Box<dyn Error>> {
 fn temperature_density() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
-    let density = MOL / METER.powi::<P3>();
+    let density = MOL / METER.powi::<3>();
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
         .density(density)
@@ -58,7 +57,7 @@ fn temperature_total_moles_volume() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
     let total_moles = MOL;
-    let volume = METER.powi::<P3>();
+    let volume = METER.powi::<3>();
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
         .volume(volume)
@@ -74,7 +73,7 @@ fn temperature_total_moles_density() -> Result<(), Box<dyn Error>> {
     let saft = PcSaft::new(propane_parameters()?.0);
     let temperature = 300.0 * KELVIN;
     let total_moles = MOL;
-    let density = MOL / METER.powi::<P3>();
+    let density = MOL / METER.powi::<3>();
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
         .density(density)
@@ -131,7 +130,7 @@ fn pressure_temperature_initial_density() -> Result<(), Box<dyn Error>> {
     let state = StateBuilder::new(&&saft)
         .temperature(temperature)
         .pressure(pressure)
-        .initial_density(MOL / METER.powi::<P3>())
+        .initial_density(MOL / METER.powi::<3>())
         .build()?;
     assert_relative_eq!(
         state.pressure(Contributions::Total),
diff --git a/crates/feos/tests/saftvrmie/critical_properties.rs b/crates/feos/tests/saftvrmie/critical_properties.rs
index 3a4a585ed..c78eb110a 100644
--- a/crates/feos/tests/saftvrmie/critical_properties.rs
+++ b/crates/feos/tests/saftvrmie/critical_properties.rs
@@ -3,12 +3,11 @@ use feos::saftvrmie::{SaftVRMie, test_utils};
 use feos_core::{SolverOptions, State};
 use quantity::*;
 use std::collections::HashMap;
-use typenum::P3;
 
 /// Critical data reported in Lafitte et al.
 pub fn critical_data() -> HashMap<&'static str, (Temperature, Pressure, MassDensity)> {
     let mut data = HashMap::new();
-    let kg_m3 = KILOGRAM / METER.powi::<P3>();
+    let kg_m3 = KILOGRAM / METER.powi::<3>();
     let mpa = MEGA * PASCAL;
     let k = KELVIN;
 
diff --git a/py-feos/Cargo.toml b/py-feos/Cargo.toml
index 9bc040b8f..9affe6b99 100644
--- a/py-feos/Cargo.toml
+++ b/py-feos/Cargo.toml
@@ -37,7 +37,6 @@ serde_json = { workspace = true }
 indexmap = { workspace = true }
 rayon = { workspace = true, optional = true }
 itertools = { workspace = true }
-typenum = { workspace = true }
 paste = { workspace = true }
 
 feos = { workspace = true }
diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs
index 20fd4ed19..426d5b099 100644
--- a/py-feos/src/dft/profile.rs
+++ b/py-feos/src/dft/profile.rs
@@ -1,5 +1,7 @@
 macro_rules! impl_profile {
     ($struct:ident, $arr:ident, $arr2:ident, $si_arr:ident, $si_arr2:ident, $py_arr2:ident, [$([$ind:expr, $ax:ident]),+]$(, $si_arr3:ident)?) => {
+
+
         #[pymethods]
         impl $struct {
             /// Calculate the residual for the given profile.
@@ -118,7 +120,7 @@ macro_rules! impl_profile {
             fn entropy_density(
                 &mut self,
                 contributions: PyContributions,
-            ) -> PyResult<typenum::Quot<Entropy<$si_arr<f64>>, Volume>> {
+            ) -> PyResult<<Entropy<$si_arr<f64>> as std::ops::Div<Volume>>::Output> {
                 Ok(self.0.profile.entropy_density(contributions.into()).map_err(PyFeosError::from)?)
             }
 
@@ -166,33 +168,33 @@ macro_rules! impl_profile {
             }
             $(
                 #[getter]
-                fn get_drho_dmu(&self) -> PyResult<typenum::Quot<Density<$si_arr3<f64>>, MolarEnergy>> {
+                fn get_drho_dmu(&self) -> PyResult<<Density<$si_arr3<f64>> as std::ops::Div<MolarEnergy>>::Output> {
                     Ok(self.0.profile.drho_dmu().map_err(PyFeosError::from)?)
                 }
             )?
 
             #[getter]
-            fn get_dn_dmu(&self) -> PyResult<typenum::Quot<Moles<DMatrix<f64>>, MolarEnergy>> {
+            fn get_dn_dmu(&self) -> PyResult<<Moles<DMatrix<f64>> as std::ops::Div<MolarEnergy>>::Output> {
                 Ok(self.0.profile.dn_dmu().map_err(PyFeosError::from)?)
             }
 
             #[getter]
-            fn get_drho_dp(&self) -> PyResult<typenum::Quot<Density<$si_arr2<f64>>, Pressure>> {
+            fn get_drho_dp(&self) -> PyResult<<Density<$si_arr2<f64>> as std::ops::Div<Pressure>>::Output> {
                 Ok(self.0.profile.drho_dp().map_err(PyFeosError::from)?)
             }
 
             #[getter]
-            fn get_dn_dp(&self) -> PyResult<typenum::Quot<Moles<DVector<f64>>, Pressure>> {
+            fn get_dn_dp(&self) -> PyResult<<Moles<DVector<f64>> as std::ops::Div<Pressure>>::Output> {
                 Ok(self.0.profile.dn_dp().map_err(PyFeosError::from)?)
             }
 
             #[getter]
-            fn get_drho_dt(&self) -> PyResult<typenum::Quot<Density<$si_arr2<f64>>, Temperature>> {
+            fn get_drho_dt(&self) -> PyResult<<Density<$si_arr2<f64>> as std::ops::Div<Temperature>>::Output> {
                 Ok(self.0.profile.drho_dt().map_err(PyFeosError::from)?)
             }
 
             #[getter]
-            fn get_dn_dt(&self) -> PyResult<typenum::Quot<Moles<DVector<f64>>, Temperature>> {
+            fn get_dn_dt(&self) -> PyResult<<Moles<DVector<f64>> as std::ops::Div<Temperature>>::Output> {
                 Ok(self.0.profile.dn_dt().map_err(PyFeosError::from)?)
             }
         }
diff --git a/py-feos/src/eos/mod.rs b/py-feos/src/eos/mod.rs
index 86b6a1447..1147ac020 100644
--- a/py-feos/src/eos/mod.rs
+++ b/py-feos/src/eos/mod.rs
@@ -7,8 +7,10 @@ use nalgebra::{DVector, DVectorView, Dyn};
 use numpy::{PyArray1, PyReadonlyArray1, ToPyArray};
 use pyo3::prelude::*;
 use quantity::*;
+use std::ops::Div;
 use std::sync::Arc;
-use typenum::Quot;
+
+type Quot<T1, T2> = <T1 as Div<T2>>::Output;
 
 mod constructors;
 #[cfg(feature = "epcsaft")]
diff --git a/py-feos/src/phase_equilibria.rs b/py-feos/src/phase_equilibria.rs
index 77f95df3d..c3786e300 100644
--- a/py-feos/src/phase_equilibria.rs
+++ b/py-feos/src/phase_equilibria.rs
@@ -1,10 +1,10 @@
 use crate::{
-    eos::{parse_molefracs, PyEquationOfState},
+    PyVerbosity,
+    eos::{PyEquationOfState, parse_molefracs},
     error::PyFeosError,
     ideal_gas::IdealGasModel,
     residual::ResidualModel,
     state::{PyContributions, PyState, PyStateVec},
-    PyVerbosity,
 };
 use feos_core::{
     Contributions, EquationOfState, PhaseDiagram, PhaseDiagramHetero, PhaseEquilibrium, ResidualDyn,
@@ -17,7 +17,6 @@ use pyo3::prelude::*;
 use quantity::*;
 use std::ops::Deref;
 use std::sync::Arc;
-use typenum::P3;
 
 /// A thermodynamic two phase equilibrium state.
 #[pyclass(name = "PhaseEquilibrium")]
@@ -1082,7 +1081,7 @@ impl PyPhaseDiagram {
             self.0
                 .liquid()
                 .density()
-                .convert_to(MOL / METER.powi::<P3>())
+                .convert_to(MOL / METER.powi::<3>())
                 .into_raw_vec_and_offset()
                 .0,
         );
@@ -1091,7 +1090,7 @@ impl PyPhaseDiagram {
             self.0
                 .vapor()
                 .density()
-                .convert_to(MOL / METER.powi::<P3>())
+                .convert_to(MOL / METER.powi::<3>())
                 .into_raw_vec_and_offset()
                 .0,
         );
@@ -1137,7 +1136,7 @@ impl PyPhaseDiagram {
                 self.0
                     .liquid()
                     .mass_density()
-                    .convert_to(KILOGRAM / METER.powi::<P3>())
+                    .convert_to(KILOGRAM / METER.powi::<3>())
                     .into_raw_vec_and_offset()
                     .0,
             );
@@ -1146,7 +1145,7 @@ impl PyPhaseDiagram {
                 self.0
                     .vapor()
                     .mass_density()
-                    .convert_to(KILOGRAM / METER.powi::<P3>())
+                    .convert_to(KILOGRAM / METER.powi::<3>())
                     .into_raw_vec_and_offset()
                     .0,
             );
diff --git a/py-feos/src/state.rs b/py-feos/src/state.rs
index 2b5492d85..e7b54211c 100644
--- a/py-feos/src/state.rs
+++ b/py-feos/src/state.rs
@@ -1,7 +1,7 @@
 use crate::eos::parse_molefracs;
 use crate::{
-    eos::PyEquationOfState, error::PyFeosError, ideal_gas::IdealGasModel, residual::ResidualModel,
-    PyVerbosity,
+    PyVerbosity, eos::PyEquationOfState, error::PyFeosError, ideal_gas::IdealGasModel,
+    residual::ResidualModel,
 };
 use feos_core::{
     Contributions, DensityInitialization, EquationOfState, FeosError, ResidualDyn, State, StateVec,
@@ -13,9 +13,10 @@ use pyo3::exceptions::{PyIndexError, PyValueError};
 use pyo3::prelude::*;
 use quantity::*;
 use std::collections::HashMap;
-use std::ops::{Deref, Neg, Sub};
+use std::ops::{Deref, Div, Neg, Sub};
 use std::sync::Arc;
-use typenum::{Quot, P3};
+
+type Quot<T1, T2> = <T1 as Div<T2>>::Output;
 
 type DpDn<T> = Quantity<T, <_Pressure as Sub<_Moles>>::Output>;
 type InvT<T> = Quantity<T, <_Temperature as Neg>::Output>;
@@ -1633,7 +1634,7 @@ impl PyStateVec {
             String::from("density"),
             states
                 .density()
-                .convert_to(MOL / METER.powi::<P3>())
+                .convert_to(MOL / METER.powi::<3>())
                 .into_raw_vec_and_offset()
                 .0,
         );
@@ -1658,7 +1659,7 @@ impl PyStateVec {
                 String::from("mass density"),
                 states
                     .mass_density()
-                    .convert_to(KILOGRAM / METER.powi::<P3>())
+                    .convert_to(KILOGRAM / METER.powi::<3>())
                     .into_raw_vec_and_offset()
                     .0,
             );