Unifying machine learning and quantum chemistry with a deep neural network for molecular wavefunctions — TEI-tagged view

K. T. SchüttpersName, M. GasteggerpersName, A. TkatchenkopersName, K.-R. MüllerpersName, R. J. MaurerpersName (2019date). “Unifying machine learning and quantum chemistry with a deep neural network for molecular wavefunctions.” Nature Communicationstitle·j.

Object of study — a deep neural network predicting molecular wavefunctionsterm·key in an atomic-orbital basis

Research question — how MLterm·mth models can capture the electronic degrees of freedom that property models miss

Contribution — proof of principle that direct MLterm·mth models of electronic structureterm·key can enhance quantum chemistryterm·key

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Unifying machine learningterm·mth and quantum chemistryterm·key with a deep neural network for molecular wavefunctionsterm·key Abstract Machine learningterm·mth advances chemistry and materials science by enabling large-scale exploration of chemical space based on quantum chemical calculations. While these models supply fast and accurate predictions of atomistic chemical properties, they do not explicitly capture the electronic degrees of freedom of a molecule, which limits their applicability for reactive chemistry and chemical analysis. Here we present a deep learningterm·key framework for the prediction of the quantum mechanical wavefunctionterm·key in a local basis of atomic orbitalsterm·key from which all other ground-state properties can be derived. This approach retains full access to the electronic structureterm·key via the wavefunctionterm·key at force-field-like efficiency and captures quantum mechanics in an analytically differentiable representation. On several examples, we demonstrate that this opens promising avenues to perform inverse designterm·key of molecular structures for targeting electronic property optimisation and a clear path towards increased synergy of machine learningterm·mth and quantum chemistryterm·key.

Introduction Machine learningterm·mth (MLterm·mth) methods reach ever deeper into quantum chemistryterm·key and materials simulation, delivering predictive models of interatomic potential energy surfaces 1 – 6 , molecular forces 7 , 8 , electron densities 9 , density functionals 10 , and molecular response properties such as polarisabilities 11 , and infrared spectra 12 . Large data sets of molecular properties calculated from quantum chemistryterm·key or measured from experiment are equally being used to construct predictive models to explore the vast chemical compound space 13 – 17 to find new sustainable catalyst materials 18 , and to design new synthetic pathways 19 . Recent research has explored the potential role of machine learningterm·mth in constructing approximate quantum chemical methods 20 , as well as predicting MP2 and coupled cluster energies from Hartree–Fock orbitals 21 , 22 . There have also been approaches that use neural networks as a basis representation of the wavefunctionterm·key 23 – 25 .

Most existing MLterm·mth models have in common that they learn from quantum chemistryterm·key to describe molecular properties as scalar, vector, or tensor fields 26 , 27 . Figure 1a shows schematically how quantum chemistryterm·key data of different electronic properties, such as energies or dipole moments, is used to construct individual MLterm·mth models for the respective properties. This allows for the efficient exploration of chemical space with respect to these properties. Yet, these MLterm·mth models do not explicitly capture the electronic degrees of freedom in molecules that lie at the heart of quantum chemistryterm·key. All chemical concepts and physical molecular properties are determined by the electronic Schrödinger equation and derive from the ground-state wavefunctionterm·key. Thus, an electronic structureterm·key MLterm·mth model that directly predicts the ground-state wavefunctionterm·key (see Fig. 1b ) would not only allow to obtain all ground-state properties, but could open avenues towards new approximate quantum chemistryterm·key methods based on an interface between MLterm·mth and quantum chemistryterm·key. Hegde and Bowen 28 have explored this idea using kernel ridge regression to predict the band structure and ballistic transmission in a limited study on straining single-species bulk systems with up to μfour atomic orbitalsterm·key. Another recent example of this scheme is the prediction of coupled-cluster singles and doubles amplitudes from MP2-derived properties by Townsend and Vogiatzis 29 . Fig. 1 Synergy of quantum chemistryterm·key and machine learningterm·mth. a Forward model: MLterm·mth predicts chemical properties based on reference calculations. If another property is required, an additional MLterm·mth model has to be trained. b Hybrid model: MLterm·mth predicts the wavefunctionterm·key. All ground state properties can be calculated and no additional MLterm·mth is required. The wavefunctionsterm·key can act as an interface between MLterm·mth and QM In this work, we develop a deep learningterm·key framework that provides an accurate MLterm·mth model of molecular electronic structureterm·key via a direct representation of the electronic Hamiltonianterm·key in a local basis representation. The model provides a seamless interface between quantum mechanics and MLterm·mth by predicting the eigenvalue spectrum and molecular orbitalsterm·dom (MOsterm·dom) of the Hamiltonianterm·key for organic molecules close to ‘chemical accuracyterm·key’ (~0.04 eV). This is achieved by training a flexible MLterm·mth model to capture the chemical environment of atoms in molecules and of pairs of atoms. Thereby, it provides access to electronic properties that are important for chemical interpretation of reactions such as charge populations, bond orders, as well as dipole and quadrupole moments without the need of specialised MLterm·mth models for each property. We demonstrate how our model retains the conceptual strength of quantum chemistryterm·key by performing an ML-driven molecular dynamicsterm·key simulation of malondialdehyde showing the evolution of the electronic structureterm·key during a proton transfer while reducing the computational cost by 2–3 orders of magnitude. As we obtain a symmetry-adapted and analytically differentiable representation of the electronic structureterm·key, we are able to optimise electronic properties, such as the HOMO-LUMO gapterm·key, in a step towards inverse designterm·key of molecular structures. Beyond that, we show that the electronic structureterm·key predicted by our approach may serve as input to further quantum chemical calculations. For example, wavefunction restartsterm·mth based on this MLterm·mth model provide a significant speed-up of the self-consistent field procedureterm·mth (SCFterm·mth) due to a reduced number of iterationsterm·mth, without loss of accuracy. The latter showcases that quantum chemistryterm·key and machine learningterm·mth can be used in tandem for future electronic structureterm·key methods.

