Machine Learning of Dynamic Electron Correlation ... - ACS Publications

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Machine Learning of Dynamic Electron Correlation Energies from Topological Atoms James L. McDonagh,† Arnaldo F. Silva,† Mark A. Vincent,‡ and Paul L. A. Popelier*,†,‡ †

Manchester Institute of Biotechnology, The University of Manchester, 131 Princess Street, Manchester M1 7DN, Great Britain School of Chemistry, The University of Manchester, Oxford Road, Manchester M13 9PL, Great Britain

‡

S Supporting Information *

ABSTRACT: We present an innovative method for predicting the dynamic electron correlation energy of an atom or a bond in a molecule utilizing topological atoms. Our approach uses the machine learning method Kriging (Gaussian Process Regression with a non-zero mean function) to predict these dynamic electron correlation energy contributions. The true energy values are calculated by partitioning the MP2 two-particle density-matrix via the Interacting Quantum Atoms (IQA) procedure. To our knowledge, this is the first time such energies have been predicted by a machine learning technique. We present here three important proof-of-concept cases: the water monomer, the water dimer, and the van der Waals complex H2···He. These cases represent the final step toward the design of a full IQA potential for molecular simulation. This final piece will enable us to consider situations in which dispersion is the dominant intermolecular interaction. The results from these examples suggest a new method by which dispersion potentials for molecular simulation can be generated. such purely mathematical constructs are not required,30 and indeed, they are not if one works with topological atoms,31,32 as we do here. These are essentially regions that naturally appear in a system’s electron density without using any parameters. They are space-filling, which means that there are no gaps between them, and they do not overlap. The latter property is important because topological atoms do not penetrate each other (see ref 33 and Figure S1 in the Supporting Information (SI)), and hence, no associated damping function is required. The space-filling nature of the topological atoms also avoids the well-known polarization catastrophe (occurring with overlapping charge clouds) and, hence, avoids the corresponding damping functions. Consequently, when using topological atoms, there is no need to work out how such damping functions (for either penetration or polarization catastrophe) change with atom type. Machine learning (ML) is a popular tool in many areas of science and has found a variety of applications in chemistry.34−40 Here, we present a novel application in which models capable of accurate predictions of dynamic electron correlation energy are employed. This means that, in exchange for a relatively small number of initial expensive correlation calculations, we can make accurate predictions of unknown correlation energies for related geometries in a fraction of the

1. INTRODUCTION Dynamic electron correlation is a quantum mechanical phenomenon with profound ramifications for intermolecular interactions.1−3 Typical quantum chemical methods commonly employed, such as DFT, generally do not completely include many of these correlation effects, and in some cases they are not included at all. To include them completely requires expensive quantum mechanical methods. However, a plethora of methods now exists to obtain Electron Correlation Energies (ECE). These methods range from empirical potentials in simulation methods4−12 to a posteriori corrections in DFT or ab initio wave function methods.13−18 In addition, more recently, many-body correction terms have been employed to supplement DFT.16,19 These methods have been applied to a wide variety of important chemical systems including the solid state,20−22 biological materials,23 and more recently nanostructures.24−26 Recent work has seen much improvement in methods16,27 to account for ECE. Although electron correlation is proportionally a small contribution energetically, it is ubiquitous and has been shown to be responsible for macroscopic phenomena, such as an explanation of the gecko’s ability to adhere to a variety of surfaces.28 Many of these proposed corrections rely upon the use of socalled damping functions to attenuate the interactions between two atoms, which would become unrealistically large without such a mathematical device. The use of these damping functions has been questioned (Section 2.11 in ref 29). Ideally, © XXXX American Chemical Society

Received: November 15, 2017

A

DOI: 10.1021/acs.jctc.7b01157 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

capable of analyzing the geometrical deviation of a molecule over a set of GAUSSIAN0968 input files and certain molecular dynamics trajectory files. 2.2. Correlation Energy. For each of TYCHE’s generated geometries the ab initio program GAUSSIAN0968 provides the correlated part of the two-particle reduced density-matrix (2PDM), arising from an MP2/uncon-6-31++G(d,p) calculation, where uncon emphasizes that the basis set is uncontracted. We stress that the 6-31++G(d,p) basis set is uncontracted and thus provides a better description of the systems considered here compared to the normal contracted form. The 2PDM, which is a four-dimensional matrix, is then normalized to the primitive atomic orbital basis, and its elements are formally denoted dcorr,prim . A detailed derivation of jklm this matrix is provided in Section 2 of the SI. The fourdimensional matrix dcorr,prim is provided as input for the IQA63 jklm method as implemented in the in-house code MORFI, a local derivative of MORPHY.69 IQA is a quantum mechanically rigorous prescription to partition the electron density, and hence the energy, of a molecule (or complex) into topological atoms.32,70,71 The electron correlation energy VAB ee,corr between two topological atoms A and B is calculated via a sixdimensional quadrature, implemented as a succession of two three-dimensional integrals carried out over the respective volumes Ω of atoms A and B:

time. This method holds the promise of a new way to derive van der Waals energies for use in simulations. In order to generate such machine learning models, the chemical system must be appropriately represented. There are many approaches to representing a chemical system for input to machine learning methods.41−44 Here, we use our own method,44 which employs the Atomic Local Frame (ALF) approach. This technique defines a local axis system, in spherical coordinates, centered on an atom of interest. This leads to a relative portrayal of the system’s geometry as viewed from the atom of interest. An ALF is generated for each atom of the system, and the respective atomic ECE is then linked with it as a target “true” value, for that given geometry. The ML models task is to find a relationship between the ALF and the ECE value for each geometry. Our ultimate aim is to apply this method to the novel force field FFLUX29 (formerly called QCTFF45), which we are currently developing. FFLUX, is a novel next-generation force field based upon a ML method called Kriging.46,47 This approach has previously showcased the ability of such Kriging models to predict, with promising accuracy, the Coulombic and exchange energies for chemical systems in a variety of molecular arrangements. Kriging models have been applied successfully on a variety of systems, including ethanol,48 (peptide-capped) alanine,49 the microhydrated sodium ion,49 N-methylacetamide (NMA) and histidine,44 the four aromatic (peptide-capped) amino acids,50 all naturally occurring amino acids,51 helical deca-alanines,52,53 water clusters,54 cholesterol,55 and carbohydrates.56 Although the roots of FFLUX lie in electrostatics,57 particularly multipolar58 electrostatics,59 nonelectrostatic energy contributions were also successfully trained for using Kriging: for example, the atomic kinetic energy,60 which is well-defined61 for topological atoms.62 All energies are determined for each atom in a molecule of interest by the Interacting Quantum Atoms63 (IQA) approach, which divides a molecule up into its constituent atoms based on the electron density’s gradient vector field. Subsequently, IQA was employed to obtain both the interatomic and the intra-atomic exchange energies, which were also successfully Kriged64 for the B3LYP functional.65 In the current article we present the final energy contribution needed in order to complete the full description of the FFLUX force field for a chemical system’s dynamics: here we show how the MP2 correlation energy can be partitioned and successfully predicted by ML methods. We will present results for both intramolecular interactions (i.e., covalent) and intermolecular interactions (i.e., noncovalent including a hydrogen bond).

NG

AB Vee,corr =

j

j=1 k=1

∫Ω

d r2

NG

j

B

AB Vee,corr =

NG

l

corr,prim ∑ ∑ kjk ∑ ∑ klmdjklm ∫

ΩA

l=1 m=1

d r1Gjk(r1 − R jk) ×

1 Glm(r2 − R lm) r12 NG

l

corr,prim corr E ΩA, jk , ΩB, lm ∑ ∑ kjk ∑ ∑ klmdjklm j=1 k=1

l=1 m=1

(1)

The fact that the product of two Gaussian functions in eq 1 is a third Gaussian function has been employed to reduce the number of Gaussians to two from the four primitive ones, corresponding to j, k, l, and m indices. Note that A and B can refer to the same atom, in which case eq 1 returns the ECE within a single atom (i.e., the intra-atomic dynamic electron correlation). In eq 1, kjk and klm are prefactors from the producing of the primitive Gaussian basis functions. The upper bounds in the summations over indices j and k, or over l and m, mean that only the lower triangle and the diagonal terms of dcorr,prim need to be considered. Full details of this method are jklm available elsewhere.30,72 2.3. Machine Learning (ML). After having partitioned the MP2 energy, which arises from each of the TYCHE generated geometries, a training set file is produced. This file contains geometries as input and energies as output. In particular, this file defines the molecule in terms of ML input variables, which are called features, and are in this case the ALF coordinates.49 Because we are only interested in training for internal geometry changes, there are 3N−6 features for nonlinear systems, where N is the number of atoms. The training set files are automatically generated using a Python script (see Section 3 of the SI) and in-house programs. The three ML methods Kriging,44 Random Forest (RF),73 and Support Vector Regression (SVR)74,75 are employed in the present study and are discussed in detail in Section 4 of the SI. Briefly, Kriging, which is the machine learning (ML) method of

2. METHODS 2.1. Sampling. In order to predict atomic correlation energies, as a function of molecular geometry, coordinates and the MP2 energy must be provided to train and test the model. The geometries of both the training sets and test sets are generated via an in-house sampling program called TYCHE.56,66 In the current work we use TYCHE’s normal mode option. The details of this process have been discussed elsewhere.56,66 Briefly, TYCHE applies a fixed temperature as an energy source and a Boltzmann weighting system over the available modes thereby populating the accessible normal modes. A large number of geometries can be efficiently generated by this method. The local program HERMES assesses the geometrical deviations and produces visualization files. HERMES is available on GitHub from ref 67 and is B


Article


molecular energy into inter- and intra-atomic contributions. To provide a sense of scale, Table 1 gives the partitioned and total ECE of an arbitrarily distorted water molecule. The SI gives an overlay of all the distorted water geometries (Figure S4).

choice as an interpolator, maps a set of input data to an output response value, via a Gaussian process conditioned by a kernel function. This method is implemented76 in our program FEREBUS. While Kriging has several variants, we have applied only one form and that is simple Kriging. In simple Kriging the output value (y) depends on a set of n input points x and is modeled by the mean of the output (μ) plus an “error” or deviation (ε) term, which relates to, and is explicitly dependent on, the geometry-specific effects (marked by x) on the ECE (or y). Second, RF is an ensemble learning method based upon the concept of decision trees. A forest of decision trees is generated, in which each tree is “grown” using a subset of the training data and features. The nodes, branching these trees, are decided on based upon minimizing an error metric, which is the Root Mean Square Error (RMSE) in this case. Finally, SVR searches for an optimal regression within a margin of acceptable error (ε) in a higher order feature space via the “kernel trick”. The RF and SVR codes are used as implemented in the R package, Caret.77 Hyperparameter information for all models is provided in Section 4 of the SI, and final Kriging hyperparameters for each model are reported in the SI in Section 5. Note that due to the different training methods and hyperparameter optimization techniques the Kriging results are not directly comparable with the RF and SVR results. We have chosen the RF and SVR algorithms, with a 10-fold cross-validation, as these provide examples of widely accessible methods. Previously, such methods have been successfully applied to property prediction.78−82

Table 1. Partitioned and Total ECE Calculated by MORFI for a Distorted Water Molecule with an MP2/uncon-6-31+ +G(d,p) Wave Functiona intra-atomic

