  Time-resolved diffuse scattering experiments have gained increasing attention
due to their potential to reveal non-equilibrium dynamics of crystal lattice
vibrations with full momentum resolution. Although progress has been made in
interpreting experimental data on the basis of one-phonon scattering,
understanding the role of individual phonons can be sometimes hindered by
multi-phonon excitations. In Ref. [arXiv:2103.10108] we have introduced a
rigorous approach for the calculation of the all-phonon inelastic scattering
intensity of solids from first-principles. In the present work, we describe our
implementation in detail and show that multi-phonon interactions are captured
efficiently by exploiting translational and time-reversal symmetries of the
crystal. We demonstrate its predictive power by calculating the diffraction
patterns of monolayer molybdenum disulfide (MoS$_2$), bulk MoS$_2$, and black
phosphorus (bP), and we obtain excellent agreement with our measurements of
thermal electron diffuse scattering. Remarkably, our results show that
multi-phonon excitations dominate in bP across multiple Brillouin zones, while
in MoS$_2$ play a less pronounced role. We expand our analysis for each system
and we examine the effect of individual atomic and interatomic vibrational
motion on the diffuse scattering signals. Our findings indicate that distinct
features are explained by the collective displacement of MoS and specific pairs
of P atoms. We further demonstrate that the special displacement method
reproduces the thermally distorted configuration which generates precisely the
all-phonon diffraction pattern. The present methodology opens the way for
high-throughput calculations of the scattering intensity in crystals and the
accurate interpretation of static and time-resolved diffuse scattering
experiments.