Results Atomic representation of molecular electronic structureterm·key In quantum chemistryterm·key, the wavefunctionterm·key associated with the electronic Hamiltonianterm·key Ĥ is typically expressed by anti-symmetrised products of single-electron functions or molecular orbitalsterm·dom. These are represented in a local atomic orbitalterm·key basis of spherical atomic functions ψ m = ∑ i c m i ϕ i with varying angular momentum. As a consequence, one can write the electronic Schrödinger equation in matrix form 1 Hc m = ϵ m Sc m , where the Hamiltonianterm·key matrix H may correspond to the Fock or Kohn–Sham matrix, depending on the chosen level of theory 30 . In both cases, the Hamiltonianterm·key and overlap matrices are defined as: 2 H i j = ϕ i Ĥ ϕ j and 3 S i j = ϕ i ∣ ϕ j . The eigenvalues ϵ m and electronic wavefunctionterm·key coefficients c m i contain the same information as H and S where the electronic eigenvalues are naturally invariant to rigid molecular rotations, translations or permutation of equivalent atoms. Unfortunately, as a function of atomic coordinates and changing molecular configurations, eigenvalues and wavefunctionterm·key coefficients are not well-behaved or smooth. State degeneracies and electronic level crossings provide a challenge to the direct prediction of eigenvalues and wavefunctionsterm·key with MLterm·mth techniques. We address this problem with a deep learningterm·key architecture that directly describes the Hamiltonianterm·key matrix in local atomic orbitalterm·key representation.

SchNOrbterm·mth deep learningterm·key framework SchNOrbterm·mth (SchNet for Orbitalsterm·mth) presents a framework that captures the electronic structureterm·key in a local representation of atomic orbitalsterm·key that is common in quantum chemistryterm·key. Figure 2a gives an overview of the proposed architecture. SchNOrbterm·mth extends the deep tensor neural network SchNet 31 to represent electronic wavefunctionsterm·key. The core idea is to construct symmetry-adapted pairwise features Ω i j l to represent the block of the Hamiltonianterm·key matrix corresponding to atoms i , j . They are written as a product of rotationally invariant ( λ = 0) and covariant ( λ > 0) components ω i j λ which ensures that – given a sufficiently large feature space – all rotational symmetries up to angular momentum l can be represented: 4 Ω i j l = ∏ λ = 0 l ω i j λ with 0 ≤ l ≤ 2 L 5 ω i j λ = p i j λ ⊗ 1 D for λ = 0 p i j λ ⊗ r i j r i j W λ for λ > 0 , here, r ij is the vector pointing from atom i to atom j , p i j λ ∈ R B are rotationally invariant coefficients and W λ ∈ R 3× D are learnable parameters projecting the features along D randomly chosen directions. This allows to rotate the different factors of Ω i j l ∈ R B ⋅ D relative to each other and further increases the flexibility of the model for D > 3. In case of λ = 0, the coefficients are independent of the directions due to rotational invariance. Fig. 2 Prediction of electronic properties with SchNOrbterm·mth. a Illustration of the network architecture. The neural network architecture consists of three steps (grey boxes) starting from initial representations of atom types and positions (top), continuing with the construction of representations of chemical environments of atoms and atom pairs (middle) before using these to predict energy and Hamiltonianterm·key matrix respectively (bottom). The left path through the network to the energy prediction E is rotationally invariant by design, while the right pass to the Hamiltonianterm·key matrix H allows for a maximum angular momentum L of predicted orbitals by employing a multiplicative construction of the basis ω ij using sequential interaction passes l = 0…2 L . The onsite and offsite blocks of the Hamiltonianterm·key matrix are treated separately. The prediction of overlap matrix S is performed analogously. b Illustration of the SchNet interaction block 32 . c Illustration of SchNorbterm·mth interaction block. The pairwise representation h i j l of atoms i , j is constructed by a factorised tensor layer f tensor from atomic representations as well as the interatomic distance. Using this, rotationally invariant interaction refinements v i m and basis coefficients p i j l are computed. d Loewdin population analysis for uracil based on the density matrix calculated from the predicted Hamiltonianterm·key and overlap matrices. e Mean abs. errors of lowest 20 orbitals (13 occupied + 7 virtual) of ethanol for Hartree–Fock and DFT@PBE. f The predicted (solid black) and reference (dashed grey) orbital energies of an ethanol molecule for DFT. Shown are the last four occupied and first four unoccupied orbitals, including HOMO and LUMO. The associated predicted and reference molecular orbitalsterm·dom are compared for four selected energy levels We obtain the coefficients p i j λ from an atomistic neural network, as shown in Fig. 2a . Starting from atom type embeddings x i 0 , rotationally invariant representations of atomistic environments x i T are computed by applying T consecutive interaction refinements. These are by construction invariant with respect to rotation, translation and permutations of atoms. This part of the architecture is equivalent to the SchNet model for atomistic predictions (see refs. 32 , 33 ). In addition, we construct representations of atom pairs i , j that will enable the prediction of the coefficients p i j λ . This is achieved by 2 L + 1 SchNOrbterm·mth interaction blocks, which compute the coefficients p i j λ with a given angular momentum λ with respect to the atomic environment of the respective atom pair ij . This corresponds to adapting the atomic orbitalterm·key interaction based on the presence and position of atomic orbitalsterm·key in the vicinity of the atom pair. As shown in Fig. 2b , the coefficient matrix depends on pair interactions mlp pair of atoms i , j as well as environment interactions mlp env of atom pairs ( i , m ) and ( n , j ) for neighbouring atoms m , n . These are crucial to enable the model to capture the orientation of the atom pair within the molecule for which pair-wise interactions of atomic environments are not sufficient.