interatomic

total

O

H1

H2

O−H1

O−H2

H1−H2

H2O

−521.0

−22.5

−25.6

15.9

18.4

−5.1

−539.8

a

All values are quoted in kJ/mol. This table provides the quantitative details of the energy partitioning of a single geometry of the water monomer. These values are provided to give an indication of the magnitude of the differences in ECE between individual components of the total ECE.

We can see clearly from Table 1 that the oxygen’s intraatomic correlation energy is the dominant term. In addition, it is clearly noticeable from these energies that this is a distorted geometry with asymmetric hydrogen intra-atomic energies and interaction energies with oxygen. However, the energyminimized structure, used as a seed for the distorted structures, was symmetric. Note that the total molecular dynamic ECE of −539.8 kJ/mol is still relatively very small (0.3%) compared to the total energy for water, which is of the order of −200,000 kJ/ mol. Considering all of the geometries, the average absolute error over the data set was 1.6 kJ/mol with a maximum absolute error of 1.7 kJ/mol (Figure S8). 3.1.3. Kriging Results. A training set of 40 geometries was created for the construction of Kriging models, which were tested on 60 geometries. Here we discuss the ML of seven energy terms: the O, H1, and H2 intra-atomic energies, the O− H1, O−H2, and H1−H2 interatomic energies, and finally the total ECE. Table 2 summarizes the statistics of the predictions of the various ECE contributions with respect to the MORFI reconstruction energies. Table 2 shows that the models are of excellent quality with all q2 values of better than 0.99 and all RMSE values less than 0.5 kJ/mol. We note that the largest maximum absolute error originates from the oxygen intra-atomic energy and the same error is found in the total ECE, but even this error represents only 1.5% of the range of the predictions. From this point onward we focus our attention on the total ECE predictions. Below we display (Figure 2A) a so-called S-curve of the total ECE, which is a plot showing the percentage (y-axis) of total ECE predictions below a certain absolute error (x-axis). Table 2 and Figure 2A display excellent predictive accuracy, with all total ECE predictions below an absolute error of 0.5 kJ/mol and 90% of these predictions below an absolute error of 0.01 kJ/mol, over an energy range of 33.9 kJ/mol. This level of accuracy is sufficient for use in a molecular force field. This is particularly promising given the training set required consisted of only 40 geometries and was able to predict 60 varying geometries. Previous results for the IQA multipolar electrostatics and exchange have required much larger training sets (at least 5−10 times) to achieve similar levels of accuracy.64 As a final test, it is useful to compare three MP2 correlation energies: (i) the original one from GAUSSIAN09 (MP2/ uncon-6-31++G(d,p)), (ii) the “reconstructed” one, obtained from MORFI, which introduces a small numerical error via the

3. RESULTS 3.1. Water Monomer. 3.1.1. Introduction. The first proofof-concept example we present is that of a single water molecule. This case study tests the ability of our ML methods to predict the atomic correlation energy contributions of a molecular system. We refer to these contributions as intramolecular correlation energies, which are rarely studied yet are important to fully comprehend a chemical system’s electronic structure. We note that dispersion energies are traditionally considered inter-molecular in nature and determined by Symmetry Adapted Perturbation Theory, but in our approach they can be intra-molecular as well. However, the current method is not limited to through-space interactions only and can even calculate correlation within an atom bound inside a molecule.83 Thus, the intramolecular correlation energy can be divided up into intra-atomic correlation energy (i.e., the correlation energy within a single atom in the molecule of interest) and interatomic correlation energy (i.e., the correlation energy between two atoms in the molecule of interest). We note that these different atoms do not have to be directly bonded to each other. This full partitioning enables us to study the effects of electron correlation on both chemical bonds and atoms individually. For clarity, we note that the term interatomic may refer to intermolecular or intramolecular interactions. 3.1.2. Energy Partitioning Results. The program TYCHE generated 100 distorted water geometries around one energy minimum. HERMES determined the O−H bonds had a distortion range of 0.2 Å (0.86−1.06 Å) with an average bond length of ∼0.95 Å, while the H−O−H angle had an average of 105.5° and a range of 80.7°. A visualization of the distortion is shown in Section 5 of the SI. A locally modified version of GAUSSIAN09 generated the dcorr,prim matrix for each of these jklm distorted geometries. Subsequently, MORFI partitioned the C


Article


Table 2. Statistics for the Predictions of the Kriging Models for All Seven ECE Components over a Set of Water Monomer Geometries Optimized Using MP2/uncon-6-31++G(d,p)a measure

O intra

H1 intra

H2 intra

O−H1 inter

O−H2 inter

H1−H2 inter

total ECE

max absolute error 99% absolute error 90% absolute error q2 r2 mean absolute error (MAE) root mean squared error (RMSE) range

0.52 0.47 0.10 0.9998 0.9999 0.05 0.12 33.90

0.07 0.04 0.01 0.9999 0.9999 0.00 0.01 9.16

0.12 0.09 0.02 0.9999 0.9999 0.01 0.03 9.03

0.23 0.20 0.03 0.9997 0.9998 0.02 0.05 11.42

0.09 0.07 0.04 0.9999 0.9999 0.01 0.02 11.65

0.03 0.02 0.00 0.9999 0.9999 0.00 0.00 5.93

0.49 0.48 0.09 0.9998 0.9999 0.05 0.12 33.67

a

All energies are quoted in kJ/mol. These statistics display the quality and predictive power of the Kriging models produced for each individual ECE term and the total ECE.

