john.chodera<EMAIL_ADDRESS>
# Espaloma-0.3.0: Machine-learned molecular mechanics force field for the
simulation of protein-ligand systems and beyond.

Kenichiro Takaba (ORCID: 0000-0002-2481-8830)  Computational and Systems
Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer
Center, New York, N.Y. 10065  Pharmaceutical Research Center, Advanced Drug
Discovery, Asahi Kasei Pharma Corporation, Shizuoka 410-2321, Japan  Iván
Pulido (ORCID: 0000-0002-7178-8136)  Computational and Systems Biology
Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center,
New York, N.Y. 10065  Mike Henry (ORCID: 0000-0002-3870-9993)  Computational
and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan
Kettering Cancer Center, New York, N.Y. 10065  Hugo MacDermott-Opeskin
(ORCID:0000-0002-7393-7457)  Computational and Systems Biology Program, Sloan
Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, N.Y.
10065  John D. Chodera (ORCID: 0000-0003-0542-119X)  Computational and Systems
Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer
Center, New York, N.Y. 10065  Yuanqing Wang111Work partly done while at
Memorial Sloan Kettering Cancer Center. (ORCID: 0000-0003-4403-2015)  New York
University, New York, N.Y. 10004  Computational and Systems Biology Program,
Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York,
N.Y. 10065

###### Abstract

