the low-fidelity potential modified by perturbations of varying $\epsilon$ and computing a numerical derivative for each position $r_i$. These modified Lennard Jones potentials and resulting average potential energies as a function of size of the perturbation $\epsilon$ are shown in Fig. \[fig:calcFD\](a) and (b), respectively. However, this brute force approach is computationally very intensive considering that we may need to sample the functional derivative at 100 values of $r_i$. Even with only three values of $\epsilon$ for each separation distance one would need to perform 300 separate MD simulations. To alleviate the computational cost of the approach we now derive a perturbative approach to calculate the functional derivative that can be computed with little overhead with respect to the nominal simulation using the low-fidelity potential, $\phi_0$.

At the heart of the calculation of the FD is the evaluation of canonical averages of the QoIs with the interatomic potential and a Gaussian perturbation. Recognizing that the modified Hamiltonian of the system can be additively decomposed into the kinetic energy, the original LJ potential energy and the functional variation (or perturbation) for these two-body potentials we can write the canonical ensemble average of quantity Q as: $$\label{eq:ensavg}
\langle Q \rangle_H = \frac{\int 