diff --git a/Cargo.toml b/Cargo.toml
index 58f415ba0..046042203 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -52,6 +52,10 @@ criterion = "0.4"
 name = "state_properties"
 harness = false
 
+[[bench]]
+name = "state_creation"
+harness = false
+
 [[bench]]
 name = "dual_numbers"
 harness = false
diff --git a/benches/README.md b/benches/README.md
index bc9ab59e5..0689e2cbe 100644
--- a/benches/README.md
+++ b/benches/README.md
@@ -13,4 +13,5 @@ cargo bench --profile=release-lto --features=pcsaft --bench=dual_numbers
 |Name|Description|Cargo features|
 |--|--|--|
 |`dual_numbers`|Helmholtz energy function evaluated using `StateHD` with different dual number types.|`pcsaft`|
-|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`|
\ No newline at end of file
+|`state_properties`|Properties of `State`. Including state creation using the natural variables of the Helmholtz energy (no density iteration).|`pcsaft`|
+|`state_creation`|Different constructors of `State` and `PhaseEquilibrium` including critical point calculations. For pure substances and mixtures.|`pcsaft`|
\ No newline at end of file
diff --git a/benches/dual_numbers.rs b/benches/dual_numbers.rs
index bf9f37f38..b38607b4b 100644
--- a/benches/dual_numbers.rs
+++ b/benches/dual_numbers.rs
@@ -48,8 +48,11 @@ fn bench_dual_numbers<E: EquationOfState>(
     group.bench_function("a_dual", |b| {
         b.iter(|| a_res((&state.eos, &state.derive1(Derivative::DV))))
     });
+    group.bench_function("a_dual2", |b| {
+        b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV))))
+    });
     group.bench_function("a_hyperdual", |b| {
-        b.iter(|| a_res((&state.eos, &state.derive2(Derivative::DV, Derivative::DV))))
+        b.iter(|| a_res((&state.eos, &state.derive2_mixed(Derivative::DV, Derivative::DV))))
     });
     group.bench_function("a_dual3", |b| {
         b.iter(|| a_res((&state.eos, &state.derive3(Derivative::DV))))
@@ -66,7 +69,7 @@ fn pcsaft(c: &mut Criterion) {
         IdentifierOption::Name,
     )
     .unwrap();
-    bench_dual_numbers(c, "methane", state_pcsaft(parameters));
+    bench_dual_numbers(c, "dual_numbers_pcsaft_methane", state_pcsaft(parameters));
 
     // water (4C, polar)
     let parameters = PcSaftParameters::from_json(
@@ -76,7 +79,7 @@ fn pcsaft(c: &mut Criterion) {
         IdentifierOption::Name,
     )
     .unwrap();
-    bench_dual_numbers(c, "water_4c_polar", state_pcsaft(parameters));
+    bench_dual_numbers(c, "dual_numbers_pcsaft_water_4c_polar", state_pcsaft(parameters));
 
     // methane, ethane, propane
     let parameters = PcSaftParameters::from_json(
@@ -86,7 +89,7 @@ fn pcsaft(c: &mut Criterion) {
         IdentifierOption::Name,
     )
     .unwrap();
-    bench_dual_numbers(c, "methane_ethane_propane", state_pcsaft(parameters));
+    bench_dual_numbers(c, "dual_numbers_pcsaft_methane_ethane_propane", state_pcsaft(parameters));
 }
 
 /// Benchmark for the PC-SAFT equation of state.
@@ -116,7 +119,7 @@ fn methane_co2_pcsaft(c: &mut Criterion) {
     let x = arr1(&[0.15, 0.85]);
     let moles = &x * 10.0 * MOL;
     let state = State::new_nvt(&eos, temperature, volume, &moles).unwrap();
-    bench_dual_numbers(c, "methane_co2", state);
+    bench_dual_numbers(c, "dual_numbers_pcsaft_methane_co2", state);
 }
 
 criterion_group!(bench, pcsaft, methane_co2_pcsaft);
diff --git a/benches/state_creation.rs b/benches/state_creation.rs
new file mode 100644
index 000000000..55f1f3a11
--- /dev/null
+++ b/benches/state_creation.rs
@@ -0,0 +1,171 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use feos::pcsaft::{PcSaft, PcSaftParameters};
+use feos_core::{
+    parameter::{IdentifierOption, Parameter},
+    Contributions, DensityInitialization, EquationOfState, PhaseEquilibrium, State,
+};
+use ndarray::{Array, Array1};
+use quantity::si::*;
+use std::sync::Arc;
+
+/// Evaluate NPT constructor
+fn npt<E: EquationOfState>(
+    (eos, t, p, n, rho0): (
+        &Arc<E>,
+        SINumber,
+        SINumber,
+        &SIArray1,
+        DensityInitialization<SIUnit>,
+    ),
+) {
+    State::new_npt(eos, t, p, n, rho0).unwrap();
+}
+
+/// Evaluate critical point constructor
+fn critical_point<E: EquationOfState>((eos, n): (&Arc<E>, Option<&SIArray1>)) {
+    State::critical_point(eos, n, None, Default::default()).unwrap();
+}
+
+/// VLE for pure substance for given temperature or pressure
+fn pure<E: EquationOfState>((eos, t_or_p): (&Arc<E>, SINumber)) {
+    PhaseEquilibrium::pure(eos, t_or_p, None, Default::default()).unwrap();
+}
+
+/// Evaluate temperature, pressure flash.
+fn tp_flash<E: EquationOfState>((eos, t, p, feed): (&Arc<E>, SINumber, SINumber, &SIArray1)) {
+    PhaseEquilibrium::tp_flash(eos, t, p, feed, None, Default::default(), None).unwrap();
+}
+
+fn bubble_point<E: EquationOfState>((eos, t, x): (&Arc<E>, SINumber, &Array1<f64>)) {
+    PhaseEquilibrium::bubble_point(
+        eos,
+        t,
+        x,
+        None,
+        None,
+        (Default::default(), Default::default()),
+    )
+    .unwrap();
+}
+
+fn dew_point<E: EquationOfState>((eos, t, y): (&Arc<E>, SINumber, &Array1<f64>)) {
+    PhaseEquilibrium::dew_point(
+        eos,
+        t,
+        y,
+        None,
+        None,
+        (Default::default(), Default::default()),
+    )
+    .unwrap();
+}
+
+fn bench_states<E: EquationOfState>(c: &mut Criterion, group_name: &str, eos: &Arc<E>) {
+    let ncomponents = eos.components();
+    let x = Array::from_elem(ncomponents, 1.0 / ncomponents as f64);
+    let n = &x * 100.0 * MOL;
+    let crit = State::critical_point(&eos, Some(&n), None, Default::default()).unwrap();
+    let vle = if ncomponents == 1 {
+        PhaseEquilibrium::pure(&eos, crit.temperature * 0.95, None, Default::default()).unwrap()
+    } else {
+        PhaseEquilibrium::tp_flash(
+            &eos,
+            crit.temperature,
+            crit.pressure(Contributions::Total) * 0.95,
+            &crit.moles,
+            None,
+            Default::default(),
+            None,
+        )
+        .unwrap()
+    };
+
+    let mut group = c.benchmark_group(group_name);
+    group.bench_function("new_npt_liquid", |b| {
+        b.iter(|| {
+            npt((
+                &eos,
+                vle.liquid().temperature,
+                vle.liquid().pressure(Contributions::Total) * 1.01,
+                &n,
+                DensityInitialization::Liquid,
+            ))
+        })
+    });
+    group.bench_function("new_npt_vapor", |b| {
+        b.iter(|| {
+            npt((
+                &eos,
+                vle.vapor().temperature,
+                vle.vapor().pressure(Contributions::Total) * 0.99,
+                &n,
+                DensityInitialization::Vapor,
+            ))
+        })
+    });
+    group.bench_function("critical_point", |b| {
+        b.iter(|| critical_point((&eos, Some(&n))))
+    });
+    if ncomponents != 1 {
+        group.bench_function("tp_flash", |b| {
+            b.iter(|| {
+                tp_flash((
+                    &eos,
+                    crit.temperature,
+                    crit.pressure(Contributions::Total) * 0.99,
+                    &n,
+                ))
+            })
+        });
+
+        group.bench_function("bubble_point", |b| {
+            b.iter(|| bubble_point((&eos, vle.liquid().temperature, &vle.liquid().molefracs)))
+        });
+
+        group.bench_function("dew_point", |b| {
+            b.iter(|| dew_point((&eos, vle.vapor().temperature, &vle.vapor().molefracs)))
+        });
+    } else {
+        group.bench_function("pure_t", |b| {
+            b.iter(|| pure((&eos, vle.vapor().temperature)))
+        });
+        group.bench_function("pure_p", |b| {
+            b.iter(|| pure((&eos, vle.vapor().pressure(Contributions::Total))))
+        });
+    }
+}
+
+fn pcsaft(c: &mut Criterion) {
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+    bench_states(c, "state_creation_pcsaft_methane", &eos);
+
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane", "ethane", "propane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+    bench_states(c, "state_creation_pcsaft_methane_ethane", &eos);
+
+    let parameters = PcSaftParameters::from_json(
+        vec!["methane", "ethane", "propane"],
+        "./parameters/pcsaft/gross2001.json",
+        None,
+        IdentifierOption::Name,
+    )
+    .unwrap();
+    let eos = Arc::new(PcSaft::new(Arc::new(parameters)));
+    bench_states(c, "state_creation_pcsaft_methane_ethane_propane", &eos);
+}
+
+criterion_group!(bench, pcsaft);
+criterion_main!(bench);
diff --git a/benches/state_properties.rs b/benches/state_properties.rs
index e9b386c2b..db4373568 100644
--- a/benches/state_properties.rs
+++ b/benches/state_properties.rs
@@ -50,7 +50,7 @@ fn properties_pcsaft(c: &mut Criterion) {
     let x = arr1(&[1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]);
     let m = &x * 100.0 * MOL;
 
-    let mut group = c.benchmark_group("methane_ethane_propane");
+    let mut group = c.benchmark_group("state_properties_pcsaft_methane_ethane_propane");
     group.bench_function("a", |b| {
         b.iter(|| property((&eos, S::helmholtz_energy, t, v, &m, Contributions::Total)))
     });
diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md
index 475f68937..d428a2d69 100644
--- a/feos-core/CHANGELOG.md
+++ b/feos-core/CHANGELOG.md
@@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - Added `Sync` and `Send` as supertraits to `EquationOfState`. [#57](https://github.com/feos-org/feos/pull/57)
+- Added `Dual2_64` dual number to `HelmholtzEnergy` trait to facilitate efficient non-mixed second order derivatives. [#94](https://github.com/feos-org/feos/pull/94)
+- Renamed `PartialDerivative::Second` to `PartialDerivative::SecondMixed`. [#94](https://github.com/feos-org/feos/pull/94)
+- Added `PartialDerivative::Second` to enable non-mixed second order partial derivatives. [#94](https://github.com/feos-org/feos/pull/94)
+- Changed `dp_dv_` and `ds_dt_` to use `Dual2_64` instead of `HyperDual64`. [#94](https://github.com/feos-org/feos/pull/94)
+- Added `get_or_insert_with_d2_64` to `Cache`. [#94](https://github.com/feos-org/feos/pull/94)
 
 ## [0.3.1] - 2022-08-25
 ### Added
diff --git a/feos-core/src/equation_of_state.rs b/feos-core/src/equation_of_state.rs
index e2efca5c7..406da7d0a 100644
--- a/feos-core/src/equation_of_state.rs
+++ b/feos-core/src/equation_of_state.rs
@@ -2,7 +2,7 @@ use crate::errors::{EosError, EosResult};
 use crate::state::StateHD;
 use crate::EosUnit;
 use ndarray::prelude::*;
-use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64};
+use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualNum, DualVec64, HyperDual, HyperDual64, Dual2_64};
 use num_traits::{One, Zero};
 use quantity::{QuantityArray1, QuantityScalar};
 use std::fmt;
@@ -28,6 +28,7 @@ pub trait HelmholtzEnergy:
     + HelmholtzEnergyDual<Dual64>
     + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
     + HelmholtzEnergyDual<HyperDual64>
+    + HelmholtzEnergyDual<Dual2_64>
     + HelmholtzEnergyDual<Dual3_64>
     + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
     + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
@@ -46,6 +47,7 @@ impl<T> HelmholtzEnergy for T where
         + HelmholtzEnergyDual<Dual64>
         + HelmholtzEnergyDual<Dual<DualVec64<3>, f64>>
         + HelmholtzEnergyDual<HyperDual64>
+        + HelmholtzEnergyDual<Dual2_64>
         + HelmholtzEnergyDual<Dual3_64>
         + HelmholtzEnergyDual<HyperDual<Dual64, f64>>
         + HelmholtzEnergyDual<HyperDual<DualVec64<2>, f64>>
@@ -99,6 +101,7 @@ pub trait IdealGasContribution:
     + IdealGasContributionDual<Dual64>
     + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
     + IdealGasContributionDual<HyperDual64>
+    + IdealGasContributionDual<Dual2_64>
     + IdealGasContributionDual<Dual3_64>
     + IdealGasContributionDual<HyperDual<Dual64, f64>>
     + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
@@ -115,6 +118,7 @@ impl<T> IdealGasContribution for T where
         + IdealGasContributionDual<Dual64>
         + IdealGasContributionDual<Dual<DualVec64<3>, f64>>
         + IdealGasContributionDual<HyperDual64>
+        + IdealGasContributionDual<Dual2_64>
         + IdealGasContributionDual<Dual3_64>
         + IdealGasContributionDual<HyperDual<Dual64, f64>>
         + IdealGasContributionDual<HyperDual<DualVec64<2>, f64>>
diff --git a/feos-core/src/python/statehd.rs b/feos-core/src/python/statehd.rs
index ab4dc2ae6..c1f70f99e 100644
--- a/feos-core/src/python/statehd.rs
+++ b/feos-core/src/python/statehd.rs
@@ -1,7 +1,9 @@
 use crate::StateHD;
 use ndarray::Array1;
-use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64};
-use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
+use num_dual::python::{
+    PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64,
+};
+use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
 use pyo3::prelude::*;
 use std::convert::From;
 
@@ -75,6 +77,16 @@ impl From<StateHD<Dual<DualVec64<3>, f64>>> for PyStateDDV3 {
     }
 }
 
+#[pyclass(name = "StateD2")]
+#[derive(Clone)]
+pub struct PyStateD2(StateHD<Dual2_64>);
+
+impl From<StateHD<Dual2_64>> for PyStateD2 {
+    fn from(s: StateHD<Dual2_64>) -> Self {
+        Self(s)
+    }
+}
+
 #[pyclass(name = "StateD3")]
 #[derive(Clone)]
 pub struct PyStateD3(StateHD<Dual3_64>);
@@ -167,6 +179,7 @@ impl_state_hd!(
     HyperDual<Dual64, f64>
 );
 impl_state_hd!(PyStateD, PyDual64, StateHD<Dual64>, Dual64);
+impl_state_hd!(PyStateD2, PyDual2_64, StateHD<Dual2_64>, Dual2_64);
 impl_state_hd!(PyStateD3, PyDual3_64, StateHD<Dual3_64>, Dual3_64);
 impl_state_hd!(
     PyStateD3D,
diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs
index e218ee60e..939e73fd8 100644
--- a/feos-core/src/python/user_defined.rs
+++ b/feos-core/src/python/user_defined.rs
@@ -1,8 +1,10 @@
 use crate::python::statehd::*;
 use crate::*;
 use ndarray::prelude::*;
-use num_dual::python::{PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64};
-use num_dual::{Dual, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
+use num_dual::python::{
+    PyDual2_64, PyDual3Dual64, PyDual3_64, PyDual64, PyHyperDual64, PyHyperDualDual64,
+};
+use num_dual::{Dual, Dual2_64, Dual3, Dual3_64, Dual64, DualVec64, HyperDual, HyperDual64};
 use numpy::convert::IntoPyArray;
 use pyo3::prelude::*;
 use quantity::python::PySIArray1;
@@ -215,6 +217,7 @@ impl_helmholtz_energy!(
     PyHyperDualDualVec64_3,
     HyperDual<DualVec64<3>, f64>
 );
+impl_helmholtz_energy!(PyStateD2, PyDual2_64, Dual2_64);
 impl_helmholtz_energy!(PyStateD3, PyDual3_64, Dual3_64);
 impl_helmholtz_energy!(PyStateD3D, PyDual3Dual64, Dual3<Dual64, f64>);
 impl_helmholtz_energy!(PyStateD3DV2, PyDual3DualVec64_2, Dual3<DualVec64<2>, f64>);
diff --git a/feos-core/src/state/cache.rs b/feos-core/src/state/cache.rs
index 3723ed7e6..daf2a3e4c 100644
--- a/feos-core/src/state/cache.rs
+++ b/feos-core/src/state/cache.rs
@@ -50,6 +50,26 @@ impl Cache {
         }
     }
 
+    pub fn get_or_insert_with_d2_64<F: FnOnce() -> Dual2_64>(
+        &mut self,
+        derivative: Derivative,
+        f: F,
+    ) -> f64 {
+        if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed(derivative, derivative)) {
+            self.hit += 1;
+            value
+        } else {
+            self.miss += 1;
+            let value = f();
+            self.map.insert(PartialDerivative::Zeroth, value.re);
+            self.map
+                .insert(PartialDerivative::First(derivative), value.v1[0]);
+            self.map
+                .insert(PartialDerivative::SecondMixed(derivative, derivative), value.v2[0]);
+            value.v2[0]
+        }
+    }
+
     pub fn get_or_insert_with_hd64<F: FnOnce() -> HyperDual64>(
         &mut self,
         derivative1: Derivative,
@@ -58,7 +78,7 @@ impl Cache {
     ) -> f64 {
         let d1 = min(derivative1, derivative2);
         let d2 = max(derivative1, derivative2);
-        if let Some(&value) = self.map.get(&PartialDerivative::Second(d1, d2)) {
+        if let Some(&value) = self.map.get(&PartialDerivative::SecondMixed(d1, d2)) {
             self.hit += 1;
             value
         } else {
@@ -69,8 +89,10 @@ impl Cache {
                 .insert(PartialDerivative::First(derivative1), value.eps1[0]);
             self.map
                 .insert(PartialDerivative::First(derivative2), value.eps2[0]);
-            self.map
-                .insert(PartialDerivative::Second(d1, d2), value.eps1eps2[(0, 0)]);
+            self.map.insert(
+                PartialDerivative::SecondMixed(d1, d2),
+                value.eps1eps2[(0, 0)],
+            );
             value.eps1eps2[(0, 0)]
         }
     }
@@ -89,8 +111,10 @@ impl Cache {
             self.map.insert(PartialDerivative::Zeroth, value.re);
             self.map
                 .insert(PartialDerivative::First(derivative), value.v1);
-            self.map
-                .insert(PartialDerivative::Second(derivative, derivative), value.v2);
+            self.map.insert(
+                PartialDerivative::SecondMixed(derivative, derivative),
+                value.v2,
+            );
             self.map
                 .insert(PartialDerivative::Third(derivative), value.v3);
             value.v3
diff --git a/feos-core/src/state/mod.rs b/feos-core/src/state/mod.rs
index 13ff9792a..db6778e5a 100644
--- a/feos-core/src/state/mod.rs
+++ b/feos-core/src/state/mod.rs
@@ -211,7 +211,8 @@ impl Derivative {
 pub(crate) enum PartialDerivative {
     Zeroth,
     First(Derivative),
-    Second(Derivative, Derivative),
+    Second(Derivative),
+    SecondMixed(Derivative, Derivative),
     Third(Derivative),
 }
 
@@ -685,15 +686,28 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         let mut v = Dual64::from(self.reduced_volume);
         let mut n = self.reduced_moles.mapv(Dual64::from);
         match derivative {
-            Derivative::DT => t.eps[0] = 1.0,
-            Derivative::DV => v.eps[0] = 1.0,
-            Derivative::DN(i) => n[i].eps[0] = 1.0,
+            Derivative::DT => t = t.derive(),
+            Derivative::DV => v = v.derive(),
+            Derivative::DN(i) => n[i] = n[i].derive(),
+        }
+        StateHD::new(t, v, n)
+    }
+
+    /// Creates a [StateHD] taking the first and second (partial) derivatives.
+    pub fn derive2(&self, derivative: Derivative) -> StateHD<Dual2_64> {
+        let mut t = Dual2_64::from(self.reduced_temperature);
+        let mut v = Dual2_64::from(self.reduced_volume);
+        let mut n = self.reduced_moles.mapv(Dual2_64::from);
+        match derivative {
+            Derivative::DT => t = t.derive(),
+            Derivative::DV => v = v.derive(),
+            Derivative::DN(i) => n[i] = n[i].derive(),
         }
         StateHD::new(t, v, n)
     }
 
     /// Creates a [StateHD] taking the first and second (partial) derivatives.
-    pub fn derive2(
+    pub fn derive2_mixed(
         &self,
         derivative1: Derivative,
         derivative2: Derivative,
@@ -702,14 +716,14 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         let mut v = HyperDual64::from(self.reduced_volume);
         let mut n = self.reduced_moles.mapv(HyperDual64::from);
         match derivative1 {
-            Derivative::DT => t.eps1[0] = 1.0,
-            Derivative::DV => v.eps1[0] = 1.0,
-            Derivative::DN(i) => n[i].eps1[0] = 1.0,
+            Derivative::DT => t = t.derive1(),
+            Derivative::DV => v = v.derive1(),
+            Derivative::DN(i) => n[i] = n[i].derive1(),
         }
         match derivative2 {
-            Derivative::DT => t.eps2[0] = 1.0,
-            Derivative::DV => v.eps2[0] = 1.0,
-            Derivative::DN(i) => n[i].eps2[0] = 1.0,
+            Derivative::DT => t = t.derive2(),
+            Derivative::DV => v = v.derive2(),
+            Derivative::DN(i) => n[i] = n[i].derive2(),
         }
         StateHD::new(t, v, n)
     }
diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs
index e84be6645..9a580926b 100644
--- a/feos-core/src/state/properties.rs
+++ b/feos-core/src/state/properties.rs
@@ -50,8 +50,13 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
                     -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).eps[0]
                         * (U::reference_energy() / v.reference())
                 }
-                PartialDerivative::Second(v1, v2) => {
-                    let new_state = self.derive2(v1, v2);
+                PartialDerivative::Second(v) => {
+                    let new_state = self.derive2(v);
+                    -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln()).v2[0]
+                        * (U::reference_energy() / (v.reference() * v.reference()))
+                }
+                PartialDerivative::SecondMixed(v1, v2) => {
+                    let new_state = self.derive2_mixed(v1, v2);
                     -(new_state.moles.sum() * new_state.temperature * new_state.volume.ln())
                         .eps1eps2[(0, 0)]
                         * (U::reference_energy() / (v1.reference() * v2.reference()))
@@ -82,8 +87,15 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
                     cache.get_or_insert_with_d64(v, &computation) * U::reference_energy()
                         / v.reference()
                 }
-                PartialDerivative::Second(v1, v2) => {
-                    let new_state = self.derive2(v1, v2);
+                PartialDerivative::Second(v) => {
+                    let new_state = self.derive2(v);
+                    let computation =
+                        || self.eos.evaluate_residual(&new_state) * new_state.temperature;
+                    cache.get_or_insert_with_d2_64(v, &computation) * U::reference_energy()
+                        / (v.reference() * v.reference())
+                }
+                PartialDerivative::SecondMixed(v1, v2) => {
+                    let new_state = self.derive2_mixed(v1, v2);
                     let computation =
                         || self.eos.evaluate_residual(&new_state) * new_state.temperature;
                     cache.get_or_insert_with_hd64(v1, v2, &computation) * U::reference_energy()
@@ -114,8 +126,14 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
                         * U::reference_energy()
                         / v.reference()
                 }
-                PartialDerivative::Second(v1, v2) => {
-                    let new_state = self.derive2(v1, v2);
+                PartialDerivative::Second(v) => {
+                    let new_state = self.derive2(v);
+                    (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).v2[0]
+                        * U::reference_energy()
+                        / (v.reference() * v.reference())
+                }
+                PartialDerivative::SecondMixed(v1, v2) => {
+                    let new_state = self.derive2_mixed(v1, v2);
                     (self.eos.ideal_gas().evaluate(&new_state) * new_state.temperature).eps1eps2
                         [(0, 0)]
                         * U::reference_energy()
@@ -190,16 +208,16 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
     }
 
     fn dp_dv_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
-        -self.get_or_compute_derivative(PartialDerivative::Second(DV, DV), evaluate)
+        -self.get_or_compute_derivative(PartialDerivative::Second(DV), evaluate)
     }
 
     fn dp_dt_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
-        -self.get_or_compute_derivative(PartialDerivative::Second(DV, DT), evaluate)
+        -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DT), evaluate)
     }
 
     fn dp_dni_(&self, evaluate: Evaluate) -> QuantityArray1<U> {
         QuantityArray::from_shape_fn(self.eos.components(), |i| {
-            -self.get_or_compute_derivative(PartialDerivative::Second(DV, DN(i)), evaluate)
+            -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DN(i)), evaluate)
         })
     }
 
@@ -209,19 +227,19 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
 
     fn dmu_dt_(&self, evaluate: Evaluate) -> QuantityArray1<U> {
         QuantityArray::from_shape_fn(self.eos.components(), |i| {
-            self.get_or_compute_derivative(PartialDerivative::Second(DT, DN(i)), evaluate)
+            self.get_or_compute_derivative(PartialDerivative::SecondMixed(DT, DN(i)), evaluate)
         })
     }
 
     fn dmu_dni_(&self, evaluate: Evaluate) -> QuantityArray2<U> {
         let n = self.eos.components();
         QuantityArray::from_shape_fn((n, n), |(i, j)| {
-            self.get_or_compute_derivative(PartialDerivative::Second(DN(i), DN(j)), evaluate)
+            self.get_or_compute_derivative(PartialDerivative::SecondMixed(DN(i), DN(j)), evaluate)
         })
     }
 
     fn ds_dt_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
-        -self.get_or_compute_derivative(PartialDerivative::Second(DT, DT), evaluate)
+        -self.get_or_compute_derivative(PartialDerivative::Second(DT), evaluate)
     }
 
     fn d2s_dt2_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs
index cc5dc099f..b458d386b 100644
--- a/feos-dft/src/functional_contribution.rs
+++ b/feos-dft/src/functional_contribution.rs
@@ -36,6 +36,7 @@ impl_helmholtz_energy!(f64);
 impl_helmholtz_energy!(Dual64);
 impl_helmholtz_energy!(Dual<DualVec64<3>, f64>);
 impl_helmholtz_energy!(HyperDual64);
+impl_helmholtz_energy!(Dual2_64);
 impl_helmholtz_energy!(Dual3_64);
 impl_helmholtz_energy!(HyperDual<Dual64, f64>);
 impl_helmholtz_energy!(HyperDual<DualVec64<2>, f64>);
@@ -76,6 +77,7 @@ pub trait FunctionalContribution:
     + FunctionalContributionDual<Dual64>
     + FunctionalContributionDual<Dual<DualVec64<3>, f64>>
     + FunctionalContributionDual<HyperDual64>
+    + FunctionalContributionDual<Dual2_64>
     + FunctionalContributionDual<Dual3_64>
     + FunctionalContributionDual<HyperDual<Dual64, f64>>
     + FunctionalContributionDual<HyperDual<DualVec64<2>, f64>>
@@ -161,6 +163,7 @@ impl<T> FunctionalContribution for T where
         + FunctionalContributionDual<Dual64>
         + FunctionalContributionDual<Dual<DualVec64<3>, f64>>
         + FunctionalContributionDual<HyperDual64>
+        + FunctionalContributionDual<Dual2_64>
         + FunctionalContributionDual<Dual3_64>
         + FunctionalContributionDual<HyperDual<Dual64, f64>>
         + FunctionalContributionDual<HyperDual<DualVec64<2>, f64>>