A general computational scheme for the (nonrelativistic) Bethe logarithm is developed, opening the route to "routine" evaluation of the leading-order quantum electrodynamics correction (QED) relevant for spectroscopic applications for small polyatomic and polyelectronic molecular systems. The implementation relies on the Schwartz method and minimization of a Hylleraas functional. In relation with electronically excited states, a projection technique is considered, which ensures positive definiteness of the functional over the entire parameter (photon momentum) range. Using this implementation, the Bethe logarithm is converged to a relative precision better than 1:103 for selected electronic states of the two-electron H2 and H3+, and the three-electron He2+ and H+H2 molecular systems. The present work focuses on nuclear configurations near the local minimum of the potential energy surface, but the computations can be repeated also for other structures.