The Hamiltonianterm·key matrix is obtained by treating on-site and off-site blocks separately. Given a basis of atomic orbitalsterm·key up to angular momentum L , we require pair-wise environments with angular momenta up to 2 L to describe all Hamiltonianterm·key blocks 6 H ~ i j = H off Ω i j l 0 ≤ l ≤ 2 L + 1 for i ≠ j H on Ω i m l m≠i 0 ≤ l ≤ 2 L + 1 for i = j , The predicted Hamiltonianterm·key is obtained through symmetrisation H = 1 2 ( H ~ + H ~ ⊺ ) . H off and H on are modelled by neural networks that are described in detail in the methods section. The overlap matrix S can be obtained in the same manner. Based on this, the orbital energies and coefficients can be calculated according to Eq. ( 1 ). The computational cost of the diagonalisationterm·mth is negligible for the molecules and basis sets we study here (<1 ms). For large basis sets and molecules, when the diagonalisationterm·mth starts to dominate the computational cost, our method requires only a single diagonalization instead of one per SCFterm·mth step. In addition to the Hamiltonianterm·key and overlap matrices, we predict the total energy separately as a sum over atom-wise energy contributions, in analogy with the conventional SchNet treatment 32 to drive the molecular dynamicsterm·key simulations.

Learning electronic structureterm·key and derived properties The proposed SchNOrbterm·mth architecture allows us to perform predictions of total energies, Hamiltonianterm·key and overlap matrices in end-to-end fashion using a combined regression loss. We train separate neural networks for several data sets of water as well as ethanol, malondialdehyde, and uracil from the MD17 dataset 7 . The reference calculations were performed with Hartree-Fock (HF) and density functional theory (DFT) with the PBE exchange correlation functional 34 . The employed Gaussian atomic orbitalterm·key bases include angular momenta up to l = 2 (d-orbitals). We augment the training data by adding rotated geometries and correspondingly rotated Hamiltonianterm·key and overlap matrices to learn the correct rotational symmetries (see Methods section). Detailed model and training settings for each data set are listed in Supplementary Table 1 .

As Supplementary Table 2 shows, the total energies could be predicted up to a mean absolute error below 2 meV for the molecules. The predictions show mean absolute errors below 8 meV for the Hamiltonianterm·key and below 1 × 10 −4 for the overlap matrices. We examine how these errors propagate to orbital energy and coefficients. Figure 2e shows mean absolute errors for energies of the lowest 20 molecular orbitalsterm·dom for ethanol reference calculations using DFT as well as HF. The errors for the DFT reference data are consistently lower. Beyond that, the occupied orbitals (1–13) are predicted with higher accuracy (<20 meV) than the virtual orbitals (~100 meV). We conjecture that the larger error for virtual orbitals arises from the fact that these are not strictly defined by the underlying data from the HF and Kohn–Sham DFT calculations. Virtual orbitals are only defined up to an arbitrary unitary transformation. Their physical interpretation is limited and, in HF and DFT theory, they do not enter in the description of ground-state properties. For the remaining data sets, the average errors of the occupied orbitals are <10 meV for water and malondialdehyde, as well as 48 meV for uracil. This is shown in detail in Supplementary Fig. 1 . The orbital coefficients are predicted with cosine similarities ≥90% (see Supplementary Fig. 2 ). Figure 2f depicts the predicted and reference orbital energies for the frontier MOsterm·dom of ethanol (solid and dotted lines, respectively), as well as the orbital shapes derived from the coefficients. Both occupied and unoccupied energy levels are reproduced with high accuracy, including the highest occupied (HOMO) and lowest unoccupied orbitals (LUMO). This trend is also reflected in the overall shape of the orbitals. Even the slightly higher deviations in the orbital energies observed for the third and fourth unoccupied orbital only result in minor deformations. The learned covariance of molecular orbitalsterm·dom for rotations of a water molecule is shown in Supplementary Fig. 3 .