grid used in its energy partitioning, and (iii) the Kriged one, obtained from FEREBUS. For convenience, we take the correlation of the first geometry in the test set as a reference. Figure S9 then shows relative energy from each calculation technique. We can see clearly that the relative energies are consistent between all three methods of calculation. In terms of chemistry, we are generally interested in the relative energy difference between configurations. Hence, it is vital for any calculation method to produce this accurately. Here we clearly demonstrate the utility of the Kriging models to meet this criterion. 3.1.4. RF and SVR Results. The same input features that were employed in the Kriging models were provided to the RF and SVR ML methods. Both RF and SVR were applied to predict the total ECE only, with the results summarized in Table 3.

Table 4. Partitioned and Total ECE Energies (kJ/mol) from MORFI for a Distorted H2···He Complex with a MP2/ uncon-6-31++G(d,p) Wave Functiona intra

RF

SVR

r ±σ RMSE ± σ

0.79 ± 0.01 3.67 ± 0.11

0.87 ± 0.02 2.85 ± 0.28

2

a

complex

H1

H2

He−H1

He−H2

H1−H2

total

−66.3

−45.7

−45.7

0.1

−0.1

22.5

−135.3

a

This table provides the quantitative details of the energy partitioning of a single geometry of the H2···He complex.

of the interactions is stabilizing, while the other is destabilizing in this geometric arrangement. However, taken together, the overall intermolecular interaction is neutral in this case although generally stabilizing. Similar observations have been made previously.30,83 For all the distorted geometries of the H 2 ···He complexes, the maximum error in MORFI’s reconstruction of the correlation energy, for a given geometry, was 1.33 kJ/mol with an average error of 1.29 kJ/mol over a range of approximately 4 kJ/mol. It is shown in Figure S13 that, unlike for the water example, the correlation energy was more consistent over the range of geometries, with over half of these geometries having an ECE in the range −136 to −137 kJ/mol. Additionally, we note that the reconstruction isis similar in accuracy to that of the H2O monomer with a reasonably constant offset of approximately 1 kJ/mol. We have additionally tested the applicability of this method to the atomization of the H2···He complex. These results are presented in Section 7.2 of the SI and show comparable accuracy levels to the current example but over a larger energy range of several hundreds of kJ/mol. 3.2.2. Kriging Results. For the H2···He van der Waals complex a total of 15 geometries were used as training data, and 35 geometries were employed as prediction points. There are once again seven IQA ECE terms in this system. Table 5 summarizes the statistics arising from the Kriging models predictions of the test set with errors measured against the MORFI energy reconstructions. The statistics in Table 5 show the same trends as in Table 2. The quality of the models is excellent as displayed by all q2 values being over 0.98 and MAE and RMSE errors being less than 0.1 kJ/mol. The S-curve in Figure 2B displays the predictive quality of this model. The maximum absolute prediction errors are largest for the total ECE energy, with the single largest max absolute error for a single component being for H1−H2 interatomic ECE. Numerically, these energy values are all notably smaller than that of the water molecule. As with the water monomer example we will now focus on the total ECE. It is clear from Table 5, Figure 2B, and Figure S14, that

Table 3. Total ECE Predictions by RF and SVR ML Methods for the Water Monomera data/measure

inter

He

The RMSE is expressed in kJ/mol.

The results in Table 3 clearly show that the total ECE is a predictable quantity in this case, for both RF and SVR approaches. In both these cases reasonable r2 coefficients are matched by fairly low RMSE values for the predictions. Additionally the standard deviation of the predictions over the 10-fold cross-validations is low, showing that the results are of a reasonable accuracy independent of the training set selected. 3.2. H2···He. The second example is that of the van der Waals complex H2···He, which includes intermolecular as well as intramolecular predictions. 3.2.1. Energy Partitioning Results. The details of the distorted geometries for H2···He are given in Section 5 of the SI. Table 4 gives the breakdown of the ECE for the first distorted geometry as an example. These values are provided to give an indication of the magnitude of the differences in ECE between individual components of the total ECE. Table 4 shows that the helium intra-atomic energy is the largest single contribution, but it is not as singularly dominant as the oxygen intra-atomic correlation energy in the water examples. The distortions produce geometries for this complex ranging from nearly linear to asymmetric T-shape (see Section 5.2 of the SI). We also note that the interaction energies between H and He show different signs. This suggests that one D


Article


Table 5. Summary Statistics for the Predictions of the Kriging Models with Respect to the Seven Generated ECE Terms (kJ/ mol) over a Set of H2···He Geometries at the MP2/uncon-6-31++G(d,p) Level of Theorya

a

measure

H1 intra

H2 intra

He intra

H1−H2 inter

He−H1 inter

He−H2 inter

total ECE

max absolute error 99% absolute error 90% absolute error q2 r2 mean absolute error (MAE) root mean squared error (RMSE) range

0.09 0.06 0.01 0.9999 0.9999 0.01 0.00 7.60

0.08 0.06 0.01 0.9999 0.9999 0.01 0.00 7.61

0.00 0.00 0.00 0.9963 0.9963 0.00 0.00 0.02

0.10 0.09 0.06 0.9998 0.9999 0.02 0.00 11.59

0.00 0.00 0.00 0.9905 0.9925 0.00 0.00 0.02

0.00 0.00 0.00 0.9962 0.9973 0.00 0.00 0.04

0.43 0.20 0.12 0.9853 0.9867 0.04 0.00 3.60

These statistics display the quality and predictive power of the Kriging models produced for each individual ECE term and the total ECE.

the accuracy of the prediction is very good for this intermolecular case. All predictions of the total ECE have an absolute error below 0.5 kJ/mol, and 90% of predictions have an absolute error below 0.12 kJ/mol. Once again this example required very few training points, 15 in total, and yet has a wide domain of applicability, predicting 35 test set points from the resultant models. In terms of relative energy, Figure S14 shows the relative energy differences from the three calculation methods explained above. Despite the increased MORFI reconstruction error the relative errors remain very consistent across the three methods although there are small deviations. The fact that the Kriging models are very close to the energy landscape computed from MORFI demonstrates the ability and utility of these models as a basis for the dispersion energy component of a molecular force field. 3.2.3. RF and SVR Results. As with the water monomer example, we provided the same features to the RF and SVR ML models, for predictions of the total ECE. The results are summarized in Table 6.

