A recently developed lower bound theory for Coulombic problems (E. Pollak, R. Martinazzo, J. Chem. Theory Comput. 2021Comput. , 17, 1535 is further developed and applied to the highly accurate calculation of the ground-state energy of two-(He, Li + , and H − ) and three-(Li) electron atoms. The method has been implemented with explicitly correlated many-particle basis sets of Gaussian type, on the basis of the highly accurate (Ritz) upper bounds they can provide with relatively small numbers of functions. The use of explicitly correlated Gaussians is developed further for computing the variances, and the necessary modifications are here discussed. The computed lower bounds are of submilli-Hartree (parts per million relative) precision and for Li represent the best lower bounds ever obtained.Although not yet as accurate as the corresponding (Ritz) upper bounds, the computed bounds are orders of magnitude tighter than those obtained with other lower bound methods, thereby demonstrating that the proposed method is viable for lower bound calculations in quantum chemistry applications. Among several aspects, the optimization of the wave function is shown to play a key role for both the optimal solution of the lower bound problem and the internal check of the theory.