We analyze the performance of gradient-free local density approximation (LDA) and gradient-dependent generalized gradient approximation (GGA) functionals in a density functional theory variational calculations based on the total energy bifunctional (E[rho(1),rho(2)]). These approximations are applied to the exchange-correlation energy and to the nonadditive component of the kinetic energy of the complex. Benchmark ab initio interaction energies taken from the literature for 25 intermolecular complexes for which the interaction energies fall into the 0.1-3.0 kcal/mol range are used as reference. At the GGA level, the interaction energies derived from E[rho(1),rho(2)] are more accurate than the Kohn-Sham ones. LDA leads to very good interaction energies for such complexes where the rho(1),rho(2) overlap is very small (Ne-Ne, Ar-Ar, for instance) but it is not satisfactory for such cases where the overlap is larger. Introduction of gradient-dependent terms into the approximate part of E[rho(1),rho(2)] improves significantly the overall accuracy of the interaction energies. Gradient-dependent functionals applied in E[rho(1),rho(2)] lead to the average error and the average absolute error of the interaction energies amounting to 0.08 kcal/mol and 0.29 kcal/mol, respectively. (C) 2003 American Institute of Physics.