The MLterm·mth model uses about 93 million parameters to predict a large Hamiltonianterm·key matrix with >100 atomic orbitalsterm·key. This size is comparable to state-of-the-art neural networks for the generation of similarly sized images 35 . Supplementary Table 6 shows the computational costs of calculating the reference data, training the network and predicting Hamiltoniansterm·key. While training of SchNOrbterm·mth took about 80 h, performing the required DFT reference calculations remains the bottleneck for obtaining a trained network, in particular for larger molecules. Our approach to predicting Hamiltonianterm·key matrices leads to accelerations of 2–3 orders of magnitude.

As SchNOrbterm·mth learns the electronic structureterm·key of molecular systems, all chemical properties that are defined as quantum mechanical operators on the wavefunctionsterm·key can be computed from the MLterm·mth prediction without the need to train a separate model. We investigate this feature by directly calculating electronic dipole and quadrupole moments from the orbital coefficients predicted by SchNOrbterm·mth, as well as the HF total energies for the ethanol molecule. The corresponding mean absolute errors are reported in Supplementary Tables 4 and 5 . The calculation of energies and forces from the coefficients requires the evaluation of the core Hamiltonianterm·key (HF) or exchange correlation terms (DFT), for which we currently resort to the ORCA code. To avoid this computational overhead and obtain highly accurate predictions for molecular dynamicsterm·key simulations with mean absolute error below 1 meV, we predict energies and forces directly as a sum of atomic contributions 31 . Regarding the electrostatic moments, excellent agreement with the electronic structureterm·key reference is observed for the majority of molecules (<0.054 D for dipoles and <0.058 D Å for quadrupoles). The only deviation from this trend is observed for uracil, where a loss function minimising only the errors of Hamiltonianterm·key and overlap matrices is too limited. The dipole moment depends strongly on the molecular electron density derived from the orbital coefficients, which are never learned directly.

Beyond that, we have studied the prediction accuracy for ethanol when using the larger def2-tzvp basis set which includes f-orbitals ( l = 3). While the predictions of the Hamiltonianterm·key and overlap matrices remain remarkably accurate with 8.3 meV and 10 −6 , respectively, the derived properties exhibit large errors, e.g., an MAE of 0.4775 eV for the orbital energies. For large numbers of orbitals, errors in the Hamiltonianterm·key can accumulate due to the diagonalisationterm·mth. This problem could be solved by improving the neural network architecture to further reduce the prediction errorterm·mth or introducing a density dependent term into the loss function, which will be explored in future investigations.

In this case, a similar accuracy as the other methods could in principle be reached upon the addition of more reference data points. The above results demonstrate the utility of combining a learned Hamiltonianterm·key with quantum operators. This makes it possible to access a wide range of chemical properties without the need for explicitly developing specialised neural network architectures.

Chemical insights from electronic deep learningterm·key Recently, a lot of research has focused on explaining predictions of MLterm·mth models 36 – 38 aiming both at the validation of the model 39 , 40 as well as the extraction of scientific insight 17 , 31 , 41 . However, these methods explain MLterm·mth predictions either in terms of the input space, atom types and positions in this case, or latent features such as local chemical potentials 31 , 42 . In quantum chemistryterm·key however, it is more common to analyse electronic properties in terms of the MOsterm·dom and properties derived from the electronic wavefunctionterm·key, which are direct output quantities of the SchNOrbterm·mth architecture.

Molecular orbitalsterm·dom encode the distribution of electrons in a molecule, thus offering direct insights into its underlying electronic structureterm·key. They form the basis for a wealth of chemical bonding analysis schemes, bridging the gap between quantum mechanics and abstract chemical concepts, such as bond orders and atomic partial charges 30 . These quantities are invaluable tools in understanding and interpreting chemical processes based on molecular reactivity and chemical bonding strength. As SchNOrbterm·mth yields the MOsterm·dom, we are able to apply population analysis to our MLterm·mth predictions. Figure 2d shows Loewdin partial atomic charges and bond orders for the uracil molecule. Loewdin charges provide a chemically intuitive measure for the electron distribution and can e.g., aid in identifying potential nucleophilic or electrophilic reaction sites in a molecule. The negatively charged carbonyl oxygens in uracil, for example, are involved in forming RNA base pairs. The corresponding bond orders provide information on the connectivity and types of bonds between atoms. In the case of uracil, the two double bonds of the carbonyl groups are easily recognisable (bond order 2.12 and 2.14, respectively). However, it is also possible to identify electron delocalisation effects in the pyrimidine ring, where the carbon double bond donates electron density to its neighbours. A population analysis for malondialdehyde, as well as population prediction errorsterm·mth for all molecules can be found in Supplementary Fig. 4 and Supplementary Table 3 .

