diff --git a/dumux/material/constraintsolvers/fluidsystemconstraintsolver.hh b/dumux/material/constraintsolvers/fluidsystemconstraintsolver.hh
index f69967655601e21d4f60184250dfcab126c7116e..29e5375275fcb9ff116b5db6b916f8551b25efdb 100644
--- a/dumux/material/constraintsolvers/fluidsystemconstraintsolver.hh
+++ b/dumux/material/constraintsolvers/fluidsystemconstraintsolver.hh
@@ -32,6 +32,7 @@
 
 #include <dumux/common/exceptions.hh>
 #include <dumux/common/valgrind.hh>
+#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
 
 namespace Dumux {
 /*!
@@ -58,49 +59,7 @@ namespace Dumux {
  * - if the setEnthalpy parameter is true, also specific enthalpies of *all* phases
  */
 template <class Scalar, class FluidSystem>
-class FluidSystemConstraintSolver
-{
-    static constexpr int numPhases = FluidSystem::numPhases;
-
-public:
-    /*!
-     * \brief @copybrief Dumux::FluidSystemConstraintSolver
-     *
-     * \param fluidState A container with the current (physical) state of the fluid
-     * \param paramCache A container for iterative calculation of fluid composition
-     * \param setViscosity Should the viscosity be set in the fluidstate?
-     * \param setEnthalpy Should the enthalpy be set in the fluidstate?
-     */
-    template <class FluidState, class ParameterCache>
-    static void solve(FluidState & fluidState,
-                      ParameterCache & paramCache,
-                      const bool setViscosity,
-                      const bool setEnthalpy)
-    {
-
-        // In this function the actual work is done.
-        // Either tables for solubility or functional relations are used therein.
-        // One way or the other, this function needs to set all the mole fractions
-        // in the fluidstate.
-        FluidSystem::calculateEquilibriumMoleFractions(fluidState, paramCache);
-
-        for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
-            Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
-            fluidState.setDensity(phaseIdx, value);
-
-            if (setViscosity) {
-                value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
-                fluidState.setViscosity(phaseIdx, value);
-            }
-
-            if (setEnthalpy) {
-                value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-                fluidState.setEnthalpy(phaseIdx, value);
-            }
-        }
-    }
-};
-
+using FluidSystemConstraintSolver DUNE_DEPRECATED_MSG("Use MiscibleMultiPhaseComposition from dumux/material/constraintsolvers/misciblemultiphasecomposition.hh")
+= MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
 } // end namespace Dumux
-
 #endif
diff --git a/dumux/material/fluidsystems/h2on2kinetic.hh b/dumux/material/fluidsystems/h2on2kinetic.hh
index 20f84ce8fa420807fa4d5f3686db613c22f05726..a84dd2100a5679db0b5ed0c2baa0d3582258ffd1 100644
--- a/dumux/material/fluidsystems/h2on2kinetic.hh
+++ b/dumux/material/fluidsystems/h2on2kinetic.hh
@@ -129,7 +129,9 @@ public:
      *        \param calcCompIdx The component for which the composition in the other phase is to be
      *               calculated.
      */
+
     template <class FluidState>
+    DUNE_DEPRECATED_MSG("FluidSystems should not compute equilibrium mole fractionos. Please use a constraintsolver e.g. the MiscibleMultiPhaseComposition instead")
     static void calculateEquilibriumMoleFractionOtherPhase(FluidState & fluidState,
                                                     const ParameterCache & paramCache,
                                                     const unsigned int referencePhaseIdx,
@@ -259,6 +261,7 @@ public:
      *        \param paramCache A container for iterative calculation of fluid composition
      */
     template <class FluidState>
+    DUNE_DEPRECATED_MSG("FluidSystems should not compute equilibrium mole fractionos. Please use a constraintsolver e.g. the MiscibleMultiPhaseComposition instead")
     static void calculateEquilibriumMoleFractions(FluidState & fluidState,
                                                   const ParameterCache & paramCache)
     {
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 25d711f61a709f05107405dd889fc5b26b98e5b1..0ecec664225ce30a0307ab8e29af20ec02ea2faf 100644
--- a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
@@ -133,6 +133,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 {
@@ -269,7 +270,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_ ) ;
@@ -418,7 +422,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