diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
index ccd6fff7d..0c615714a 100644
--- a/feos-core/src/cubic.rs
+++ b/feos-core/src/cubic.rs
@@ -13,6 +13,7 @@ use num_dual::DualNum;
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
 use std::fmt;
+use std::marker::PhantomData;
 use std::sync::Arc;
 
 const KB_A3: f64 = 13806490.0;
@@ -150,8 +151,8 @@ struct PengRobinsonContribution {
     parameters: Arc<PengRobinsonParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for PengRobinsonContribution {
-    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PhantomData<D>, D> for PengRobinsonContribution {
+    fn helmholtz_energy(&self, state: &StateHD<D>, _: &PhantomData<D>) -> D {
         // temperature dependent a parameter
         let p = &self.parameters;
         let x = &state.molefracs;
@@ -188,13 +189,13 @@ pub struct PengRobinson {
     /// Parameters
     parameters: Arc<PengRobinsonParameters>,
     /// Non-ideal contributions to the Helmholtz energy
-    contributions: Vec<Box<dyn HelmholtzEnergy>>,
+    contributions: Vec<Box<dyn HelmholtzEnergy<Self>>>,
 }
 
 impl PengRobinson {
     /// Create a new equation of state from a set of parameters.
     pub fn new(parameters: Arc<PengRobinsonParameters>) -> Self {
-        let contributions: Vec<Box<dyn HelmholtzEnergy>> =
+        let contributions: Vec<Box<dyn HelmholtzEnergy<Self>>> =
             vec![Box::new(PengRobinsonContribution {
                 parameters: parameters.clone(),
             })];
@@ -222,12 +223,19 @@ impl Components for PengRobinson {
 }
 
 impl Residual for PengRobinson {
+    type Properties<D> = PhantomData<D>;
+    type Inner = Self;
+
+    fn properties<D: DualNum<f64>>(&self, _: D) -> PhantomData<D> {
+        PhantomData
+    }
+
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
         let b = (moles * &self.parameters.b).sum() / moles.sum();
         0.9 / b
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self>>] {
         &self.contributions
     }
 
diff --git a/feos-core/src/equation_of_state/helmholtz_energy.rs b/feos-core/src/equation_of_state/helmholtz_energy.rs
index 11d2c2703..af438b4c1 100644
--- a/feos-core/src/equation_of_state/helmholtz_energy.rs
+++ b/feos-core/src/equation_of_state/helmholtz_energy.rs
@@ -1,4 +1,4 @@
-use crate::StateHD;
+use crate::{Residual, StateHD};
 use num_dual::*;
 use std::fmt;
 
@@ -9,49 +9,52 @@ use std::fmt;
 /// the specific types in the supertraits of [HelmholtzEnergy]
 /// so that the implementor can be used as a Helmholtz energy
 /// contribution in the equation of state.
-pub trait HelmholtzEnergyDual<D: DualNum<f64>> {
+pub trait HelmholtzEnergyDual<P, D: DualNum<f64>> {
     /// The Helmholtz energy contribution $\beta A$ of a given state in reduced units.
-    fn helmholtz_energy(&self, state: &StateHD<D>) -> D;
+    fn helmholtz_energy(&self, state: &StateHD<D>, properties: &P) -> D;
 }
 
 /// Object safe version of the [HelmholtzEnergyDual] trait.
 ///
 /// The trait is implemented automatically for every struct that implements
 /// the supertraits.
-pub trait HelmholtzEnergy:
-    HelmholtzEnergyDual<f64>
-    + HelmholtzEnergyDual<Dual64>
-    + HelmholtzEnergyDual<Dual<DualSVec64<3>, f64>>
-    + HelmholtzEnergyDual<HyperDual64>
-    + HelmholtzEnergyDual<Dual2_64>
-    + HelmholtzEnergyDual<Dual3_64>
-    + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
-    + HelmholtzEnergyDual<HyperDual<DualSVec64<2>, f64>>
-    + HelmholtzEnergyDual<HyperDual<DualSVec64<3>, f64>>
-    + HelmholtzEnergyDual<Dual2<Dual64, f64>>
-    + HelmholtzEnergyDual<Dual3<Dual64, f64>>
-    + HelmholtzEnergyDual<Dual3<DualSVec64<2>, f64>>
-    + HelmholtzEnergyDual<Dual3<DualSVec64<3>, f64>>
+pub trait HelmholtzEnergy<P: Residual>: HelmholtzEnergyDual<P::Properties<f64>, f64>
+    + HelmholtzEnergyDual<P::Properties<Dual64>, Dual64>
+    + HelmholtzEnergyDual<P::Properties<Dual<DualSVec64<3>, f64>>, Dual<DualSVec64<3>, f64>>
+    + HelmholtzEnergyDual<P::Properties<HyperDual64>, HyperDual64>
+    + HelmholtzEnergyDual<P::Properties<Dual2_64>, Dual2_64>
+    + HelmholtzEnergyDual<P::Properties<Dual3_64>, Dual3_64>
+    + HelmholtzEnergyDual<P::Properties<HyperDual<Dual64, f64>>, HyperDual<Dual64, f64>>
+    + HelmholtzEnergyDual<P::Properties<HyperDual<DualSVec64<2>, f64>>, HyperDual<DualSVec64<2>, f64>>
+    + HelmholtzEnergyDual<P::Properties<HyperDual<DualSVec64<3>, f64>>, HyperDual<DualSVec64<3>, f64>>
+    + HelmholtzEnergyDual<P::Properties<Dual2<Dual64, f64>>, Dual2<Dual64, f64>>
+    + HelmholtzEnergyDual<P::Properties<Dual3<Dual64, f64>>, Dual3<Dual64, f64>>
+    + HelmholtzEnergyDual<P::Properties<Dual3<DualSVec64<2>, f64>>, Dual3<DualSVec64<2>, f64>>
+    + HelmholtzEnergyDual<P::Properties<Dual3<DualSVec64<3>, f64>>, Dual3<DualSVec64<3>, f64>>
     + fmt::Display
     + Send
     + Sync
 {
 }
 
-impl<T> HelmholtzEnergy for T where
-    T: HelmholtzEnergyDual<f64>
-        + HelmholtzEnergyDual<Dual64>
-        + HelmholtzEnergyDual<Dual<DualSVec64<3>, f64>>
-        + HelmholtzEnergyDual<HyperDual64>
-        + HelmholtzEnergyDual<Dual2_64>
-        + HelmholtzEnergyDual<Dual3_64>
-        + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
-        + HelmholtzEnergyDual<HyperDual<DualSVec64<2>, f64>>
-        + HelmholtzEnergyDual<HyperDual<DualSVec64<3>, f64>>
-        + HelmholtzEnergyDual<Dual2<Dual64, f64>>
-        + HelmholtzEnergyDual<Dual3<Dual64, f64>>
-        + HelmholtzEnergyDual<Dual3<DualSVec64<2>, f64>>
-        + HelmholtzEnergyDual<Dual3<DualSVec64<3>, f64>>
+impl<T, P: Residual> HelmholtzEnergy<P> for T where
+    T: HelmholtzEnergyDual<P::Properties<f64>, f64>
+        + HelmholtzEnergyDual<P::Properties<Dual64>, Dual64>
+        + HelmholtzEnergyDual<P::Properties<Dual<DualSVec64<3>, f64>>, Dual<DualSVec64<3>, f64>>
+        + HelmholtzEnergyDual<P::Properties<HyperDual64>, HyperDual64>
+        + HelmholtzEnergyDual<P::Properties<Dual2_64>, Dual2_64>
+        + HelmholtzEnergyDual<P::Properties<Dual3_64>, Dual3_64>
+        + HelmholtzEnergyDual<P::Properties<HyperDual<Dual64, f64>>, HyperDual<Dual64, f64>>
+        + HelmholtzEnergyDual<
+            P::Properties<HyperDual<DualSVec64<2>, f64>>,
+            HyperDual<DualSVec64<2>, f64>,
+        > + HelmholtzEnergyDual<
+            P::Properties<HyperDual<DualSVec64<3>, f64>>,
+            HyperDual<DualSVec64<3>, f64>,
+        > + HelmholtzEnergyDual<P::Properties<Dual2<Dual64, f64>>, Dual2<Dual64, f64>>
+        + HelmholtzEnergyDual<P::Properties<Dual3<Dual64, f64>>, Dual3<Dual64, f64>>
+        + HelmholtzEnergyDual<P::Properties<Dual3<DualSVec64<2>, f64>>, Dual3<DualSVec64<2>, f64>>
+        + HelmholtzEnergyDual<P::Properties<Dual3<DualSVec64<3>, f64>>, Dual3<DualSVec64<3>, f64>>
         + fmt::Display
         + Send
         + Sync
diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs
index a1ac42a64..bca430600 100644
--- a/feos-core/src/equation_of_state/mod.rs
+++ b/feos-core/src/equation_of_state/mod.rs
@@ -79,11 +79,20 @@ impl<I: IdealGas, R: Components + Sync + Send> IdealGas for EquationOfState<I, R
 }
 
 impl<I: IdealGas, R: Residual> Residual for EquationOfState<I, R> {
+    type Properties<D> = R::Properties<D>;
+    type Inner = R::Inner;
+
+    fn properties<D: num_dual::DualNum<f64>>(
+        &self,
+        temperature: D,
+    ) -> <Self::Inner as Residual>::Properties<D> {
+        self.residual.properties(temperature)
+    }
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
         self.residual.compute_max_density(moles)
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self::Inner>>] {
         self.residual.contributions()
     }
 
diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs
index e8e2819e2..82b8d64b4 100644
--- a/feos-core/src/equation_of_state/residual.rs
+++ b/feos-core/src/equation_of_state/residual.rs
@@ -5,10 +5,20 @@ use crate::{EosError, EosResult};
 use ndarray::prelude::*;
 use num_dual::*;
 use num_traits::{One, Zero};
+use std::marker::PhantomData;
 use std::ops::Div;
 
 /// A reisdual Helmholtz energy model.
 pub trait Residual: Components + Send + Sync {
+    type Properties<D>;
+    /// For wrapper structs, the inner type. In basic cases `Self`.
+    type Inner: Residual;
+
+    fn properties<D: DualNum<f64>>(
+        &self,
+        temperature: D,
+    ) -> <Self::Inner as Residual>::Properties<D>;
+
     /// Return the maximum density in Angstrom^-3.
     ///
     /// This value is used as an estimate for a liquid phase for phase
@@ -19,7 +29,7 @@ pub trait Residual: Components + Send + Sync {
 
     /// Return a slice of the individual contributions (excluding the ideal gas)
     /// of the equation of state.
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>];
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self::Inner>>];
 
     /// Molar weight of all components.
     ///
@@ -29,11 +39,13 @@ pub trait Residual: Components + Send + Sync {
     /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
     fn evaluate_residual<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D
     where
-        dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
+        dyn HelmholtzEnergy<Self::Inner>:
+            HelmholtzEnergyDual<<Self::Inner as Residual>::Properties<D>, D>,
     {
+        let properties = self.properties(state.temperature);
         self.contributions()
             .iter()
-            .map(|c| c.helmholtz_energy(state))
+            .map(|c| c.helmholtz_energy(state, &properties))
             .sum()
     }
 
@@ -44,11 +56,13 @@ pub trait Residual: Components + Send + Sync {
         state: &StateHD<D>,
     ) -> Vec<(String, D)>
     where
-        dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
+        dyn HelmholtzEnergy<Self::Inner>:
+            HelmholtzEnergyDual<<Self::Inner as Residual>::Properties<D>, D>,
     {
+        let properties = self.properties(state.temperature);
         self.contributions()
             .iter()
-            .map(|c| (c.to_string(), c.helmholtz_energy(state)))
+            .map(|c| (c.to_string(), c.helmholtz_energy(state, &properties)))
             .collect()
     }
 
@@ -193,11 +207,18 @@ impl Components for NoResidual {
 }
 
 impl Residual for NoResidual {
+    type Properties<D> = PhantomData<D>;
+    type Inner = Self;
+
+    fn properties<D: DualNum<f64>>(&self, _: D) -> PhantomData<D> {
+        PhantomData
+    }
+
     fn compute_max_density(&self, _: &Array1<f64>) -> f64 {
         1.0
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self>>] {
         &[]
     }
 
diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs
index 913406f10..d95e7f2cd 100644
--- a/feos-dft/src/functional.rs
+++ b/feos-dft/src/functional.rs
@@ -15,6 +15,7 @@ use petgraph::graph::{Graph, UnGraph};
 use petgraph::visit::EdgeRef;
 use petgraph::Directed;
 use std::borrow::Cow;
+use std::marker::PhantomData;
 use std::ops::{Deref, MulAssign};
 use std::sync::Arc;
 
@@ -95,11 +96,18 @@ impl<F: HelmholtzEnergyFunctional> Components for DFT<F> {
 }
 
 impl<F: HelmholtzEnergyFunctional> Residual for DFT<F> {
+    type Properties<D> = PhantomData<D>;
+    type Inner = Self;
+
+    fn properties<D: DualNum<f64>>(&self, _: D) -> PhantomData<D> {
+        PhantomData
+    }
+
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
         self.0.compute_max_density(moles)
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self>>] {
         unreachable!()
     }
 
@@ -109,14 +117,16 @@ impl<F: HelmholtzEnergyFunctional> Residual for DFT<F> {
 
     fn evaluate_residual<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D
     where
-        dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
+        dyn HelmholtzEnergy<Self>: HelmholtzEnergyDual<PhantomData<D>, D>,
     {
         self.0
             .contributions()
             .iter()
-            .map(|c| (c as &dyn HelmholtzEnergy).helmholtz_energy(state))
+            .map(|c| (c as &dyn HelmholtzEnergy<Self>).helmholtz_energy(state, &PhantomData))
             .sum::<D>()
-            + self.ideal_chain_contribution().helmholtz_energy(state)
+            + self
+                .ideal_chain_contribution()
+                .helmholtz_energy(state, &PhantomData)
     }
 
     fn evaluate_residual_contributions<D: DualNum<f64> + Copy>(
@@ -124,7 +134,7 @@ impl<F: HelmholtzEnergyFunctional> Residual for DFT<F> {
         state: &StateHD<D>,
     ) -> Vec<(String, D)>
     where
-        dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
+        dyn HelmholtzEnergy<Self>: HelmholtzEnergyDual<PhantomData<D>, D>,
     {
         let mut res: Vec<(String, D)> = self
             .0
@@ -133,13 +143,14 @@ impl<F: HelmholtzEnergyFunctional> Residual for DFT<F> {
             .map(|c| {
                 (
                     c.to_string(),
-                    (c as &dyn HelmholtzEnergy).helmholtz_energy(state),
+                    (c as &dyn HelmholtzEnergy<Self>).helmholtz_energy(state, &PhantomData),
                 )
             })
             .collect();
         res.push((
             self.ideal_chain_contribution().to_string(),
-            self.ideal_chain_contribution().helmholtz_energy(state),
+            self.ideal_chain_contribution()
+                .helmholtz_energy(state, &PhantomData),
         ));
         res
     }
diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs
index 89bbf8499..02f6916b2 100644
--- a/feos-dft/src/functional_contribution.rs
+++ b/feos-dft/src/functional_contribution.rs
@@ -5,11 +5,18 @@ use ndarray::RemoveAxis;
 use num_dual::*;
 use num_traits::{One, Zero};
 use std::fmt::Display;
+use std::marker::PhantomData;
 
 macro_rules! impl_helmholtz_energy {
     ($number:ty) => {
-        impl HelmholtzEnergyDual<$number> for Box<dyn FunctionalContribution> {
-            fn helmholtz_energy(&self, state: &StateHD<$number>) -> $number {
+        impl HelmholtzEnergyDual<PhantomData<$number>, $number>
+            for Box<dyn FunctionalContribution>
+        {
+            fn helmholtz_energy(
+                &self,
+                state: &StateHD<$number>,
+                _: &PhantomData<$number>,
+            ) -> $number {
                 // calculate weight functions
                 let weight_functions = self.weight_functions(state.temperature);
 
diff --git a/feos-dft/src/ideal_chain_contribution.rs b/feos-dft/src/ideal_chain_contribution.rs
index d42e68a00..a38e6e544 100644
--- a/feos-dft/src/ideal_chain_contribution.rs
+++ b/feos-dft/src/ideal_chain_contribution.rs
@@ -3,6 +3,7 @@ use feos_core::{EosResult, HelmholtzEnergyDual, StateHD};
 use ndarray::*;
 use num_dual::DualNum;
 use std::fmt;
+use std::marker::PhantomData;
 
 #[derive(Clone)]
 pub struct IdealChainContribution {
@@ -19,8 +20,8 @@ impl IdealChainContribution {
     }
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for IdealChainContribution {
-    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PhantomData<D>, D> for IdealChainContribution {
+    fn helmholtz_energy(&self, state: &StateHD<D>, _: &PhantomData<D>) -> D {
         let segments = self.component_index.len();
         if self.component_index[segments - 1] + 1 != segments {
             return D::zero();
diff --git a/src/association/mod.rs b/src/association/mod.rs
index c3c74dff0..4c8468e90 100644
--- a/src/association/mod.rs
+++ b/src/association/mod.rs
@@ -332,7 +332,7 @@ impl<P: HardSphereProperties> Association<P> {
     }
 }
 
-impl<D: DualNum<f64> + Copy + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual<D>
+impl<D: DualNum<f64> + Copy + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDual<Self, D>
     for Association<P>
 {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
diff --git a/src/hard_sphere/mod.rs b/src/hard_sphere/mod.rs
index 62b494795..786b2c5f8 100644
--- a/src/hard_sphere/mod.rs
+++ b/src/hard_sphere/mod.rs
@@ -116,7 +116,9 @@ impl<P> HardSphere<P> {
     }
 }
 
-impl<D: DualNum<f64> + Copy, P: HardSphereProperties> HelmholtzEnergyDual<D> for HardSphere<P> {
+impl<D: DualNum<f64> + Copy, P: HardSphereProperties> HelmholtzEnergyDual<Self, D>
+    for HardSphere<P>
+{
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         let p = &self.parameters;
         let zeta = p.zeta(state.temperature, &state.partial_density, [0, 1, 2, 3]);
diff --git a/src/pcsaft/eos/dispersion.rs b/src/pcsaft/eos/dispersion.rs
index 416778685..75d67ce9a 100644
--- a/src/pcsaft/eos/dispersion.rs
+++ b/src/pcsaft/eos/dispersion.rs
@@ -1,4 +1,4 @@
-use super::PcSaftParameters;
+use super::{PcSaft, PcSaftParameters};
 use crate::hard_sphere::HardSphereProperties;
 use feos_core::{HelmholtzEnergyDual, StateHD};
 use num_dual::DualNum;
@@ -65,7 +65,7 @@ pub struct Dispersion {
     pub parameters: Arc<PcSaftParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for Dispersion {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PcSaft, D> for Dispersion {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         // auxiliary variables
         let n = self.parameters.m.len();
diff --git a/src/pcsaft/eos/hard_chain.rs b/src/pcsaft/eos/hard_chain.rs
index a046c893c..87adb1167 100644
--- a/src/pcsaft/eos/hard_chain.rs
+++ b/src/pcsaft/eos/hard_chain.rs
@@ -1,4 +1,4 @@
-use super::PcSaftParameters;
+use super::{PcSaft, PcSaftParameters};
 use crate::hard_sphere::HardSphereProperties;
 use feos_core::{HelmholtzEnergyDual, StateHD};
 use ndarray::Array;
@@ -10,7 +10,7 @@ pub struct HardChain {
     pub parameters: Arc<PcSaftParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for HardChain {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PcSaft, D> for HardChain {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         let p = &self.parameters;
         let d = self.parameters.hs_diameter(state.temperature);
diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs
index af2ad08d7..358477b28 100644
--- a/src/pcsaft/eos/mod.rs
+++ b/src/pcsaft/eos/mod.rs
@@ -9,6 +9,7 @@ use feos_core::{
 use ndarray::Array1;
 use std::f64::consts::{FRAC_PI_6, PI};
 use std::fmt;
+use std::marker::PhantomData;
 use std::sync::Arc;
 use typenum::P2;
 
@@ -44,7 +45,7 @@ impl Default for PcSaftOptions {
 pub struct PcSaft {
     parameters: Arc<PcSaftParameters>,
     options: PcSaftOptions,
-    contributions: Vec<Box<dyn HelmholtzEnergy>>,
+    contributions: Vec<Box<dyn HelmholtzEnergy<Self>>>,
 }
 
 impl PcSaft {
@@ -108,13 +109,20 @@ impl Components for PcSaft {
 }
 
 impl Residual for PcSaft {
+    type Properties<D> = PhantomData<D>;
+    type Inner = Self;
+
+    fn properties<D: num_dual::DualNum<f64>>(&self, temperature: D) -> PhantomData<D> {
+        todo!()
+    }
+
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
         self.options.max_eta * moles.sum()
             / (FRAC_PI_6 * &self.parameters.m * self.parameters.sigma.mapv(|v| v.powi(3)) * moles)
                 .sum()
     }
 
-    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy>] {
+    fn contributions(&self) -> &[Box<dyn HelmholtzEnergy<Self>>] {
         &self.contributions
     }
 
diff --git a/src/pcsaft/eos/polar.rs b/src/pcsaft/eos/polar.rs
index 9efe1c428..b6129afc0 100644
--- a/src/pcsaft/eos/polar.rs
+++ b/src/pcsaft/eos/polar.rs
@@ -1,4 +1,4 @@
-use super::PcSaftParameters;
+use super::{PcSaft, PcSaftParameters};
 use crate::hard_sphere::HardSphereProperties;
 use feos_core::{HelmholtzEnergyDual, StateHD};
 use ndarray::prelude::*;
@@ -165,7 +165,7 @@ pub struct Dipole {
     pub parameters: Arc<PcSaftParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for Dipole {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PcSaft, D> for Dipole {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Dipole);
         let p = &self.parameters;
@@ -242,7 +242,7 @@ pub struct Quadrupole {
     pub parameters: Arc<PcSaftParameters>,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for Quadrupole {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PcSaft, D> for Quadrupole {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         let m = MeanSegmentNumbers::new(&self.parameters, Multipole::Quadrupole);
         let p = &self.parameters;
@@ -329,7 +329,7 @@ pub struct DipoleQuadrupole {
     pub variant: DQVariants,
 }
 
-impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<D> for DipoleQuadrupole {
+impl<D: DualNum<f64> + Copy> HelmholtzEnergyDual<PcSaft, D> for DipoleQuadrupole {
     fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
         let p = &self.parameters;