The SchNOrbterm·mth architecture enables an accurate prediction of the electronic structureterm·key across molecular configuration space, which provides for rich chemical interpretation during molecular reaction dynamics. Figure 3a shows an excerpt of a molecular dynamicsterm·key simulation of malondialdehyde that was driven by atomic forces predicted using SchNOrbterm·mth. It depicts the proton transfer together with the relevant MOsterm·dom and the electronic density. Supplementary Video 1 shows a side-by-side comparison between the predicted and reference HOMO-2 orbital during this excerpt of the trajectory. The density paints an intuitive picture of the reaction as it migrates along with the hydrogen. This exchange of electron density during proton transfer is also reflected in the orbitals. Their dynamical rearrangement indicates an alternation between single and double bonds. The latter effect is hard to recognise based on the density alone and demonstrates the wealth of information encoded in the molecular wavefunctionsterm·key. Fig. 3 Proton transfer in malondialdehyde. a Excerpt of the MD trajectory showing the proton transfer, the electron density as well as the relevant MOsterm·dom HOMO-2 and HOMO-3 for three configurations (I, II, III). b Forces exerted by the MOsterm·dom on the transferred proton for configurations I and II. c Density of states broadened across the proton transfer trajectory. MO energies of the equilibrium structure are indicated by grey dashed lines. The inset shows a zoom of HOMO-2 and HOMO-3 Figure 3b depicts the forces the different MOsterm·dom exert onto the hydrogen atom exchanged during the proton transfer. All forces are projected onto the reaction coordinate, where positive values correspond to a force driving the proton towards the product state. In the initial configuration I, most forces lead to attraction of the hydrogen atom to the right oxygen. In the intermediate configuration II, orbital rearrangement results in a situation where the majority of orbital force contributions on the hydrogen atom become minimal, representing mostly non-bonding character between oxygens and hydrogen. One exception is MO 13, depicted in the inset of Fig. 3b . Due to a minor deviation from a symmetric O–H–O arrangement, the orbital represents a one-sided O–H bond, exerting forces that promote the reaction. The intrinsic fluctuations during the proton transfer molecular dynamicsterm·key are captured by the MOsterm·dom as can be seen in Fig. 3c . This shows the distribution of orbital energies encountered during the reaction. As would be expected, both HOMO-2 and HOMO-3 (inset, orange and blue respectively), which strongly participate in the proton transfer, show significantly broadened peaks due to strong energy variations in the dynamics. This example nicely shows the chemically intuitive interpretation that can be obtained by the electronic structureterm·key prediction of SchNOrbterm·mth.

Deep learning-enhanced quantum chemistryterm·key An essential paradigm of chemistry is that the molecular structure defines chemical properties. Inverse chemical design turns this paradigm on its head by enabling property-driven chemical structure exploration. The SchNOrbterm·mth framework constitutes a suitable tool to enable inverse chemical design due to its analytic representation of electronic structureterm·key in terms of the atomic positions. We can therefore obtain analytic derivatives with respect to the atomic positions, which provide the ability to optimise electronic properties. Figure 4a shows the minimisation and maximisation of the HOMO-LUMO gapterm·key ϵ gap of malondialdehyde as an example. We perform gradient descent and ascent from a randomly selected configuration r ref until convergence at r min and r max , respectively. We are able to identify structures which minimise and maximise the gap from its initial 3.15 to 2.68 eV at r min and 3.59 eV at r max . While in this proof of concept these changes were predominantly caused by local deformations in the carbon–carbon bonds indicated in Fig. 4a , they present an encouraging prospect how electronic surrogate models such as SchNOrbterm·mth can contribute to computational chemical design using more sophisticated optimisation methods, such as alchemical derivatives 43 or reinforcement learning 44 . Fig. 4 Applications of SchNOrbterm·mth. a Optimisation of the HOMO-LUMO gapterm·key. HOMO and LUMO with energy levels are shown for a randomly drawn configuration of the malonaldehyde dataset (centre) as well as for configurations that were obtained from minimising or maximising the HOMO-LUMO gapterm·key prediction using SchNOrbterm·mth (left and right, respectively). For the optimised configurations, the difference of the orbitals are shown in green (increase) and violet (decrease). The dominant geometrical change is indicated by the black arrows. b The predicted MO coefficients for the uracil configurations from the test set are used as a wavefunctionterm·key guess to obtain accurate solutions from DFT at a reduced number of self-consistent-field (SCFterm·mth) iterations. This reduces the required SCFterm·mth iterations by an average of 77% using a Newton solver. In terms of runtime, it is more efficient to use SOSCF, even though this saves only 15% of iterations for uracil MLterm·mth applications for electronic structureterm·key methods have usually been one-directional, i.e., MLterm·mth models are trained to predict the outputs of calculations. On the other hand, models in the spirit of Fig. 1b , such as SchNOrbterm·mth, offer the prospect of providing a deeper integration with quantum chemistryterm·key methods by substituting parts of the electronic structureterm·key calculation. SchNOrbterm·mth directly predicts wavefunctionsterm·key based on quantum chemistryterm·key data, which in turn, can serve as input for further quantum chemical calculations. For example, in the context of HF or DFT calculations, the relevant equations are solved via a self-consistent field approach (SCFterm·mth) that determines a set of MOsterm·dom. The convergence with respect to SCFterm·mth iteration steps largely determines the computational speed of an electronic structureterm·key calculation and strongly depends on the quality of the initialisation for the wavefunctionterm·key. The coefficients predicted by SchNOrbterm·mth can serve as such an initialisation of SCFterm·mth calculations. To this end, we generated wavefunctionterm·key files for the ORCA quantum chemistryterm·key package 45 from the predicted SchNOrbterm·mth coefficients, which were then used to initialise SCFterm·mth calculations. Figure 4b depicts the SCFterm·mth convergence for three sets of computations on the uracil molecule: using the standard initialisation techniques of quantum chemistryterm·key codes, and the SchNorbterm·mth coefficients with or without a second order solver. Nominally, only small improvements are observed using SchNorbterm·mth coefficients in combination with a conventional SCFterm·mth solver. This is due to the various strategies employed in electronic structureterm·key codes in order to provide a numerically robust SCFterm·mth procedure. By performing SCFterm·mth calculations with a second order solver, which would not converge using a less accurate starting point than our SchNorbterm·mth MO coefficients, the efficiency of our combined MLterm·mth and second order SCFterm·mth approach becomes apparent. Convergence is obtained in only a fraction of the original iterations, reducing the number of cycles by ~77%. Similarly, Supplementary Fig. 5 shows the reduction of SCFterm·mth iteration by ~73% for malondialdehyde. However, since second-order optimisation steps are more costly, it is more time-efficient to perform conventional SOSCF which reduces the convergence time by 13% and 16% for uracil and malondialdehyde, respectively.

