An ab initio method is developed for variational grand-canonical molecular electronic structure of open systems based on the Gibbs–Peierls–Boguliobov inequality. We describe the theory and a practical method for performing the calculations within standard quantum chemistry codes using Gaussian basis sets. The computational effort scales similarly to the ground-state Hartree–Fock method. The quality of the approximation is studied on a hydrogen molecule by comparing to the exact Gibbs free energy, computed using full configuration-interaction calculations. We find the approximation quite accurate, with errors similar to those of the Hartree–Fock method for ground-state zero-temperature calculations. A further demonstration is given of the temperature effects on the bending potential curve for water. Some future directions and applications of the method are discussed. Several appendices give the mathematical and algorithmic details of the method.