to optimizing the model parameters to minimize this loss, summed over all DFT calculations in the dataset. We select the characteristic length scale $\ell_0 = 0.01$ to balance energy and force training terms.

Standard optimization techniques require calculating the gradient of $L$ with respect to all model parameters $W$ (modern ML potentials might have $\sim 10^5$ trainable model parameters). Because forces $\hat{\mathbf f}_{j}$ appear in $L$, calculating $\mathrm{d} L / \mathrm{d} W$ involves *second* derivatives of the ANI energy output, i.e., $\mathrm{d}^2 \hat E / \mathrm{d} W \mathrm{d} \mathbf r_j$. Fortunately, direct calculation of all such second partial derivatives is not required. Frameworks such as TensorFlow or PyTorch support iterated backpropagation, thereby enabling efficient calculation of the full gradient $\mathrm{d} L / \mathrm{d} W$ [@Yao2017ThePhysics; @Lubbers2018; @Zhang2018DeepMechanics]. However, we train our models in Neurochem, a C++/CUDA implementation which is many times faster than comparable PyTorch implementations. The source code transformation required to perform backpropagation *twice* would be very difficult to perform manually within a hand-written CUDA/C++ code.

To train to force data, we have devised an efficient finite differencing scheme to approximate $\mathrm{d} L / \mathrm{d} W$. Our scheme involves just two first-order backpropagation calculations, $\mathrm{d} \hat{E}_i[\mathbf r^\pm] / \mathrm{d} 