diff --git a/feos-core/CHANGELOG.md b/feos-core/CHANGELOG.md
index 0b3e7d654..4f3bd109b 100644
--- a/feos-core/CHANGELOG.md
+++ b/feos-core/CHANGELOG.md
@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 ### Added
 - Added new functions `isenthalpic_compressibility`, `thermal_expansivity` and `grueneisen_parameter` to `State`. [#154](https://github.com/feos-org/feos/pull/154)
+- Readded `PhaseEquilibrium::new_npt` to the public interface in Rust and Python.  [#164](https://github.com/feos-org/feos/pull/164)
 
 ### Changed
 - Changed constructors of `Parameter` trait to return `Result`s. [#161](https://github.com/feos-org/feos/pull/161)
diff --git a/feos-core/src/phase_equilibria/mod.rs b/feos-core/src/phase_equilibria/mod.rs
index 7c020e352..7574ce458 100644
--- a/feos-core/src/phase_equilibria/mod.rs
+++ b/feos-core/src/phase_equilibria/mod.rs
@@ -195,7 +195,13 @@ impl<E: EquationOfState> PhaseEquilibrium<E, 2> {
         Self([vapor, liquid])
     }
 
-    pub(super) fn new_npt(
+    /// Creates a new PhaseEquilibrium that contains two states at the
+    /// specified temperature, pressure and moles.
+    ///
+    /// The constructor can be used in custom phase equilibrium solvers or,
+    /// e.g., to generate initial guesses for an actual VLE solver.
+    /// In general, the two states generated are NOT in an equilibrium.
+    pub fn new_npt(
         eos: &Arc<E>,
         temperature: SINumber,
         pressure: SINumber,
diff --git a/feos-core/src/python/phase_equilibria.rs b/feos-core/src/python/phase_equilibria.rs
index 72951802d..4e4950ded 100644
--- a/feos-core/src/python/phase_equilibria.rs
+++ b/feos-core/src/python/phase_equilibria.rs
@@ -227,6 +227,46 @@ macro_rules! impl_phase_equilibrium {
                 )?))
             }
 
+            /// Creates a new PhaseEquilibrium that contains two states at the
+            /// specified temperature, pressure and moles.
+            ///
+            /// The constructor can be used in custom phase equilibrium solvers or,
+            /// e.g., to generate initial guesses for an actual VLE solver.
+            /// In general, the two states generated are NOT in an equilibrium.
+            ///
+            /// Parameters
+            /// ----------
+            /// eos : EquationOfState
+            ///     The equation of state.
+            /// temperature : SINumber
+            ///     The system temperature.
+            /// pressure : SINumber
+            ///     The system pressure.
+            /// vapor_moles : SIArray1
+            ///     Amount of substance of the vapor phase.
+            /// liquid_moles : SIArray1
+            ///     Amount of substance of the liquid phase.
+            ///
+            /// Returns
+            /// -------
+            /// PhaseEquilibrium
+            #[staticmethod]
+            pub fn new_npt(
+                eos: $py_eos,
+                temperature: PySINumber,
+                pressure: PySINumber,
+                vapor_moles: &PySIArray1,
+                liquid_moles: &PySIArray1
+            ) -> PyResult<Self> {
+                Ok(Self(PhaseEquilibrium::new_npt(
+                    &eos.0,
+                    temperature.into(),
+                    pressure.into(),
+                    vapor_moles,
+                    liquid_moles
+                )?))
+            }
+
             #[getter]
             fn get_vapor(&self) -> PyState {
                 PyState(self.0.vapor().clone())