Molecular mechanics (MM) force fields—the models that characterize the energy
landscape of molecular systems via simple pairwise and polynomial terms—have
traditionally relied on human expert-curated, inflexible, and poorly
extensible discrete chemical parameter assignment rules, namely atom or
valence types. Recently, there has been significant interest in using graph
neural networks to replace this process, while enabling the parametrization
scheme to be learned in an end-to-end differentiable manner directly from
quantum chemical calculations or condensed-phase data. In this paper, we
extend the Espaloma [75] end-to-end differentiable force field construction
approach by incorporating both energy and force fitting directly to quantum
chemical data into the training process. Building on the OpenMM SPICE dataset
[23], we curate a dataset containing chemical spaces highly relevant to the
broad interest of biomolecular modeling, covering small molecules, proteins,
and RNA. The resulting force field, espaloma 0.3.0, self-consistently
parametrizes these diverse biomolecular species, accurately predicts quantum
chemical energies and forces, and maintains stable quantum chemical energy-
minimized geometries. Surprisingly, this simple approach produces highly
accurate protein-ligand binding free energies when self-consistently
parametrizing protein and ligand. This approach—capable of fitting new force
fields to large quantum chemical datasets in one GPU-day—shows significant
promise as a path forward for building systematically more accurate force
fields that can be easily extended to new chemical domains of interest. The
espaloma 0.3.0 force field is available for use directly or within OpenMM via
the open-source Espaloma package [https://github.com/choderalab/espaloma], and
both the code and datasets for constructing this force field are openly
available [ https://github.com/choderalab/refit-espaloma].

Molecular mechanics force fields [15, 29] have long been a workhorse for
biomolecular simulation and computer-aided drug design [45, 66]. These
models—characterizing the conformational energetics and noncovalent
interactions—have become instrumental to a wide variety of in silico
biomolecular modeling tasks, including enumeration of putative bioactive
conformations [14], hit identification via virtual screening [5], prediction
of membrane permeability [71], simulations of biomolecular dynamics [61], and
estimation of protein:ligand binding free energies via alchemical free energy
calculations [54].

Traditional molecular mechanics force fields employ a set of rules termed atom
typing to group atoms to share sets of parameters based on their elemental
identity, hybridization, aromaticity, and number and properties of neighbors,
subsequently assigning tabulated parameters based on the atom types composing
the valence term [56, 15, 29]. This practice manifests several pathologies:
Firstly, atoms of distinct chemical environments might share the same type and
thereby parameters, leading to suboptimal chemical resolution; reducing errors
requires the introduction of more fine-grained atom type descriptions,
resulting in an explosion in the number of necessary atom types and parameter
assignment table entries, especially for torsions [56]. Secondly, atom typing
rules are traditionally constructed by human experts, whose design optimality
cannot be quantitatively assessed; also, the evolution of such methodologies
would constantly require engineer time and is not sustainable. Finally,
discrepancies in these assignments might arise among different implementations
in chemoinformatics toolkits [56].

More concerningly, biomolecular systems relevant to drug discovery are
frequently highly heterogeneous. Atom type-based force fields are frequently
developed separately for each class of biomolecules, enormously increasing the
complexity of ensuring these force fields remain compatible, and risking poor
accuracy when multiple classes of biomolecules must interact. As such,
extension or expansion to new classes of biomolecules or chemical spaces
becomes a time-consuming ordeal. For instance, the recent AmberTools 23
release [8] recommends combining independently developed force fields to
simulate systems containing proteins [53], DNA [82], RNA [60], water [31],
monovalent [36, 37] and divalent [48, 47, 49] counterions, lipids [28],
carbohydrates [42], glycoconjugates [17, 18], small molecules [73, 72], post-
translational modifications [38], and nucleic acid modifications [70].

Recently, Espaloma [75] (extendable surrogate potential optimized by message-
passing) emerges as a novel approach to parametrize molecular mechanics (MM)
force fields. Rather than relying on human experts, Espaloma (Figure 1)
employs a graph neural network model (conditioned on neural network parameters
$\Phi_{\mathtt{NN}}$) to construct continuous vectorial representations of
atoms based on their chemical environments (Stage 1), combining them in a
symmetry-preserving manner, to produce invariant latent-space representations
of atom, bond, angle, and torsion chemical environments (Stage 2). These atom,
bond, angle, and torsion representations can then be fed into neural modules
to produce molecular mechanics parameters $\Phi_{\mathtt{FF}}$ directly (Stage
3), rather than requiring lookup tables to assign parameters—enabling smooth
interpolation between known chemical environments. Atomic partial charges can
also be generated by this approach, with recent applications using a geometry-
independent charge equilibration approach to rapidly generate AM1-BCC [33, 34]
quality charges [74, 77]. This framework also enables rapid experimentation
with new classes of molecular force fields, such as the addition of terms to
represent point polarizability or higher-order couplings between valence terms
in class II force fields [15, 29].

Constructing a complete force field with Espaloma is straightforward: Because
the entire process of parameter assignment is end-to-end differentiable, the
Espaloma neural network parameters ($\Phi_{\mathtt{NN}}$) can be optimized
directly using a standard machine learning framework to fit quantum chemical
and/or experimental data directly [75]. Operating on chemical graphs, Espaloma
can rapidly and self-consistently generate parameters for any biomolecules
with elemental coverage in its training set without the need to combine
distinct force fields. Earlier work has demonstrated how both proteins and
small molecules could rapidly and self-consistently be parameterized in under
a second using a model jointly trained on a limited set of quantum chemical
data covering small molecules and amino acids [75].

Here, we demonstrate how this approach can scale to foundational quantum
chemical datasets that provide extensive coverage of heterogeneous
biomolecular space. We present a new generation Espaloma
model—espaloma-0.3.0—trained on a broad chemical space that includes small
druglike molecules, peptides, nucleic acids, and fragments of a broad variety
of biomolecular species. We fit both energies and forces from a rich quantum
chemical dataset containing over 15,000 unique molecules and 1M conformations,
which augment both the foundational SPICE quantum chemical dataset [23] and
open quantum chemical datasets provided by the Open Force Field Initiative
[63] hosted through QCArchive [68]. To ensure the model remains well-behaved
when trained on very large datasets, we fine-tune Espaloma’s loss function to
incorporate signifantly more regularization to enhance its stability. With
extensive testing and benchmarking, we show that espaloma-0.3.0 is highly
capable in popular MD modeling tasks, including modeling protein-ligand
binding.

The main contributions of this paper are summarized as follows:

  * •

Chemically and conformationally diverse datasets of interest to biomolecular
modeling are curated to enable the construction of molecular mechanics (MM)
force fields that can be applied to small molecules, proteins, RNA, and
beyond.

  * •

A training strategy that uses both energies and forces along with
regularization is developed to enable Espaloma to parametrize force fields
that demonstrate robustness and good agreement with quantum chemical energy
surfaces.

  * •

The newly constructed force field, espaloma-0.3.0, produces not only accurate
quantum conformational energetics, but also accurate protein-ligand binding
free energy estimates vis-à-vis experimental data; this holds true for when
espaloma-0.3.0 is used to parametrize both protein and small molecule, or just
the small molecule (when combined with a protein force field espaloma-0.3.0 is
designed to be compatible with). The accuracy and correlation are on par with
popular legacy force fields produced by decades of expert-guided work.

##  1 Espaloma provides a flexible, end-to-end differentiable framework for
assigning molecular mechanics (MM) parameters using graph neural networks
(GNNs)

Figure 1:  Espaloma is an end-to-end differentiable molecular mechanics
parameter assignment scheme for arbitrary organic molecules. espaloma
(extendable surrogate potential optimized by message-passing) is a modular
approach for directly computing molecular mechanics force field parameters
$\Phi_{\mathrm{FF}}$ from a chemical graph $\mathcal{G}$ such as a small
molecule or biopolymer via a process that is fully differentiable in the model
parameters $\Phi_{\mathrm{NN}}$. In Stage 1, a graph neural network is used to
generate continuous latent atom embeddings describing local chemical
environments from the chemical graph. In Stage 2, these atom embeddings are
transformed into feature vectors that preserve appropriate symmetries for
atom, bond, angle, and proper/improper torsion inference via Janossy pooling
[57]. In Stage 3, molecular mechanics parameters are directly predicted from
these feature vectors using feed-forward neural nets. This parameter
assignment process is performed once per molecular species, allowing the
potential energy to be rapidly computed using standard molecular mechanics or
molecular dynamics frameworks thereafter. The collection of parameters
$\Phi_{\mathrm{NN}}$ describing the espaloma model can be considered as the
equivalent complete specification of a traditional molecular mechanics force
field such as GAFF [73, 72]/AM1-BCC [33, 34] in that it encodes the equivalent
of traditional typing rules, parameter assignment tables, and even partial
charge models. Figure reproduced from the arXiv preprint of Wang et al., [75]
under license http://arxiv.org/licenses/nonexclusive-distrib/1.0/.

Espaloma [75] (Figure 1) operates analogously to a legacy force field typing
scheme to generate molecular mechanics (MM) force field parameters
$\Phi_{\text{FF}}$ from a molecular graph $\mathcal{G}$ and graph neural
network (NN) parametrized by neural model parameters $\Phi_{\text{NN}}$,

$\Phi_{\text{FF}}\leftarrow\texttt{espaloma}(\mathcal{G},\Phi_{\text{NN}}).$
(1)

The resulting parameters $\Phi_{\text{FF}}$ can then be subsequently used in a
standard molecular mechanics package to compute the MM energy and forces for
any conformation, as with a standard MM force field.

Espaloma parametrizes molecular systems in three sequential stages (Figure 1):

Stage 1: Graph neural networks generate a continuous vectorial atom embedding,
replacing legacy discrete atom typing. First, using chemoinformatics toolkits
such as RDKit [64], the molecular system is abstracted as a graph, with nodes
denoting atoms and edges covalent bonds. Espaloma uses graph neural networks
[81, 76, 41, 3, 21, 79, 74, 78, 22, 26] as a replacement for rule-based
chemical environment perception [56] to operate on this graph. These neural
architectures learn useful representations of atomic chemical environments
from simple input features by updating and pooling embedding vectors via
message-passing iterations to neighboring atoms [26]. As such, symmetries in
chemical graphs (chemical equivalencies) are inherently preserved, while a
rich, continuous, and differentiably learnable representation of the atomic
environment is derived.

Stage 2: Symmetry-preserving pooling generates continuous bond, angle, and
torsion embeddings, replacing discrete types. The representations determined
by GNNs in Stage 1 are used to come up with bond, angle, and torsion
representations in a symmetry-preserving manner, where the relevant equivalent
atom permutations are listed and summed up via Janossy pooling [57].

Stage 3: Neural parametrization of atoms, bonds, angles, and torsions replaces
tabulated parameter assignment. In the final stage, feed-forward neural
networks learn the mapping from these symmetry-preserving invariant atom,
bond, angle, and torsion encodings to MM parameters $\Phi_{\text{FF}}$ that
reflect the specific chemical environments appropriate for these terms. Each
distinct parameter class (such as atom parameters, bond parameters, angle
parameters, and torsion parameters) is assigned by a separate neural network,
making this stage fully modular. This stage is analogous to the final table
lookup step in traditional force field construction but with significant added
flexibility arising from the continuous embedding that captures the chemical
environment specific to the potential energy term being assigned. The output
of Stage 3, and the entire pipeline, is a set of force field parameters
$\Phi_{\mathtt{FF}}$ uniquely determined by the neural network conditioned on
its associated weights $\Phi_{\mathtt{NN}}$.

Overall, the Espaloma framework is end-to-end differentiable—the error in
energy (or the function thereof, such as forces) can be backpropagated to
optimize the force field parameters $\Phi_{\mathtt{FF}}$, and thereby neural
network parameters $\Phi_{\mathtt{NN}}$ that govern how they are produced from
the input molecule. Stage 3 is especially modular and flexible: New force
field terms that act on atoms, bonds, angles, torsions, or combinations
thereof can easily be added and the entire force field refit starting from
either an existing $\Phi_{\mathtt{NN}}$ or training from scratch. In this way,
Espaloma provides a rapid and flexible approach to experimenting with
different potential functions (such as the addition of point polarizability or
exploration of alternative functional forms) or retraining with augmented
training datasets.

##  2 An extensive quantum chemical dataset was curated to provide coverage of
small molecules, proteins, nucleic acids, and other biomolecules

While the initial Espaloma model described in Wang et al., [75] had extremely
limited chemical space coverage, it demonstrated the potential to provide
self-consistent coverage of multiple classes of biomolecules. In Wang et al.,
[75], the quantum chemical training set was limited to 792 small molecules and
1528 di- and tri-peptides, providing a useful proof of concept that was
competitive with other noncommercial small molecule and protein force fields,
but far from providing comprehensive coverage of chemical space relevant to
drug discovery and biomolecular modeling.

To demonstrate that Espaloma can scale to provide high-quality coverage on
large regions of biomolecular space when systematically parameterized on a
large quantum chemical dataset, we constructed a significantly expanded
dataset with the intention of producing a force field that was broadly
applicable to biomolecular modeling and computer-aided drug discovery (Table
1). This quantum chemical dataset is is built from several components that
provide complementary coverage of relevant biomolecular chemistries: From the
foundational SPICE dataset [23], we extracted a large set of druglike small
molecules selected from PubChem [50], dipeptides and their common protonation
and tuatomeric variants, and diverse molecular fragments providing broad
coverage of biomolecules from the DES370K dataset [19]; from the OpenFF 1.x
("Parsley") [63] and 2.x ("Sage") [7] datasets, we extracted optimization and
torsion-drive datasts for diverse small molecules; a diverse set of dipeptide,
tripeptides, disulfide-bridged, bioactive, and cyclic peptides from the
PepConf dataset [62]; a peptide torsion scan set generated by the Open Force
Field Consortium for the OpenFF 3.x ("Rosemary") force field [9]; and a new
set of RNA nucleosides, trinucleotides, and diverse experimental RNA fragments
from the Nucleic Acid Database [13] and RNA Structure Atlas [58] to extend
coverage to this important and growing class of drug targets. Complete details
of the dataset construction and composition are given in SI Section B.

All quantum chemical energies were computed with the Open Force Field (OpenFF)
standard level of quantum chemical theory (B3LYP-D3BJ/DZVP) [63, 7], generated
with the open source psi4 quantum chemistry package [69] using the QCArchive
[68] QCFractal infrastructure via OpenFF QCSubmit [32] workflows. The final
dataset used to construct and assess the espaloma-0.3.0 model contains 17,427
unique molecules and 1,188,317 conformations in total. This dataset is
publicly available for download (see SI Section B).

Dataset Name | Category | Mols | Confs | Split |  | espaloma-0.3.0
---
Energy RMSE (kcal/mol)
Force RMSE (kcal/mol $\cdot$ ${\text{Å}}^{-1}$)
|  Baseline Force Field (Test molecules)
---
Energy RMSE (kcal/mol)
Force RMSE (kcal/mol $\cdot$ ${\text{Å}}^{-1}$)
|  |  |  |  | Train (80%) | Test (10%) | gaff-2.11 | openff-2.0.0 | openff-2.1.0 | ff14SB/RNA.OL3
| SPICE-Pubchem [23, 39]
---
(Dataset)
| Small
---
molecule
14110 | 608436 | 80:10:10 |  | $2.21_{2.20}^{2.22}$
---
$6.18_{6.17}^{6.18}$
| $2.52_{2.48}^{2.56}$
---
$6.87_{6.85}^{6.89}$
| $4.64_{4.60}^{4.69}$
---
$14.38_{14.35}^{14.40}$
| $4.43_{4.38}^{4.47}$
---
$13.94_{13.92}^{13.96}$
| $4.71_{4.66}^{4.75}$
---
$15.53_{15.51}^{15.55}$
| SPICE-DES-Monomers [23, 20]
---
(Dataset)
| Small
---
molecule
369 | 18435 | 80:10:10 |  | $1.50_{1.46}^{1.53}$
---
$5.94_{5.91}^{5.97}$
| $1.68_{1.56}^{1.80}$
---
$6.05_{5.97}^{6.14}$
| $2.12_{1.97}^{2.30}$
---
$10.74_{10.61}^{10.88}$
| $2.63_{2.43}^{2.84}$
---
$12.07_{11.94}^{12.21}$
| $2.71_{2.50}^{2.95}$
---
$12.59_{12.46}^{12.71}$
| Gen2-Opt [1]
---
(OptimizationDataset)
| Small
---
molecule
1024 | 244989 | 80:10:10 |  | $1.66_{1.63}^{1.70}$
---
$3.65_{3.64}^{3.67}$
| $2.05_{1.95}^{2.18}$
---
$4.53_{4.49}^{4.58}$
| $3.36_{3.28}^{3.46}$
---
$10.94_{10.91}^{10.96}$
| $2.68_{2.59}^{2.78}$
---
$10.97_{10.95}^{11.00}$
| $2.77_{2.69}^{2.86}$
---
$11.31_{11.28}^{11.33}$
| Gen2-Torsion [1]
---
(TorsionDriveDataset)
| Small
---
molecule
729 | 25832 | 80:10:10 |  | $2.71_{2.66}^{2.76}$
---
$4.26_{4.24}^{4.28}$
| $2.40_{2.34}^{2.46}$
---
$4.35_{4.33}^{4.37}$
| $2.94_{2.85}^{3.05}$
---
$10.37_{10.32}^{10.41}$
| $2.49_{2.44}^{2.55}$
---
$10.75_{10.71}^{10.80}$
| $2.77_{2.70}^{2.84}$
---
$10.62_{10.58}^{10.66}$
| SPICE-Dipeptide [23]
---
(Dataset)
Peptide | 677 | 26279 | 80:10:10 |  | $3.27_{3.23}^{3.30}$
---
$7.99_{7.97}^{8.00}$
| $3.15_{3.05}^{3.25}$
---
$7.85_{7.80}^{7.90}$
| $4.28_{4.15}^{4.42}$
---
$11.98_{11.92}^{12.05}$
| $4.15_{4.01}^{4.29}$
---
$12.01_{11.95}^{12.06}$
| $4.35_{4.22}^{4.50}$
---
$11.61_{11.56}^{11.66}$
| $4.44_{4.29}^{4.57}$
---
$11.82_{11.77}^{11.88}$
| Pepconf-Opt [2, 62]
---
(OptimizationDataset)
Peptide | 557 | 166291 | 80:10:10 |  | $2.78_{2.73}^{2.84}$
---
$3.50_{3.48}^{3.51}$
| $3.09_{2.95}^{3.27}$
---
$3.76_{3.72}^{3.81}$
| $5.16_{5.05}^{5.27}$
---
$10.15_{10.12}^{10.17}$
| $4.33_{4.19}^{4.48}$
---
$9.61_{9.58}^{9.63}$
| $4.62_{4.50}^{4.75}$
---
$9.79_{9.77}^{9.81}$
| $4.38_{4.28}^{4.49}$
---
$10.13_{10.11}^{10.15}$
| Protein-Torsion
---
(TorsionDriveDataset)
Peptide | 62 | 48999 | 80:10:10 |  | $2.51_{2.47}^{2.55}$
---
$4.08_{4.06}^{4.10}$
| $2.08_{2.01}^{2.14}$
---
$3.62_{3.57}^{3.67}$
| $3.22_{3.08}^{3.35}$
---
$8.01_{7.97}^{8.05}$
| $2.97_{2.85}^{3.10}$
---
$8.77_{8.73}^{8.81}$
| $3.37_{3.22}^{3.54}$
---
$9.02_{8.98}^{9.06}$
| $3.19_{3.06}^{3.31}$
---
$9.29_{9.24}^{9.32}$
| RNA-Diverse
---
(Dataset)
RNA | 64 | 3703 | 80:10:10 |  | $4.22_{4.08}^{4.37}$
---
$4.41_{4.40}^{4.42}$
| $4.06_{3.69}^{4.40}$
---
$4.35_{4.32}^{4.38}$
| $5.72_{5.27}^{6.16}$
---
$16.94_{16.82}^{17.06}$
| $5.76_{5.29}^{6.23}$
---
$18.19_{18.02}^{18.38}$
| $6.44_{5.96}^{6.99}$
---
$19.68_{19.51}^{19.84}$
| $5.97_{5.54}^{6.44}$
---
$19.15_{18.97}^{19.33}$
| RNA-Trinucleotide
---
(Dataset)
RNA | 64 | 35811 | 0:0:100 |  |  | $3.86_{3.77}^{3.96}$
---
$4.24_{4.23}^{4.25}$
| $5.93_{5.88}^{5.99}$
---
$17.09_{17.08}^{17.10}$
| $5.92_{5.87}^{5.97}$
---
$18.70_{18.68}^{18.72}$
| $6.44_{6.38}^{6.49}$
---
$20.03_{20.01}^{20.04}$
| $6.05_{6.00}^{6.11}$
---
$19.70_{19.69}^{19.72}$
| RNA-Nucleoside
---
(Dataset)
RNA | 4 | 9542 | 100:0:0 |  | $1.33_{1.28}^{1.39}$
---
$4.17_{4.15}^{4.18}$
|  | $2.72_{2.67}^{2.77}$
---
$13.80_{13.77}^{13.84}$
| $2.30_{2.26}^{2.34}$
---
$14.02_{13.98}^{14.05}$
| $2.00_{1.97}^{2.03}$
---
$15.41_{15.38}^{15.45}$
| $2.48_{2.44}^{2.51}$
---
$14.18_{14.14}^{14.21}$

Table 1: espaloma-0.3.0 can directly fit quantum chemical potential energies
and forces more accurately than the baseline force fields. Espaloma was fit to
quantum chemical potential energies and forces from various QM datasets
sourced from QCArchive [68], covering a broad chemical space that includes
small molecules, peptides, and RNA molecules (see SI Section B). The entire
dataset consists of 17,427 unique molecules and 1,188,317 conformations. These
datasets were extracted from three different QCArchive workflows:
BasicDataset, OptimizationDataset, and TorsionDriveDataset. Except for RNA-
Trinucleotide and RNA-nucleoside dataset, all datasets were partitioned into
train, validate, and test sets in an 80:10:10 ratio based on molecules. Due to
the lack of chemically diversity, only the RNA-diverse dataset was partitioned
similarly to the other datasets, whereas the entire RNA-nucleoside dataset was
considered during training while the RNA-trinucleotide dataset was used solely
for the test set. The number of molecules and total conformations for each
dataset is annotated in the table. We report the root mean square error (RMSE)
on the training and test sets, along with the performance of other force
fields as baselines on the test set. Note that the reported baseline force
field RMSE for RNA-Nucleoside uses the same train dataset to report
esaploma-0.3.0 RMSE. All statistics are computed with predicted and reference
energies centered to have a zero mean for each molecule similar to the
previous work [75]. The 95% confidence intervals, annotated in the results,
were calculated by bootstrapping molecular conformers with replacement using
1000 replicates.

##  3 espaloma-0.3.0 accurately predicts quantum chemical energies and forces

Leveraging the datasets generated in Section 2, we refit the Espaloma force
field. Similar to the original implementation, since the van der Waals
parameters affect the physical property prediction, which is computationally
challenging to optimize, we only optimize the valence parameters while using
openff-2.0.0 [6]—we leave the tuning of van der Waals parameters for future
study. For partial charges, following the protocol of Wang et al., 2019b [76],
we predict the electronegativity and hardness of atoms before using charge
equilibration [27] to predict atomic partial charges while preserving the
total charge exactly. We improve the training pipeline of the original
Espaloma force field to increase its stability and data efficiency:

  * •

Quantum mechanics forces are incorporated into training.

  * •

Strong L2 regularization is applied on proper and improper torsion force
constants.

  * •

We truncate the improper torsion terms to include only $n=1,2$ periodicities.

To generate the espaloma-0.3.0 force field, we randomly shuffle the datasets
and split each into train, test, and validation sets (80%, 10%, and 10%,
respectively) based on chemical species (unique isomeric SMILES strings).
Since the MM force field is incapable of reproducing quantum chemical heats of
formation, which are reflected as an additive offset in quantum chemical
energy targets, we shift the reference energy to have zero mean; note that
when deployed, the absolute value of MM energy is not physically meaningful
and traditional MM force fields are never used to simulate bond-breaking
events.

As shown in Table 1, espaloma-0.3.0 significantly outperforms all baseline
force fields in predicting QM energies and forces, demonstrating the ability
of espaloma-0.3.0 to recapitulate the QM energy landscape, while legacy force
fields widely popular in the field of MD modeling yield considerable energy
errors and huge force errors (on average twice to thrice of that of
espaloma-0.3.0) vis-à-vis QM calculations. The performance superiority holds
true across diverse chemical landscapes (as discussed in Section 2)—small
molecules and biopolymers and conformational landscape—MD trajectories,
optimization trajectories, and torsion scans, suggesting general utility of
espaloma-0.3.0 in a wide array of chemical and biochemical modeling tasks, as
evidenced in Section 5.

##  4 espaloma-0.3.0 preserves quantum mechanics energy minima

Figure 2:  Espaloma-derived molecular mechanics force fields preserve the
location of quantum chemical minima. Kernel density estimate (KDE) plot
(upper) and empirical cumulative density function (CDF) of aligned root mean
square deviation (RMSD) between quantum chemical minima and MM-optimized
minima using parameters from either espaloma-generated molecular mechanics or
legacy force fields. This qualitative experiment shows that the energy minima
determined by espaloma-0.3.0 is close to the QM minima.

We examined whether the espaloma-0.3.0 force field described in Section 3 not
only quantitatively reproduced quantum chemical equilibrium conformational
energetics, but was able to also qualitatively preserve quantum chemical local
minima. To assess this, we initiated MM energy minimizations from QM-optimized
conformations from the FDA Optimization Dataset in QCArchive [32].
Specifically, we selected the minimum quantum chemical energy snapshot for
each molecule in the dataset and minimized with espaloma- derived force field,
openff-2.0.0 [6], openff-2.1.0 [4], and gaff-2.11 MM force fields using an
L-BFGS optimizer with a 0.1 kJ/mol convergence tolerance. Figure 2 shows a
kernel density estimate (KDE, upper) and an empirical cumulative density
function (CDF, lower) of the aligned rooted mean square distance (RMSD) from
the quantum chemical minimum produced by each force field. We noticed that
espaloma-0.3.0 preserves quantum mechanics energy minima better (“not
optimized away”) compared to legacy force fields.

Figure 3:  espaloma-0.3.0 can be used for accurate protein-ligand alchemical
free energy calculations. Protein-ligand (PL) alchemical free energy
calculations were calculated for TYK2 (10 ns), CDK2 (10 ns), MCL1 (15 ns), P38
(20 ns) using a curated PL-benchmark dataset which comprises 76 ligands in
total. Here, we used Perses 0.10.1 relative free energy calculation
infrastructure [65], based on OpenMM 8.0.0 [24], to assess the accuracy of
espaloma-0.3.0 and openff-2.1.0 for comparison. (a) The openff-2.1.0 with
protein parametrized with Amber ff14SB force field [53] achieves an absolute
free energy ($\Delta G$) RMSE of 1.01 [95% CI: 0.73, 1.33] kcal/mol. (b) The
espaloma-0.3.0 for predicting valence parameters and partial charges of small
molecules combined with Amber ff14SB force field for proteins achieves an
absolute free energy ($\Delta G$) RMSE of 1.13 [95% CI: 0.86, 1.47] kcal/mol.
(c) Parametrizing small molecule and protein self-consistently with
espaloma-0.3.0 achieves absolute free energy ($\Delta G$) RMSE of 1.02 [95%
CI: 0.74, 1.37] kcal/mol which is comparable to those obtained by (a) and (b).
All systems were solvated with TIP3P water and neutralized with 300 mM NaCl
salt using Joung and Cheatham monovalent counterions. The light and dark gray
regions depict the -0.5 kcal/mol and 1.0 kcal/mol confidence bounds.

|  |  |  |  |  | ff14SB+openff-2.1.0
---|---|---|---|---|---|---
|  |  |  |  |  | Relative | Absolute
System | PDB ID | Compds | Edges |  $\Delta\Delta$G(kcal/mol) | ns/replica | RMSE | MUE | RMSE | MUE | R2 | Spearman $\rho$
Tyk2 | 4GIH | 13 | 12 | 3.47 | 10 | $0.54_{0.36}^{0.71}$ | $0.45_{0.28}^{0.62}$ | $0.50_{0.36}^{0.64}$ | $0.42_{0.27}^{0.57}$ | $0.80_{0.53}^{0.93}$ | $0.89_{0.75}^{0.96}$
Cdk2 | 1H1Q | 10 | 9 | 2.78 | 10 | $1.43_{1.04}^{1.75}$ | $1.29_{0.80}^{1.67}$ | $0.74_{0.50}^{0.93}$ | $0.63_{0.41}^{0.86}$ | $0.48_{0.13}^{0.85}$ | $0.69_{0.30}^{0.92}$
Mcl1 | 4HW3 | 25 | 24 | 4.19 | 15 | $1.50_{0.83}^{2.12}$ | $1.02_{0.63}^{1.55}$ | $1.36_{0.77}^{2.01}$ | $0.97_{0.66}^{1.41}$ | $0.50_{0.35}^{0.73}$ | $0.71_{0.57}^{0.86}$
P38 | 3FLY | 28 | 27 | 3.81 | 20 | $1.06_{0.81}^{1.30}$ | $0.87_{0.65}^{1.09}$ | $0.90_{0.60}^{1.19}$ | $0.69_{0.50}^{0.92}$ | $0.57_{0.38}^{0.78}$ | $0.76_{0.63}^{0.89}$
|  |  |  |  |  | ff14SB+espaloma-0.3.0
|  |  |  |  |  | Relative | Absolute
System | PDB ID | Compds | Edges |  $\Delta\Delta$G(kcal/mol) | ns/replica | RMSE | MUE | RMSE | MUE | R2 | Spearman $\rho$
Tyk2 | 4GIH | 13 | 12 | 3.47 | 10 | $0.70_{0.34}^{0.98}$ | $0.52_{0.28}^{0.80}$ | $0.48_{0.29}^{0.65}$ | $0.37_{0.23}^{0.55}$ | $0.79_{0.49}^{0.95}$ | $0.89_{0.71}^{0.97}$
Cdk2 | 1H1Q | 10 | 9 | 2.78 | 10 | $1.15_{0.85}^{1.44}$ | $1.05_{0.73}^{1.36}$ | $0.56_{0.32}^{0.74}$ | $0.46_{0.27}^{0.66}$ | $0.63_{0.27}^{0.92}$ | $0.80_{0.53}^{0.96}$
Mcl1 | 4HW3 | 25 | 24 | 4.19 | 15 | $1.38_{0.90}^{1.96}$ | $1.06_{0.76}^{1.44}$ | $1.51_{0.90}^{2.15}$ | $1.08_{0.74}^{1.56}$ | $0.60_{0.42}^{0.80}$ | $0.77_{0.63}^{0.90}$
P38 | 3FLY | 28 | 27 | 3.81 | 20 | $1.03_{0.81}^{1.26}$ | $0.82_{0.59}^{1.05}$ | $1.10_{0.86}^{1.32}$ | $0.88_{0.63}^{1.13}$ | $0.38_{0.11}^{0.64}$ | $0.62_{0.34}^{0.80}$
|  |  |  |  |  | espaloma-0.3.0
|  |  |  |  |  | Relative | Absolute
System | PDB ID | Compds | Edges |  $\Delta\Delta$G(kcal/mol) | ns/replica | RMSE | MUE | RMSE | MUE | R2 | Spearman $\rho$
Tyk2 | 4GIH | 13 | 12 | 3.47 | 10 | $0.67_{0.45}^{0.87}$ | $0.56_{0.35}^{0.76}$ | $0.46_{0.33}^{0.58}$ | $0.40_{0.28}^{0.53}$ | $0.81_{0.64}^{0.94}$ | $0.90_{0.79}^{0.97}$
Cdk2 | 1H1Q | 10 | 9 | 2.78 | 10 | $0.84_{0.58}^{1.05}$ | $0.75_{0.51}^{0.99}$ | $0.63_{0.48}^{0.76}$ | $0.58_{0.41}^{0.74}$ | $0.47_{0.14}^{0.82}$ | $0.68_{0.41}^{0.90}$
Mcl1 | 4HW3 | 25 | 24 | 4.19 | 15 | $1.44_{0.96}^{1.99}$ | $1.10_{0.76}^{1.50}$ | $1.40_{0.78}^{2.09}$ | $1.00_{0.67}^{1.43}$ | $0.56_{0.40}^{0.78}$ | $0.75_{0.63}^{0.88}$
P38 | 3FLY | 28 | 27 | 3.81 | 20 | $1.02_{0.77}^{1.24}$ | $0.79_{0.56}^{1.04}$ | $0.91_{0.68}^{1.13}$ | $0.75_{0.57}^{0.95}$ | $0.47_{0.24}^{0.68}$ | $0.68_{0.49}^{0.82}$

Table 2: Protein-ligand alchemical free energy calculation benchmarks show
espaloma-0.3.0 achieves high accuracy that is competitive to well-established
force fields. Here, we report several different metrics to assess the
performance of the protein-ligand binding benchmark results including root
mean square error (RMSE), mean unsigned error (MUE), the square of the
correlation coefficient (R2), and the Spearman’s rank correlation coefficient
($\rho$) along with 95% CI for each metric. The initial PDB ID, number of
compounds, number of edges, the binding affinity range ($\Delta\Delta$G), and
the simulation time per replica are reported in the table.

##  5 espaloma-0.3.0 supports accurate alchemical free energy calculations.

To assess whether the small molecule parameters and charges generated by
espaloma-0.3.0 achieves competitive performance with a well-established force
field, openff-2.1.0 [4], we used the perses 0.10.1 relative alchemical free
energy calculation infrastructure [65] to compare the performance on a curated
protein-ligand benchmark dataset. We work with four well-studied protein-
ligand binding systems:

  * •

Tyk2 (PDB: 4GIH) [51], a non-receptor tyrosine-protine kinase, has therapeutic
significance in inflammatory bowel diseases (IBD); this particularly popular
system has good convergence;

  * •

Cdk2 (PDB: 1H1Q) [16], a cyclin-dependent kinase, is involved in molecular
pathology of cancer and is, therefore, a popular target for structure-based
drug design; we use this system, complex with cyclin A, to test the capability
of parametrizing multiple protein subunits;

  * •

P38 (PDB: 3FLY) [43] is a mitogen-activated protein (MAP) kinase which is a
central componenet in signaling networks in mammalian cell types; this target
is another well-studied system, but is expected to be more challenging
compared to Tyk2 and Cdk2 because of the larger ligand transformations and
exploration of structure-activity relationships with multiple R-groups from
different scaffold positions;

  * •

Mcl1 (PDB: 4HW3) [25] (myeloid cell leukemia 1) is a member of the Bcl-2
family of proteins, is overexpressed in various cancers and promotes aberrant
survival of tumor cells; this target entails all ligands with a net charge of
-1 and includes scaffold hopping, thus chosen to test the capability of
simulating free energy calculations for charged ligands and scaffold hopping;

These target systems have been adopted from Open Force Field [11] protein-
ligand benchmark where the protein structures and ligand poses, as well as the
ligand transformation networks were manually curated (see SI Section D).
Within each system, we benchmark two approaches of parametrization:

  * •

Ligand-only. We first apply espaloma-0.3.0 to only the ligand region of a
protein-ligand system, to provide head-to-head comparison with a well-
established small molecule force field openff-2.1.0 [4]; the protein
parameters were assigned by Amber ff14SB force field [53].

  * •

Protein-ligand self-consistently. We apply espaloma-0.3.0 to both the ligand
as well as the protein region of the system; this is to test the capability of
espaloma-0.3.0 to entirely replace the force field parametrization pipeline.

All systems were solvated with TIP3P water [35] and neutralized with the Joung
and Cheatham monovalent counterions [36].

In Figure 3 and Table 2, we illustrate that espaloma-0.3.0 has comparable
protein-ligand binding free energy performance with the combination of Amber
ff14SB [53] and openff-2.1.0 [4] force fields—the reported error and
correlation statistics have overlapping confidence interval. In almost all
cases, the performance of using espaloma-0.3.0 on both the protein and ligand
system are slightly higher than using espaloma-0.3.0 on ligand only and
resorting to Amber ff14SB [53] for protein parametrization.

We also compared the protein-ligand binding free energy calculations using the
first generation espaloma model, espaloma-0.2.2 [75], which was trained on a
limited quantum chemical dataset compared to espaloma-0.3.0. The system was
prepared similarily to those described above. In SI Figure 4, espaloma-0.2.2
significantly underperforms compared to espaloma-0.3.0 for the Cdk2 system due
a large outlier. espaloma-0.2.2 also demonstrates lesser performance on the
Tyk2 system, as illustrated in SI Figure 5. Importantly, the protein-ligand
binding free energy calculations were unstable for Mcl1 and P38, with many of
the ligand transformations being suspended during the simulation. These
results indicate that espaloma-0.3.0, trained on an extensive quantum chemical
dataset and with an improved training strategy, has resulted in the
development of a superior espaloma model.

##  6 Conclusion

Here, we demonstrate that the new espaloma-0.3.0 force field can capture the
quantitative and qualitative behavior of quantum mechanics on a wide variety
of chemical landscapes. As a result, it produces accurate absolute and
relative binding free energies estimates. We hope this work will lay the
foundations to inspire the design of new generations of machine learning-
empowered molecular mechanics force fields.

###  6.1 Chemically and conformationally diverse quantum chemical datasets
curated to construct biomolecular MM force fields

Though Wang et al., [75] laid the foundation of machine learning molecular
mechanics (MM) force field construction, the example force field released with
its implementation is demonstrative rather than impactful because of the
limited diversity it was trained out. In this paper, we tackle this problem
squarely by assembling a high-quality dataset covering chemical spaces and
conformational regions of interest to biomolecular simulation. We release this
dataset with our implementation in the hope that this will encourage the
community to further optimize MM force fields.

###  6.2 espaloma-0.3.0 recapitulates the energy landscape of quantum
mechanics

It has been shown, in this paper in Table 1 and previously [11] that legacy
force fields typically have considerable disagreement with quantum mechanics
(QM) calculations in terms of the force and energy predictions. With carefully
crafted training and regularization strategies, we show that espaloma-0.3.0 is
much closer to the QM energy landscapes and qualitatively behave similarly to
QM as well.

###  6.3 espaloma-0.3.0 accurately predicts protein-ligand binding free energy
with further improvement possible via direct optimization with regard to
experimental free energy

The field of computational chemistry has relied on decades-old force fields
(with limited continuous evolvement) for MD simulation and free energy
calculations. In this paper, we show, by merely fitting quantum mechanics
data, we can curate a force field just as powerful in a few GPU hours, as
opposed to countless engineer years. For future study, we plan to further
optimize force fields under the Espaloma framework directly incorporating
experimental binding free energy targets.

##  7 Disclosures

J.D.C. is a current member of the Scientific Advisory Board of OpenEye
Scientific Software, Redesign Science, Ventus Therapeutics, and Interline
Therapeutics, and has equity interests in Redesign Science and Interline
Therapeutics. The Chodera laboratory receives or has received funding from
multiple sources, including the National Institutes of Health, the National
Science Foundation, the Parker Institute for Cancer Immunotherapy, Relay
Therapeutics, Entasis Therapeutics, Silicon Therapeutics, EMD Serono (Merck
KGaA), AstraZeneca, Vir Biotechnology, Bayer, XtalPi, Interline Therapeutics,
the Molecular Sciences Software Institute, the Starr Cancer Consortium, the
Open Force Field Consortium, Cycle for Survival, a Louis V. Gerstner Young
Investigator Award, and the Sloan Kettering Institute. A complete funding
history for the Chodera lab can be found at http://choderalab.org/funding.
Y.W. has limited financial interests in Flagship Pioneering, Inc.

##  8 Acknowledgements

The authors thank OpenEye Scientific Software for providing a free academic
license to the OpenEye Toolkits. Y.W. acknolwedges support from the Schmidt
Science Fellowship, in partnership with the Rhodes Trust, and the Simons
Center for Computational Physical Chemistry at New York University. J.D.C.
acknowledges support from NIH grant P30 CA008748, NIH grant R01 GM132386, NIH
grant R01 GM121505, and the Sloan Kettering Institute.

## References

  * gen, [2020] (2020).  Openff sandbox gen2 optimization dataset.
  * pep, [2020] (2020).  Pepconf dataset.
  * Battaglia et al., [2018] Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., et al. (2018).  Relational inductive biases, deep learning, and graph networks.  arXiv preprint arXiv:1806.01261.
  * Behara et al., [2023] Behara, P. K., Gokey, T., Cavender, C., Horton, J., Wang, L., Jang, H., Wagner, J., Cole, D., Bayly, C., and Mobley, D. (2023).  openforcefield/openff-forcefields.
  * Bender et al., [2021] Bender, B. J., Gahbauer, S., Luttens, A., Lyu, J., Webb, C. M., Stein, R. M., Fink, E. A., Balius, T. E., Carlsson, J., Irwin, J. J., et al. (2021).  A practical guide to large-scale docking.  Nature protocols, 16(10):4799–4832.
  * Boothroyd et al., [2022] Boothroyd, S., Behara, P. K., Madin, O., Hahn, D., Jang, H., Gapsys, V., Wagner, J., Horton, J., Dotson, D., Thompson, M., et al. (2022).  Development and benchmarking of open force field 2.0.0—the sage small molecule force field.
  * Boothroyd et al., [2023] Boothroyd, S., Behara, P. K., Madin, O. C., Hahn, D. F., Jang, H., Gapsys, V., Wagner, J. R., Horton, J. T., Dotson, D. L., Thompson, M. W., et al. (2023).  Development and benchmarking of open force field 2.0.0: The Sage small molecule force field.  Journal of Chemical Theory and Computation.
  * Case et al., [2023] Case, D., Aktulga, H., Belfon, K., Ben-Shalom, I., Berryman, J., Brozell, S., Cerutti, D., Cheatham, T., III, Cisneros, G., Cruzeiro, V., Darden, T., Forouzesh, N., Giambasu, G., Giese, T., Gilson, M., Gohlke, H., Goetz, A., Harris, J., Izadi, S., Izmailov, S., Kasavajhala, K., Kaymak, M., King, E., Kovalenko, A., Kurtzman, T., Lee, T., Li, P., Lin, C., Liu, J., Luchko, T., Luo, R., Machado, M., Man, V., Manathunga, M., Merz, K., Miao, Y., Mikhailovskii, O., Monard, G., Nguyen, H., O’Hearn, K., Onufriev, A., Pan, F., Pantano, S., Qi, R., Rahnamoun, A., Roe, D., Roitberg, A., Sagui, C., Schott-Verdugo, S., Shajan, A., Shen, J., Simmerling, C., Skrynnikov, N., Smith, J., Swails, J., Walker, R., Wang, J., Wang, J., Wei, H., Wu, X., Wu, Y., Xiong, Y., Xue, Y., York, D., Zhao, S., Zhu, Q., , and (2023), P. K. (2023).  Amber 2023.
  * Cavender et al., [2023] Cavender, C. E., Behara, P. K., Boothroyd, S., Dotson, D. L., Horton, J. T., Mitchell, J. A., Pulido, I. J., Thompson, M. W., Wagner, J., Wang, L., Chodera, J. D., Cole, D. J., Mobley, D. L., Shirts, M. R., and Gilson, M. K. (2023).  Development and benchmarking of an open, self- consistent force field for proteins and small molecules from the Open Force Field Initiative.  Funding from National Institutes of Health grant GM132386.
  * Chodera et al., [2021] Chodera, J., Rizzi, A., Naden, L., Beauchamp, K., Grinaway, P., Fass, J., Wade, A., Rustenburg, B., Ross, G. A., Krämer, A., Macdonald, H. B., Rodríguez-Guerra, J., dominicrufa, Simmonett, A., Swenson, D. W. H., hb0402, Henry, M., Roet, S., and Silveira, A. (2021).  Choderalab/openmmtools: 0.20.3 Bugfix Release.  Zenodo.
  * Chodera et al., [2020] Chodera, J., Wiewiora, R., Stern, C., and peastman (2020).  openmm/openmm-forcefields: Fix GAFF AM1-BCC charging bug for some molecules.
  * Chodera and Shirts, [2011] Chodera, J. D. and Shirts, M. R. (2011).  Replica exchange and expanded ensemble simulations as gibbs sampling: Simple improvements for enhanced mixing.  The Journal of chemical physics, 135(19):194110.
  * Coimbatore Narayanan et al., [2014] Coimbatore Narayanan, B., Westbrook, J., Ghosh, S., Petrov, A. I., Sweeney, B., Zirbel, C. L., Leontis, N. B., and Berman, H. M. (2014).  The nucleic acid database: new features and capabilities.  Nucleic acids research, 42(D1):D114–D122.
  * Coutsias et al., [2016] Coutsias, E. A., Lexa, K. W., Wester, M. J., Pollock, S. N., and Jacobson, M. P. (2016).  Exhaustive conformational sampling of complex fused ring macrocycles using inverse kinematics.  Journal of chemical theory and computation, 12(9):4674–4687.
  * Dauber-Osguthorpe and Hagler, [2019] Dauber-Osguthorpe, P. and Hagler, A. T. (2019).  Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there?  Journal of computer-aided molecular design, 33(2):133–203.
  * Davies et al., [2002] Davies, T. G., Bentley, J., Arris, C. E., Boyle, F. T., Curtin, N. J., Endicott, J. A., Gibson, A. E., Golding, B. T., Griffin, R. J., Hardcastle, I. R., et al. (2002).  Structure-based design of a potent purine-based cyclin-dependent kinase inhibitor.  Nature structural biology, 9(10):745–749.
  * DeMarco and Woods, [2009] DeMarco, M. L. and Woods, R. J. (2009).  Atomic-resolution conformational analysis of the gm3 ganglioside in a lipid bilayer and its implications for ganglioside–protein recognition at membrane surfaces.  Glycobiology, 19(4):344–355.
  * DeMarco et al., [2010] DeMarco, M. L., Woods, R. J., Prestegard, J. H., and Tian, F. (2010).  Presentation of membrane-anchored glycosphingolipids determined from molecular dynamics simulations and nmr paramagnetic relaxation rate enhancement.  Journal of the American Chemical Society, 132(4):1334–1338.
  * [19] Donchev, A. G., Taube, A. G., Decolvenaere, E., Hargus, C., McGibbon, R. T., Law, K.-H., Gregersen, B. A., Li, J.-L., Palmo, K., Siva, K., et al. (2021a).  Quantum chemical benchmark databases of gold-standard dimer interaction energies.  Scientific data, 8(1):55.
  * [20] Donchev, A. G., Taube, A. G., Decolvenaere, E., Hargus, C., McGibbon, R. T., Law, K.-H., Gregersen, B. A., Li, J.-L., Palmo, K., Siva, K., et al. (2021b).  Quantum chemical benchmark databases of gold-standard dimer interaction energies.  Scientific data, 8(1):55.
  * Du et al., [2018] Du, J., Zhang, S., Wu, G., Moura, J. M. F., and Kar, S. (2018).  Topology Adaptive Graph Convolutional Networks.  arXiv:1710.10370 [cs, stat].
  * Duvenaud et al., [2015] Duvenaud, D. K., Maclaurin, D., Iparraguirre, J., Bombarell, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P. (2015).  Convolutional networks on graphs for learning molecular fingerprints.  In Advances in neural information processing systems, pages 2224–2232.
  * Eastman et al., [2023] Eastman, P., Behara, P. K., Dotson, D. L., Galvelis, R., Herr, J. E., Horton, J. T., Mao, Y., Chodera, J. D., Pritchard, B. P., Wang, Y., et al. (2023).  Spice, a dataset of drug-like molecules and peptides for training machine learning potentials.  Scientific Data, 10(1):11.
  * Eastman et al., [2017] Eastman, P., Swails, J., Chodera, J. D., McGibbon, R. T., Zhao, Y., Beauchamp, K. A., Wang, L.-P., Simmonett, A. C., Harrigan, M. P., Stern, C. D., et al. (2017).  Openmm 7: Rapid development of high performance algorithms for molecular dynamics.  PLoS computational biology, 13(7):e1005659.
  * Friberg et al., [2013] Friberg, A., Vigil, D., Zhao, B., Daniels, R. N., Burke, J. P., Garcia-Barrantes, P. M., Camper, D., Chauder, B. A., Lee, T., Olejniczak, E. T., et al. (2013).  Discovery of potent myeloid cell leukemia 1 (mcl-1) inhibitors using fragment-based methods and structure-based design.  Journal of medicinal chemistry, 56(1):15–30.
  * Gilmer et al., [2017] Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. (2017).  Neural message passing for quantum chemistry.  In International conference on machine learning, pages 1263–1272. PMLR.
  * Gilson et al., [2003] Gilson, M. K., Gilson, H. S., and Potter, M. J. (2003).  Fast assignment of accurate partial atomic charges: an electronegativity equalization method that accounts for alternate resonance forms.  Journal of chemical information and computer sciences, 43(6):1982–1997.
  * Gould et al., [2018] Gould, I., Skjevik, A., Dickson, C., Madej, B., and Walker, R. (2018).  Lipid17: A comprehensive amber force field for the simulation of zwitterionic and anionic lipids.  Manuscript in preparation.
  * Hagler, [2019] Hagler, A. T. (2019).  Force field development phase ii: Relaxation of physics-based criteria… or inclusion of more rigorous physics into the representation of molecular energetics.  Journal of computer-aided molecular design, 33(2):205–264.
  * Hamilton et al., [2017] Hamilton, W., Ying, Z., and Leskovec, J. (2017).  Inductive representation learning on large graphs.  In Advances in neural information processing systems, pages 1024–1034.
  * Horn et al., [2004] Horn, H. W., Swope, W. C., Pitera, J. W., Madura, J. D., Dick, T. J., Hura, G. L., and Head-Gordon, T. (2004).  Development of an improved four-site water model for biomolecular simulations: Tip4p-ew.  The Journal of chemical physics, 120(20):9665–9678.
  * Horton, [2022] Horton, J. (2022).  openforcefield/openff-qcsubmit: 0.3.1.
  * Jakalian et al., [2000] Jakalian, A., Bush, B. L., Jack, D. B., and Bayly, C. I. (2000).  Fast, efficient generation of high-quality atomic charges. am1-bcc model: I. method.  Journal of computational chemistry, 21(2):132–146.
  * Jakalian et al., [2002] Jakalian, A., Jack, D. B., and Bayly, C. I. (2002).  Fast, efficient generation of high-quality atomic charges. am1-bcc model: Ii. parameterization and validation.  Journal of computational chemistry, 23(16):1623–1641.
  * Jorgensen et al., [1983] Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983).  Comparison of simple potential functions for simulating liquid water.  The Journal of chemical physics, 79(2):926–935.
  * Joung and Cheatham III, [2008] Joung, I. S. and Cheatham III, T. E. (2008).  Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations.  The journal of physical chemistry B, 112(30):9020–9041.
  * Joung and Cheatham III, [2009] Joung, I. S. and Cheatham III, T. E. (2009).  Molecular dynamics simulations of the dynamic and energetic properties of alkali and halide ions using water-model-specific ion parameters.  The Journal of Physical Chemistry B, 113(40):13279–13290.
  * Khoury et al., [2013] Khoury, G. A., Thompson, J. P., Smadbeck, J., Kieslich, C. A., and Floudas, C. A. (2013).  Forcefield_ptm: Ab initio charge and amber forcefield parameters for frequently occurring post-translational modifications.  Journal of chemical theory and computation, 9(12):5653–5674.
  * Kim et al., [2023] Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., Yu, B., et al. (2023).  Pubchem 2023 update.  Nucleic Acids Research, 51(D1):D1373–D1380.
  * Kingma and Ba, [2014] Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kipf and Welling, [2016] Kipf, T. N. and Welling, M. (2016).  Semi-supervised classification with graph convolutional networks.  CoRR, abs/1609.02907.
  * Kirschner et al., [2008] Kirschner, K. N., Yongye, A. B., Tschampel, S. M., González-Outeiriño, J., Daniels, C. R., Foley, B. L., and Woods, R. J. (2008).  Glycam06: a generalizable biomolecular force field. carbohydrates.  Journal of computational chemistry, 29(4):622–655.
  * Labute and Ebert, [2021] Labute, P. and Ebert, M. (2021).  Optimizing simulations protocols for relative free energy calculations.  In Free Energy Methods in Drug Discovery: Current State and Future Directions, pages 227–245. ACS Publications.
  * Landrum et al., [2023] Landrum, G., Tosco, P., Kelleyand, B., Ric, Cosgrove, D., sriniker, gedeck, Vianello, R., NadineSchneider, Kawashima, E., N, D., Jones, G., Dalke, A., Cole, B., Swain, M., Turk, S., AlexanderSavelyev, Vaucher, A., Wójcikowski, M., Take, I., Probst, D., Ujihara, K., Scalfani, V. F., guillaume godin, Lehtivarjo, J., Pahl, A., Walker, R., Berenger, F., jasondbiggs, and strets123 (2023).  rdkit/rdkit: 2023_03_2 (Q1 2023) Release.
  * Leach, [2001] Leach, A. R. (2001).  Molecular modelling: principles and applications.  Pearson education.
  * Leimkuhler and Matthews, [2016] Leimkuhler, B. and Matthews, C. (2016).  Efficient molecular dynamics using geodesic integration and solvent–solute splitting.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2189):20160138.
  * Li and Merz Jr, [2014] Li, P. and Merz Jr, K. M. (2014).  Taking into account the ion-induced dipole interaction in the nonbonded model of ions.  Journal of chemical theory and computation, 10(1):289–297.
  * Li et al., [2013] Li, P., Roberts, B. P., Chakravorty, D. K., and Merz Jr, K. M. (2013).  Rational design of particle mesh ewald compatible lennard-jones parameters for+ 2 metal cations in explicit solvent.  Journal of chemical theory and computation, 9(6):2733–2748.
  * Li et al., [2015] Li, P., Song, L. F., and Merz Jr, K. M. (2015).  Parameterization of highly charged metal ions using the 12-6-4 lj-type nonbonded model in explicit water.  The Journal of Physical Chemistry B, 119(3):883–895.
  * Li et al., [2010] Li, Q., Cheng, T., Wang, Y., and Bryant, S. H. (2010).  Pubchem as a public resource for drug discovery.  Drug discovery today, 15(23-24):1052–1057.
  * Liang et al., [2013] Liang, J., Tsui, V., Van Abbema, A., Bao, L., Barrett, K., Beresini, M., Berezhkovskiy, L., Blair, W. S., Chang, C., Driscoll, J., et al. (2013).  Lead identification of novel and selective tyk2 inhibitors.  European journal of medicinal chemistry, 67:175–187.
  * Macdonald et al., [2022] Macdonald, H. B., dfhahn, Henry, M., Chodera, J., Dotson, D., Glass, W., and Pulido, I. (2022).  openforcefield/openff-arsenic: v0.2.1.
  * Maier et al., [2015] Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015).  ff14sb: improving the accuracy of protein side chain and backbone parameters from ff99sb.  Journal of chemical theory and computation, 11(8):3696–3713.
  * Mey et al., [2020] Mey, A. S., Allen, B. K., Macdonald, H. E. B., Chodera, J. D., Hahn, D. F., Kuhn, M., Michel, J., Mobley, D. L., Naden, L. N., Prasad, S., et al. (2020).  Best practices for alchemical free energy calculations [article v1. 0].  Living journal of computational molecular science, 2(1).
  * [55] Mobley, D. L., Bannan, C. C., Rizzi, A., Bayly, C. I., Chodera, J. D., Lim, V. T., Lim, N. M., Beauchamp, K. A., Shirts, M. R., Gilson, M. K., et al. (2018a).  Open force field consortium: Escaping atom types using direct chemical perception with smirnoff v0. 1.  BioRxiv, page 286542.
  * [56] Mobley, D. L., Bannan, C. C., Rizzi, A., Bayly, C. I., Chodera, J. D., Lim, V. T., Lim, N. M., Beauchamp, K. A., Slochower, D. R., Shirts, M. R., et al. (2018b).  Escaping atom types in force fields using direct chemical perception.  Journal of chemical theory and computation, 14(11):6076–6092.
  * Murphy et al., [2018] Murphy, R. L., Srinivasan, B., Rao, V. A., and Ribeiro, B. (2018).  Janossy pooling: Learning deep permutation-invariant functions for variable-size inputs.  CoRR, abs/1811.01900.
  * Parlea et al., [2016] Parlea, L. G., Sweeney, B. A., Hosseini-Asanjan, M., Zirbel, C. L., and Leontis, N. B. (2016).  The rna 3d motif atlas: Computational methods for extraction, organization and evaluation of rna motifs.  Methods, 103:99–119.
  * Paszke et al., [2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. (2019).  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc.
  * Pérez et al., [2007] Pérez, A., Marchán, I., Svozil, D., Sponer, J., Cheatham III, T. E., Laughton, C. A., and Orozco, M. (2007).  Refinement of the amber force field for nucleic acids: improving the description of $\alpha$/$\gamma$ conformers.  Biophysical journal, 92(11):3817–3829.
  * Prasad et al., [2018] Prasad, S., Mobley, D., Braun, E., Mayes, H., Monroe, J., Zuckerman, D., et al. (2018).  Best practices for foundations in molecular simulations [article v1. 0].  Living Journal of Computational Molecular Science, 1:1–28.
  * Prasad et al., [2019] Prasad, V. K., Otero-de La-Roza, A., and DiLabio, G. A. (2019).  Pepconf, a diverse data set of peptide conformational energies.  Scientific data, 6(1):1–9.
  * Qiu et al., [2021] Qiu, Y., Smith, D. G., Boothroyd, S., Jang, H., Hahn, D. F., Wagner, J., Bannan, C. C., Gokey, T., Lim, V. T., Stern, C. D., et al. (2021).  Development and benchmarking of open force field v1.0.0—the Parsley small-molecule force field.  Journal of chemical theory and computation, 17(10):6262–6280.
  * RDKit, online, [2013] RDKit, online (2013).  RDKit: Open-source cheminformatics.  http://www.rdkit.org.  [Online; accessed 11-April-2013].
  * Rufa et al., [2022] Rufa, D. A., Zhang, I., Bruce Macdonald, H. E., Grinaway, P. B., Pulido, I., Henry, M. M., Rodríguez-Guerra, J., Wittmann, M., Albanese, S. K., Glass, W. G., Silveira, A., Schaller, D., Naden, L. N., and Chodera, J. D. (2022).  Perses.  Please cite this software using these metadata.
  * Schlick, [2010] Schlick, T. (2010).  Molecular modeling and simulation: an interdisciplinary guide, volume 2.  Springer.
  * Shirts and Chodera, [2008] Shirts, M. R. and Chodera, J. D. (2008).  Statistically optimal analysis of samples from multiple equilibrium states.  The Journal of chemical physics, 129(12):124105.
  * Smith et al., [2021] Smith, D. G., Altarawy, D., Burns, L. A., Welborn, M., Naden, L. N., Ward, L., Ellis, S., Pritchard, B. P., and Crawford, T. D. (2021).  The molssi qcarchive project: An open-source platform to compute, organize, and share quantum chemistry data.  Wiley Interdisciplinary Reviews: Computational Molecular Science, 11(2):e1491.
  * Smith et al., [2020] Smith, D. G., Burns, L. A., Simmonett, A. C., Parrish, R. M., Schieber, M. C., Galvelis, R., Kraus, P., Kruse, H., Di Remigio, R., Alenaizan, A., et al. (2020).  Psi4 1.4: Open-source software for high-throughput quantum chemistry.  The Journal of chemical physics, 152(18).
  * Tan et al., [2018] Tan, D., Piana, S., Dirks, R. M., and Shaw, D. E. (2018).  Rna force field with accuracy comparable to state-of-the-art protein force fields.  Proceedings of the National Academy of Sciences, 115(7):E1346–E1355.
  * Tse et al., [2019] Tse, C. H., Comer, J., Sang Chu, S. K., Wang, Y., and Chipot, C. (2019).  Affordable membrane permeability calculations: permeation of short-chain alcohols through pure-lipid bilayers and a mammalian cell membrane.  Journal of Chemical Theory and Computation, 15(5):2913–2924.
  * Wang et al., [2006] Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006).  Automatic atom type and bond type perception in molecular mechanical calculations.  Journal of molecular graphics and modelling, 25(2):247–260.
  * Wang et al., [2004] Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004).  Development and testing of a general amber force field.  Journal of computational chemistry, 25(9):1157–1174.
  * [74] Wang, M., Zheng, D., Ye, Z., Gan, Q., Li, M., Song, X., Zhou, J., Ma, C., Yu, L., Gai, Y., et al. (2019a).  Deep graph library: A graph-centric, highly-performant package for graph neural networks.  arXiv preprint arXiv:1909.01315.
  * Wang et al., [2022] Wang, Y., Fass, J., Kaminow, B., Herr, J. E., Rufa, D., Zhang, I., Pulido, I., Henry, M., Bruce Macdonald, H. E., Takaba, K., and Chodera, J. D. (2022).  End-to-end differentiable construction of molecular mechanics force fields.  Chem. Sci., 13:12016–12033.
  * [76] Wang, Y., Fass, J., Stern, C. D., Luo, K., and Chodera, J. (2019b).  Graph nets for partial charge prediction.  arXiv preprint arXiv:1909.07903.
  * Wang et al., [2023] Wang, Y., Pulido, I., Takaba, K., Kaminow, B., Scheen, J., Wang, L., and Chodera, J. D. (2023).  Espalomacharge: Machine learning-enabled ultra-fast partial charge assignment.  arXiv preprint arXiv:2302.06758.
  * [78] Wang, Y., Sun, Y., Liu, Z., Sarma, S. E., Bronstein, M. M., and Solomon, J. M. (2019c).  Dynamic graph cnn for learning on point clouds.  Acm Transactions On Graphics (tog), 38(5):1–12.
  * Wu et al., [2019] Wu, F., Zhang, T., Souza Jr, A. H. d., Fifty, C., Yu, T., and Weinberger, K. Q. (2019).  Simplifying graph convolutional networks.  arXiv preprint arXiv:1902.07153.
  * Xu, [2019] Xu, H. (2019).  Optimal measurement network of pairwise differences.  Journal of Chemical Information and Modeling, 59(11):4720–4728.
  * Xu et al., [2018] Xu, K., Hu, W., Leskovec, J., and Jegelka, S. (2018).  How powerful are graph neural networks?  arXiv preprint arXiv:1810.00826.
  * Zgarbová et al., [2015] Zgarbová, M., Sponer, J., Otyepka, M., Cheatham III, T. E., Galindo-Murillo, R., and Jurecka, P. (2015).  Refinement of the sugar–phosphate backbone torsion beta for amber force fields improves the description of z-and b-dna.  Journal of chemical theory and computation, 11(12):5723–5736.

