  We describe a method for imposing the correct electron-nucleus (e-n) cusp in
molecular orbitals expanded as a linear combination of (cuspless) Gaussian
basis functions. Enforcing the e-n cusp in trial wave functions is an important
asset in quantum Monte Carlo calculations as it significantly reduces the
variance of the local energy during the Monte Carlo sampling. In the method
presented here, the Gaussian basis set is augmented with a small number of
Slater basis functions. Note that, unlike other e-n cusp correction schemes,
the presence of the Slater function is not limited to the vicinity of the
nuclei. Both the coefficients of these cuspless Gaussian and cusp-correcting
Slater basis functions may be self-consistently optimized by diagonalization of
an orbital-dependent effective Fock operator. Illustrative examples are
reported for atoms (\ce{H}, \ce{He} and \ce{Ne}) as well as for a small
molecular system (\ce{BeH2}). For the simple case of the \ce{He} atom, we
observe that, with respect to the cuspless version, the variance is reduced by
one order of magnitude by applying our cusp-corrected scheme.
