  In this work we develop and implement a method for calculation of the Bethe
logarithm for many-electron atoms. This quantity is required to evaluate the
leading-order quantum electrodynamics correction to the energy and properties
of atomic and molecular systems beyond the Dirac theory (the Lamb shift). The
proposed formalism is based on the mean-field representation of the
ground-state electronic wavefunction and of the response functions required in
the Schwartz method [C. Schwartz, Phys. Rev. {\bf 123}, 1700 (1961)]. We
discuss difficulties encountered in the calculations with the emphasis on the
specific basis set requirements in the vicinity of the atomic nucleus. This
problem is circumvented by introducing a modified basis set of exponential
functions which are able to accurately represent the gradient of hydrogen-like
orbitals. The Bethe logarithm is computed for ground electronic states of atoms
from hydrogen to magnesium and, additionally, for argon. Whenever possible, the
results are compared with the available reference data from the literature. In
general, the mean-field approximation introduces a surprisingly small error in
the calculated values, suggesting that the electron correlation effects are of
minor importance in determination of the Bethe logarithm. Finally, we propose a
robust scheme to evaluate the Lamb shift for arbitrary light molecular systems
at little computational cost. As an illustration, the method is used to
calculate Lamb shifts of the vibrational levels of the nitrogen molecule.
