Compositional multiphase flow in subsurface porous media is becoming increasingly attractive due to issues related with enhanced oil recovery, greenhouse effect and global warming, and the urgent need for development in unconventional oil/gas reservoirs. One key effort prior to construct the mathematical model governing the compositional multiphase flow is to determine the phase compositions of the fluid mixture, and then calculate other related physical properties. In this paper, recent progress on phase equilibrium calculations in subsurface reservoirs have been reviewed and concluded with authors’ own analysis. Phase equilibrium calculation is the main approach to perform such calculation, which could be conducted using two different types of flash calculation algorithms: The NPT flash and NVT flash. NPT flash calculations are proposed early, well developed within the last few decades and now become the most commonly used method. However, it fails to remain the physical meanings in the solution as a cubic equation, derived from equation of state, is often needed to solve. Alternatively, NVT flash can handle the phase equilibrium calculations as well, without the pressure known a priori. Recently, Diffuse Interface Models, which were proved to keep a high consistency with thermodynamic laws, have been introduced in the phase calculation, incorporating the realistic equation of state (EOS), e.g. Peng-Robinson EOS. In NVT flash, Helmholtz free energy is minimized instead of Gibbs free energy used in NPT flash, and this energy density is treated with convex-concave splitting technique. A semi-implicit numerical scheme is designed to process the dynamic model, which ensures the thermodynamic stability and then preserve the fast convergence property. A positive definite coefficient matrix is designed to meet the Onsager Reciprocal Principle so as to keep the entropy increasing property in the presence of capillary pressure, which is required by the thermodynamic laws. The robustness of the proposed algorithm is verified via two numerical examples, one of which has up to seven components. In the complex fluid mixture, special phenomena could be capture from the global minimum of TPD functions as well as the phase envelope resulted from the phase equilibrium calculations. It can be found that the boundary between the single-phase and vapor–liquid phase regions will move in the presence of capillary pressure, and then the area of each region will change accordingly. Some remarks have been concluded at the end, as well as suggestions on potential topics for future studies.