Machine learning of accurate energy-conserving molecular force fields

Abstract

Using conservation of energy-a fundamental property of closed classical and quantum mechanical systems-we develop an efficient gradient-domain machine learning (GDML) approach to construct accurate molecular force fields using a restricted number of samples from ab initio molecular dynamics (AIMD) trajectories. The GDML implementation is able to reproduce global potential energy surfaces of intermediate-sized molecules with an accuracy of 0.3 kcal mol^-1 for energies and 1 kcal mol^-1 Å̊^-1 for atomic forces using only 1000 conformational geometries for training. We demonstrate this accuracy for AIMD trajectories of molecules, including benzene, toluene, naphthalene, ethanol, uracil, and aspirin. The challenge of constructing conservative force fields is accomplished in our work by learning in a Hilbert space of vector-valued functions that obey the law of energy conservation. The GDML approach enables quantitative molecular dynamics simulations for molecules at a fraction of cost of explicit AIMD calculations, thereby allowing the construction of efficient force fields with the accuracy and transferability of high-level ab initio methods.

Keywords: energy conservation; force field; gradient field; kernel regression; machine learning; molecular dynamics; path integrals; potential-energy surface.

Fig. 1. The construction of ML models: First, reference data from an MD trajectory are sampled.

(A) The geometry of each molecule is encoded in a descriptor. This representation introduces elementary transformational invariances of energy and constitutes the first part of the prior. A kernel function then relates all descriptors to form the kernel matrix—the second part of the prior. The kernel function encodes similarity between data points. Our particular choice makes only weak assumptions: It limits the frequency spectrum of the resulting model and adds the energy conservation constraint. Hess, Hessian. (C) These general priors are sufficient to reproduce good estimates from a restricted number of force samples. (B) A comparable energy model is not able to reproduce the PES to the same level of detail.

Fig. 2. Modeling the true vector field (leftmost subfigure) based on a small number of vector samples

With GDML, a conservative vector field estimate ${\widehat{\mathit{f}}}_{\mathrm{F}}$ is obtained directly. A naïve estimator ${\displaystyle {\widehat{\mathit{f}}}_{\mathrm{F}}^{-}}$ with independent predictions for each element of the output vector is not capable of imposing energy conservation constraints. We perform a Helmholtz decomposition of this nonconservative vector field to show the error component that violates the law of energy conservation. This is the portion of the overall prediction error that was avoided with GDML because of the addition of the energy conservation constraint.

Fig. 3. Efficiency of the GDML predictor versus a model that has been trained on energies.

(A) Required number of samples for a force prediction performance of MAE (1 kcal mol⁻¹ Å⁻¹) with the energy-based model (gray) and GDML (blue). The energy-based model was not able to achieve the targeted performance with the maximum number of 63,000 samples for aspirin. (B) Force prediction errors for the converged models (same number of partial derivative samples and energy samples). (C) Energy prediction errors for the converged models. All reported prediction errors have been estimated via cross-validation.

Fig. 4. Results of classical and PIMD simulations.

The recently developed estimators based on perturbation theory were used to evaluate structural and electronic observables (30). (A) Comparison of the interatomic distance distributions, $\mathit{h}(\mathit{r})={\langle \frac{2}{\mathit{N}(\mathit{N}-1)}{\displaystyle {\sum }_{\mathit{i}<\mathit{j}}^{\mathit{N}}}\mathrm{\delta }(\mathit{r}-||\overrightarrow{{\mathit{r}}_{\mathit{i}}}-\overrightarrow{{\mathit{r}}_{\mathit{j}}}||)\rangle }_{\mathit{P},\mathit{t}}$, obtained from GDML (blue line) and DFT (dashed red line) with classical MD (main frame), and PIMD (inset). a.u., arbitrary units. (B) Probability distribution of the dihedral angles (corresponding to carboxylic acid and ester functional groups) using a 20 ps time interval from a total PIMD trajectory of 200 ps.