Figure 1. Water dimer atom labeling.

water dimer described at the MP2/uncon-6-31++G(d,p) level. As we observed for the water monomer, the oxygen’s intraatomic correlation energy is the largest contribution to the total ECE. The second largest contributions come from hydrogen’s intra-atomic ECE and the O−H interatomic terms. This is followed by the intramolecular H···H interactions and the intermolecular O···O interaction. Notice that there is no symmetry in either the inter- or intra-atomic energies due to the distorted nature of the geometry generated by TYCHE, although the minimum of the dimer has Cs symmetry. Considering all of the geometries, the maximum absolute error in MORFI’s reconstruction of the total ECE is 4.0 kJ/mol with an average absolute error of 3.6 kJ/mol. Figure S23 presents the MORFI reconstruction of the total MP2 energy, along with the MP2/uncon-6-31++G(d,p) energy. 3.3.3. Kriging Results. A training set of 70 geometries was used for the construction of the Kriging models, which were tested for the prediction of 30 geometries. The number of features used for each model was 3N−6 = 3 × 6 − 6 = 12, which is notably larger than for the water monomer and H2··· He complex (3 features each). We present the results for all 21 IQA ECE components in Table S4. Note that the O1···O4 (O···O intermolecular) and the O1···H6 (O···H intermolecular) terms are directly involved in the hydrogen bond.84 Table S4 shows that the quality of the models is similar to those of the water monomer and the H2···He complex with the q2 values still being very high as most values are larger than 0.98, with the only exception being H6···O4 interaction. For this interaction in particular the maximum absolute error is 0.5 kJ/mol. Note that all RMSE values are below 0.2 kJ/mol. The models for the oxygen atoms and the total ECE are still very strong, with q2 values being larger than 0.99. We now focus the discussion on the total ECE, as for the water monomer and H2···He before. Figure 2C confirms the accuracy seen in Table S4 for the prediction of the total ECE. All predictions have a maximum absolute error less than 1.0 kJ/mol, which is below

Table 6. Total ECE Predictions by the RF and SVR ML Methods for the H2···He van der Waals Complexa

a

data/measure

RF

SVR

r2 ± σ RMSE ± σ

0.78 ± 0.02 0.32 ± 0.01

0.7 ± 0.1 0.39 ± 0.06


Again the total ECE is a predictable quantity by the RF and SVR methods. Both achieve reasonable r2 values and low RMSEs. We note the standard deviations over the 10-fold cross-validation are once again low, showing a consistency over different choices of training and test set. 3.3. H2O···H2O. 3.3.1. Introduction. The third case study is the water dimer. We employed TYCHE to generate 100 distorted water dimer geometries from a reference optimized geometry. To provide a sense of the magnitude of the displacements generated by TYCHE, the following statistics were determined by HERMES with reference to the atom labeling in Figure 1. The O−H bonds had an average bond length of 0.97 Å with distortion ranges for the bonds H2−O1, H3−O1, H5−O4, and H6−O4 being 0.21, 0.27, 0.24, and 0.18 Å, respectively. The two water molecules (i.e., H2−O1−H3 and H5−O4−H6), respectively, had average bond angles of 105.48°, over a range of 32.72°, and 104.61° over a range of 37.14°. 3.3.2. Energy Partitioning Results. Table S3 gives the magnitude of the ECE for an arbitrary example of a distorted E


Article


Figure 2. S-curves for the (A) water monomer, (B) H2···He complex, (C) water dimer, and a (D) comparison plot. The S-curves show the percentage of Kriging predictions with a given maximum error compared to the MORFI reconstruction.

ology has recently been applied to larger water clusters.84 It may be possible, through refinements to the numerical integration grid, to further improve the results by reducing the MORFI reconstruction error. This will require further investigations. 3.3.4. RF and SVR Results. The same input features, as for the Kriging ML, were provided for the RF and SVR ML methods, but only the total ECE was predicted. The results are summarized in Table 7.

4.0 kJ/mol (commonly taken as chemical accuracy), and 90% of the predictions have an absolute error below 0.25 kJ/mol. This level of accuracy is appropriate for application to molecular force field development. A training set of 70 geometries predicted 30 geometries with an excellent level of accuracy. Even though the results for the water monomer were slightly better, the number of features necessary to describe the water dimer is four times larger (12 versus 3), demonstrating the increased conformational freedom of the water dimer complex. It is therefore unsurprising that the accuracy is slightly reduced. In terms of the relative energies of the predictions, Figure S24 shows the relative energy for the test set geometries calculated in the three possible ways (i.e., original MP2 energy, MORFI-reconstructed energy and FEREBUS predicted energy). It is clear that the MORFI and MP2/uncon-6-31+ +G(d,p) methods track each other very closely in terms of relative energies. It is also clear that the Kriging predictions track closely to the other method’s correlation energies and, hence, reproduce the relative energies with excellent accuracy. This test clearly shows the use of such models in a molecular force field that then enables an accurate representation of the relative energetic differences between conformations, as demonstrated for all three systems. We have additionally calculated the dissociation energy of the water dimer at the same level of theory that our Kriging model was developed at, which results in a dissociation energy of 26.7 kJ/mol. This dissociation energy was obtained as the difference between the water dimer and two isolated water molecules, at the MP2/ uncon-6-31++G(d,p) level. This energy partitioning method-

Table 7. Prediction of the Total ECE by RF and SVR for the Water Dimera

a

measure

RF

SVR

