diff --git a/src/cubic.rs b/src/cubic.rs
index 8062487..744e93d 100644
--- a/src/cubic.rs
+++ b/src/cubic.rs
@@ -17,6 +17,7 @@ use num_dual::DualNum;
 use quantity::si::{SIArray1, SIUnit};
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
+use std::fmt;
 use std::rc::Rc;
 
 const KB_A3: f64 = 13806490.0;
@@ -165,12 +166,51 @@ impl Parameter for PengRobinsonParameters {
     }
 }
 
+struct PengRobinsonContribution {
+    parameters: Rc<PengRobinsonParameters>,
+}
+
+impl<D: DualNum<f64>> HelmholtzEnergyDual<D> for PengRobinsonContribution {
+    fn helmholtz_energy(&self, state: &StateHD<D>) -> D {
+        // temperature dependent a parameter
+        let p = &self.parameters;
+        let x = &state.molefracs;
+        let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
+            .mapv(|x| x.powi(2))
+            * &p.a;
+
+        // Mixing rules
+        let mut ak_mix = D::zero();
+        for i in 0..ak.len() {
+            for j in 0..ak.len() {
+                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
+            }
+        }
+        let b = (x * &p.b).sum();
+
+        // Helmholtz energy
+        let n = state.moles.sum();
+        let v = state.volume;
+        n * ((v / (v - b * n)).ln()
+            - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
+                * ((v * (SQRT_2 - 1.0) + b * n) / (v * (SQRT_2 + 1.0) - b * n)).ln())
+    }
+}
+
+impl fmt::Display for PengRobinsonContribution {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "Peng Robinson")
+    }
+}
+
 /// A simple version of the Peng-Robinson equation of state.
 pub struct PengRobinson {
     /// Parameters
     parameters: Rc<PengRobinsonParameters>,
     /// Ideal gas contributions to the Helmholtz energy
     ideal_gas: Joback,
+    /// Non-ideal contributions to the Helmholtz energy
+    contributions: Vec<Box<dyn HelmholtzEnergy>>,
 }
 
 impl PengRobinson {
@@ -180,9 +220,14 @@ impl PengRobinson {
             || Joback::default(parameters.tc.len()),
             |j| Joback::new(j.clone()),
         );
+        let contributions: Vec<Box<dyn HelmholtzEnergy>> =
+            vec![Box::new(PengRobinsonContribution {
+                parameters: parameters.clone(),
+            })];
         Self {
             parameters,
             ideal_gas,
+            contributions,
         }
     }
 }
@@ -202,42 +247,7 @@ impl EquationOfState for PengRobinson {
     }
 
     fn residual(&self) -> &[Box<dyn HelmholtzEnergy>] {
-        unreachable!()
-    }
-
-    fn evaluate_residual<D: DualNum<f64>>(&self, state: &StateHD<D>) -> D {
-        // temperature dependent a parameter
-        let p = &self.parameters;
-        let x = &state.molefracs;
-        let ak = (&p.tc.mapv(|tc| (D::one() - (state.temperature / tc).sqrt())) * &p.kappa + 1.0)
-            .mapv(|x| x.powi(2))
-            * &p.a;
-
-        // Mixing rules
-        let mut ak_mix = D::zero();
-        for i in 0..ak.len() {
-            for j in 0..ak.len() {
-                ak_mix += (ak[i] * ak[j]).sqrt() * (x[i] * x[j] * (1.0 - p.k_ij[(i, j)]));
-            }
-        }
-        let b = (x * &p.b).sum();
-
-        // Helmholtz energy
-        let n = state.moles.sum();
-        let v = state.volume;
-        n * ((v / (v - b * n)).ln()
-            - ak_mix / (b * SQRT_2 * 2.0 * state.temperature)
-                * ((v * (SQRT_2 - 1.0) + b * n) / (v * (SQRT_2 + 1.0) - b * n)).ln())
-    }
-
-    fn evaluate_residual_contributions<D: DualNum<f64>>(
-        &self,
-        state: &StateHD<D>,
-    ) -> Vec<(String, D)>
-    where
-        dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
-    {
-        vec![("Peng-Robinson".into(), self.evaluate_residual(state))]
+        &self.contributions
     }
 
     fn ideal_gas(&self) -> &dyn IdealGasContribution {
diff --git a/src/equation_of_state.rs b/src/equation_of_state.rs
index 9d3e86f..6218747 100644
--- a/src/equation_of_state.rs
+++ b/src/equation_of_state.rs
@@ -165,9 +165,6 @@ pub trait EquationOfState {
     fn residual(&self) -> &[Box<dyn HelmholtzEnergy>];
 
     /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
-    ///
-    /// For simple equations of state (see e.g. `PengRobinson`) it might be
-    /// easier to overwrite this function instead of implementing `residual`.
     fn evaluate_residual<D: DualNum<f64>>(&self, state: &StateHD<D>) -> D
     where
         dyn HelmholtzEnergy: HelmholtzEnergyDual<D>,
@@ -179,10 +176,7 @@ pub trait EquationOfState {
     }
 
     /// Evaluate the reduced Helmholtz energy of each individual contribution
-    /// and return them together with a string representatino of the contribution.
-    ///
-    /// If `evaluate_residual` is implemented instead of `residual`, this function
-    /// also needs to be overwritten to avoid panics.
+    /// and return them together with a string representation of the contribution.
     fn evaluate_residual_contributions<D: DualNum<f64>>(
         &self,
         state: &StateHD<D>,