  Molecular-orbital-based machine learning (MOB-ML) enables the prediction of
accurate correlation energies at the cost of obtaining molecular orbitals.
Here, we present the derivation, implementation, and numerical demonstration of
MOB-ML analytical nuclear gradients which are formulated in a general
Lagrangian framework to enforce orthogonality, localization, and Brillouin
constraints on the molecular orbitals. The MOB-ML gradient framework is general
with respect to the regression technique (e.g., Gaussian process regression or
neural networks) and the MOB feature design. We show that MOB-ML gradients are
highly accurate compared to other ML methods on the ISO17 data set while only
being trained on energies for hundreds of molecules compared to energies and
gradients for hundreds of thousands of molecules for the other ML methods. The
MOB-ML gradients are also shown to yield accurate optimized structures, at a
computational cost for the gradient evaluation that is comparable to
Hartree-Fock theory or hybrid DFT.
