  Harmonic calculations based on density-functional theory are generally the
method of choice for the description of phonon spectra of metals and
insulators. The inclusion of anharmonic effects is, however, delicate as it
relies on perturbation theory requiring a considerable amount of computer time,
fast increasing with the cell size. Furthermore, perturbation theory breaks
down when the harmonic solution is dynamically unstable or the anharmonic
correction of the phonon energies is larger than the harmonic frequencies
themselves.We present a stochastic implementation of the self-consistent
harmonic approximation valid to treat anharmonicity at any temperature in the
non-perturbative regime. The method is based on the minimization of the free
energy with respect to a trial density matrix described by an arbitrary
harmonic Hamiltonian. The minimization is performed with respect to all the
free parameters in the trial harmonic Hamiltonian, namely, equilibrium
positions, phonon frequencies and polarization vectors. The gradient of the
free energy is calculated following a stochastic procedure. The method can be
used to calculate thermodynamic properties, dynamical properties and anharmonic
corrections to the Eliashberg function of the electron-phonon coupling. The
scaling with the system size is greatly improved with respect to perturbation
theory. The validity of the method is demonstrated in the strongly anharmonic
palladium and platinum hydrides. In both cases we predict a strong anharmonic
correction to the harmonic phonon spectra, far beyond the perturbative limit.
In palladium hydrides we calculate thermodynamic properties beyond the
quasiharmonic approximation, while in PtH we demonstrate that the high
superconducting critical temperatures at 100 GPa predicted in previous
calculations based on the harmonic approximation are strongly suppressed when
anharmonic effects are included.
