PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments and Partial Charges

TL;DR

PhysNet is a deep neural network architecture that accurately predicts energies, forces, dipole moments, and partial charges of chemical systems, demonstrating state-of-the-art performance and good generalization to larger molecules and complex phenomena.

Contribution

This work introduces PhysNet, a novel neural network architecture tailored for chemical property prediction, with improved accuracy and generalization capabilities over existing models.

Findings

Achieves state-of-the-art results on QM9, MD17, and ISO17 benchmarks.

Successfully generalizes from small peptides to larger proteins like Ala10.

Explicit electrostatics inclusion is crucial for correct asymptotic potential energy surface description.

Abstract

In recent years, machine learning (ML) methods have become increasingly popular in computational chemistry. After being trained on appropriate ab initio reference data, these methods allow to accurately predict the properties of chemical systems, circumventing the need for explicitly solving the electronic Schr\"odinger equation. Because of their computational efficiency and scalability to large datasets, deep neural networks (DNNs) are a particularly promising ML algorithm for chemical applications. This work introduces PhysNet, a DNN architecture designed for predicting energies, forces and dipole moments of chemical systems. PhysNet achieves state-of-the-art performance on the QM9, MD17 and ISO17 benchmarks. Further, two new datasets are generated in order to probe the performance of ML models for describing chemical reactions, long-range interactions, and condensed phase systems. It is shown that explicitly including electrostatics in energy predictions is crucial for a qualitatively correct description of the asymptotic regions of a potential energy surface (PES). PhysNet models trained on a systematically constructed set of small peptide fragments (at most eight heavy atoms) are able to generalize to considerably larger proteins like deca-alanine (Ala$_{10}$): The optimized geometry of helical Ala$_{10}$ predicted by PhysNet is virtually identical to ab initio results (RMSD = 0.21 \r{A}). By running unbiased molecular dynamics (MD) simulations of Ala$_{10}$ on the PhysNet-PES in gas phase, it is found that instead of a helical structure, Ala$_{10}$ folds into a wreath-shaped configuration, which is more stable than the helical form by 0.46 kcal mol$^{-1}$ according to the reference ab initio calculations.

Tables

Table 1. Table 1: Hyperparameters of all models used in this work (unless specified otherwise).

hyperparameter	value	significance
$F$	$128$	dimensionality of feature space
$K$	$64$	number of radial basis functions
$N_{\mathrm{module}}$	$5$	number of stacked modular building blocks
$N_{\mathrm{residual}}^{\mathrm{atomic}}$	$2$	number of residual blocks for atom-wise refinements
$N_{\mathrm{residual}}^{\mathrm{interaction}}$	$3$	number of residual blocks for refinements of proto-message
$N_{\mathrm{residual}}^{\mathrm{output}}$	$1$	number of residual blocks in output blocks
$r_{\mathrm{cut}}$	10 Å	cut-off radius for interactions in the neural network

Table 2. Table 2: Number of structures for each species present in the S N 2 reactions dataset.

