  We introduce a novel approach for a fully quantum description of coupled
electron-ion systems from first principles. It combines the variational quantum
Monte Carlo (QMC) solution of the electronic part with the path integral (PI)
formalism for the quantum nuclear dynamics. On the one hand, the PI molecular
dynamics includes nuclear quantum effects by adding a set of fictitious
classical particles (beads) aimed at reproducing nuclear quantum fluctuations
via a harmonic kinetic term. On the other hand, variational QMC can provide
Born-Oppenheimer (BO) potential energy surfaces with a precision comparable to
the most advanced post Hartree-Fock approaches, and with a favorable scaling
with the system size. To deal with the intrinsic QMC noise, we generalize the
PI molecular dynamics using a Langevin thermostat correlated according to the
covariance matrix of QMC nuclear forces. The variational parameters of the QMC
wave function are evolved during the nuclear dynamics, such that the BO
potential energy surface is unbiased. Statistical errors on the wave function
parameters are reduced by resorting to bead grouping average, which we show to
be accurate and well controlled. Our general algorithm relies on a Trotter
breakup between the dynamics driven by ionic forces and the one set by the
harmonic interbead couplings. The latter is exactly integrated even in presence
of the Langevin thermostat, thanks to the mapping onto an Ornstein-Uhlenbeck
process. This framework turns out to be very efficient also in the case of
deterministic ionic forces. The new implementation is validated on the Zundel
ion by direct comparison with standard PI Langevin dynamics calculations made
with a coupled cluster potential energy surface. Nuclear quantum effects are
confirmed to be dominant over thermal effects well beyond room temperature
giving the excess proton an increased mobility by quantum tunneling.
