In this study, the PC-SAFT equation of state is used for vapor-liquid equilibrium calculations using as independent variables the mixture composition, density and temperature. The method is based on unconstrained minimisation of the Helmholtz Free energy via a combination of the successive substitution iteration and Newton-Raphson minimisation methods with line-search; the positive definiteness of the Hessian is guaranteed by a modified Cholesky decomposition. The algorithm consists of two stages; initially, the mixture is assumed to be a single-phase and its stability is assessed; in case of being found unstable, a second stage of phase splitting (flash) takes place, in which the pressure of the fluid and compositions of both the liquid and vapor phases are calculated. The reliability of two different methods presented in the existing literature, (i) using mole numbers and (ii) using the logarithm of the equilibrium constants as iterative variables, is evaluated in terms of both iterations and computational time needed to reach convergence, for seven test cases. These include both single and multicomponent Diesel fuel surrogates, known to give incomplete density information when using pressure and temperature as independent variables. Results show that iterating with the logarithm of the equilibrium constants also reproduces this issue, while it requires a smaller number of iterations than using with mole numbers as independent variables. However, the total computational time needed for the latter case is vastly inferior. Pressure and vapor volume fraction fields are discussed for a range of temperatures and densities, apart from the number of iterations needed during the flash calculation stage. A performance comparison is obtained against the Peng-Robinson equation of state, showing similar number of iterations required but a computational time increasing with the number of components. While for a single component PC-SAFT needs around 3 times more CPU time, for 4 components it is 6 times and for a mixture of 8 components it increases up to 14 times. Finally, the method is demonstrated to converge unconditionally for all conditions tested.