r2 ± σ RMSE ± σ

0.87 ± 0.01 4.44 ± 0.1

0.91 ± 0.01 3.89 ± 0.26


The results in Table 7 show that the total ECEs predictability by RF and SVR has slightly diminished compared to the water monomer. For both RF and SVR the r2 coefficient and RMSE are good, with RF RMSE slightly over the gold standard 4 kJ/ mol, which is considered chemical accuracy. Both have seen a small reduction in predictive accuracy compared to the water monomer. We note that the standard deviations of the predictions have remained fair, suggesting the methods are consistent in their predictive accuracy with respect to the training and test sets selected. As we noted previously for the Kriging models, one would expect a reduction in the accuracy on going from the water F


Article


an opportunity to consider ML of correlation energies from other ab initio methods.85,86

monomer to the water dimer, as the system size, and thus chemical complexity, increases. This is also seen for the RF and SVR methods. The raw results, in terms of RMSE and r2, show that Kriging has maintained a higher level of predictive power. We note that previous work has enabled the parameters of the Kriging method to be highly optimized to the current task and thus a direct comparison of Kriging, RF, and SVR should not be made based on these data.

■

CONCLUSION The total MP2/uncon-6-31++G(d,p) dynamic electron correlation energy is a readily predictable quantity for small molecular systems by the Kriging machine learning method, which provides excellent results. Kriging can predict, for all the test geometries, their total ECE within 0.5 kJ/mol of the actual MORFI reconstruction energy. Additionally, the training set sizes required to achieve this level of accuracy are considerably smaller (order of tens), compared to those required for the Coulombic and exchange energy counterparts (hundreds or even thousands). These results also represent the last proof-ofconcept step for obtaining the energy components of the FFLUX force field. In its current state, FFLUX already contains well-defined Coulombic and exchange energy terms. It is also of note that the correlation energy is calculated from the same wave function as the other Coulombic and exchange energy terms. From this viewpoint, the method presents a highly coherent way to incorporate ECE, requiring no additional functions or corrections, which also means that we do not regard any interactions to be “special cases”. In a wider context, these results suggest a novel method to include ECE into molecular simulations. In the future, this work may also be applicable to generate energy corrections for other methods, such as DFT.

■

