We present a new method to solve general systems of equations containing complementarity conditions, with a special focus on those arising in the thermodynamics of multicomponent multiphase mixtures at equilibrium. Indeed, the unified formulation introduced by Lauser et al. [Adv. Water Res. 34 (2011), 957-966] has recently emerged as a promising way to automatically handle the appearance and disappearance of phases in porous media compositional multiphase flows. From a mathematical viewpoint and after discretization in space and time, this leads to a system consisting of algebraic equations and nonlinear complementarity equations. Such a system exhibit serious convergence difficulties for the existing semismooth and smoothing methods. This observation led us to design a new strategy called NPIPM (Non-Parametric Interior-Point Method ). Inspired from interior-point methods in optimization, the technique we propose avoids any parameter management while ensuring good theoretical convergence results. These are validated by extensive numerical tests, in which we compare NPIPM to the Newton-min method.