It should be noted, that this combined approach does not introduce any approximations into the electronic structureterm·key method itself and yields exactly the same results as the full computation. Another example of integration of the SchNOrbterm·mth deep learningterm·key framework with quantum chemistryterm·key, as shown in Fig. 1b , is the use of predicted wavefunctionsterm·key and MO energies based on Hartree–Fock as starting point for post-Hartree–Fock correlation methods such as Møller–Plesset perturbation theory (MP2). Supplementary Table 5 presents the mean absolute error of an MP2 calculation for ethanol based on wavefunctionsterm·key predicted from SchNOrbterm·mth. The associated prediction errorterm·mth for the test set is 83 meV. Compared to the overall HF and MP2 energy, the relative error of SchNOrbterm·mth amounts to 0.01% and 0.06%, respectively. For the MP2 correlation energy, we observe a deviation of 17%, the reason of which is inclusion of virtual orbitals in the calculation of the MP2 integrals. However, even in this case, the total error only amounts to a deviation of 93 meV.

Discussion The SchNOrbterm·mth framework provides an analytical expression for the electronic wavefunctionsterm·key in a local atomic orbitalterm·key representation as a function of molecular composition and atom positions. While previous approaches have predicted Hamiltoniansterm·key of single-species bulk materials in a small basis set for limited geometric deformations 28 , SchNorbterm·mth has been shown to enable the accurate predictions of molecular Hamiltoniansterm·key in a basis of >100 atomic orbitalsterm·key up to angular momentum l = 2 for a much larger configuration space obtained from molecular dynamicsterm·key simulations. As a consequence, the model provides access to atomic derivatives of wavefunctionsterm·key, which include molecular orbitalterm·dom energy derivatives, Hamiltonianterm·key derivatives, which can be used to approximate nonadiabatic couplings 46 , as well as higher order derivatives that describe the electron-nuclear response of the molecule. Thus, the SchNOrbterm·mth framework preserves the benefits of interatomic potentials while enabling access to the electronic structureterm·key as predicted by quantum chemistryterm·key methods.

SchNOrbterm·mth opens up completely new applications to ML-enhanced molecular simulation. This includes the construction of interatomic potentials with electronic properties that can facilitate efficient photochemical simulations during surface hopping trajectory dynamics or Ehrenfest-type mean-field simulations, but also enables the development of new ML-enhanced approaches to inverse molecular design via electronic property optimisation.

The SchNOrbterm·mth neural network architecture demonstrates that an accurate prediction of electronic structureterm·key is feasible. However, with increasing number of atomic orbitalsterm·key, the diagonalisationterm·mth of the Hamiltonianterm·key leads to the accumulation of prediction errorsterm·mth and becomes the bottleneck of the prediction. Thus, for larger molecules and basis sets, the accuracy of SchNOrbterm·mth will have to be further improved. More research on the architecture is also required in order to reduce the required amount of training data and parameters, e.g., by adding more prior knowledge to the model. An import step into this direction is to encode the full rotational symmetries of the basis into the architecture, replacing the current data augmentation scheme. Alternatively, on the basis of the SchNOrbterm·mth framework, intelligent preprocessing of quantum chemistryterm·key data in the form of effective Hamiltoniansterm·key or optimised minimal basis representations 47 can be developed in the future. Such preprocessing will also pave the way towards the prediction of the electronic structureterm·key based on post-HF correlated wavefunctionterm·key methods and post-DFT quasiparticle methods.

