Training a force field for proteins and small molecules from scratch

Abstract

Force fields for molecular dynamics are usually developed manually, limiting their transferability and making systematic exploration of functional forms challenging. We developed a graph neural network that assigns all force field parameters for diverse molecules using continuous atom typing. The freely-available model, called Garnet, was trained on quantum mechanical, condensed phase and protein nuclear magnetic resonance data without the use of existing parameters. The resulting force field shows comparable performance to current force fields on small molecules, folded proteins, protein complexes and disordered proteins. It shows similar results to popular approaches for relative binding free energy predictions across a range of targets. Assessing different functional forms shows that the double exponential potential is a flexible and accurate alternative to the Lennard-Jones potential. Garnet provides a platform for automated, reproducible force field discovery that brings the benefits of machine learning to classical force fields.

Figures (15)

• Figure 1: The Garnet model. (A) An overview of the model architecture, which is inspired by Espaloma Takaba2024Wang2022. NN refers to a fully-connected neural network. The double exponential $\alpha$ and $\beta$ global parameters are also trained by the model. (B) The order of batches during training. Calculating gradients to train on GB3 NMR data with ensemble reweighting is slow, so gradients are re-used and new gradients are only computed after 100 batches. A random subset of 200 of the available frames is used to calculate the gradients each time. (C) The double exponential $\sigma$ and $\varepsilon$ parameters for atoms in the SPICE test set calculated with the trained model. $\sigma$ and $\varepsilon$ are similar to the corresponding parameters in the Lennard-Jones potential with $\sigma$ being proportional to the van der Waals radius. By construction, $\sigma$ is allowed to range from 0.05 nm to 0.5 nm and $\varepsilon$ is allowed to range from 0.02 kJ/mol to 1.5 kJ/mol.
• Figure 1: Exploring the embedding space of continuous atom types. Ubiquitin was parameterised with Garnet and the 64-dimensional embedding of each atom was calculated. t-SNE was used to reduce this to two dimensions for visualisation. The corresponding Amber14SB atom types found in different regions are shown, along with water atoms (Owat and Hwat).
• Figure 2: Energy minimisation of small molecules in the OpenFF Industry Benchmark dataset Horton2025. (A) Violin plots for the different structural deviation metrics. The whisker range is three times the median absolute deviation, clamped to zero when needed. Points falling outside this range are considered outliers, are not included in the violin density estimation, and are shown as small "x" markers. The black circular marker denotes the median of each distribution. (B) An example of molecules well captured (top) and poorly captured (bottom) by the Garnet force field, with the QM structure in green and the MM structure in blue. (C) Violin plots for the energy deviation metric, $\Delta\Delta\mathrm{E}$, showing how well energy differences near the QM minimum are captured by the force fields. Values near zero represent accurate energy wells.
• Figure 2: Partial charge distribution in Garnet. The partial charge of each atom in ubiquitin (1,231 atoms) is assigned with Garnet or Amber14SB Case2025; the distribution of the absolute values of the partial charges is shown. For Garnet, the mean is 0.33 and the median is 0.19. For Amber14SB, the mean is 0.25 and the median is 0.11.
• Figure 3: Performance of Garnet on simulating folded proteins. (A) RMSD to the native structure for the simulated proteins, each replica represented in a different line style. The RMSD is smoothed by taking the mean over a window of values extending 10 snapshots either side. (B) RMSE for the hydrogen-bond related scalar couplings of GB3 and Ubq. (C) ANE for the scalar couplings of the four studied proteins. In all cases the error bars estimate a 99.9% confidence interval, computed through bootstrapping. HEWL lacked data for the couplings related to backbone torsions. GB3 was used during model training.