The phase and chemical equilibrium problem is of paramount importance for predicting fluid phase behavior for a very large number of separation process applications. The ubiquity of the flash calculation in chemical engineering is just one example of its prevalence. Process simulators need to be able to reliably and efficiently predict the correct number of phases that will exist at equilibrium, and the distribution of components within those phases. However, the equation based and local optimisation approaches which are in common use can provide no theoretical guarantee that the true equilibrium solution will be obtained. For conditions of constant temperature and pressure, a necessary and sufficient condition for the true equilibrium solution is that (i) the total Gibbs free energy of the system be at its global minimum, or (ii) the minimum of the tangent plane distance function be nonnegative for all phase models used to represent the system. In this contribution, the goal is to obtain equilibrium solutions corresponding to a global minimum of the Gibbs free energy as efficiently as possible, for cases where the liquid phase can be modeled by the NRTL UNIQUAC, UNIFAC, Wilson, modified Wilson and ASOG equations. Vapor phases whose behavior can be described as ideal can also be handled. In achieving this goal, there are two distinct problems of relevance: (i) the minimisation of the Gibbs free energy, denoted (G), and (ii) the minimisation of the tangent plane distance function, or the tangent plane stability criterion, denoted (S). The approach builds on previous contributions by the authors. McDonald and Floudas (Comput. chem. Engng., 19(11), 1995) use the Global optimization, GOP, algorithm of Floudas and Visweswaran (1990, 1993) to obtain global solutions for (G) when employing the NRTL equation. For the UNIQUAC equation, McDonald and Floudas (JOGO, 5:205, 1993) show how the Gibbs function can be transformed into the difference of two convex functions, so that a branch and bound deterministic global optimisation algorithm can be used to solve (G) globally. McDonald and Floudas (AIChE J., 41(7):1798, 1995) demonstrate how global solutions can be obtained for (S) when using the NRTL or UNIQUAC equations. McDonald and Floudas (I&EC Res., 34.1674, 1995) show how a similar branch and bound algorithm can also be applied for the UNIFAC, Wilson, modified Wilson and ASOG models to solve both (G) and (S). The key contribution is the presentation of an algorithm called GLOPEQ (Global OPtimization for the Phase and chemical EQuilibrium problem) which is theoretically guaranteed to converge to the global (sic equilibrium) solution for the activity coefficient models indicated, no matter the starting point. GLOPEQ solves (G) and (S) in tandem, using (G) to generate candidate equilibrium solutions which can then be verified for thermodynamic stability by solving (S). Two key features of the combined algorithm are that (i) as much information as is possible is obtained from local searches, and (ii) it is preferable to verify a globally stable equilibrium solution using the tangent plane criterion, as this problem contains fewer variables than the minimisation of the Gibbs free energy. Both these features improve computational performance greatly. The GLOPEQ algorithm uses the tangent plane criterion to establish if the Gibbs free energy of a current equilibrium solution can be reduced. However, the main difference is that deterministic global optimisation is used to solve the subproblems (G) and (S) as opposed to local optimisation techniques, thus allowing the above guarantees to be made, unlike all other previous approaches.