This work serves as a first proof of principle that direct MLterm·mth models of electronic structureterm·key based on quantum chemistryterm·key can be constructed and used to enhance further quantum chemistryterm·key calculations. We have presented an immediate consequence of this by reducing the number of DFT-SCF iterations with wavefunctionsterm·key predicted via SchNOrbterm·mth. The presented model delivers derived electronic properties that can be formulated as quantum mechanical expectation values. This provides an important step towards a full integration of MLterm·mth and quantum chemistryterm·key into the scientific discovery cycle.

Methods Reference data Reference configurations were sampled at random from the MD17 dataset 7 for each molecule. The number of selected configurations per molecule is given in Supplementary Table 1 . All reference calculations were carried out with the ORCA quantum chemistryterm·key code 45 using the def2-SVP basis set 48 . Integration grid levels of 4 and 5 were employed during SCFterm·mth iterations and the final computation of properties, respectively. Unless stated otherwise, the default ORCA SCFterm·mth procedure was used, which is based on the Pulay method 49 . For the remaining cases, the Newton–Raphson procedure implemented in ORCA was employed as a second order SCFterm·mth solver. SCFterm·mth convergence criteria were set to VeryTight. DFT calculations were carried out using the PBE functional 34 . For ethanol, additional HF computations were performed.

Molecular dynamicsterm·key simulations for malondialdehyde were carried out with SchNetPack 50 . The equations of motions were integrated using a timestep of 0.5 fs. Simulation temperatures were kept at 300 K with a Langevin thermostat 51 employing a time constant of 100 fs. Trajectories were propagated for a total of 50 ps, of which the first 10 ps were discarded.

Details on the neural network architecture In the following we describe the neural network depicted in Fig. 2 in detail. We use shifted softplus activation functions 7 ssp ( x ) = ln 1 2 e x + 1 2 throughout the architecture. Linear layers are written as 8 linear ( x ) = W ⊺ x + b with input x ∈ R n in , weights W ∈ R n in × n out and bias b ∈ R n out . Fully-connected neural networks with one hidden layer are written as 9 mlp ( x ) = W 2 ⊺ ssp W 1 ⊺ x + b 1 + b 2 with weights W 1 ∈ R n in × n hidden and W 2 ∈ R n hidden × n out and biases b 1 , b 2 accordingly. Model parameters are shared within layers across atoms and interactions, but never across layers. We omit layer indices for clarity.

The representations of atomic environments are constructed with the neural network structure as in SchNet. In the following, we summarise this first part of the model. For further details, please refer to SchüttpersName et al. 32 . First, each atom is assigned an initial element-specific embedding 10 x i 0 = a Z i ∈ R B , where Z i is the nuclear charge and B is the number of atom-wise features. In this work, we use B = 1000 for all models. The representations are refined using SchNet interaction layers (Fig. 2b ). The main component is a continuous-filter convolutional layer (cfconv) 11 cfconv ( ( x 1 , r 1 ) , … , ( x n , r n ) ) = ∑ j ≠ i x j ∘ W filter ( r i j ) i = 1 … n which takes a spatial filter W filter : R → R B 12 W filter ( r i j ) = mlp ( g ( r i j ) ) f cutoff ( r i j ) with 13 g ( r i j ) = exp ( − γ ( r i j − k Δ μ ) 2 ) 0 ≤ k ≤ r c ∕ Δ μ 14 f cutoff ( r i j ) = 0.5 × 1 + cos π r r c r < r c 0 r ≥ r c , where r c is the cutoff radius and Δ μ is the grid spacing of the radial basis function expansion of interatomic distance r ij . While this adds spatial information to the environment representations for each feature separately, the crosstalk between features is performed atom-wise by fully-connected layers (linear and mlp atom in Fig. 2b ) to obtain the refinements v i t , where t is the current interaction iteration. The refined atom representations are then 15 x i t = x i t − 1 + v i t . These representations of atomic environments are employed by SchNet to predict chemical properties via atom-wise contributions. However, in order to extend this scheme to the prediction of the Hamiltonianterm·key, we need to construct representations of pair-wise environments in a second interaction phase.

The Hamiltonianterm·key matrix is of the form 16 H = H 11 ⋯ H 1 j ⋯ H 1 n ⋮ ⋱ ⋮ ⋮ H i 1 ⋯ H i j ⋯ H i n ⋮ ⋮ ⋱ ⋮ H n 1 ⋯ H n j ⋯ H n n where a matrix block H i j ∈ R n ao , i × n ao , j depends on the atoms i , j within their chemical environment as well as on the choice of n ao, i and n ao, j atomic orbitalsterm·key, respectively. Therefore, SchNOrbterm·mth builds representations of these embedded atom pairs based on the previously constructed representations of atomic environments. This is achieved through the SchNOrbterm·mth interaction module (see Fig. 2c ).

