appendix. Use of the coordination number simply allows to further reduce computation and avoid modifying the MD code. Performing the calculation directly within the MD force loop with atomic positions would be equivalent excluding slight discrepancies from discretization.

Using Eq. \[eq:FD2\] with the first term in the difference calculated using Eq. \[eq:perturbative\] we compute the functional derivative with respect to Gaussian perturbations centered at $r_i$ with width $\sigma = 0.1$ Å and height $\epsilon=\pm0.00075$ and $\pm0.00375$ eV. The perturbation width was chosen to minimize (localize) the perturbation while retaining sufficient sampling. The heights were similarly minimized, due to the difficulty of converging exponentially weighted averages, while ensuring a measurable effect from the perturbations. The numerical derivative is then evaluated by computing the slope of the QoI with respect to $\epsilon$ (Fig. \[fig:calcFD\](b)) at a set of interatomic separations ranging from zero past the potential cutoff in increments of $0.05$ Å.

The thick blue lines in Figure \[fig:FD\] show the functional derivatives obtained in this manner for all physical conditions and both QoIs. These curves were averaged from multiple calculations of the FD with independent randomized samplings from the total of 64 ns of simulation time. The functional derivative displays significant 