species	count
${\left[{\mathrm{FCH}}_3\mathrm{Cl}\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{FCH}}_3\mathrm{Br}\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{FCH}}_{3}I\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{ClCH}}_3\mathrm{Br}\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{ClCH}}_{3}I\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{BrCH}}_{3}I\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{44\hspace{0.17em}501}$
${\left[{\mathrm{FCH}}_{3}F\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{24\hspace{0.17em}801}$
${\left[{\mathrm{ClCH}}_3\mathrm{Cl}\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{24\hspace{0.17em}801}$
${\left[{\mathrm{BrCH}}_3\mathrm{Br}\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{24\hspace{0.17em}801}$
${\left[{\mathrm{ICH}}_{3}I\right]}^{-}$ $]{}_{32}{}^{+-}$	$\mathrm{24\hspace{0.17em}801}$
${\mathrm{CH}}_{3}F$ $]{}_{32}{}^{+-}$	$3500$
${\mathrm{CH}}_3\mathrm{Cl}$ $]{}_{32}{}^{+-}$	$3500$
${\mathrm{CH}}_3\mathrm{Br}$ $]{}_{32}{}^{+-}$	$3500$

Table 3. Table 3: Mean absolute errors in kcal mol -1 for energy predictions on the QM9 dataset for several models reported in the literature and different training set sizes. Results for PhysNet are averaged over five independent runs, whereas PhysNet-ens5 is the performance of an ensemble of five PhysNet models. Best results are shown in bold.

$N_{\mathrm{train}}$ + $N_{\mathrm{valid}}$	enn-s2s-ens5	DTNN	SchNet	HIP-NN	PhysNet	PhysNet-ens5
$\mathrm{110\hspace{0.17em}426}$	0.33	–	0.31	0.26	0.19	0.14
$\mathrm{100\hspace{0.17em}000}$	–	0.84	0.34	0.26	0.19	0.14
$\mathrm{50\hspace{0.17em}000}$	–	0.94	0.49	0.35	0.30	0.24

Table 4. Table 4: Mean absolute errors for predictions of energy (in kcal mol -1 ) and forces (in kcal mol -1  Å -1 ) for molecules in the MD17 dataset for several models reported in the literature. All models utilize a combined 50k structures for training and validation, but using different reference data: DTNN and HIP-NN are trained on energies only, GDML is trained on forces only, and SchNet and PhysNet are trained on energies and forces. Results for PhysNet are averaged over five independent runs, whereas PhysNet-ens5 is the performance of an ensemble of five PhysNet models. The best results in each category are shown in bold.

		DTNN	HIP-NN	GDML	SchNet	PhysNet	PhysNet-ens5
Aspirin	energy	—	—	0.13	0.12	0.12	0.12
	forces	—	—	0.02	0.33	0.06	0.04
Benzene	energy	0.04	0.06	0.07	0.07	0.07	0.07
	forces	—	—	0.24	0.17	0.15	0.14
Ethanol	energy	—	—	0.05	0.05	0.05	0.05
	forces	—	—	0.09	0.05	0.03	0.02
Malonaldehyde	energy	0.19	0.09	0.08	0.08	0.07	0.07
	forces	—	—	0.09	0.08	0.04	0.03
Naphthalene	energy	—	—	0.12	0.11	0.12	0.12
	forces	—	—	0.03	0.11	0.04	0.03
Salicylic acid	energy	0.41	0.20	0.11	0.10	0.11	0.11
	forces	—	—	0.03	0.19	0.04	0.03
Toluene	energy	0.18	0.14	0.09	0.09	0.10	0.10
	forces	—	—	0.05	0.09	0.03	0.03
Uracil	energy	—	—	0.11	0.10	0.10	0.10
	forces	—	—	0.03	0.11	0.03	0.03
total count of best in class		1	0	6	6	6	11

Table 5. Table 5: Mean absolute errors for predictions of energy (in kcal mol -1 ) and forces (in kcal mol -1  Å -1 ) for the two variants of the ISO17 benchmark. PhysNet is compared with the performance of SchNet. Results for PhysNet are averaged over five independent runs, whereas PhysNet-ens5 is the performance of an ensemble of five PhysNet models. Best results are shown in bold.

		SchNet	PhysNet	PhysNet-ens5
known molecules /	energy	0.36	0.10	0.10
unknown conformations	forces	1.00	0.12	0.08
unknown molecules /	energy	2.40	2.94	2.86
unknown conformations	forces	2.18	1.38	1.13

Table 6. Table 6: Mean absolute errors for predictions of energy (in kcal mol -1 ), forces (in kcal mol -1  Å -1 ) and dipole moments (in D) for the S N 2 reactions dataset for PhysNet with and without long-range augmentation (see Eqs.  11  and  12 ). Results for PhysNet are averaged over five independent runs, whereas PhysNet-ens5 is the performance of an ensemble of five models. The best results in each category are shown in bold. Note that because of the partial charge correction scheme (see Eq. 14 ), total charge is always predicted exactly. However, for completeness, the error for the prediction of total charge (in e 𝑒 e ) using the uncorrected partial charges is also reported.

	PhysNet	PhysNet-ens5	PhysNet	PhysNet-ens5
	without long-range	without long-range	with long-range	with long-range
energy	0.071	0.070	$0.009$	$0.009$
forces	0.035	0.032	$0.012$	$0.009$
dipole	—	—	$0.0044$	$0.0042$
charge	—	—	$\mathrm{0.000\hspace{0.17em}23}$	$\mathrm{0.000\hspace{0.17em}19}$

Table 7. Table 7: Mean absolute errors for predictions of energy and energy per atom (in kcal mol -1 ), forces (in kcal mol -1 ) Å -1 and dipole moments (in D) for the solvated protein fragments dataset. Results for PhysNet are averaged over five independent runs, whereas PhysNet-ens5 is the performance of an ensemble of five PhysNet models. The best results in each category are shown in bold. Note that because of the charge correction scheme (see Eq. 14 ), total charge is always predicted exactly. However, for completeness, the error for the prediction of total charge (in e 𝑒 e ) using the uncorrected partial charges is also reported.

	PhysNet	PhysNet-ens5
energy	$1.03$	$0.95$
energy/atom	$0.054$	$0.050$
forces	$0.88$	$0.72$
dipole	$0.060$	$0.054$
charge	$0.004$	$0.003$

Equations

$$\mathbf{y}=\mathbf{Wx}+\mathbf{b}$$

$$\mathbf{y}={\mathbf{W}}_2\sigma \left({\mathbf{W}}_1\mathbf{x}+{\mathbf{b}}_1\right)+{\mathbf{b}}_2$$

$${\mathbf{x}}_i^0={\mathbf{e}}_{Z_i}$$

$${\mathbf{x}}_i^{l+2}={\mathbf{x}}_i^l+{\mathbf{W}}^{l+1}\sigma \left({\mathbf{W}}^l\sigma ({\mathbf{x}}_i^l)+{\mathbf{b}}^l\right)+{\mathbf{b}}^{l+1}$$

$${\mathbf{x}}_i^{l+1}={\mathbf{u}}^l\circ {\mathbf{x}}_i^l+{\mathbf{W}}^l\sigma ({\mathbf{v}}_i^l)+{\mathbf{b}}^l$$

$${\tilde{\mathbf{v}}}_i^l=\sigma \left({{\mathbf{W}}_{\mathbf{I}}}^l\sigma ({\mathbf{x}}_i^l)+{{\mathbf{b}}_{\mathbf{I}}}^l\right)+\sum _{j\ne i}{\mathbf{G}}^l\mathbf{g}(r_{ij})\circ \sigma \left({{\mathbf{W}}_{\mathbf{J}}}^l\sigma ({\mathbf{x}}_j^l)+{{\mathbf{b}}_{\mathbf{J}}}^l\right)$$

$$g_k(r_{ij})=\varphi (r_{ij})\cdot \exp \left(-{\beta }_k{\left(\exp (-r_{ij})-{\mu }_k\right)}^2\right)$$

$$\varphi (r_{ij})=\{\begin{array}{ll}1-6{\left(\frac{r_{ij}}{r_{\mathrm{cut}}}\right)}^5+15{\left(\frac{r_{ij}}{r_{\mathrm{cut}}}\right)}^4-10{\left(\frac{r_{ij}}{r_{\mathrm{cut}}}\right)}^3& r_{ij}<r_{\mathrm{cut}}\\ 0& r_{ij}\ge r_{\mathrm{cut}}\end{array}$$

$${\mathbf{y}}_i^m={\mathbf{W}}_{\mathrm{out}}^m\sigma ({\mathbf{x}}_i^l)+{\mathbf{b}}_{\mathrm{out}}^m$$

$${\mathbf{y}}_i={\mathbf{s}}_{Z_i}\circ \left(\sum _{m=1}^{N_{\mathrm{module}}}{\mathbf{y}}_i^m\right)+{\mathbf{c}}_{Z_i}$$

$$E=\sum _{i=1}^N{E}_i$$

$$E=\sum _{i=1}^N{E}_i+k_e\sum _{i=1}^N\sum _{j>i}^N{\tilde{q}}_i{\tilde{q}}_j\chi (r_{ij})+E_{D3}$$

$$\chi (r_{ij})=\varphi (2r_{ij})\frac{1}{\sqrt{r_{ij}^2+1}}+\left(1-\varphi (2r_{ij})\right)\frac{1}{r_{ij}}$$

$${\tilde{q}}_i=q_i-\frac{1}{N}\left(\sum _{j=1}^N{q}_j-Q\right)$$

$$\mathbf{p}=\sum _{i=1}^N{\tilde{q}}_i{\mathbf{r}}_i$$

$$\begin{array}{r}\mathcal{L}=w_E\left|E-E^{\mathrm{ref}}\right|+\frac{w_F}{3N}\sum _{i=1}^N\sum _{\alpha =1}^3\left|-\frac{\partial E}{\partial r_{i,\alpha }}-F_{i,\alpha }^{\mathrm{ref}}\right|+{\mathcal{L}}_{\mathrm{nh}}\end{array}$$

$$\begin{array}{r}\mathcal{L}=w_E\left|E-E^{\mathrm{ref}}\right|+\frac{w_F}{3N}\sum _{i=1}^N\sum _{\alpha =1}^3\left|-\frac{\partial E}{\partial r_{i,\alpha }}-F_{i,\alpha }^{\mathrm{ref}}\right|+w_Q\left|\sum _{i=1}^N{q}_i-Q^{\mathrm{ref}}\right|\\ +\frac{w_p}{3}\sum _{\alpha =1}^3\left|\sum _{i=1}^N{q}_i{r}_{i,\alpha }-p_{\alpha }^{\mathrm{ref}}\right|+{\mathcal{L}}_{\mathrm{nh}}\end{array}$$

$${\mathcal{L}}_{\mathrm{nh}}=\frac{{\lambda }_{\mathrm{nh}}}{N}\sum _{i=1}^N\sum _{m=2}^{N_{\mathrm{module}}}\frac{(E_i^m{)}^2}{(E_i^m{)}^2+(E_i^{m-1}{)}^2}$$

$${\mathcal{L}}_{\mathrm{nh}}=\frac{{\lambda }_{\mathrm{nh}}}{N}\sum _{i=1}^N\sum _{m=2}^{N_{\mathrm{module}}}\frac{1}{2}\left[\frac{(E_i^m{)}^2}{(E_i^m{)}^2+(E_i^{m-1}{)}^2}+\frac{(q_i^m{)}^2}{(q_i^m{)}^2+(q_i^{m-1}{)}^2}\right]$$

Full text

K****eywords PhysNet $\cdot$ Computational chemistry  $\cdot$ Neural network  $\cdot$ Machine learning  $\cdot$ Potential energy surface  $\cdot$ Molecular dynamics

1 Introduction

As was stated by Dirac already in 1929,1 the Schrödinger equation (SE) in principle contains all that is necessary to describe the whole of chemistry. Unfortunately, the SE can only be solved in closed form for the simplest systems, hence computational and numerical methods have been devised to find approximate solutions. However, even with these approximations, solving the electronic SE is computationally demanding and, depending on the accuracy required and the approximations used, is only tractable for a limited number of atoms.2

For this reason, machine learning (ML) methods have become increasingly popular in recent years in order to circumvent the solution of the SE altogether. Such approaches give a computer the ability to learn patterns in data without being explicitly programmed3 and have been used in the past to estimate properties of unknown compounds or structures after being trained on a reference dataset.4; 5; 6; 7 Of particular interest in this context is the energy $E$, as the forces derived from it can be used to drive molecular dynamics (MD) simulations. Since the Hamiltonian of a chemical system is uniquely determined by the external potential, which in turn depends on a set of nuclear charges $\{Z_{i}\}$ and atomic positions $\{\mathbf{r}_{i}\}$, all information necessary to determine $E$ is contained in $\{Z_{i},\mathbf{r}_{i}\}$. Hence, there must exist an exact mapping $f:\{Z_{i},\mathbf{r}_{i}\}\mapsto E$, which is usually referred to as a potential energy surface (PES).

Artificial neural networks (NNs)8; 9; 10; 11; 12; 13; 14 are a particular class of ML algorithms proven to be general function approximators15; 16 and thus ideally suited to learn a representation of the PES. For small systems, PESs based on NNs have been designed in the spirit of a many-body expansion,17; 18; 19 but these approaches scale poorly for large systems, because they typically involve a large number of individual NNs (one for each term in the many-body expansion). An alternative approach, known as high-dimensional neural network (HDNN),20 decomposes the total energy of a chemical system into atomic contributions and uses a single NN (or one for each element) for the energy prediction. It relies on the chemically intuitive assumption that the contribution of an atom to the total energy depends mainly on its local environment.

Two variants of HDNNs can be distinguished: One is “descriptor-based”, also referred to as Behler-Parrinello networks,20 for which the environment of an atom is encoded in a hand-crafted descriptor,21; 22; 23; 24 which is used as input of a multi-layer feed-forward NN. Examples for this kind of approach are ANI25 and TensorMol.26 In the second “message-passing”27 variant, nuclear charges and Cartesian coordinates are used as input and a deep neural network (DNN) is used to exchange information between individual atoms, such that a meaningful representation of their chemical environments is learned directly from the data. This approach was first introduced by the DTNN28 and has since been refined in other DNN architectures, for example SchNet29 or HIP-NN.30 While both types of HDNN perform well, the message-passing variant was found to be able to adapt better to the training data and usually achieves a better performance.31

This work introduces PhysNet, a HDNN of the message-passing type designed based on physical principles. It is shown that PhysNet improves upon or matches state-of-the-art performance on the QM9,32 MD17,33 and ISO1729 benchmark data sets. Further, two new benchmark datasets are presented: The first set contains structures probing the PES of SN2 reactions of methyl halides with halide anions. Capturing the correct long-range physics is particularly challenging for this dataset due to the presence of strong charge-dipole interactions. By comparing different variants of PhysNet, it is demonstrated that the explicit inclusion of electrostatics in the energy prediction is important for correctly describing asymptotic regions of the PES.

The second dataset contains structures for small protein-like compounds (consisting of at most eight heavy atoms) and clusters with water molecules. Available benchmark datasets cover conformational and chemical degrees of freedom, but they usually do not probe many-body intermolecular interactions, which are important in the description of condensed phase systems. Due to their biological importance, proteins in aqueous solution are a particularly relevant system of this kind. However, even small proteins contain hundreds of atoms, which makes ab initio reference calculations for them prohibitively expensive. Fortunately, it is possible to construct a predictive ML model for large molecules by training it only on smaller molecules that are structurally similar.34 These so-called “amons” can be readily constructed by considering a large molecule as chemical graph, generating all possible connected subgraphs with a fixed number of heavy atoms and saturating the resulting structures with hydrogen atoms.34 Due to the fact that most proteins are comprised of only $20$ different amino acids, many bonding patterns are shared between proteins and a comparatively small number of amons is sufficient to cover all possibilities.

It is shown that a PhysNet model trained on this data is able to accurately predict interaction energies of common sidechain-sidechain and backbone-backbone interactions in proteins. The model also generalizes to considerably larger molecules than it was trained on: When an ensemble35 of PhysNet models is used to optimize the geometry of helical deca-alanine (Ala10), the results are almost indistinguishable from a structure optimized at the reference DFT level of theory ($\mathrm{RMSD}=0.21$ Å). By running unbiased MD simulations on the PhysNet-PES, a “wreath-shaped” configuration of Ala10 is found, which, according to ab initio calculations, is as stable as the helical form in gas phase (the wreath-shaped form is even slightly lower in energy by $0.46$ kcal mol*-1*).

In section 2, the PhysNet architecture and the process used for generating the reference data for SN2 reactions and solvated protein fragments is described in detail. The performance of PhysNet on these datasets and other commonly used quantum chemical benchmarks is reported in section 3. Further, the generalization capabilities of the model are explored by applying it to Ala10. Finally, the results are discussed and summarized in section 4.

2 Methods

In section 2.1, the building blocks of PhysNet and its complete architecture are described. Further, the procedure used for training PhysNet on reference data (how the neural network parameters are fitted) is described in section 2.2. Because currently available benchmark datasets do not cover chemical reactions or many-body intermolecular interactions (as exhibited by condensed phase systems), two additional sets of reference data for SN2 reactions and solvated protein fragments were constructed and their generation is detailed in section 2.3.

2.1 Neural network

The basic building block of every fully-connected NN are so-called dense layers. They take an input vector $\mathbf{x}\in\mathbb{R}^{n_{\rm in}}$ and output a vector $\mathbf{y}\in\mathbb{R}^{n_{\rm out}}$ based on the transformation

[TABLE]

where $\mathbf{W}\in\mathbb{R}^{n_{\rm out}\times n_{\rm in}}$ and $\mathbf{b}\in\mathbb{R}^{n_{\rm out}}$ are parameters and $n_{\rm in}$ and $n_{\rm out}$ denote the dimensionality of input and output vectors, respectively. Note that a single dense layer can only represent a linear transformation from input to output. In order to model arbitrary non-linear relationships, at least two dense layers need to be stacked and combined with a (non-linear) activation function $\sigma$, i. e.

[TABLE]

Throughout this work, the notation $\sigma(\mathbf{x})$ means that a scalar function $\sigma$ is applied to a vector $\mathbf{x}$ entrywise, i.e. $\sigma(\mathbf{x})=[\sigma(x_{1})\ \cdots\ \sigma(x_{n_{\rm in}})]^{\mathsf{T}}$.

Two dense layers combined according to Eq. 2 can already approximate any (non-linear) mapping between input $\mathbf{x}$ and output $\mathbf{y}$, provided that the first (“hidden”) layer is “wide” enough (contains sufficiently many neurons) and an appropriate activation function $\sigma$ is used.15 Many different choices for $\sigma$ are possible,36 for example, popular activation functions include $\tanh(x)$15 or $\max\left(0,x\right)$.37 In this work, the shifted softplus function29 given by $\sigma(x)=\log\left(e^{x}+1\right)-\log\left(2\right)$ is used (see Fig. 1).

While shallow neural networks composed of just two dense layers (Eq. 2) are already capable of approximating arbitrary functions,36 deep neural networks composed of more than two layers were shown to be exponentially more parameter efficient in their approximation capability.38

In this work, several reusable building blocks are combined in a modular fashion to form a deep neural network architecture, which predicts atomic contributions to properties (such as energy) of a chemical system composed of $N$ atoms based on atomic features $\mathbf{x}_{i}\in\mathbb{R}^{F}$ (here, $F$ denotes the dimensionality of the feature space). The features simultaneously encode information about nuclear charge $Z$ and local atomic environment of each atom $i$ and are constructed by iteratively refining an initial representation depending solely on $Z_{i}$ through coupling with the feature vectors $\mathbf{x}_{j}$ of all atoms $j\neq i$ within a cut-off radius $r_{\rm cut}$.

In the following, the individual building blocks of the neural network and important underlying concepts are described in more detail. The complete architecture is schematically represented in Fig. 2. Note that for the remainder of this section, superscripts $l$ are used to denote features or parameters of layer $l$. All models are implemented in the TensorFlow39 framework.

Embedding layer

An embedding is a mapping from a discrete object to a vector of real numbers. For example, word embeddings41 find wide spread use in the field of natural language processing, where semantically similar words are mapped such that they appear close to each other in the embedding space. In this work, atomic numbers $Z\in\mathbb{N}$ are mapped to embeddings $\mathbf{e}_{Z}\in\mathbb{R}^{F}$, where the entries of the embedding vectors $\mathbf{e}_{Z}$ are parameters. The embedding layer initializes the atomic features of an atom with nuclear charge $Z$ to the corresponding embedding vector $\mathbf{e}_{Z}$ (Eq. 3).

[TABLE]

The output of the embedding layer is passed to a stack of $N_{\rm module}$ modules sharing the same composition, but not parameters.

Module

Each module in the stack (except for the first one, which receives its input from the embedding layer) takes the output of the previous module and couples the features $\mathbf{x}_{i}$ of each atom $i$ with the features $\mathbf{x}_{j}$ of all atoms $j$ within the cut-off distance $r_{\rm cut}$ through an interaction block. The features are then further refined atom-wise through $N^{\rm atomic}_{\rm residual}$ residual blocks. Subsequently, the computation splits into two branches: One branch passes the atomic features onwards to the next module in the stack (if present) without further modification, whereas the other branch passes the features to an output block, which computes the module’s contribution to the final prediction. The split into two branches helps to decouple the feature representations passed between modules from the prediction task at hand.

Residual block

The ability of a neural network to model arbitrary functions should always increase, or at least remain the same, when the depth (i.e. the amount of dense layers stacked on top of each other) is increased, as additional layers could in principle always reduce to the identity mapping and should therefore never decrease the performance. However, this is not observed in practice: As neural networks get deeper, they become increasingly difficult to train because of the vanishing gradients problem,42 which leads to a degradation of their performance. In order to alleviate this, it was proposed to add shortcut connections to the neural network architecture that skip one or several layers, creating a so-called residual block.43

Since their first introduction, the design of residual blocks was further refined to allow completely unhindered gradient flow through all layers of a neural network.40 It was shown that stacking these so-called pre-activation residual blocks allows successfully training neural networks more than $1000$ layers deep.40 In this work, pre-activation residual blocks are used extensively to refine the atomic features according to

[TABLE]

where $\mathbf{x}_{i}^{l}$ and $\mathbf{x}_{i}^{l+2}$ are input and output features, respectively, and $\mathbf{W}^{l},\mathbf{W}^{l+1}\in\mathbb{R}^{F\times F}$ and $\mathbf{b}^{l},\mathbf{b}^{l+1}\in\mathbb{R}^{F}$ are parameters.

Interaction block

In order to predict properties that depend on the environment of an atom $i$, it is important to model interactions with its surrounding atoms $j$ in a chemically and physically meaningful manner. In doing so, it is crucial that all known invariances of the property of interest are respected: For example, the energy of a molecular system is known to be invariant with respect to translation, rotation and permutation of equivalent atoms.21 Further, since it is a valid assumption that most (but not all) chemical interactions are inherently short-ranged,24 it is meaningful to introduce a cut-off radius $r_{cut}$, such that only interactions within the local environment of an atom are considered. Apart from encoding chemical knowledge directly in the modelling of interactions, this approach has the important computational advantage of making predictions scale linearly with system size.

Based on these design principles, the feature vector $\mathbf{x}$ of an atom is refined by interacting with its local environment through a “message”27 $\mathbf{v}\in\mathbb{R}^{F}$ according to

[TABLE]

where $\mathbf{u}^{l},\mathbf{b}^{l}\in\mathbb{R}^{F}$ and $\mathbf{W}^{l}\in\mathbb{R}^{F\times F}$ are parameters and ‘$\circ$’ denotes the Hadamard (entrywise) product (see Fig. 2C). The gating vector $\mathbf{u}$, inspired by the gated recurrent unit,44 allows individual entries of the feature vector to be damped or reinforced during the update. The final message $\mathbf{v}$ used in Eq. 5 is obtained by passing a proto-message ${\tilde{\mathbf{v}}}$ through $N_{\rm residual}^{\rm interaction}$ residual blocks (Eq.  4), where ${\tilde{\mathbf{v}}}$ is given by

[TABLE]

and $\mathbf{W_{I}}^{l},\mathbf{W_{J}}^{l}\in\mathbb{R}^{F\times F}$, $\mathbf{b_{I}}^{l},\mathbf{b_{J}}^{l}\in\mathbb{R}^{F}$ and $\mathbf{G}^{l}\in\mathbb{R}^{F\times K}$ are parameters and $r_{ij}$ denotes the Euclidean distance between atoms $i$ and $j$. The vector $\mathbf{g}(r_{ij})=\left[g_{1}(r_{ij})\ \cdots\ g_{K}(r_{ij})\right]^{\mathsf{T}}$ is composed of the values of $K$ radial basis functions of the form

[TABLE]

where $\mu_{k},\beta_{k}\in\mathbb{R}_{>0}$ are parameters that specify centre and width of $g_{k}(r_{ij})$, respectively, and $\phi(r_{ij})$ is a smooth cut-off function given by45

[TABLE]

that ensures continuous behaviour when an atom enters or leaves the cut-off sphere. The vector $\mathbf{G}^{l}\mathbf{g}(r_{ij})\in\mathbb{R}^{F}$ takes the role of a learnable attention mask46 that selects different features based on the pairwise distance $r_{ij}$ between atoms. Note that the Gaussian in Eq. 7 takes $\exp(-r_{ij})$ instead of $r_{ij}$ as its argument, which biases attention masks towards a functional form that decays exponentially with $r_{ij}$ (see Fig. 3). Such a bias is meaningful for a chemical system, as it entails the physical knowledge that bound state wave functions in two-body systems decay exponentially. Since only pairwise distances are used in Eq. 6, the output of an interaction block is automatically translationally and rotationally invariant, while the commutative property of summation ensures permutational invariance.47

Output block

Each output block passes the atomic features through $N^{\rm output}_{\rm residual}$ additional residual blocks and computes the output $\mathbf{y}_{i}^{m}\in\mathbb{R}^{n_{\rm out}}$ of module $m$ for atom $i$ by a linear transformation of the activated features according to

[TABLE]

where $\mathbf{W}_{\rm out}^{m}\in\mathbb{R}^{F\times n_{\rm out}}$ and $\mathbf{b}_{\rm out}^{m}\in\mathbb{R}^{n_{\rm out}}$ are parameters. How many entries the output vector $\mathbf{y}_{i}^{m}$ has depends on how many atomic properties are predicted at once. In this work, two variants are considered: The first version predicts atomic energy contributions along with atomic partial charges (i.e. $\mathbf{y}_{i}^{m}=\left[E_{i}^{m}\ q_{i}^{m}\right]^{\mathsf{T}}$ and $n_{\rm out}=2$), whereas the second version predicts just atomic energy contributions (i.e. $\mathbf{y}_{i}^{m}=\left[E_{i}^{m}\right]$ and $n_{\rm out}=1$). In principle, other properties could be predicted as well.

Final prediction

The final atomic properties $\mathbf{y}_{i}$ are obtained by summing the contributions of the individual modules according to

[TABLE]

where $\mathbf{s}_{Z},\mathbf{c}_{Z}\in\mathbb{R}^{n_{\rm out}}$ are learnable element-specific scale and shift parameters depending on the nuclear charge $Z_{i}$ of atom $i$. The scaling and shifting of the output decouples the values of other parameters from the numeric range of target properties, which depends mainly on the chosen system of units. Element-specific (instead of global) parameters are motivated by a previous observation indicating that atomic properties of distinct elements can span vastly different ranges.24

The final prediction for the total energy of a system of interest composed of $N$ atoms is obtained by summation of the atomic energy contributions $E_{i}$:

[TABLE]

A potential shortcoming of Eq. 11 is the fact that all long-range interactions contributing to $E$ beyond the cut-off radius $r_{\rm cut}$ cannot be properly accounted for. As long as $r_{\rm cut}$ is chosen sufficiently large, this is not an issue. However, in order to account for electrostatic interactions which decay with the inverse of the distance, a large cut-off would be necessary and reduce the computational efficiency of the model. Fortunately, their functional form is known and they can be explicitly included when computing $E$. Other types of long-range interactions for which the functional form is also known analytically, for example dispersion corrections like DFT-D3,48 can be included as well.

Because of the shortcomings of Eq. 11, for a PhysNet model that also predicts atomic partial charges, $E$ is calculated by

[TABLE]

instead, where $E_{\rm D3}$ is the DFT-D3 dispersion correction,48 $k_{e}$ is Coulomb’s constant, $\tilde{q}_{i}$ and $\tilde{q}_{j}$ are corrected partial charges (see Eq. 14) of atoms $i$ and $j$, and $\chi(r_{ij})$ is given by

[TABLE]

where $\phi(r_{ij})$ is the cut-off function given by Eq. 8. Here, $\chi(r_{ij})$ smoothly interpolates between the correct $r_{ij}^{-1}$ dependence of Coulomb’s law at long-range and a damped term at small distances to avoid the singularity at $r_{ij}=0$ (see Fig. 4). The corrected partial charges $\tilde{q}_{i}$ are obtained from the partial charges $q_{i}$ predicted by the neural network (Eq. 9) according to

[TABLE]

where $Q$ is the total charge of the system. As neural networks are a purely numerical algorithm, it is not guaranteed a priori that the sum of all predicted atomic partial charges $q_{i}$ is equal to the total charge $Q$ (although the result is usually very close when the neural network is properly trained), so a correction scheme like Eq. 14 is necessary to guarantee charge conservation.

It should be pointed out that the summation over all atom pairs when evaluating the long-range interactions in Eq. 12 makes the evaluation of $E$ scale quadratically with system size. Fortunately, various schemes to recover linear scaling are described in the literature, e.g. Ewald summation49 or cut-off methods50, and can be applied without modification.

Note that the concept of augmenting a neural network potential with long range interactions is not novel and was first proposed in [\NAT@swatrue\NAT@parfalse\NAT@citexnum[][]artrith2011high]. However, most previous works use a separately trained neural network to predict atomic partial charges, in contrast to the present work, which uses a single neural network to learn both, atomic energy contributions and partial charges (see Eqs. 9 and 10). Aside from computational advantages (only a single neural network needs to be trained and evaluated), shared feature representations in such “multi-task learning” are believed to increase the generalization capabilities (transferability) of a model.

Apart from allowing the computation of long-range electrostatic interactions, the partial charges $\tilde{q}_{i}$ can also be used to predict the electric dipole moment $\mathbf{p}$ of a structure according to

[TABLE]

where $\mathbf{r}_{i}$ are the Cartesian coordinates of atom $i$. The ability to predict $\mathbf{p}$ is useful for example for the calculation of infrared spectra.

Finally, analytical derivatives of $E$ (see Eqs. 11 and 12) with respect to the Cartesian coordinates $\{\mathbf{r}_{1},\dots,\mathbf{r}_{N}\}$ of the atoms, for example in order to derive the forces $\mathbf{F}_{i}$ acting on each atom $i$, are readily obtained by reverse mode automatic differentiation.

Hyperparameters

The PhysNet architecture can be tuned by hyperparameters that control width and depth of the neural network (see Fig. 2). While it would be possible to optimize hyperparameters for individual learning tasks, for example via a grid search, it was found that this is not necessary for good performance. For simplicity, all models used in this work share the same architecture with the hyperparameters summarized in Table 1, unless specified otherwise.

2.2 Training

Before the training of a neural network starts, its parameters need to be initialized. All entries of the embedding vectors $\mathbf{e}_{Z}$ are initialized with random values uniformly distributed between $-\sqrt{3}$ and $\sqrt{3}$ (such that they have unit expected variance) and the weight matrices $\mathbf{W}$,$\mathbf{W_{I}}$ and $\mathbf{W_{J}}$ are initialized to random orthogonal matrices with entries scaled such that their variance corresponds to the value recommended in the Glorot initialization scheme. The entries of all bias vectors $\mathbf{b}$, $\mathbf{b_{I}}$, $\mathbf{b_{J}}$, $\mathbf{b}_{\rm out}$ and matrices $\mathbf{G}$ and $\mathbf{W}_{\rm out}$ are initialized to zero, whereas the entries of the gating vectors $\mathbf{u}$ are initialized to one. The centres $\mu_{k}$ of the of the radial basis functions are set to $K$ equally spaced values between $\exp(-r_{\rm cut})$ and $1$ and their widths $\beta_{k}$ are initialized to $\left(2K^{-1}\left(1-\exp(-r_{\rm cut})\right)\right)^{-2}$ (see Fig. 3).

After initialization, the parameters of the neural network are optimized by minimizing a loss function $\mathcal{L}$ using AMSGrad with a learning rate of $10^{-3}$ (other hyperparameters of the optimizer are set to the default values recommended in [\NAT@swatrue\NAT@parfalse\NAT@citexnum[][]reddi2018convergence]) and a batch size of $32$ randomly chosen reference structures. For predicting energies without long-range augmentation (Eq. 11) the loss function is

[TABLE]

whereas when energies and charges are predicted (Eq. 11), the loss function is

[TABLE]

Here, $E^{\rm ref}$ and $Q^{\rm ref}$ are reference energy and total charge, $p_{\alpha}^{\rm ref}$ are the Cartesian components of the reference dipole moment $\mathbf{p}^{\rm ref}$, $F_{i,\alpha}^{\rm ref}$ are the Cartesian components of the reference force $\mathbf{F}_{i}^{\rm ref}$ acting on atom $i$ and $r_{i,\alpha}$ is the $\alpha$th Cartesian coordinate of atom $i$. The energy prediction $E$ is given either by Eq. 11 or Eq. 12, depending on which variant of PhysNet is used.

The weighting hyperparameters $w_{E}$, $w_{F}$, $w_{Q}$ and $w_{p}$ determine the relative contribution of the individual error terms to the loss term. Note that the numeric ranges of the error terms (and therefore their contribution to $\mathcal{L}$) also depend on the chosen system of units. For simplicity, weighting hyperparameters are not optimized for individual learning tasks and instead always set to $w_{E}=w_{Q}=w_{p}=1$ and $w_{F}=10^{2}$ (when all quantities are measured in atomic units). The higher relative weight of force errors is motivated by the fact that forces alone determine the dynamics of a chemical system and accurate force predictions are therefore most important for MD simulations. Note that for datasets where any of the reference quantities used in Eqs. 16 and 17 are not available, the corresponding weight is set to zero.

The term $\mathcal{L}_{\rm nh}$ is a “non-hierarchicality penalty”, inspired by a similar regularization method introduced in [\NAT@swatrue\NAT@parfalse\NAT@citexnum[][]lubbers2018hierarchical], given either by

[TABLE]

or

[TABLE]

depending on which variant of PhysNet is used, and $\lambda_{\rm nh}$ is the corresponding regularization hyperparameter. The $\mathcal{L}_{\rm nh}$ term penalizes when the predictions of individual modules do not decay with increasing depth in the hierarchy. Since deeper feature representations of atoms capture increasingly higher-order interactions, such a regularization is motivated by the fact that higher-order terms in many-body expansions of the energy are known to decay rapidly in magnitude. For simplicity, $\lambda_{\rm nh}$ is not tuned for individual learning tasks and instead always set to $10^{-2}$.

During training, an exponential moving average of all parameter values is kept using a decay rate of $0.999$. Overfitting is prevented using early stopping: After every epoch (one pass over all reference structures in the training set), the loss function is evaluated on a validation set of reference structures using the parameters obtained from the exponential moving average. After training, the model that performed best on the validation set is selected. Since the validation set is used indirectly during the training procedure, the performance of the final models (see section 3) is always measured on a separate test set.

Note that only true quantum mechanical observables, such as total energy, forces or dipole moments, are used as reference when training the neural network (see Eqs. 16 and 17). While it would also be possible to train directly on atomic energies and partial charges obtained using a decomposition method, such schemes are essentially arbitrary and it is unclear whether the corresponding decompositions are meaningful. Further, it is not always guaranteed that the quantities obtained from such methods vary smoothly when the molecular geometry changes, which makes it difficult for a neural network to learn them. By only relying on quantum mechanical observables, the model automatically learns to perform a smooth decomposition in a data-driven way.

2.3 Dataset generation

In the following, the generation of two new benchmark datasets is described, which probe chemical reactivity, long-range electrostatics, and many-body intermolecular interactions.

SN2 reactions

The SN2 reactions dataset probes chemical reactions of the kind X*-* + H3C–Y $\rightarrow$ X–CH3 + Y*-* and contains structures for all possible combinations of X,Y $\in$ {F, Cl, Br, I}. The reactions of methyl halides with halide anions are a prototypical examples for chemical reactions and have been studied extensively both experimentally and theoretically. It consists of different geometries for the high-energy transition regions, ion-dipole bound state complexes and long-range ($>10$ Å) interactions of CH3X molecules with Y*-* ions. The dataset also includes various structures for several smaller molecules that can be formed in fragmentation reactions, such as CH3X, HX, CHX or CH2X*-* with X $\in$ {F, Cl, Br, I}, as well as geometries for H2, CH2, CH${}_{3}^{+}$ and XY interhalogen compounds for all possible combinations of X,Y $\in$ {F, Cl, Br, I}. In total, the dataset provides reference energies, forces, and dipole moments for $452\,709$ structures calculated at the DSD-BLYP-D3(BJ)/def2-TZVP level of theory.

Different conformations for each species present in the SN2 reactions dataset were sampled by running MD simulations at a temperature of $5000$ K with a time step of $0.1$ fs using the Atomic Simulation Environment (ASE). The necessary forces were obtained with the semi-empirical PM7 method implemented in MOPAC2016. Structures were saved every $10$ steps and for each of them, reference energies, forces and dipole moments were calculated at the DSD-BLYP-D3(BJ)/def2-TZVP level of theory using the ORCA 4.0.1 code. The DSD-BLYP functional is one of the best performing double hybrid methods in the GMTKN55 benchmark. All MD simulations were started from the PM7-optimized geometries. For the reaction complexes $[\mathrm{XCH}_{3}\mathrm{Y}]^{-}$, MD simulations were randomly started either from the respective van der Waals complexes, or, in order to better sample the transition regions, from the transition state (calculated with PM7) of the respective reaction. Further, long-range interactions were sampled by choosing a random conformation from the CH3X MD simulations and randomly placing an ion Y- in the vicinity of the CH3X molecule such that its distance to any other atom is at most $16$ Å. Table 2 lists the number of conformations for each species.

Solvated protein fragments

The solvated protein fragments dataset probes many-body intermolecular interactions between “protein fragments” and water molecules, which are important for the description of many biologically relevant condensed phase systems. It contains structures for all possible “amons” (hydrogen-saturated covalently bonded fragments) of up to eight heavy atoms (C, N, O, S) that can be derived from chemical graphs of proteins containing the $20$ natural amino acids connected via peptide bonds or disulfide bridges. For amino acids that can occur in different charge states due to (de-)protonation (i.e. carboxylic acids that can be negatively charged or amines that can be positively charged), all possible structures with up to a total charge of $\pm 2e$ are included. These structures are augmented with solvated variants containing a varying number of water molecules such that the total number of heavy atoms does not exceed $21$. The dataset also contains randomly sampled dimer interactions of protein fragments, as well as structures of pure water with up to $40$ molecules. For all included structures, several conformations sampled with MD simulations at $1000$ K are included. In total, the dataset contains reference energies, forces and dipole moments for $2\,731\,180$ structures calculated at the revPBE-D3(BJ)/def2-TZVP level of theory. On average, the structures contain $21$ atoms (with a maximum of $120$ atoms) and consist of 63% hydrogen, 19% carbon, 12% oxygen, 5% nitrogen, and 1% sulfur atoms.

The different structures in the dataset were constructed as follows: All amons with up to eight heavy atoms (C, N, O, S) were constructed according to the method described in [\NAT@swatrue\NAT@parfalse\NAT@citexnum[][]huang2017chemical] for all possible chemical graphs of proteins containing the $20$ natural amino acids connected via peptide bonds or disulfide bridges. For amons derived from amino acids that can occur in different charge states due to (de-)protonation, variants for all charge states are included. This results in $2307$ different molecules. In order to sample interactions with solvent molecules, the amon structures were augmented by randomly placing up to $20$ water molecules in their vicinity, such that the total number of heavy atoms does not exceed $21$. This results in $29\,991$ additional structures. Further, other intermolecular interactions were sampled by generating all possible dimers from amons with up to $3$ heavy atoms, resulting in $867$ possible combinations. Important interactions between different amino acids were included by adding sidechain-sidechain and backbone-backbone complexes from the BioFragment Database ($3480$ structures). Further, interactions in pure water were sampled by constructing water clusters with up to $21$ molecules. Each water cluster is complemented by a variant with an additional proton, as well as with a variant lacking one proton in order to sample the different possible charge states of water. This results in $24$ additional structures.

All structures were optimized using the semi-empirical PM7 method implemented in MOPAC2016. Starting from the optimized geometry, $100$ different conformations for each structure were sampled by running MD simulations (at the same level of theory) at a temperature of $1000$ K with a time step of $0.1$ fs using the ASE. Structures were saved every $10$ steps and for each of them, reference energies, forces and dipole moments were calculated at the revPBE-D3(BJ)/def2-TZVP level of theory using the ORCA 4.0.1 code. The revPBE functional is one of the best performing GGA functionals in the GMTKN55 benchmark.

While this initial data already covers many different chemical situations, it is not guaranteed that the contained structures cover chemical and configurational space sufficiently well to account for all situations that might be relevant in MD simulations. For this reasons, the initial dataset was iteratively augmented using an adaptive sampling method: An ensemble of three PhysNet models trained (see section 2.2) on the initial dataset is used to run MD simulations and all structures for which their predictions deviate by more than a threshold value (here $1$ kcal mol*-1*) are saved. Discrepancies in the predictions are a strong indicator that the dataset used for training the models does not contain sufficient information for a particular conformation. For each structure saved in this process, energies, forces and dipole moments were calculated with the reference ab initio method and added to the dataset. Afterwards, the models were retrained and the sampling process repeated. In total, the dataset was adaptively augmented in this way for four times, after which significant deviations between the predictions were found to be rare. Finally, energies, forces and dipole moments were calculated with the reference method for $10\,000$ structures of $40$ water molecules in a spherical arrangement (obtained by running MD simulations with PM7, see above) to include training examples similar to bulk phase water. The final dataset contains data for $2\,731\,180$ structures.

3 Results

In this section, PhysNet is applied to various quantum-chemical datasets that all probe different aspects of chemical space (i.e. chemical and/or conformational degrees of freedom). Apart from the well-established benchmarks QM9, MD17, and ISO17, the model is also applied to the two new datasets introduced in section 2.3.

QM9

The QM9 dataset is a widely used benchmark for the prediction of several properties of molecules in equilibrium. It consists of geometric, energetic, electronic, and thermodynamic properties for $\approx$134k small organic molecules made up of H, C, O, N, and F atoms. These molecules correspond to the subset of all species with up to nine heavy atoms (C, O, N, and F) out of the GDB-17 chemical universe database. All properties are calculated at the B3LYP/6-31G(2df,p) level of theory. About 3k molecules within QM9 fail a geometric consistency check or are difficult to converge and are commonly removed from the dataset. Since the QM9 contains only equilibrium geometries, the benchmark probes just chemical (but no conformational) degrees of freedom.

Table 3 compares the performance for predicting the energy on the pruned QM9 dataset ($\approx$131k structures) of various models published in the literature with PhysNet. Results for PhysNet are averaged over five independent runs using the same training set. The performance of PhysNet can be further improved by bagging: While a single PhysNet model already improves upon the state-of-the-art, an ensemble of five PhysNet models (PhysNet-ens5) brings the error down even further.

MD17

The MD17 dataset is a collection of structures, energies and forces from ab initio MD simulations of eight small organic molecules. All trajectories are calculated at a temperature of $500$ K and a resolution of $0.5$ fs using the PBE+vdW-TS electronic structure method. The datasets range in size from 150k to almost 1M conformations and cover energy ranges between $20$ to $48$ kcal mol*-1* and force components between $266$ to $570$ kcal mol*-1* Å*-1*. The task is to predict energies (and forces) using a separate model for each molecule. Since each task is limited to a single molecule, the MD17 benchmark probes only conformational (an no chemical) degrees of freedom.

Table 4 compares the performance of various ML models published in the literature on the MD17 benchmark with PhysNet. Results for PhysNet are averaged over five independent runs using the same training set and their ensemble prediction is also reported (PhysNet-ens5). When comparing different models, it should be noted that they use different subsets of the available data for training: DTNN and HIP-NN are trained on energies only, GDML is trained on forces only, and SchNet and PhysNet are trained on energies and forces. PhysNet-ens5 matches or even improves upon state-of-the-art performance for all molecules in at least one category (energy or forces).

ISO17

The ISO17 dataset consists of short MD trajectories of $127$ isomeric molecules with the composition C7O2H10 drawn randomly from the largest set of isomers in QM9. Each trajectory samples 5k conformations at a resolution of $1$ fs. In total, the dataset contains 635k structures, for which energies and forces, calculated at the PBE+vdW-TS level of theory, are reported. The task is to predict energies (and forces) for two different scenarios: In the first variant (known molecules / unknown conformations), the training set contains $\approx$80% of all molecules and conformations (400k structures for training 4k structures for validation) and the task is to predict the remaining $\approx$20% of conformations for the same subset of molecules present in the training set (101k structures). Thus, the first variant tests the generalization capabilities of the model for unknown conformations of previously seen molecules. In the second, more challenging variant (unknown molecules / unknown conformations), the training set remains the same, but the task is to predict all 5k conformations of the $\approx$20% of molecules not present in the training set (130k structures). Here, generalization capabilities of the model are tested for unknown conformations of unknown molecules. Both variants of the ISO17 benchmark probe chemical and conformational degrees of freedom.

Table 5 compares the performance of SchNet to the average performance of five PhysNet models, as well as their ensemble prediction (PhysNet-ens5) for the two variants of the ISO17 benchmark. While for the first variant of the benchmark (known molecules / unknown conformations), PhysNet achieves state-of-the-art performance on both energies and forces, for the second variant (unknown molecules / unknown conformations), PhysNet improves upon SchNet only for force predictions, but performs slightly worse for energies. This is likely due to the higher relative weight of force data during the training process (see section 2.2) and it is possible that the results on energy could be improved by tuning the corresponding weighting hyperparameters.

SN2 reactions

This benchmark dataset was newly generated for the present work, for a detailed description of the dataset, see section 2.3. The task is to predict energies, forces and dipole moments using a single model for all structures contained in the dataset, testing the generalization capabilities of the model across chemical and conformational degrees of freedom, chemical reactions, and challenging long-range intermolecular interactions.

Table 6 lists the performance of PhysNet with and without explicit long-range electrostatic interactions (see Eqs. 12 and 11). All models were trained on 400k structures with 5k structures used for validation. The results in each case are averaged over five independent runs and the performance of ensembles of five PhysNet models is also reported. Because of the partial charge correction scheme (see Eq. 14), the total charge is always predicted exactly. However, for completeness, the error for the prediction of total charge (in $e$) using the uncorrected partial charges is also given.

The model without explicit inclusion of long-range interactions (Eq. 11) performs significantly worse. This is to be expected, as for this dataset, ion-dipole interactions, which decay with the square of the distance, play an important role for determining the overall energy. As their influence extends well beyond the cut-off distance $r_{\rm cut}$ (here $10$ Å), a model without long-range augmentation cannot properly account for them. This effect is also seen in Fig. 5, which shows minimum energy paths (MEPs) for all SN2 reactions of the kind X*-* + H3C–Y $\rightarrow$ X–CH3 + Y*-* covered in the dataset (all possible combinations X–Y with X,Y $\in$ {F, Cl, Br, I}) along the reaction coordinate defined by the distance difference $r_{\rm CY}-r_{\rm CX}$. While PhysNet including explicit long-range interactions (Eq. 12) is able to reproduce the reference energies accurately across the whole range of values of the reaction coordinate, the NN without long-range interactions (Eq. 11) shows qualitatively wrong asymptotic behaviour (see Figs. 5 and 6). A correct description of the asymptotics is crucial for quantitative predictions of reaction rates with MD simulations, as errors can strongly influence the maximum impact parameter for collisions at which a reaction is still possible. Including the long-range behaviour in the functional form of an ML model has been used previously in the construction of PESs with kernel ridge regression and has the additional advantage that less reference data is needed in asymptotic regions. Note that the MEPs calculated with the reference ab initio method are not included in the training data.

Solvated protein fragments

For a detailed description of the dataset, see section 2.3. The task is to predict all properties (energy, forces and dipole moments) using a single model, which tests the generalization capabilities across chemical and conformational degrees of freedom in gas and solution phase, proton transfer reactions and challenging many-body intermolecular interactions.

Table 7 lists the performance of five PhysNet models (including long-range interactions, see Eq. 12) and their ensemble prediction (PhysNet-ens5). All models were trained on 2560k structures with 100k structures used for validation. Because the solvated protein fragments dataset contains structures with widely different numbers of atoms (up to 120), the MAE for energy predictions per atom is also reported.

Since non-covalent interactions play a crucial rule for the structure of large systems like proteins, PhysNet-ens5 was also used to predict interaction energies for sidechain-sidechain interactions (SSIs, $3380$ structures) and backbone-backbone interactions (BBIs, $100$ structures) in the BioFragment Database and compared to values calculated at the reference revPBE-D3(BJ)/def2-TZVP level of theory. For each case, interaction energies were determined by subtracting monomer energies from the energy of the complex. The predictions of PhysNet-ens5 correlate well with the reference values ($R^{2}>0.99$, see Fig. 7) and have mean absolute errors of 0.28 kcal mol*-1* and 0.21 kcal mol*-1* for the SSI and BBI complexes, respectively. Note that although complexes from the BioFragment Database are included in the solvated protein fragments dataset (see section 2.3), the reference data contains only total energies and models were therefore never directly trained to reproduce interaction energies. Despite this fact, PhysNet is able to learn a meaningful decomposition of the total energy into intramolecular and intermolecular contributions and predict interaction energies accurately.

In order to test whether predictions also generalize to larger molecules, deca-alanine (Ala10), which is a widely used model system to study protein folding dynamics, is considered as a test case. Starting from a previously published helical structure of Ala10 (capped with an acetylated N-terminus and amidated C-terminus), its geometry was optimized with the BFGS algorithm using PhysNet-ens5, as well as revPBE-D3(BJ)/def2-TZVP to determine the necessary energy gradients. The energies (relative to free atoms) of the optimized structures are $-11339.49$ kcal mol*-1* and $-11317.05$ kcal mol*-1* for the ab initio method and PhysNet-ens5, respectively, which corresponds to a relative prediction error of about $0.20$%. Although PhysNet-ens5 predicts the optimized structure to be about $0.207$ kcal mol*-1* atom*-1* less stable than the ab initio method, both optimized geometries are structurally almost indistinguishable ($\mathrm{RMSD}=0.21$ Å, see Fig. 8A). This result is remarkable, considering that the “protein fragments” in the solvated proteins dataset contain at most eight heavy atoms (see section 2.3), whereas Ala10 consists of $54$ heavy atoms ($109$ atoms in total).

As a final test, the folding of Ala10 was investigated by running unbiased Langevin dynamics with the ASE at a temperature of $300$ K and using a time step of $0.1$ fs. The necessary forces were obtained from the predictions of PhysNet-ens5. Starting from the optimized structure of stretched Ala10, the simulation was run for a total of $400\,000$ time steps ($40$ ps). After about $30$ ps of simulation, Ala10 folds into a wreath-shaped structure (see Fig. 8B), in which it remains for the remainder of the simulation.

In order to determine whether the PES explored during the dynamics is representative of the PES computed using the reference method, the energy of $20$ structures sampled at $2$ ps intervals along the trajectory was evaluated with PhysNet-ens5 and revPBE-D3(BJ)/def2-TZVP. On average, the prediction error for these structures is $0.233$ kcal mol*-1* atom*-1* (0.23% relative error), with minimum and maximum errors of $0.072$ kcal mol*-1* atom*-1* (0.07% relative error) and $0.405$ kcal mol*-1* atom*-1* (0.39% relative error), respectively. Finally, to determine whether the wreath-shaped conformation obtained at the end of the trajectory is a local minimum on the Ala10 PES, its geometry was optimized with BFGS using PhysNet-ens5, as well as revPBE-D3(BJ)/def2-TZVP to determine the necessary energy gradients. The energies (relative to free atoms) of the optimized structures are $-11339.95$ kcal mol*-1* and $-11337.07$ kcal mol*-1* for the ab initio method and the PhysNet-ens5, respectively (which corresponds to a relative error of about $0.03$%). While both optimized geometries are structurally similar ($\mathrm{RMSD}=0.52$ Å, see Fig. 8B), the PhysNet-ens5 predicts the wreath-shaped geometry to be more stable than the helical form by about $0.184$ kcal mol*-1* atom*-1*, whereas according to the ab initio method, both structures have almost the same energy (the wreath-shaped geometry is still more stable, but only by about $0.004$ kcal mol*-1* atom*-1*). The RMSD of Ala10 with respect to the optimized wreath-shaped structure along the trajectory is shown in Fig. 9.

One of the attractive future prospects of ML-based PESs is the possibility that they become accurate and computationally efficient alternatives to empirical force fields (FFs). Such FFs are parametrized, empirical functions to describe chemical and biological systems and their parameters are fitted to a combination of experimental (e.g. X-ray structures, infrared spectroscopy, or hydration free energies) and computer-generated data (e.g. partial charges). From this perspective, it is instructive to consider the question whether an empirical FF yields similar results as PhysNet-ens5 when compared to ab initio values. For this, the energy of 20 structures sampled along the trajectory (see above) was evaluated with the CHARMM program using the CHARMM36 all atom FF. In addition, the helical and wreath-shaped structures were (re)optimized using the FF and compared to those obtained from PhysNet-ens5 and revPBE-D3(BJ)/def2-TZVP. For both FF-optimized structures, the RMSD is $0.29$ Å with respect to the structures obtained from revPBE-D3(BJ)/def2-TZVP, which is comparable to RMSDs of $0.21$ Å and $0.52$ Å for the helical and wreath-shaped structures obtained using PhysNet-ens5 (see above). To allow a meaningful comparison between FF and ab initio energies, all values are taken relative to those of the optimized wreath-shaped structure (the minimum energy according to revPBE-D3(BJ)/def2-TZVP). On average, the FF energies differ by $0.430$ kcal mol*-1* atom*-1* from the revPBE-D3(BJ)/def2-TZVP reference values, which is approximately a factor of $2$ larger than energies from PhysNet-ens5 ($0.233$ kcal mol*-1* atom*-1*). Given the faster evaluation time of an empirical FF (by a factor of $\approx$10–100), this performance appears to be acceptable. It should be noted that the standard CHARMM36-FF does not contain refinements such as multipolar interactions or polarization effects, which can increase the accuracy of FFs, but also make them less computationally efficient.

4 Discussion and Conclusion

In the present work, the PhysNet neural network architecture was introduced and tested on several common quantum-chemical benchmark datasets. It matches or improves state-of-the-art performance of machine learning models for all tested benchmarks, in some cases decreasing the error of previously published models between 50–90%. While kernel-based approaches like GDML (or its successor sGDML) achieve similar or sometimes even better performance on some benchmarks like MD17, kernel-based methods have the disadvantage that the computational cost of evaluating them scales linearly with training set size. The cost of evaluating PhysNet on the other hand is independent of the amount of data used for training and only depends on the chosen architecture. For this reason, neural network based models are better suited for constructing PESs that are transferable between many different chemical systems and require large amounts of training data (see for example the solvated protein fragments benchmark).

Two new datasets were introduced that address chemical situations not covered in other published datasets, namely chemical reactivity and many-body intermolecular interactions important for condensed phase systems. It was demonstrated that incorporating physical knowledge into the model (e.g. electrostatic contributions) can be crucial for its predictive quality.

While optimized helical structures of Ala10 using the PhysNet-PES (trained on the solvated protein fragments dataset) and the reference ab initio method are almost identical (see Fig. 8A), the relative error in the energy prediction of PhysNet-ens5 is about an order of magnitude larger than for the wreath-shaped structure of Ala10 (see Figure 8B). A possible explanation for this discrepancy could be the large dipole moment of helical protein structures due to the cumulative effect of the individual dipole moments of carbonyl groups aligned along the helix axis. The electric field associated with a large dipole moment likely leads to strong polarization effects, which potentially influence the total energy substantially. While polarization effects can be captured implicitly by PhysNet due to its ability to assign environment-dependent partial charges to atoms, it is likely that the structures included in the training data do not contain sufficient information to describe the cumulative polarization effects of multiple aligned dipole moments. A larger dataset of reference structures including helical motifs would likely be needed for a proper description of such phenomena.

In summary, PhysNet is able to accurately predict energies and forces for a wide range of structures across chemical and conformational degrees of freedom and different datasets. For SN2 reactions of methyl halides with halide anions, it was shown that including long-range electrostatic interactions explicitly in the model significantly improves the qualitative shape of the predicted PES close to and beyond the cut-off radius. Further, it was shown that PhysNet can distinguish between intra- and intermolecular contributions in SSIs and BBIs of proteins in a meaningful manner. When trained on a large set of small reference structures, the PhysNet model is able to generalize to larger structures like Ala10 with similar structural motifs. This result suggests that with a systematically constructed set of small reference structures, it is possible to build a transferable model applicable to a wide range of chemical systems. However, some large-scale effects, for example strong electric fields due to multiple aligned microscopic dipole moments, might not be properly accounted for when training only on small molecules.

Acknowledgments

The authors acknowledge financial support from the Swiss National Science Foundation (NCCR-MUST and Grant No. 200021-7117810) and the University of Basel.

Dataset and code availability

The SN2 reactions and the solvated protein fragments datasets are available from www.quantum-machine.org or directly from [\NAT@swatrue\NAT@parfalse\NAT@citexnum[][]sn2_reactions_dataset,solvated_protein_fragments_dataset]. The code for training PhysNet models is freely available from https://github.com/MeuwlyGroup/PhysNet.