First, a raw representation of atom pairs is obtained using a factorised tensor layer 31 , 52 : 17 h i j λ = ftensor ( x i λ , x j λ , r i j ) = ssp ( linear 2 [ linear 1 ( x i λ ) ∘ linear 1 ( x j λ ) ∘ W filter ( r i j ) ] ) . The layers linear 1 : R B ↦ R B map the atom representations to the factors, while the filter-generating network W filter : R → R B is defined analogously to Eq. ( 12 ) and directly maps to the factor space. In analogy to how the SchNet interactions are used to build atomic environments, SchNOrbterm·mth interactions are applied several times, where each instance further refines the atomic environments of the atom pairs with additive corrections 18 x i T + λ = x i T + λ − 1 + v i T + λ 19 v i T + λ = mlp atom ∑ j ≠ i h i j l as well as constructs pair-wise features: 20 p i j λ = p i j λ , pair + ∑ m ≠ i p m j λ , env + ∑ n ≠ j p i n λ , env 21 p i j λ , pair = mlp pair ( h i j λ ) 22 p i j λ , env = mlp env ( h i j λ ) where mlp pair models the direct interactions of atoms i , j while mlp env models the interactions of the pair with neighbouring atoms. As described above, the atom pair coefficients are used to form a basis set ω i j λ = p i j 0 ⊗ 1 D for λ = 0 p i j λ ⊗ r i j r i j W for λ > 0 where λ corresponds to the angular momentum channel and W ∈ R 3× D are learnable parameters to project along D directions. For all results in this work, we used D = 4. For interactions between s-orbitals, we consider the special case λ = 0 where the features along all directions are equal due to rotational invariance. At this point, ω i j 0 is rotationally invariant and ω i j λ > 0 is covariant. On this basis, we obtain features with higher angular momenta using: Ω i j l = ∏ λ = 0 l ω i j λ with 0 ≤ l ≤ 2 L , where features Ω i j l possess angular momentum l . The SchNOrbterm·mth representations of atom pairs embedded in their chemical environment, that were constructed from the previously constructed SchNet atom-wise features, will serve in a next step to predict the corresponding blocks of the ground-state Hamiltonianterm·key.

Finally, we assemble the Hamiltonianterm·key and overlap matrices to be predicted. Each atom pair block is predicted from the corresponding features Ω i j l : H ~ i j = H off Ω i j l 0 ≤ l ≤ 2 L + 1 for i ≠ j H on Ω i m l m≠i 0 ≤ l ≤ 2 L + 1 for i = j , where we restrict the network to linear layers in order to conserve the angular momenta: H off ( ⋅ ) = ∑ l linear off l Ω i j l [ : n ao , i , : n ao , j ] H on ( ⋅ ) = ∑ j , l linear on l Ω i j l [ : n ao , i , : n ao , i ] with linear off , linear on : R 2 L + 1 × D → R n ao,max × n ao,max , i.e., mapping to the maximal number of atomic orbitalsterm·key in the data. Then, a mask is applied to the matrix block to yield only H ~ i j ∈ R n ao , i × n ao , j . Finally, we symmetrise the prediced Hamiltonianterm·key: 23 H = 1 2 ( H ~ + H ~ ⊺ ) The overlap matrix is obtained similarly with blocks S off ⋅ = linear S,on ( Ω i j ) [ : n ao , i , : n ao , j ] S i i = S Z i . The prediction of the total energy is obtained analogously to SchNet as a sum over atom-wise energy contributions: E = ∑ i mlp E ( x i ) .

Data augmentation While SchNOrbterm·mth constructs features Ω i j l and Ω i i l with angular momenta such that the Hamiltonianterm·key matrix can be represented as a linear combination of those, it does not encode the full rotational symmetry a priori. However, this can be learned by SchNOrbterm·mth assuming the training data reflects enough rotations of a molecule. To save computing power, we reduce the amount of reference calculations by randomly rotating configurations before each training epoch using Wigner D rotation matrices 53 . Given a randomly sampled rotor R , the applied transformations are 24 r ~ i = D ( 1 ) ( R ) r i 25 F ~ i = D ( 1 ) ( R ) F i 26 H ~ μ ν = D ( l μ ) ( R ) H μ ν D ( l ν ) ( R ) 27 S ~ μ ν = D ( l μ ) ( R ) S μ ν D ( l ν ) ( R ) for atom positions r i , atomic forces F i , Hamiltonianterm·key matrix H , and overlap S .

Neural network training For the training, we used a combined loss to train on energies E , atomic forces F , Hamiltonianterm·key H and overlap matrices S simultaneously: ℓ ( H ~ , S ~ , Ẽ ) , ( H , S , E , F ) = H − H ~ F 2 + S − S ~ F 2 + ρ E − Ẽ 2 28 + 1 − ρ n atoms ∑ i = 0 n atoms F i − − ∂ Ẽ ∂ r i 2 where the variables marked with a tilde refer to the corresponding predictions and ρ determines the trade-off between total energies and forces. The neural networks were trained with stochastic gradient descent using the ADAM optimiser 54 . We reduced the learning rate using a decay factor of 0.8 after t patience epochs without improvement of the validation loss. The training is stopped at lr ≤ 5 × 10 −6 . The mini-batch sizes, patience and data set sizes are listed in Supplementary Table 1 . Afterwards, the model with lowest validation error is selected for testing.

Supplementary information Supplementary Information Peer Review File Description of Additional Supplementary Files Supplementary Video 1