diff --git a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
index 09c4f656bebde10e631eea5e31b9977e7f2336e3..16cae68b460eab5bd0cd373ec80bd636a41d624d 100644
--- a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
+++ b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
@@ -34,8 +34,6 @@
 #include <dumux/porousmediumflow/volumevariables.hh>
 
 #include <dumux/material/fluidstates/nonequilibrium.hh>
-#include <dumux/material/constraintsolvers/fluidsystemcomputefromreferencephase.hh>
-#include <dumux/material/constraintsolvers/fluidsystemconstraintsolver.hh>
 #include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
 
 namespace Dumux
@@ -304,8 +302,7 @@ class NonEquilibriumVolumeVariablesImplementation<TypeTag, true/*enableChemicalN
     using AwnSurfaceParams = typename  AwnSurface::Params;
 
     using DimLessNum = DimensionlessNumbers<Scalar>;
-    //TODO: write a good constraint solver which does that right
-    using ConstraintSolver = FluidSystemConstraintSolver<Scalar, FluidSystem>;
+    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
 public:
     /*!
      * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
@@ -533,7 +530,7 @@ class NonEquilibriumVolumeVariablesImplementation<TypeTag, true/*enableChemicalN
     using AnsSurface = typename GET_PROP_TYPE(TypeTag, AnsSurface);
     using AnsSurfaceParams = typename  AnsSurface::Params;
 
-    using ConstraintSolver = Dumux::FluidSystemConstraintSolver<Scalar, FluidSystem>;
+    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
 
 
 public:
diff --git a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
index 6205c76f098ebd45f3a0b3459a598955e07eaeae..67ed1aed79f0f343ad3797c049e2351598cf4b01 100644
--- a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
@@ -132,6 +132,7 @@ class EvaporationAtmosphereProblem: public PorousMediumFlowProblem<TypeTag>
     enum { numEnergyEqSolid = GET_PROP_VALUE(TypeTag, NumEnergyEqSolid)};
 
     static constexpr bool enableChemicalNonEquilibrium = GET_PROP_VALUE(TypeTag, EnableChemicalNonEquilibrium);
+    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
 
     // formulations
     enum {
@@ -268,7 +269,10 @@ public:
         fluidState.setTemperature(wPhaseIdx, TInitial_ ); // this value is a good one, TInject does not work
 
         // This solves the system of equations defining x=x(p,T)
-        FluidSystem::calculateEquilibriumMoleFractions(fluidState, dummyCache) ;
+        ConstraintSolver::solve(fluidState,
+                                dummyCache,
+                                /*setViscosity=*/false,
+                                /*setEnthalpy=*/false) ;
 
         // Now let's make the air phase less than fully saturated with water
         fluidState.setMoleFraction(nPhaseIdx, wCompIdx, fluidState.moleFraction(nPhaseIdx, wCompIdx)*percentOfEquil_ ) ;
@@ -417,7 +421,10 @@ private:
 
          // This solves the system of equations defining x=x(p,T)
         ParameterCache dummyCache;
-        FluidSystem::calculateEquilibriumMoleFractions(equilibriumFluidState, dummyCache) ;
+        ConstraintSolver::solve(equilibriumFluidState,
+                                dummyCache,
+                                /*setViscosity=*/false,
+                                /*setEnthalpy=*/false) ;
 
         FluidState dryFluidState(equilibriumFluidState);
         // Now let's make the air phase less than fully saturated with vapor