DISCUSSION These results provide a promising start to the exploration of Kriging predicted ECE. Previous work has shown that the task of accurately predicting electrostatic and exchange energies requires larger training data sets.64 The correlation energy partitioning is significantly more costly in computational time than the Coulombic and exchange energy analogues, but a much-reduced training set size can compensate for this additional cost. In other words, a smaller number of geometries are needed in order to generate a widely applicable model. Furthermore, the Kriging models are of excellent quality (q2 ∼ 0.99), and as the system size increases, it appears to retain this behavior for the total ECEs investigated here. However, further work is required to test that this continues to be the case for larger systems and to improve the reconstruction errors from MORFI. The fact that fewer training points are required suggests that ECE has less of a geometric dependence when contrasted with the other energy components. Hence, there is an advantage in knowing the partitioned ECE separately as the alternative approach training for a total atomic energy (i.e., including ECE, Coulomb, exchange and intra-atomic energy) would waste computational resources, as the most expensive energy contribution (ECE) requires the smallest training set. A second advantage of separate knowledge83 of ECE terms is access to the chemical energetics (i.e., dispersion) that QCT and hence FFLUX offers. The systems presented here, although small in size, represent important case studies. Thus, our study of the water monomer and dimer is the first step toward modeling the ECE effects for condensed phases. In addition, our consideration of the van der Waals complex H2···He is important, as it enables us to focus the method on predicting van der Waals dispersion energy. The dispersion energy is a key component of the ECE with important physical consequences, such as inert gas condensation, hydrocarbon liquid, and solid structure as well as some aspects of protein conformation. The predictive quality of the Kriging models is clearly shown by their S-curves (Figure 2). These plots display the ability of all of the Kriging models presented in this work to predict the total ECE correctly to within a few kJ/mol of the exact value over a range of test sets. All of the predictions are well within 4 kJ/mol (commonly asserted as chemical accuracy) of the MORFI reconstruction energies. For the RF and SVR models, the predictions made are of a good quality for all systems. The majority of the models achieves a prediction with an RMSE of approximately 4 kJ/mol or less, which is regularly cited as chemical accuracy. This is another demonstration of the opportunity in using ML methods from standard packages to learn quantum chemical energies and predict them. The Kriging models have maintained an exceptionally high accuracy level, which enables us to suggest a new path toward correlation energy corrections for simulation methods. We also note that other recent work has seen IQA partitioning applied to CCSD energies, providing

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b01157. Figures S1−S49 and Tables S1−S20 including two examples of the topological atom, explanation of the correlated part of the two particle density matrix, Python scripts used to set up training set files for ML, data on the distortion of the systems considered, ML results for RF and SVR approaches, and data on all three systems studied with different basis sets and geometries (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +44 161 3064511. E-mail: [email protected]. ORCID

James L. McDonagh: 0000-0002-2323-6898 Paul L. A. Popelier: 0000-0001-9053-1363 Funding

We acknowledge the EPSRC for funding through the award of an Established Career Fellowship to P.L.A.P. (grant EP/ K005472). A.F.S. thanks the Brazilian government’s FAPESP for the award of postdoctoral position, numbers 2014/21241-9 and 2015/22247-3. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. S. J. Davie and Dr. N. di Pasquale for their support in generating training data files for the FEREBUS Kriging engine, Dr. S. Cardamone for assistance with TYCHE, and Dr. P. I. Maxwell for assistance with the workflow process for potential generation. G


Article


■

(39) Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-R. Sci. Adv. 2017, 3, e1603015. (40) McDonagh, J. L.; Nath, N.; De Ferrari, L.; van Mourik, T. J. Chem. Inf. Model. 2014, 54, 844. (41) Handley, C. M.; Hawe, G. I.; Kell, D. B.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2009, 11, 6365. (42) Bartok, A. P.; Payne, M. C.; Kondor, R.; Csanyi, G. Phys. Rev. Lett. 2010, 104, 136403. (43) Rupp, M.; Tkatchenko, A.; Mueller, K.-B.; von Lilienfeld, O. A. Phys. Rev. Lett. 2012, 108, 058301. (44) Kandathil, S. M.; Fletcher, T. L.; Yuan, Y.; Knowles, J.; Popelier, P. L. A. J. Comput. Chem. 2013, 34, 1850. (45) Popelier, P. L. A. Int. J. Quantum Chem. 2015, 115, 1005. (46) Cressie, N. Statistics for Spatial Data; Wiley: New York, USA, 1993. (47) Rasmussen, C. E.; Williams, C. K. I. Gaussian Processes for Machine Learning; The MIT Press: Cambridge, USA, 2006. (48) Mills, M. J. L.; Popelier, P. L. A. Comput. Theor. Chem. 2011, 975, 42. (49) Mills, M. J. L.; Popelier, P. L. A. Theor. Chem. Acc. 2012, 131, 1137. (50) Fletcher, T. L.; Davie, S. J.; Popelier, P. L. A. J. Chem. Theory Comput. 2014, 10, 3708. (51) Fletcher, T. L.; Popelier, P. L. A. J. Chem. Theory Comput. 2016, 12, 2742. (52) Fletcher, T. L.; Popelier, P. L. A. Theor. Chem. Acc. 2015, 134, 135. (53) Fletcher, T. L.; Popelier, P. L. A. J. Comput. Chem. 2017, 38, 336. (54) Davie, S. J.; Maxwell, P. I.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2017, 19, 20941. (55) Fletcher, T. L.; Popelier, P. L. A. Chem. Phys. Lett. 2016, 659, 10. (56) Cardamone, S.; Popelier, P. L. A. J. Comput. Chem. 2015, 36, 2361. (57) Joubert, L.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2002, 4, 4353. (58) Popelier, P. L. A.; Stone, A. J. Mol. Phys. 1994, 82, 411. (59) Cardamone, S.; Hughes, T. J.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2014, 16, 10367. (60) Fletcher, T. L.; Kandathil, S. M.; Popelier, P. L. A. Theor. Chem. Acc. 2014, 133, 1499. (61) Bader, R. F. W.; Beddall, P. M. J. Chem. Phys. 1972, 56, 3320. (62) Popelier, P. L. A.; Aicken, F. M. J. Am. Chem. Soc. 2003, 125, 1284. (63) Blanco, M. A.; Martín Pendás, A.; Francisco, E. J. Chem. Theory Comput. 2005, 1, 1096. (64) Maxwell, P.; di Pasquale, N.; Cardamone, S.; Popelier, P. L. A. Theor. Chem. Acc. 2016, 135, 195. (65) Maxwell, P.; Martín Pendás, A.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2016, 18, 20986. (66) Hughes, T. J.; Cardamone, S.; Popelier, P. L. A. J. Comput. Chem. 2015, 36, 1844. (67) McDonagh, J.; Davie, S.; Coates, R. HERMES; 2017. (68) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. GAUSSIAN09, Revision B.01; Gaussian, Inc.: Wallingford, CT, USA, 2009.

REFERENCES

(1) Holroyd, L. F.; van Mourik, T. Chem. Phys. Lett. 2007, 442, 42. (2) Thirman, J.; Head-Gordon, M. J. Phys. Chem. Lett. 2014, 5, 1380. (3) Holroyd, L. F.; van Mourik, T. Theor. Chem. Acc. 2014, 133, 1431. (4) Price, S. L.; Leslie, M.; Welch, G. W. A.; Habgood, M.; Price, L. S.; Karamertzanis, P. G.; Day, G. M. Phys. Chem. Chem. Phys. 2010, 12, 8478. (5) Freindorf, M.; Shao, Y.; Furlani, T. R.; Kong, J. J. Comput. Chem. 2005, 26, 1270. (6) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (7) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (8) Johnson, E. R.; Becke, A. D. J. Chem. Phys. 2006, 124, 174104. (9) Tkatchenko, A.; Scheffler, M. Phys. Rev. Lett. 2009, 102, 073005. (10) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (11) Jurečka, P.; Č erný, J.; Hobza, P.; Salahub, D. R. J. Comput. Chem. 2007, 28, 555. (12) Behler, J. J. Chem. Phys. 2016, 145, 170901. (13) Grimme, S. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 211. (14) Tkatchenko, A.; DiStasio, R. A., Jr; Head-Gordon, M.; Scheffler, M. J. Chem. Phys. 2009, 131, 094106. (15) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. Phys. Rev. Lett. 2004, 93, 153004. (16) Reilly, A. M.; Tkatchenko, A. Chemical Science 2015, 6, 3289. (17) Dobson, J. F.; White, A.; Rubio, A. Phys. Rev. Lett. 2006, 96, 073201. (18) Müller, F.; Jacobs, K. Adv.Colloid Interface Sci. 2012, 179, 107. (19) Tkatchenko, A.; DiStasio, R. A., Jr; Car, R.; Scheffler, M. Phys. Rev. Lett. 2012, 108, 236402. (20) Marom, N.; DiStasio, R. A. J.; Atalla, V.; Levchenko, S.; Reilly, A. M.; Chelikowsky, J. R.; Leiserowitz, L.; Tkatchenko, A. Angew. Chem., Int. Ed. 2013, 52, 6629. (21) Reilly, A. M.; Tkatchenko, A. J. Phys. Chem. Lett. 2013, 4, 1028. (22) Stradi, D.; Barja, S.; Díaz, C.; Garnica, M.; Borca, B.; Hinarejos, J. J.; Sánchez-Portal, D.; Alcamí, M.; Arnau, A.; Vázquez de Parga, A. L.; Miranda, R.; Martín, F. Phys. Rev. Lett. 2011, 106, 186102. (23) Zhang, H.-M.; Chen, S.-L. J. Chem. Theory Comput. 2015, 11, 2525. (24) Ambrosetti, A.; Ferri, N.; DiStasio, R. A., Jr.; Tkatchenko, A. Science 2016, 351, 1171. (25) Tkatchenko, A. Adv. Funct. Mater. 2015, 25, 2054. (26) DiStasio, R. A., Jr.; von Lilienfeld, O. A.; Tkatchenko, A. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 14791. (27) Pyzer-Knapp, E. O.; Thompson, H. P. G.; Day, G. M. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 477. (28) Autumn, K.; Liang, Y. A.; Hsieh, S. T.; Zesch, W.; Chan, W. P.; Kenny, T. W.; Fearing, R.; Full, R. J. Nature 2000, 405, 681. (29) Popelier, P. L. A. Phys. Scr. 2016, 91, 033007. (30) McDonagh, J. L.; Vincent, M. A.; Popelier, P. L. A. Chem. Phys. Lett. 2016, 662, 228. (31) Bader, R. F. W. Atoms in Molecules. A Quantum Theory; Oxford Univ. Press: Oxford, Great Britain, 1990. (32) Popelier, P. L. A. In The Nature of the Chemical Bond Revisited; Frenking, G., Shaik, S., Eds.; Wiley-VCH: 2014; Chapter 8, p 271. (33) Rafat, M.; Popelier, P. L. A. J. Comput. Chem. 2007, 28, 2602. (34) Pyzer-Knapp, E. O.; Simm, G. N.; Aspuru Guzik, A. Mater. Horiz. 2016, 3, 226. (35) Salahinejad, M.; Le, T. C.; Winkler, D. A. Mol. Pharmaceutics 2013, 10, 2757. (36) Dral, P. O.; von Lilienfeld, O. A.; Thiel, W. J. Chem. Theory Comput. 2015, 11, 2120. (37) Ramakrishnan, R.; Dral, P. O.; Rupp, M.; von Lilienfeld, O. A. J. Chem. Theory Comput. 2015, 11, 2087. (38) Hansen, K.; Montavon, G.; Biegler, F.; Fazli, S.; Rupp, M.; Scheffler, M.; von Lilienfeld, O. A.; Tkatchenko, T.; Mueller, K.-R. J. Chem. Theory Comput. 2013, 9, 3404. H


Article

Journal of Chemical Theory and Computation (69) Popelier, P. L. A. Comput. Phys. Commun. 1996, 93, 212. (70) Popelier, P. L. A. In The Chemical Bond - 100 years old and getting stronger; Mingos, M., Ed.; Springer: Switzerland, 2016; p 71. (71) Popelier, P. L. A. In Challenges and Advances in Computational Chemistry and Physics dedicated to ″Applications of Topological Methods in Molecular Chemistry″; Chauvin, R., Lepetit, C., Alikhani, E., Silvi, B., Eds.; Springer: Switzerland, 2016; p 23. (72) Vincent, M. A.; Silva, A. F.; McDonagh, J. L.; Popelier, P. L. A. Int. J. Quantum Chem. 2018, e25519. (73) Breiman, L. Mach. Learning 2001, 45, 5. (74) Cortes, C.; Vapnik, V. Machine Learning 1995, 20, 273. (75) Drucker, H.; Burges, C. J. C.; Kaufman, L.; Smola, A. J.; Vapnik, V. N. Advances in Neural Information Processing Systems 1997, 9, 155. (76) Di Pasquale, N.; Davie, S. J.; Popelier, P. L. A. J. Chem. Theory Comput. 2016, 12, 1499. (77) Kuhn, M. J.Stat.Software 2008, 28, 1. (78) McDonagh, J. L.; van Mourik, T.; Mitchell, J. B. O. Mol. Inf. 2015, 34, 715. (79) McDonagh, J. L.; van Mourik, T.; Mitchell, J. B. O. Mol. Inf. 2016, 35, 538. (80) McDonagh, J.; Palmer, D. S.; van Mourik, T.; Mitchell, J. B. O. J. Chem. Inf. Model. 2016, 56, 2162. (81) Hughes, L. D.; Palmer, D. S.; Nigsch, F.; Mitchell, J. B. J. Chem. Inf. Model. 2008, 48, 220. (82) Nath, N.; de Ferrari, L.; McDonagh, J. L. Github, 2016. https:// github.com/Jammyzx1/R-ML-RF-SVM-PLS (accessed 9/11/16). (83) McDonagh, J. L.; Silva, A. F.; Vincent, M. A.; Popelier, P. L. A. J. Phys. Chem. Lett. 2017, 8, 1937. (84) Silva, A. F.; Vincent, M. A.; McDonagh, J. L.; Popelier, P. L. A. ChemPhysChem 2017, DOI: 10.1002/cphc.201700890. (85) Chávez-Calvillo, R.; García-Revilla, M.; Francisco, E.; Martín Pendás, A.; Rocha-Rinza, T. Comput. Theor. Chem. 2015, 1053, 90. (86) Holguin-Gallego, F. J.; Chavez-Calvillo, R.; Garcia-Revilla, M.; Francisco, E.; Martin Pendas, A.; Rocha-Rinza, T. J. Comput. Chem. 2016, 37, 1753.

I


Machine Learning of Dynamic Electron Correlation ... - ACS Publications

Recommend Documents