Appendix: Espaloma-0.3.0: Machine-learned molecular mechanics force field for
the simulation of protein-ligand systems and beyond.

##  Appendix A Code availability

The data and Python code used to produce the results discussed in this paper
is distributed open-source under an MIT license and is available at
https://github.com/choderalab/espaloma-0.3.0-manuscript.

Core dependencies include a modified version of Espaloma 0.3.0 release [75]
(https://github.com/choderalab/espaloma/tree/4c6155b72d00ce0190b3cb551e7e59f0adc33a56),
PyTorch 1.1.2 [59], Deep Graph Library 0.9.0 [76], and Open Force Field
Toolkit 0.10.6 [55], to refit and evaluate the espaloma model. A modified
version of Openmmforcefields 0.11.0 [11]
(https://github.com/kntkb/openmmforcefields/tree/6d2c3dcd33d9800a32032d28b6b2dca92f348a43)
was used to run all the relative alchemical protein-ligand binding free energy
calculations with Perses 0.10.1 infrastructure [65]. Espaloma 0.2.4 release
and a modified version of Espaloma 0.3.0 was used to parametrize small
molecules with espaloma-0.2.2 and espaloma-0.3.0, respectively. A modified
version of Perses 0.10.1
(https://github.com/kntkb/perses/tree/0d069fc1cf31b8cce1ae7a1482c3fa46bc1382d2)
was used to self-consistently parametrize both small molecules and proteins
with espaloma-0.3.0. A modified version of cinnabar 0.3.0 [52]
(https://github.com/kntkb/cinnabar/tree/de7bc6623fb25d75848aa1c9f538b77cd02a4b01)
was used to support arbitrary tick frequency when plotting $\Delta$G and
$\Delta$$\Delta$G plots.

##  Appendix B MolSII QCArchive QM datasets

The Python code used to download the quantum mechanics (QM) datasets from the
MolSSI QCArchive [32] is available at https://github.com/choderalab/download-
qca-datasets. The QM datasets utilized in this study were obtained from
various workflows implemented in the QCArchive ecosystem, including Dataset,
OptimizationDataset, and TorsionDriveDataset generated at the B3LYP-D3BJ/DZVP
level of theory. This level of theory was chosen to maintain consistency with
the Open Force Field initiative [63], and it is expected to balance
computational efficiency and accuracy in reproducing conformations generated
by higher-level theories.

The QM datasets in Table 1 are composed of the following datasets deposited in
QCArhive and annotated based on their respective categories.

##### Small molecules

  * •

SPICE-Pubchem 222Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-08-QMDataset-pubchem-set1-single-
points 333Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-08-QMDataset-pubchem-set2-single-
points 444Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-09-QMDataset-pubchem-set3-single-
points 555Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-09-QMDataset-pubchem-set4-single-
points 666Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-09-QMDataset-pubchem-set5-single-
points 777Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-09-QMDataset-pubchem-set6-single-
points is a collection of Dataset that contains a comprehensive and diverse
collection of small, drug-like molecules obtained from Pubchem. It includes
atoms within the range of 3 to 50, including hydrogens, and encompasses the
elements of Br, C, Cl, F, H, I, N, O, P, and S.

  * •

SPICE-DES-Monomers 888Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-15-QMDataset-DES-monomers-single-
points is a Dataset, sourced from DES370K [20], consists of small molecules
(up to 22 atoms) chosen to cover a wide range of chemical space, including the
elements of Br, C, Cl, F, H, I, N, O, P, and S.

  * •

Gen2-Opt 999Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-20-OpenFF-Gen-2-Optimization-
Set-1-Roche 101010Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-20-OpenFF-Gen-2-Optimization-
Set-2-Coverage 111111Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-20-OpenFF-Gen-2-Optimization-
Set-3-Pfizer-Discrepancy 121212Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2020-03-20-OpenFF-
Gen-2-Optimization-Set-4-eMolecules-Discrepancy 131313Source:
https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-20-OpenFF-Gen-2-Optimization-
Set-5-Bayer is a collection of OptimizationDataset that contains drug-like
molecules used for the parametrization of the OpenFF 1.2.0 ("Parsley") [63]
small molecule force field developed by the Open Force Field Consortium. This
dataset is part of the datasets used to generate the first generation espaloma
force field, espaloma-0.2.2.

  * •

Gen2-Torsion 141414Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-Set-1-Roche
151515Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-23-OpenFF-Gen-2-Torsion-
Set-1-Roche-2 161616Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-
Set-2-Coverage 171717Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-23-OpenFF-Gen-2-Torsion-
Set-2-Coverage-2 181818Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-
Set-3-Pfizer-Discrepancy 191919Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2020-03-23-OpenFF-Gen-2-Torsion-
Set-3-Pfizer-Discrepancy-2 202020Source:
https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-
Set-4-eMolecules-Discrepancy 212121Source:
https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-23-OpenFF-Gen-2-Torsion-
Set-4-eMolecules-Discrepancy-2 222222Source:
https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-Set-5-Bayer
232323Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-26-OpenFF-Gen-2-Torsion-
Set-5-Bayer-2 242424Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-03-12-OpenFF-Gen-2-Torsion-
Set-6-supplemental 252525Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2020-03-26-OpenFF-Gen-2-Torsion-
Set-6-supplemental-2 is a collection TorsionDriveDataset that contains torsion
scans of drug-like molecules which is part of the dataset used for the
parameterization of the OpenFF 2.0.0 ("Sage") [7] small molecule force field
developed by the Open Force Field Consortium.

##### Peptides

  * •

SPICE-Dipeptide 262626Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-08-QMDataset-Dipeptide-single-
points is a Dataset that contains a broad coverage of the possible dipeptides
capped with ACE and NME groups formed by the 20 natural amino acids and their
common protonation variants. This includes two forms of CYS (neutral or
negatively charged), two forms of GLU (neutral or negatively charged), two
forms of ASP (neutral or negatively charged), two forms of LYS (neutral or
positively charged), and three forms of HIS (neutral forms with a hydrogen on
either ND1 or NE2, and a positively charged form with hydrogens on both).

  * •

Pepconf-Opt 272727Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2020-10-26-PEPCONF-Optimization is a
OptimizationDataset that contains short peptides, including capped, cyclic,
and disulfide-bonded peptides originally sourced from Prasad et al. [62] and
regenerated by the Open Force Field Consortium. In this study, the default-dlc
QC specification was utilized, differing from the one used in the first
generation espaloma force field (espaloma-0.2.2) [75], leading to improved
self-consistent field convergence.

  * •

Protein-torsion 282828Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2021-11-18-OpenFF-Protein-
Dipeptide-2D-TorsionDrive 292929Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2022-02-10-OpenFF-Protein-
Capped-1-mer-Sidechains 303030Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2022-05-30-OpenFF-Protein-
Capped-3-mer-Backbones 313131Source: https://github.com/openforcefield/qca-
dataset-submission/tree/master/submissions/2023-02-06-OpenFF-Protein-
Capped-3-mer-Omega is a collection of TorsionDriveDataset that contains
various torsion scans of polypeptides generated by the Open Force Field
Consortium for the OpenFF 3.x ("Rosemary") force field [9]. These torsion
scans cover $\chi_{1}$ and $\chi_{2}$ angles in the rotatable side chains, as
well as $\phi$, $\psi$, and $\omega$ angles in the backbones.

##### RNA

  * •

RNA-Diverse 323232Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2022-07-07-RNA-basepair-triplebase-single-
points is a Dataset that contains comprehensive and diverse collection of
experimental RNA structures. It includes 138 base pair structures and 295 base
triple structures sourced from the Nucleic Acid Database [13]. Additionally,
the dataset contains 4,056 representative trinucleotide structures obtained
from the RNA Structure Atlas website [58], where the experimentally observed
internal and hairpin loop motifs, as well as junction loops of representative
sets of RNA 3D Structures with an X-ray resolution cutoff of 2.5 Å, were
segmented into all possible trinucleotide permutations, resulting in 64 unique
molecules. These trinucleotide structures are capped with O5’ hydroxyl groups
at the 5’ end and clustered to select the representative structures. For the
espaloma refitting experiment, only the trinucleotides were utilized.

  * •

RNA-Trinucleotide 333333Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2022-10-21-RNA-trinucleotide-single-points
is a Dataset that provides a broader and more diverse structural coverage of
trinucleotides compared to the RNA-Diverse dataset.

  * •

RNA-Nucleoside 343434Source: https://github.com/openforcefield/qca-dataset-
submission/tree/master/submissions/2023-03-09-RNA-nucleoside-single-points is
a Dataset that comprises a comprehensive and diverse collection of nucleosides
(adenosine, guanosine, cytidine, and uridine) without O5’ hydroxyl atoms.
These nucleosides are generated using 500K implicit solvent MD and torsion
scanning on N-glycosidic bond ($\chi$ torsion) that connects the base and
sugar, resulting in diverse sugar pucker conformations and extensive coverage
of $\chi$ torsions.

##  Appendix C Espaloma refitting experiment

The Python code used to refit and evaluate espaloma-0.3.0 is available at
https://github.com/choderalab/refit-espaloma. It should be noted that
espaloma-0.3.0 is no longer compatible with espaloma-0.2.x models and vice
versa.

###  C.1 Data preparation

The QM datasets obtained from the QCArchive [68] in SI Section B were
preprocessed prior to the refitting experiment. Molecules with a gap between
the minimum and maximum energy larger than 0.1 Hartree (62.5 kcal/mol) were
excluded. Since the van der Waals parameters affect the physical property
prediction, which is computationally challenging to optimize, we focus on
optimizing the valence parameters and use openff-2.0.0 force field [7] for the
van der Waals parameters. AM1BCC-ELF10 partial charges were pre-computed using
the OpenEye toolkit as reference charges. These charges were then used to
predict the atomic partial charges based on the predicted electronegativity
and hardness of atoms, following the same protocol described in the earlier
works by Wang et al. [75]. To ensure that each molecule was represented only
once, duplicate molecules across different datasets were merged, ensuring that
unique molecules were distributed among the train, validate, or test dataset.

###  C.2 Machine learning experimental details

####  C.2.1 Input features

One of the improvements made from the previous espaloma model [75] is the
exclusion of resonance-sensitive features, such as valences and formal
charges, in order to improve the handling of molecules with atomic resonance,
such as guanidinium and carboxylic acid. In this study, the input features of
the atoms included the one-hot encoded element, as well as the hybridization,
aromaticity, ring membership of sizes 3 to 8, atom mass, and the degree of the
atoms, which is defined as the number of directly-bonded neighbors, all
assigned using the RDKit 2023-03-4 release package [44].

####  C.2.2 Hyperparameter optimization

The hyperparameters were briefly optimized utilizing a subset of data from SI
Section B, which included OpenFF Gen2-Opt, SPICE-Dipeptide, and RNA-Diverse
datasets. The data was partitioned into train : validate : test sets in a 40 :
30 : 30 ratio. During the training process, energy and force matching were
applied, along with partial charge fitting using the charge equilibrium
approach, as in [75, 77].

$\mathcal{L}=W_{\mathtt{energy}}\mathcal{L}_{\mathtt{energy}}+W_{\mathtt{force}}\mathcal{L}_{\mathtt{force}}+W_{\mathtt{charge}}\mathcal{L}_{\mathtt{charge}}$
(2)

Following the protocol specified in Wang et al. [75], we utilized GraphSAGE
[30] as the graph neural network model, the Adam optimizer [40], and the
Rectified Linear Unit (ReLU) activation function, while maintaining the energy
and charge loss weights to 1 and 1e-3, respectively, throughout the
optimization experiment. The hyperparameters subject to optimization included
the batch size (32, 64, 128, 256), the depth of the graph neural network (2,
3, 4, 5), the depth of the Janossy pooling network (2, 3, 4, 5), the learning
rates (1e-3, 1e-4, 5e-5, 1e-5), the number of units per layer (64, 128, 256,
512), and the force weights (1, 1e-1, 1e-2, 1e-3, 1e-4) via grid search on the
validation set, and trained for 3000 epochs for each optimization experiment.

As a result, the optimal configuration was determined as follows: For the atom
embedding stage (Stage1), three GraphSAGE layers with 512 units and ReLU
activation function were employed. For the symmetry preserving pooling stage
(Stage2) and the readout stage (Stage3), four feed-forward layers with 512
units and ReLU activation, a learning rate of 1e-4, and a force loss weight of
1.

####  C.2.3 Production run

The datasets from SI Section B were partitioned into train, validate, and test
sets with a distribution of 80:10:10 ratio, respectively, with few exceptions.
Notably, the entire RNA-Nucleoside dataset was exclusively utilized for the
train set, while the entire RNA-Trinucleoside dataset was allocated for the
test set. This partitioning scheme was designed to incorporate diverse
molecular structures and enable a comprehensive evaluation of the performance
of the espaloma model.

It should be noted that the espaloma model, trained with the hyperparameters
described above, reproduced poor torsion profiles, which led to out-of-plane
geometries after molecular minimization using the trained espaloma model. We
found that this problem could be remedied by truncating the improper torsion
terms to only $n=1,2$ periodicities, instead of $n=1,...,6$ as in the original
method [75], and by utilizing regularization for the proper and improper
torsion force constants. Additionally, we discovered that increasing the loss
weight on the charge term improved the prediction of partial charges and
better recapitulated the electrostatic potentials.

Regarding these findings, the final espaloma model was trained with the
following loss function with all weights set to 1:

$\mathcal{L}=W_{\mathtt{energy}}*\mathcal{L}_{\mathtt{energy}}+W_{\mathtt{force}}*\mathcal{L}_{\mathtt{force}}+W_{\mathtt{proper}}*\mathcal{L}_{\mathtt{proper}}+W_{\mathtt{improper}}*\mathcal{L}_{\mathtt{improper}}$
(3)

To prevent overfitting and ensure optimal model performance, we applied
dropouts to the atom embedding stage (Stage 1) and symmetry-preserving stage
(Stage 2), as well as implemented an early stopping mechanism. After 800
epochs, the joint root mean square error (RMSE) loss, which incorporates both
energies and forces, was monitored using the validation set. This approach
allowed us to identify the point at which further training no longer improved
the model’s generalization capability.

##  Appendix D Protein-ligand benchmark dataset

The protein-ligand benchmark dataset can be found at
https://github.com/kntkb/protein-ligand-benchmark-custom. It consists of 4
target systems (Tyk2, Cdk2, P38, and Mcl1) and a total of 76 ligands. This
dataset was curated from the openforcefield/protein-ligand-benchmark
repository (https://github.com/openforcefield/protein-ligand-
benchmark/tree/d3387602bbeb0167abf00dfb81753d8936775dd2). The protein
structures and ligand poses, as well as the ligand transformations, were
manually curated, while the experimental results were adopted from the
original repository. The protein and ligand structures were prepared using
Maestro from Schrodinger 2022-2.

The PDB structure of a protein-ligand complex was imported and processed using
the default settings of prepwizard, along with additional options including
filling in missing side chains and loops using Prime, capping termini, and
deleting waters beyond 5.0 Å from het groups. The tautomer states of the
ligand complexed with the protein were manually inspected, and the most
reasonable state was chosen from a human perspective. For the protein
residues, the protonation and tautomer states were optimized using the default
settings of H-bond assignment. Subsequently, a restrained minimization was
performed using the OPLS4 force field, with an RMSD convergence threshold of
0.3 Å for the heavy atoms. The exported protein structure from the complex
served as the initial protein structure, and X-ray water molecules were
exported if necessary, such as buried water molecules in the binding pocket.

For the ligand poses, a flexible ligand alignment approach was applied with
respect to the PDB ligand pose found in the protein-ligand complex structure.
The default settings of ligprep were used to generate all possible ligand
tautomer states, which were then visually inspected to choose the most
reasonable state. Subsequently, ligand alignment was performed by aligning all
ligands to the PDB ligand pose found in the protein-ligand complex structure,
using the Ligand Alignment module in Maestro with Bemis-Murcko scaffold or
maximum common scaffold constrain. The ligand poses were manually adjusted,
taking into account the binding site environment, which involved rotating
ligand torsions and minimizing selected atoms to alleviate severe atom clashes
and obtain better initial poses.

Finally, the ligand transformation networks were defined manually by human
experts, creating a outward radial map with the simplest ligand in the center.
In the case of P38 and Mcl1, R-group substituent from multiple scaffold
positions and scaffold hopping were observed. In such cases, ligand
transformations were grouped into categories to resemble different structure-
activity relationship purposes while maintaining a simplified ligand
transformation network.

##  Appendix E Alchemical free energy calculations using protein-ligand
benchmark dataset

The Python code used to perform the alchemical protein-ligand binding free
energy benchmark experiment is available at https://github.com/choderalab/pl-
benchmark-espaloma-experiment. We utilized the Perses 0.10.1 relative
alchemical free energy calculation infrastructure [65], which is based on
OpenMM 8.0.0 [24], openmmtools 0.22.1 [10], and a modified version of
openmmforcefields 0.11.0 package [11]
(https://github.com/kntkb/openmmforcefields/tree/6d2c3dcd33d9800a32032d28b6b2dca92f348a43)
to support espaloam-0.3.0.

All systems were solvated with TIP3P water [35] with 9.0 Å buffer around the
protein, and the system was neutralized with the Joung and Cheatham monovalent
counterions [36] with 300 mM NaCl salt concentration. The protein was
parametrized with Amber ff14SB force field [53], and the small molecules were
parametrized with openff-2.1.0 [4], espaloma-0.3.0, or espaloma-0.2.2 [75].
Additionally, the protein-ligand was self-consistently parametrized with
espaloam-0.3.0, and a modified version of Perses 0.10.1
(https://github.com/kntkb/perses/tree/0d069fc1cf31b8cce1ae7a1482c3fa46bc1382d2)
was used to perform the protein-ligand binding free energy calculations.

Alchemical free energy calculations were simulated with replica exchange among
Hamiltonians with Gibbs sampling [12]. All simulations were performed with 12
alchemical states for 10 ns/replica for Tyk2 and Cdk2, 15 ns/replica for Mcl1,
and 20 ns/replica for P38, with replica exchange attempts made every 1 ps. The
simulations were performed at 300 K and 1 atm using a Monte Carlo Barostat and
Langevin BAOAB integrator [46] with a collision rate of 1/ps. Bonds to
hydrogen were constrained, and hydrogen atom masses were set to 3.0 amu by
transferring the masses connected to the heavy atoms, allowing for simulations
with a 4 fs timestep.

Atom mappings were generated from the provided geometries in the benchmark
set. Atoms within 0.5 Å of the transforming ligand pairs were detected as
valid mapping atoms using the use_given_geometries functionality in Perses.

PyMBAR 3.1.1 [67] was used to compute the relative free energy, while absolute
free energies up to an additive constant were estimated using a least-squares
estimation strategy [80] using a modified version of OpenFE cinnabar 0.3.0
package [52]
(https://github.com/kntkb/cinnabar/tree/de7bc6623fb25d75848aa1c9f538b77cd02a4b01).
Both experimental and calculated absolute free energies were shifted to their
respective means before computing the statistics.

Figure 4:  Training espaloma-0.3.0 on an extensive quantum chemical dataset
significantly improves protein-ligand binding affinity calculations on the
Cdk2 system. (a) We show the X-ray structure used for all free energy
calculation as well as the 2D structures of all ligands in the Cdk2 protein-
ligand benchmark dataset. An outward radial map with ligand #1 in the center
was used for the alchemical ligand transformations. We used the Perses 0.10.1
relative free energy calculation infrastructure [65] to calculate the relative
free energy and assess the performance using (b) espaloma-0.2.2 [75], (c)
espaloma-0.3.0, and (d) openff-2.1.0 [4]. Amber ff14SB force field [53] used
to parametrize the protein for all cases. espaloma-0.2.2 and espaloma-0.3.0
achieves an absolute free energy ($\Delta G$) RMSE of 1.45 [95% CI: 0.55,
2.33] and 0.56 [95% CI: 0.34, 0.73], respectively, indicating that
espaloma-0.3.0 trained on extensive quantum chemical dataset significantly
improved protein-ligand binding affinity calculations on the Cdk2 system.
Figure 5:  Training espaloma-0.3.0 on an expanded quantum chemical dataset
improves protein-ligand binding affinity on the Tyk2 system. (a) We show the
X-ray structure used for all free energy calculation as well as the 2D
structures of all ligands in the Tyk2 protein-ligand benchmark dataset. An
outward radial map with ligand #1 in the center was used for the alchemical
ligand transformations. We used the Perses 0.10.1 relative free energy
calculation infrastructure [65] to calculate the relative free energy and
assess the performance using (b) espaloma-0.2.2 [75], (c) espaloma-0.3.0, and
(d) openff-2.1.0 [4]. Amber ff14SB force field [53] used to parametrize the
protein for all cases. Notably, espaloma-0.2.2 failed to simulate the
alchemical transformation of ligand #1 to ligand #2; hence one ligand is not
reported in (b). espaloma-0.2.2 and espaloma-0.3.0 achieves an absolute free
energy ($\Delta G$) RMSE of 0.73 [95% CI: 0.34, 1.02] and 0.48 [95% CI: 0.28,
0.64], respectively, suggesting that espaloma-0.3.0 tends to show improved
performance over espaloma-0.2.2.