We review the Bardeen-Cooper-Schrieffer mean-field theory emphasizing its origin as a variational approximation for the grand potential. This is done by using the Bogoliubov inequality as the starting point. Then we write the mean-field grand potential as an explicit function of the one-particle density matrix, which turns out to be a natural generalization of the Mermin functional. This result opens the way for the application to superconducting systems of the linear scaling methods developed in the context of electronic structure theory. Finally, we show that computing the superfluid weight from the derivatives of the mean-field grand potential naturally leads to the generalized random phase approximation. Our results showcase the advantage of a density matrix-based approach and are potentially interesting for the study of disordered superconductors and superconductors with large unit cell, such as twisted bilayer graphene and other moiré materials.