Exponential Relationships Capturing Atomistic ... - ACS Publications

Subscriber access provided by Fudan University

Article

Exponential Relationships capturing atomistic Short-Range Repulsion from the Interacting Quantum Atoms (IQA) Method Alex L. Wilson, and Paul L. A. Popelier J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b10295 • Publication Date (Web): 15 Nov 2016 Downloaded from http://pubs.acs.org on November 18, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Exponential Relationships Capturing Atomistic Short-Range Repulsion from the Interacting Quantum Atoms (IQA) Method

Alex L. Wilson and Paul L.A. Popelier*

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

Abstract A topological atom is a quantum object with a well-defined intra-atomic energy, which includes kinetic energy, Coulomb and exchange energy. In the context of intermolecular interactions, this intra-atomic energy is calculated from supermolecular wave functions, by using the topological partitioning. This partitioning is parameter-free and invokes only the electron density to obtain the topological atoms. In this work no perturbation theory is used; instead, a single wave function describes the behaviour of all van der Waals complexes studied. As the monomers approach each other, frontier atoms deform, which can be monitored through a change in their shape and volume. Here we show that the corresponding atomic deformation energy is very well described by an exponential function, which matches the well-known Buckingham repulsive potential. Moreover, we recover a combination rule that enables the interatomic repulsion energy between topological atoms A and B to be expressed as a function of the interatomic repulsion energy between A and A on one hand, and between B and B on the other hand. As a result a link is established between quantum topological atomic energies and classical well-known interatomic repulsive potentials.

Paul Popelier: [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction Short-range repulsion between atoms prevents all matter from falling into infinitely attractive potential energy wells. However, of all the intermolecular forces, short-range repulsion is one of the least well studied because it is expected to become rapidly negligible. The idea that this interaction is shortrange is upheld by common simplistic potentials used in many modern force fields, such as the LennardJones (or 6-12) potential. This classical potential gives unrealistic representations of this repulsive energy because its mathematical shape is perpetuated more by computational efficiency rather than accuracy1. Calculations of this repulsive interaction, in this work and that of Smith2, have unequivocally proven that this interaction acts over longer ranges than the Lennard-Jones potential would suggest. Indeed, Smith commented that short-range repulsion plays a major role in determining atomic and molecular spacing in solids, liquids and gases 2. Alongside the electron-electron and nucleus-nucleus Coulombic repulsions, it is a major component of the steric interaction used to describe many chemical effects. The importance of the repulsive interaction is often understated because, at medium range, it is dominated by the attractive van der Waals (vdW) term. However, this dominance is only prominent at low temperature, while repulsion dominates at high temperature3. Some evidence for this understatement of the repulsive interaction stems from the fact that many molecular dynamics simulations of liquid water have been carried out without any repulsive interaction for hydrogen atoms. Still, non-hydrogen elements are described by explicit Lennard-Jones repulsive parameters and their listing is a fixture of popular (biomolecular) force fields. The origin of short-range repulsion is explained by the Pauli exclusion principle. The same four quantum numbers cannot be shared by any two electrons in an atom, thus an orbital must contain a maximum of two electrons with opposite spin. The Heitler-London antisymmetrization of the electronic wave function means that when two atoms are compressed, electrons must move into higher energy orbitals. This compression thus creates an energetic penalty, which is observed as a repulsive force, causing actual deformation of the electron density4. In the current study we will quantify this deformation energetically, using an atomic energy partitioning method explained later. Short-range repulsion naturally follows from this partitioning because it gauges atomic compression through an intra-atomic energy. This quantum mechanical explanation supports the notion that repulsion is only important when atoms are very close. Hence, widely used repulsive potentials tend to be steep, often without theoretical justification. The overall purpose of the current work is to extract a traditional, repulsive potential from atomically partitioned reduced density matrices. We aim to recover both a familiar mathematical shape (and its parameters) as well as a combination rule. The latter provides parameters for interactions between dissimilar atoms (i.e. different elements) from parameters for interactions between similar atoms (i.e. 2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


same elements). The energy decomposition scheme used in this work is one based on the Quantum Theory of Atoms in Molecules (QTAIM)5-7. This theory partitions8-10 any system (single molecule, van der Waals complex, crystalline matter etc.) to create bounded volumes of electron density about each atom, called topological atoms, without invoking parameters or an external reference state. In its original form11, QTAIM showed that the potential energy of a topological atom is trivially related to the kinetic energy of that atom. The latter is easier to calculate than the former, hence the straightforward relationship between them is most welcome. However, this relationship only exists if the molecule (i.e. the partitioned system) is in a stationary point on the potential energy surface, that is, all the forces on the nuclei vanish. Much later12, this condition of stationarity was eliminated by an algorithm that enabled the calculation of the potential energy of an atom (both intra- and interatomic), independently of the atom’s kinetic energy. This idea was taken up by a method called Interacting Quantum Atoms (IQA)13, which is the topological energy partitioning behind the current work. The fact that the topological atom has well defined kinetic energy14 makes the topological energy partitioning an attractive basis to build a force field on, especially one with a strong link15 to reduced density matrices. This force field is now called FFLUX16 and focuses on biomolecular simulation. FFLUX has an architecture very different to that of traditional force fields; FFLUX does not contain any of the traditional energy expressions (e.g. Hooke potential, Coulomb point charges, Lennard-Jones). Instead, FFLUX describes all structure and dynamics, including all chemical effects, as an interplay between four primary energy contributions: intra-atomic (i) self-energy, and inter-atomic (ii) Coulomb, (iii) exchange and (iv) dynamic correlation. The first primary energy contribution contains the kinetic energy of the atom, its intra-atomic Coulomb, exchange energy and DFT-based correlation energy. The intra-atomic self-energy is at the heart of the current study. We will show how an atom’s deformation, under the influence of other approaching atoms, generates a meaningful change in the intra-atomic self-energy. We will directly observe the deformation of atoms in a van der Waals complex and quantify the corresponding intermolecular deformation energy in response to a change in geometry, as previously observed in separate work studying intramolecular cases of repulsion by Dillen17 and Francisco et al.13. In the present work, atomic deformation energies will be successfully fitted to a well-known exponential function, thereby recovering a traditional repulsive potential from IQA. We also determine the range of this repulsion, challenging the traditional viewpoint that this energetic term dies off rapidly as a function of distance. Moreover, a successful combination rule will be selected from a number of possibilities, including two new ones proposed later. In summary, the force field FFLUX proves to be connected with a traditional repulsive potential energy function, despite its completely different architecture.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 31

2. Background and Methodology 2.1 Classical Repulsive Potentials The Lennard-Jones potential, shown in eqn (1), is used in popular force fields such as AMBER, CHARMM and OPLS because its simplicity lends itself to rapid computation of both the dispersion and short-range repulsion between any two chemical species. This interatomic potential can take the following shape for two atoms at an internuclear distance r,

 σ 12  σ 6  U Lennard − Jones = 4ε   −     r   r  

(1)

where ε is the depth of the potential well, and σ the (finite) distance at which the potential vanishes. The main issue with this potential is the lack of theoretical justification for the repulsive term, which cannot be derived to manifestly have a r-12 dependence. Not long after the introduction of the LennardJones potential, Born and Mayer proposed18 that the repulsion between atoms has an approximately exponential dependence on distance. Repulsion in a molecular complex can be interpreted as due to the overlap of the wavefunctions of the monomers, which leads to exponential behaviour. Inspired by this ansatz of interpenetrating closed-shell electron clouds, Buckingham proposed19 an alternative to the Lennard-Jones potential, known as the exp-6 potential, or

U Buckingham = Ae − Br −

C r6

(2)

Equation 2 expresses a more realistic exponential term to model repulsion and is used extensively in molecular dynamics simulations. However, the use of the Lennard-Jones potential persists in simulations due to the computationally convenient fact that r12 is simply the square of r6. But then again, powers other than 12 (for example 14) are more accurate20.

2.2 Interacting Quantum Atoms (IQA) The Interacting Quantum Atoms (IQA) methodology 13 is an energy decomposition scheme based on topological atoms, which have finite volumes and particular shapes determined by their environment. IQA has been successfully applied to a wide range of chemical systems

21-23

and offers an alternative to the

24

older methodology of Natural Energy Decomposition Analysis (NEDA) . IQA describes the total energy of ஺ ஺஻ the system as a sum of atomic self-energies ‫ܧ‬௦௘௟௙ and interaction energies ‫ܧ‬௜௡௧௘௥ , as opposed to

traditional bonded and non-bonded contributions. We note here that the topological partitioning makes no fundamental a priori distinction between bonded and non-bonded interactions: topological atoms coexist anywhere they appear in a system, as well-defined finite volumes, whether deep inside a molecule or at its edge, or whether as an ion in a crystal. Topological atoms form a space-filling agglomerate of bubble4 ACS Paragon Plus Environment

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


like volumes, that is, without gaps or overlap. They visually express how matter falls apart into atoms. Of course, molecules form (within the system they are part of) as subsets of atoms. Thus, molecules also have sharp and well-defined boundaries: they consist of topological atoms with large inter-atomic exchange energies. Conceptually, molecules distinguish themselves in this way, but through a different definition of topological atoms; they are a direct consequence of electron density patterns, regardless of whether this density is intra-molecular or inter-molecular. Equation 3 shows how the system’s energy is fully described by only intra-atomic (mono-) and interatomic (pairwise) energy contributions, n −1

n

n

A

A B< A

A AB ETotal − IQA = ∑ Eself + ∑∑ Vinter

(3)

where n is the number of atoms of the total system. The self and inter-atomic components are further decomposed as follows. Firstly, the inter-atomic interaction energy decomposition intuitively falls apart as shown in eqn (4), AB Vinter = VnnAB + VneAB + VenAB + VeeAB

(4)

஺஻ ஺஻ ܸ௡௡ accounts for nuclear-nuclear interactions, while ܸ௡௘ accounts for the interaction between the ஺஻ nucleus of atom A and the electrons of atom B, and ܸ௘௡ accounts for interaction between the nucleus of

atom B and the electrons of atom A. ܸ௘௘஺஻ accounts for electron-electron interactions and can be written as the sum of Coulomb and exchange-correlation interactions. AB VeeAB = VCoulomb + VxcAB

(5)

We have now separated out the exchange-correlation interaction from the four classical electrostatic interactions (summarized as ܸ௖௟஺஻ ), giving, AB Vinter = VclAB + VxcAB

(6)

where ܸ௫௖஺஻ describes exchange-correlation energy, which can be seen as a measure of covalency numerically dominated by the exchange part of ܸ௫௖஺஻ . This energy can be evaluated between any two atoms, and hence includes through-space interactions. Secondly, the decomposition of the (mono)atomic self-energy in eqn (3) is equally intuitive, A Eself = VneAA + VeeAA + T A

(7)

஺஺ where ܶ ஺ is the atomic kinetic energy, ܸ௡௘ represents intra-atomic electron-nuclear interactions while

ܸ௘௘஺஺ represents intra-atomic electron-electron interactions. This atomic self-energy attains large values, of the order of tens of Hartrees or hundreds of thousands of kJ/mol, for second row elements. These values are not chemically useful because they refer to the process of constructing an atom (inside a system) 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 31

A starting from its bare nucleus and its individual electrons, all at infinite separation. Obviously, Eself has a

usefulness when summed as above to produce a total energy but in the context of this work more meaningful information comes from changes in atomic self-energy. For example, if we think of an individual atom as a nucleus surrounded by an electron density cloud approaching another atom, the two clouds will be compressed with a corresponding change in volume (deformation), as discussed above. Such a deformation of atomic volume carries an associated energy penalty, which is described by changes in the self-energy of each atom. This deformation energy can be thought of as the IQA description of short-range repulsion. As another large and non-chemical value, the self-energy of a free atom serves as a natural reference value. We can easily calculate (eqn 8) the amount to which each atom is deformed by subtracting, from the value of the atom inside the system, the value of the free atom, A A EDeformation = Eself , in

system

A − Eself , free

(8)

The short-range repulsion energy between two atoms will then simply be the sum of the deformation energy of each atom, AB A B EDeformation = EDeformation + EDeformation

(9)

The validity of eqn 9 will emerge from the numerical evidence presented in the Results Section. More precisely, we will show that it makes sense to regard an “interatomic” repulsion as a sum of intra-atomic effects. This is perhaps a subtle point but conceptually important because it reinterprets the nature of repulsion between two atoms. It is not so common to think of intermolecular forces in terms of deforming atomic shapes in the context of traditional force fields, or even next-generation force fields such as SIBFA25 or AMOEBA26. The deformation of atomic shapes is a direct consequence of the nature of the topological partitioning, but this view is not confined to this type of partitioning. Indeed, short-range repulsion and deformation have been synonymous in the literature for decades 2,4,27-29. It has been shown several times30-34 that short-range repulsion is well represented by an exponential function. In this paper we therefore set out to determine whether IQA calculations of the deformation energy between two species are also exponential. To this end, we will first investigate interactions between noble gas atoms as well as interactions between small, identical molecules. In this first Results Section we will test how well a fitted exponential function of the Buckingham type represents the calculated energies. In the second Results Section, we determine whether simple mixing rules can be employed to capture the behaviour of the repulsion between two different species.


Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3 Computational Details All ab initio calculations were performed using GAUSSIAN0935 to obtain the wavefunctions of the atoms and molecules under study, at the B3LYP/aug-cc-pVTZ level of theory. Each reference molecule was first optimized at this level of theory before the wavefunctions were calculated. The optimized reference geometries were used to construct complexes at a number of relevant distances and the wavefunctions were calculated for each complex at each distance. The constructed complexes were not reoptimized as our aim was to control the orientation of approach of each molecule relative to the other. For example, each complex was arranged so that each molecule approached the other in a linear fashion, maximizing the contact between the atoms of interest. At the same time it was assured that the approach between two target atoms under study was uncluttered. To make this clear, Figure 1 shows an example of an approach between two nitrogens, each one appearing in a different ammonia molecule.

Fig.1 Staggered arrangement of two optimized ammonia molecules.

Whenever substituents such as hydrogen atoms are involved, we employ a staggered conformation (with respect to the hydrogen atoms) to allow the two nitrogen atoms to approach each other in the least sterically hindered way. This allows us to calculate the ‘purest’ repulsion between the atoms of interest in a complex, without having to account for angular dependence or substituent effects. The IQA analysis was performed by version 16.01.09 of the program AIMALL36. The choice of the level of theory was restricted because the only KS-DFT methods currently supported by AIMAll are LSDA, B3LYP and M06-2X. The extension of IQA to operate in the presence of B3LYP has been proven to give chemically meaningful results37, and also provides for the first time energy contributions that add up to a total energy that is very close to the original (unpartitioned) energy of the complex (always less than 0.15 kJ/mol). AIMAll’s highest quality ‘superhigh_leb’ grid was used throughout, as this provided the smallest integration errors, thereby giving the most accurate fits to the Buckingham potentials. IQA sees atoms as malleable portions of electron density, each with a finite volume and hence spatial extension. The fact that the atomic volumes are element-dependent influences the window of sensible separations between two approaching species. As a result it is not wise to perform all calculations 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

at a fixed range of separations for all complexes, such as the interval of 2 Å to 4 Å in 0.2 Å steps for example. Indeed, two hydrogen atoms, for example, with their nuclei placed 4 Å apart will have very small volumes whereas two sulphur atomic centres placed at the same distance apart will have much greater volumes. We therefore cannot compare the two systems at the same nuclear separation because the sulphur system will experience, by design, more repulsion than the hydrogen system as the outer electron density of the sulphur atoms will be much closer together in the sulphur system. Indeed, Smith commented2 that ‘it is somewhat surprising to find potentials for large and small atoms being compared at equal values of r, where the intensity of the disruptive forces acting on the outer electron shells at a fixed distance such as 1ao or 2ao must be extremely different as we pass through a sequence such as the rare gases from He to Rn’. To correct for element-dependent atomic volumes, we work with separation intervals that are relative rather than absolute and one-size-fits-all. In particular, we invoke, as a reference distance, the sum of the van der Waals radii of both atoms involved in the repulsive interaction. Note that these radii are those of the free atoms. For these radii we adopted the values proposed38 by Bondi, given their popularity. This sum will be different for each pair of atoms, and is taken as the informed reference value against which the separation range is set up, as a percentage. The separation corresponding to the sum of van der Waals radii is set to 100%. We then both decrease and increase the separation between the two species, away from the value at the initial separation. We decrease the separation between the two species by 30% in 5% increments and increase the separation between the two species, again by 30% and in 5% increments. We found it useful to include an additional two 2.5% increments at the shortest range in order to better map the rapidly increasing repulsion energy. Because the van der Waals radii depend on the atomic volume, this approach allows us to scale the separations based on the electron density and concomitant atomic volume, in an adaptive and responding manner. This treatment allows a direct comparison between different species, in terms of a compression or expansion percentage of the sum of the van der Waals radii. The value of 30% was chosen because at greater separations the wavefunctions were very occasionally not accurate enough far away from the nuclei, such that their topology was not stable and meaningful, causing AIMAll to fail. If the separation was decreased by more than 30% then substituents such as hydrogen start to affect the deformation energy of the atoms of interest. This is described in more detail in Section 3.1.2. As a result, ±30% was found to be an ideal range to model the deformation energy. Finally we mention that the reference point for “infinite” monomer separation was not the given complex at very large separation, described by a single wavefunction. As mentioned above, the occasional wavefunction inaccuracy at very large separation prohibited us from taking this approach. Hence, we took as a reference simply the two isolated monomers, and the respective atomic energies they provide.


Page 8 of 31

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Results 3.1 Fitting Exponential Potentials to IQA-Calculated Repulsion 3.1.1 Short-Range Repulsion between Noble Gas Atoms Following the method described above, the deformation energy between pairs of noble gas atoms was calculated at a range of separations based on compression and expansion relative to the sum of the van der Waals radii. Figure 2 shows an example of how a varying separation between two neon atoms deforms the electron density of each atom, thus changing the magnitude of the repulsive energy.

Separation relative to

Ne --- Ne

He --- Ar

the vdW separation 70%

100%

130%

Fig. 2 Distortion of malleable atomic volumes due to changing separation in noble gas van der Waals (vdW) systems.

The deformation energy was calculated for the range of noble gas dimers listed in table 1, and plotted against separation in Fig. 3. The plots in Fig. 3 appear to show exponential behaviour in the style of the first part of the Buckingham potential functional form (eqn 2). Therefore we fitted an exponential curve using the Buckingham potential (eqn 10) to each experimentally calculated curve, obtaining the parameters (A and B, not to be confused with atoms A and B, of course) which provide the best fit, and evaluated the root mean square (RMS) of the difference as an indicator of the quality of the fit. We also 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

repeated the above procedure, but this time adding an additional term to the Buckingham potential (eqn 11) to see if this would improve the RMS of each fit. The results are summarized in Table 1.

AB EDeformation = Ae − Br

AB EDeformation = Ae− Br +Cr

(10) 2

(11)

Fig.3 IQA deformation energy for 6 noble gas systems at a range of separations, each expressed as a ratio of the internuclear distance to the sum of van der Waals radii (e.g. 1.2 = 120%).

System

RMS [Eqn 10]

RMS [Eqn 11]

He-He

0.41

0.31

Ne-Ne

0.36

0.23

Ar-Ar

2.37

0.84

He-Ne

0.20

0.19

He-Ar

0.63

0.45

Ne-Ar

1.39

0.49

Average

0.89

0.42

Table 1. Energy errors (rms, kJ/mol) for each exponential fit to 6 systematically deformed noble gas systems.


Page 10 of 31

Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For eqn 10, an average RMS of less than 1 kJ/mol across a range greater than 100 kJ/mol proves that IQA-calculated repulsion energy is well described by an exponential function. Figure 4 confirms this statement, both as a superposition of energies in Fig.4a (original versus fitted) and by explicitly showing the discrepancy between the two (in Fig.4b). The results in Table 1 also indicate that adding another term to the exponential of the Buckingham potential results in a fit that is more than twice as good. However, both potentials provide an excellent fit, showing that short-range repulsion between two noble gas atoms calculated using the IQA methodology exhibits exponential behaviour. Therefore the traditional LennardJones (r-12) view of repulsion is not correct for noble gas atoms. Section 3.1.2 explores whether this exponential behaviour is universal for a range of small molecules and compares the quality of the fit from both exponential equations in each case. Also, we will directly compare Lennard-Jones to the IQAcalculated repulsion.

Fig.4a

Original

IQA-calculated

deformation Fig.4b Energy difference between original and

energies versus predicted energies using eqn 10.

fitted energy at each separation.

Separations in Figure 4 are relative to the sum of the (van der Waals) radii, as described previously in Section 2.3. To convert this representation of separation into ‘r’ (Å) for use in equations 10 and 11, one simply multiplies this separation by the sum of the radii for the two interacting atoms to achieve an ‘r’ value. The unit of the B parameter is therefore Å-1. Using Figure 4 as an example, if one wanted to predict the He-Ne deformation energy using equation 10 at a separation of 0.8 (relative to the sum of the radii), then multiplying 0.8 by the sum of the radii of He and Ne (using the Bondi values referenced in Section 2.3) gives a separation (Å), which can then be used directly in equation 10 as ‘r’.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

3.1.2 Short-Range Repulsion between Identical Small Molecules In this section we examine the repulsive interaction between small and identical molecules as they approach each other and form a dimer. Figure 5 shows the compression and expansion for two out of the eleven systems studied. As mentioned previously, the hydrogen atoms are staggered such as to avoid unwanted steric clashes affecting the repulsion between the atoms of interest. This point will be discussed further below. The eleven systems can be categorized into two types: linear X2-X2 and linear XHn-XHn (n = 1 or 2) repulsion. The repulsion between monomers in each dimer is plotted in Fig. 6, for all 11 systems. Separation relative to

H2O --- OH2

H3N --- NH3


100%

130%

Fig. 5 Distortion of malleable atomic volumes due to changing separation in two dimers.

In the same manner as for the noble gas atoms, we fit exponentials to each IQA-calculated repulsive energy curve using eqns 10 and 11. We then determine the quality of each fit from the value of the RMS of the energy difference at each point of separation. The results are summarized in Table 2.


Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 6a IQA deformation energy for 6 small X2-X2 Fig. 6b IQA deformation energy for 5 small XHndimers at a range of separations.

XHn dimers at a range of separations.

Table 2 again reveals that the extra term in the exponential’s argument (i.e. Cr 2 in eqn 11) results in a dramatic increase in the quality of each fit. In terms of the average RMS values, the extra term doubles the quality of each fit but in some cases such as S2-S2, we see more than an eight-fold improvement. An average RMS value of approximately 1 kJ/mol for eqn 11 is excellent. The average values of about double this value, for eqn 10, is still respectable given the energy range of the order of hundreds of kJ/mol for each system. This result conclusively proves that IQA-calculated repulsion (i.e. deformation energy) behaves exponentially. It is also worth noting that sometimes the quality of fit was compromised by negative deformation energies, a concept discussed below.

System

RMS [Eqn 10]

RMS [Eqn 11]

System

RMS [Eqn 10]

RMS [Eqn 11]

H2-H2

0.28

0.08

NH3-NH3

2.09

1.39

N2-N2

0.88

0.70

OH2-OH2

2.93

1.75

O2-O2

2.61

2.38

FH-FH

1.92

0.44

F2-F2

0.54

0.38

SH2-SH2

2.38

1.02

S2-S2

8.66

1.00

ClH-ClH

1.26

0.71

Cl2-Cl2

2.38

1.47

Average

2.56

1.00

Average

2.12

1.06



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2. Energy errors (rms, kJ/mol) for each exponential fit to 6 systematically deformed systems of small dimers (of identical monomers).

The current work seeks to demonstrate that quantum-mechanically computed short-range repulsion between atoms is fundamentally exponential in nature and can therefore be well represented by a fitted Buckingham potential, whereas the Lennard-Jones potential is unrealistic. Figure 7 makes this point clear by directly comparing IQA-calculated repulsion, between two nitrogen atoms in two identical ammonia molecules, to that of a parameterized Lennard-Jones potential. We plotted only the repulsive part of the standard Lennard-Jones potential (eqn 1), having used standard GAFF39 parameter values and the standard Lorentz-Berthelot mixing rules40,41 to generate the parameterized Lennard-Jones repulsive potential. We also use a standard curve-fitting procedure to vary the Lennard-Jones and Buckingham parameter values to obtain the best possible fits to the calculated data for each of the two potentials.

Fig.7 Original and fitted energy profiles of repulsion between two nitrogen atoms each in two approaching ammonia molecules.

Clearly, the shape of neither the parameterized or fitted Lennard-Jones potentials is appropriate for the description of IQA-calculated repulsion, especially when compared to the Buckingham potential, which fits almost exactly onto the calculated data. The Lennard-Jones gradient is too steep at short-range. Therefore at any very small separation, i.e. less than a 60% compression of the van der Waals radii, Figure 14 ACS Paragon Plus Environment

Page 14 of 31

Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7 would show that the repulsion will be over-predicted, but at greater separations the repulsion drops off far too quickly. For example, at 70 % compression, Lennard-Jones is under-predicting the energy by ~200 kJ/mol. In reality, however, short-range repulsion is still significant at distances greater than those suggested by a Lennard-Jones potential. For example, at 85 % compression, the actual separation between nitrogen nuclei is significant, at 2.64 Å. At this distance, the Lennard-Jones predicted repulsion is ~8 kJ/mol, but the actual calculated repulsion is 108 kJ/mol, which is therefore still a substantial interaction. This finding also impacts on hydrogen-bonding, which is a significant contributor to stability in many complex biomolecules, with bonding interactions typically on the order of tens of kJ/mol with separations of less than 2.00 Å 42. It is therefore fair to say that so-called “short-range” repulsion is a much more important interaction than is currently perceived. In summary, atoms belonging to two different and approaching molecules detect each other’s “repulsive presence” earlier than suggested by a Lennard-Jones potential. The full consequence of this finding on drug-ligand docking43 and molecular dynamics simulations will unfold as FFLUX develops to cope with such applications. FFLUX contains many machine learning models, each one storing how the energy (and multipole moments) of a given atom changes, in response to a variation in its atomic environment. The atomic energy refers to both its interatomic contribution, and more relevant for this study, also its intra-atomic energy. One can imagine that an approximate (but computationally much faster) description of repulsive (deformation) energy would be needed in certain simulations. The machine learning models of repulsive (deformation) energy could then be replaced by Buckingham potentials, which fit excellently to IQA-calculated repulsion, that are fully compatible with IQA and hence topological (atomic) multipole moments. This consistency in energy terms is important for reliability. We noticed that in a small number of systems, the deformation energy can be negative at certain separations. Figure 8 shows the most extreme case of IQA-calculated repulsion between oxygen atoms in two staggered water molecules, where the deformation energy reaches -3 kJ/mol. Although initially dismissed as an artefact due to the basis set size, an increase in basis set size did not cause the negative values to disappear. To confirm the technical validity of this result, every calculation for all systems was performed with the maximum size quadrature grid using the AIMAll program, leading to integration errors of maximum 10-5 a.u. while most are around 10-6 a.u. . We therefore must conclude that this is a real effect although it is too small to affect the quality of fit of each exponential. In real terms, a negative repulsion is an attraction, implying a very small energetic benefit occurring at large separation. The effect is so small that it will most likely be negligible when other interactions (such as electrostatics, exchange and possibly even dispersion) are taken into consideration. However, to examine this effect further, Figure 9 provides a decomposition of the total deformation energy between the two oxygen atoms. Combining eqn 7 with eqns 8 and 9 yields analogues of the deformation energy for VneAA , VeeAA and T A , that is, by subtracting the respective “free” values from the “in system” values. 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 8 Example of the occurrence (en-ellipsed) of negative deformation energies.

Fig. 9 Decomposition of the total deformation energy between two oxygen atoms in two staggered water molecules.

Figure 9 shows that the exponential shape of the repulsive curve originates from the exponential shape of the kinetic energy because VneOO and VeeOO largely cancel each other out, in this case. As explained just above, one should keep in mind that we actually plot the difference between VneOO (in system, i.e. water dimer) and VneOO (free, single water), and similarly for VeeOO and T O . Therefore, the fact that T O becomes slightly negative at long range results in the total deformation energy becoming negative. The work of Dillen17 shows a very similar interplay between TA, VneAA and VeeAA versus interatomic distance r for the pure H-H interaction in tetracyclododecane. To examine whether or not this is a coincidence, Figure 16 ACS Paragon Plus Environment

Page 16 of 31

Page 17 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10 plots the decomposition of the total deformation energy between two nitrogen atoms in two staggered ammonia molecules.

Fig. 10 Decomposition of the total deformation energy between two nitrogen atoms in two staggered ammonia molecules.

Comparing Figures 9 and 10 shows that the three terms behave similarly. In this case, as T N does not go below zero, the total deformation energy between the two nitrogen atoms does not reach sub-zero values either. We chose a maximum compression value of 70 % across all systems because at closer range (i.e. < 70 %), the hydrogen atoms begin to interfere with the repulsion between the nitrogen atoms and we start to lose the exponential behaviour. This effect occurs for the more sterically hindered systems, i.e. those with multiple hydrogen substituents such as NH3-NH3. Figure 11 shows the IQA-calculated repulsion between two nitrogen atoms in two different ammonia molecules, but now including separations smaller than 70%, down to 50 %. We fitted an exponential using eqn 9 to all points in the interval 70-130 % (RMS = 2.09 kJ/mol, see Table 2) and then extrapolated this curve to closer range. The fitted-extrapolated curve (Fig. 11, red) shows a clear departure from the original IQA-calculated curve (black). Therefore, given our interest in the pure N-N interaction, we chose a cut-off of 70% and kept this cut-off consistent across all systems.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig.11 Comparison between the IQA-calculated repulsion and the fitted exponential between nitrogen atoms in two different ammonia molecules, which is extrapolated to closer separations (below 70%).

3.1.3 Short-Range Repulsion between Mixed-Monomer Small Molecules In the previous Section, we looked at the repulsion between 11 dimers of identical small monomers. After concluding that the deformation energy between atoms appearing in identical molecules is exponential, we now investigate the same question for complexes of mixed (i.e. different) monomers. The very good exponential fits already obtained for the mixed noble gas systems suggest that for the mixedmonomer complexes studied here, there will also be no departure from the exponential behaviour. Instead of calculating the deformation energy (i.e. repulsion) between all possible mixed systems (of which there are 55), we focus on the interaction between water and every other molecule (of which there are 10), as well as a small random assortment of other interactions (another 8 complexes). Figure 12 exemplifies how some mixed systems are deformed.


Page 18 of 31

Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Separation relative to

H2O --- NH3

H2O --- H2


100%

130%

Fig. 12 Distortion of malleable atomic volumes due to changing separation in mixed-monomer van der Waals complexes.

Table 3 shows that for each mixed system exponential curves with a low RMS relative to the scale of deformation energy (hundreds of kJ/mol) can be fitted. Again, the quality of the fit can be more than doubled if eqn 11 is used instead.

System

RMS [Eqn 10]

RMS [Eqn 11]

OH2-H2

1.90

0.17

OH2-N2

2.30

OH2-NH3

System

RMS [Eqn 10]

RMS [Eqn 11]

NH3-H2

1.73

0.62

0.58

NH3-FH

1.24

0.64

2.48

1.55

N2-H2

0.58

0.41

OH2-O2

2.65

0.54

N2-S2

0.71

0.63

OH2-FH

0.72

0.45

N2-SH2

1.10

0.93

OH2-ClH

2.98

0.64

H2-F2

0.43

0.22

OH2-Cl2

4.23

0.66

H2-FH

1.56

0.26



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

OH2-F2

2.98

0.36

OH2-S2

4.66

0.71

OH2-SH2

4.28

0.74

Average

2.92

0.64

Page 20 of 31

F2-FH

2.55

0.34

Average

1.24

0.51

Table 3. Energy errors (rms, kJ/mol) for each exponential fit to 18 systematically deformed systems of mixed monomers (of which 10 are water-containing complexes, left column). The study of mixed systems allows us to determine whether the relative deformation energies of each system make chemical sense. For example, atoms that are traditionally recognized as being ‘small and hard’ may be more resistant to deformation and thus their repulsive energy curves will lie lower down due to reduced deformation energy between atoms. Conversely, ‘large and soft’ atoms would be more prone to deformation. Figure 13 tells us that these trends are largely obeyed. For example, the largest repulsions are found between oxygen and sulphur in the OH2-SH2 system, and between nitrogen and sulphur in the N2-SH2 system. On the other hand, the smallest repulsion is found between hydrogen and fluorine in the H2-F2 complex.

Fig. 13a IQA deformation energy between the oxygen atom in a water molecule and the heaviest atom in any of the 10 small molecules water intearcts with, at the 70-130% range of separations.


Page 21 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 13b IQA deformation energy between the heaviest atom in each molecule and the heaviest atom in the second molecule in the complex, for 8 small molecule complexes, at the 70-130% range of separations.

3.2 Mixing Rules We now know that for a wide variety of small molecules and noble gas atoms, the form of shortrange repulsion between atoms is modelled well by eqn 10 and even better by eqn 11. We are now in a strong position to find out if one can quickly determine the repulsion between two atoms in different molecules in the knowledge of the deformation energy of only the ‘like’ interactions. For example, to determine the deformation energy of an AB interaction from the AA and BB interactions only, where A and B cover a variety of atoms (elements). If this idea works then one does not have to compute the AB interaction directly, from scratch. It is highly desirable to be able to rapidly and reliably estimate accurate quantum mechanical calculations with relatively simplistic formulae that take a very short amount of time to compute. For this purpose, we will use mixing rules and investigate here their validity. Mixing rules are widely used in force fields but their approximate nature is widely recognized. However, this disadvantage is outweighed by the fact that the rules eliminate the need for a direct calculation of quantum mechanical energies. Moreover, these rules are typically very simple, and hence rapid to evaluate. Due to the popularity of the Lennard-Jones potential, many mixing rules were designed to give values for Lennard-Jones parameters. For example, the most used mixing rules are the LorentzBerthelot rules40,41, despite the fact that they yield significant deviations from experiment44. Other Lennard-Jones potential mixing rules that are considered to be more accurate include the Waldman21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

Hagler45 and Kong46 rules. When considering the Buckingham-style exponential potentials, commonly used rules are the Gilbert-Smith2,47,48 and Böhm-Ahlrichs rules49. Table 4 lists each set of mixing rules that were tested.

Equation 1. Gilbert-Smith

Mixing Rules

E DEF = Aij e

− Bij r 1  B jj Bii  ( Aii Bii ) ( Ajj B jj )  Aij =  Bij 1

Bij = 2 2. Böhm-Ahlrichs

3. Arithmetic Average

E DEF = Aij e

− Bij r

E DEF = Aij e

− Bij r

E DEF = Aij e

Bii + B jj

 Aij =  Aii A jj   BB Bij = 2 ii jj Bii + B jj 1 Bii

Aij = Bij =

4. Generalized Average

Bii B jj

− Bij r

Bij

1 B jj

2   

Aii + Ajj 2 Bii + B jj 2 1

 An + Anjj Aij =  ii  2  n = 0.4

n  

 Biin + B njj Bij =   2  n=2

n  

1

5. Arithmetic Average

E DEF = Aij e

− Bij r + Cij r 2

Aij = Bij = Cij =

Aii + A jj 2 Bii + B jj 2 Cii + C jj 2


Bij

2   

Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Generalized Average

E DEF = Aij e

− Bij r + Cij r 2

1

 Aiin + Anjj Aij =   2  n = 3.2

n  

 Biin + B njj Bij =   2  n = 1.35

n  

1

 ( C + s ) n + ( C + s )n ii jj Cij =   2  n = 1.35 s = 10

1

n  −s  

Table 4. List of mixing rules for each type of exponential equation.

The majority of modern force fields utilize parameterization for rapid calculation of otherwise time consuming calculations. Parameterization is built upon the idea that we calculate AA and BB values, and store them in databases so that AB values may be computed rapidly. So indeed, AA and BB parameter values must be calculated but then they can be stored so that any AB combination is possible to compute. However, to avoid having such an enormous database of AA and BB parameters, the transferable nature of atoms and functional groups in different systems can be exploited to create an atom type. Investigating atomic transferability is not touched upon in the present work but will likely form the basis of a future publication. Table 5 compares all six mixing rules to see which give the most accurate fits to the quantummechanical IQA-calculated deformation energies of the mixed systems. When constructing generalized averages (Rules 4 and 6), values of n were fitted to give the lowest overall value of the average RMS over all systems. Rule 6 introduces a technical matter that needs a little discussion. We note that the values for Cii or Cjj can be positive or negative. In the current work, the negative values are always larger than -10, a number that is associated with the choice of molecules studied here. Negative values of Cii or Cjj pose a potential practical problem because they need to be raised to a power, which mathematically can give complex numbers and multiple-valued answers. This issue can be illustrated by the simple case of n=2 in Rule 6. In that case, the calculation of Cij first squares Cii and Cjj and then takes the square root of their resulting sum, which erases the original sign of Cii and Cjj. This issue can be tackled as follows: add a “shift” (denoted s in Table 4) to both Cii and Cjj such that they are positive and then perform the operation of raising powers and taking roots. In this work this shift is set to 10, for the reason explained above. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 31

The well-known mixing rules of Gilbert and Smith2,47,48, and of Böhm and Ahlrichs49 were tested primarily on simple noble gas interactions. As a result, from the RMS values in Table 5, we see that these rules give excellent fits to the actual deformation energies for the noble gas systems. However, as we move on to larger systems, both rules deteriorate in quality. As Böhm and Ahlrichs remarked, problems arise as the difference between Bii and Bjj increases49, that is, for bigger systems. Using the average RMS of all systems as an indicator of the quality of a rule, we conclude that Rule 4, our own generalized average, is most representative of the actual calculated deformation energies.

System

Rule 1

Rule 2

Rule 3

Rule 4

Rule 5

Rule 6

RMS

RMS

RMS

RMS

RMS

RMS

He-Ne

0.38

0.38

3.87

1.94

0.56

5.15

He-Ar

1.47

1.46

4.12

1.47

2.12

13.31

Ne-Ar

1.42

1.42

1.94

2.42

1.49

0.80

OH2-H2

34.73

34.49

23.20

16.26

31.34

30.07

OH2-N2

17.43

17.25

16.02

18.32

18.47

11.06

OH2-NH3

2.57

2.62

2.73

4.06

2.71

3.60

OH2-O2

5.98

8.79

75.66

16.77

8.06

11.59

OH2-FH

14.07

13.99

28.40

3.87

11.61

1.25

OH2-ClH

3.15

3.61

44.28

15.71

5.28

8.18

OH2-Cl2

15.53

15.49

7.03

8.15

14.24

4.89

OH2-F2

3.13

3.12

1.64

2.24

3.46

43.99

OH2-S2

8.52

8.03

5.77

10.49

10.63

6.82

OH2-SH2

3.75

3.57

6.63

1.00

3.28

8.81

NH3-H2

38.01

37.97

17.16

20.19

36.06

29.72

NH3-FH

27.63

27.63

25.96

26.92

27.66

12.53

N2-H2

35.82

38.90

69.60

31.98

31.48

33.01

N2-S2

1.53

1.58

1.32

1.38

1.34

1.55


Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


N2-SH2

39.61

39.30

28.14

37.11

40.04

39.00

H2-F2

17.55

17.38

18.75

20.21

19.40

19.69

H2-FH

7.51

4.87

46.01

15.70

16.30

10.28

F2-FH

34.25

33.71

34.51

39.22

37.72

21.58

14.95

15.03

22.04

14.07

15.39

15.09

13.81

13.97

21.30

11.80

13.19

12.65

Average RMS Standard Deviation

Table 5. Comparison of the performance of each mixing rule across all systems, using RMS (in kJ/mol) as an indicator of the quality of fit of the predicted compared to the computed deformation energies.

Interestingly, although adding an extra term to the exponential potential function gives curves that fit the actual computed deformation energies better, using this extra term is not advantageous when considering mixing rules. For an extra parameter to be computationally worthwhile, it would need to provide a significant improvement in the overall average RMS value, which is not the case. The best rule is actually one of the simplest, which suits Occam’s razor. Our generalized average is simpler than the Gilbert-Smith and Böhm-Ahlrichs rules but predicts deformation energies more accurately across all systems by about 1 kJ/mol. This finding can be explained by the fact that the older rules were developed and tested on noble gas atomic interactions, which are naturally very simplistic and are not representative of many small molecule interactions. Finally we note that this study has focused on one energy term (Edef) and its link to classical repulsion, which is only a part of the total energy. Of course other energy terms such as Vxc are also present but not studied here. It is known that the atomic interaction lines are governed by Vxc, which has been proven to be related to the presence of a bond critical point. Via the additive nature of the IQA energy contributions, the Vxc contribution (or others) can dominate the repulsive deformation energy. In summary, there is no conflict in the simultaneous presence of attracting (Vxc) and repelling energy contributions (Edef).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4. Conclusion The topological energy-decomposition scheme called IQA is employed here to calculate deformation energies between pairs of atoms in small molecules. These deformation energies equate to and are thus synonymous with intermolecular short-range repulsion. We show that this repulsion is much better represented by exponential functions such as the Buckingham potential compared to the r-12 contribution in the Lennard-Jones potential. In particular, inclusion of an additional term in Buckingham potential improves the quality of every fit. Quantum-mechanical calculations on small molecules are rapid, but calculations on larger systems become increasingly computationally expensive. Instead, mixing rules are commonly used to give rapid approximations of quantum mechanical values based on parameters determined for ‘like’ systems. That is, determining an approximation for an AB interaction from previous calculations of AA and BB interactions. The quality of a mixing rule is thus determined by how well the approximation fits to the actual calculated data. We find that common rules such as Gilbert-Smith and Böhm-Ahlrichs are excellent for the simplistic noble gas systems they were developed on, but are less good when expanding to small molecules. Overall, our fitted generalized average, using the form of the potential without the additional term, performs the best on average across all systems tested. Overall, three important lessons were learnt: (i) IQA-repulsion resembles much better the exponential potential, which is known to be more accurate than the Lennard-Jones one, (ii) topological atoms detect each other’s presence (in terms of repulsion) over longer distances than expected and (iii) the interatomic nature of the repulsion between two atoms can be re-interpreted as the sum of mutual and concerted mono-atomic deformations. As a final point, we acknowledge that values of A, B and C in equations 10 and 11 may have a physical meaning. For example, we measure here the deformation energy when two atoms are compressed against each other and so the A, B and C values are probably related to physical quantities such as atomic hardness. Although we do not touch upon this topic here, it could be the focus of a future study. The raw data for the values of A, B and C can be found in the Supporting Information.

Associated Content Supporting Information List of the values of the parameters A, B and C appearing in eqns 10 and 11. This material is available free of charge via the Internet at http://pubs.acs.org.


Page 26 of 31

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgments Gratitude is also expressed to the EPSRC for the award of an Established Career Fellowship to one of us (P.L.A.P.) and to the BBSRC for a PhD studentship. References (1) Soper, A. K. Joint structure refinement of x-ray and neutron diffraction data on disordered materials: application to liquid water. J. Physics: Cond.Matter 2007, 19, 335206. (2) Smith, F. T. Atomic Distortion and the Combining Rule for Repulsive Potentials. Physical Review A 1972, 5, 1708. (3) Nedostup, V. I. Direct Methods for Derivation of Intermolecular Interaction Energies from Thermal Data. Reviews on Thermal Properties of Materials 1983, 1, 73. (4) Mitoraj, M. P.; Michalak, A.; Ziegler, T. A Combined Charge and Energy Decomposition Scheme for Bond Analysis. J.Chem.Theor.Comput. 2009, 5, 962. (5) Bader, R. F. W. Atoms in Molecules. A Quantum Theory.; Oxford Univ. Press: Oxford, Great Britain, 1990. (6) Bader, R. F. W. A Quantum-Theory of Molecular-Structure and Its Applications. Chem.Rev. 1991, 91, 893. (7) Popelier, P. L. A. Atoms in Molecules. An Introduction.; Pearson Education: London, Great Britain, 2000. (8) Popelier, P. L. A. In The Nature of the Chemical Bond Revisited; Frenking, G., Shaik, S., Eds.; Wiley-VCH, Chapter 8: 2014, p 271. (9) 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. (10) Popelier, P. L. A. In The Chemical Bond - 100 years old and getting stronger; Mingos, M., Ed.; Springer: Switzerland, 2016, p 71. (11) Bader, R. F. W.; Beddall, P. M. Virial Field Relationship for Molecular Charge Distributions and the Spatial Partitioning of Molecular Properties. J.Chem.Phys. 1972, 56, 3320. (12) Popelier, P. L. A.; Kosov, D. S. Atom-atom partitioning of intramolecular and intermolecular Coulomb energy. J.Chem.Phys. 2001, 114, 6539. (13) Blanco, M. A.; Martin Pendas, A.; Francisco, E. Interacting Quantum Atoms: A Correlated Energy Decomposition Scheme Based on the Quantum Theory of Atoms in Molecules. J.Chem.Theor.Comput. 2005, 1, 1096. (14) Fletcher, T. L.; Kandathil, S. M.; Popelier, P. L. A. The prediction of atomic kinetic energies from coordinates of surrounding atoms using kriging machine learning. Theor.Chem.Acc. 2014, 133, 1499:1. 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(15) Popelier, P. L. A. QCTFF: On the Construction of a Novel Protein Force Field. Int.J.Quant.Chem. 2015, 115, 1005. (16) Popelier, P. L. A. Molecular Simulation by Knowledgeable Quantum Atoms. Phys.Scr. 2016, 91, 033007. (17) Dillen , J. Congested Molecules. Where is the Steric Repulsion? An Analysis of the Electron Density by the Method of Interacting Quantum Atoms. Int.J.Quant.Chem. 2013, 113, 2143. (18) Born, M.; Mayer, J. E. Zur Gittertheorie der Ionenkristalle. Z.Phys. 1932, 75, 1. (19)

Buckingham, R. A.

The Classical Equation of State of Gaseous Helium, Neon and Argon.

Proc.Roy.Soc.A 1938, 168, 264. (20) Halgren, T. Merck Molecular Force Field (MMFF). J.Am.Chem.Soc. 1992, 114, 7827. (21) Eskandari, K.; Van Alsenoy, C. Hydrogen–Hydrogen Interaction in Planar Biphenyl: A Theoretical Study Based on the Interacting Quantum Atoms and Hirshfeld Atomic Energy Partitioning Methods. J.Comput.Chem. 2014, 35, 1883. (22) Jarzembska, K. N.; Dominiak, P. M. New version of the theoretical databank of transferable aspherical pseudoatoms, UBDB2011 - towards nucleic acid modelling. Acta Crystallographica Section A 2012, 68, 139. (23) Martin Pendas, A.; Francisco, E.; Blanco, M. A. Binding Energies of First Row Diatomics in the Light of the Interacting Quantum Atoms Approach. J.Phys.Chem. A 2006, 110, 12864. (24) Glendening, E. F.; Streitwieser, A. Natural energy decomposition analysis: An energy partitioning procedure for molecular Interactions with application to weak hydrogen bonding, strong ionic, and moderate donor-acceptor interactions. J.Chem.Phys. 1994, 100, 2900. (25)

Gresh, N.; Cisneros, G. A.; Darden, T. A.; Piquemal, J.-P.

Anisotropic, Polarizable Molecular

Mechanics Studies of Inter- and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy. J.Chem.Theory Comput. 2007, 3, 1960. (26) Ponder, J. W.; Wu, C.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A. J.; Head-Gordon, M.et al. Current Status of the AMOEBA Polarizable Force Field. J.Phys.Chem.B. 2010, 114, 2549. (27) Kong, C. L. Atomic distortion and the repulsive interactions of the noble gas atoms. J.Chem.Phys. 1973, 59, 968. (28) Carlton, T. S.; Gittins, C. M. Internuclear distance and repulsive potential between noble-gas atoms: Combining rules based on repulsive force. The Journal of Chemical Physics 1991, 94, 2614. (29) Zhou, N.; Lu, Z.; Wu, Q.; Zhang, Y. Improved parameterization of interatomic potentials for rare gas dimers with density-based energy decomposition analysis. The Journal of Chemical Physics 2014, 140, 214117. (30) Nyeland, C.; Toennies, J. P. Modelling of repulsive potentials from atom charge density distributions: interactions of inert gas atoms. Chemical Physics Letters 1986, 127, 172. 28 ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) Hoinkis, J.; Ahlrichs, R.; Böhm, H. J. A simple treatment of intermolecular interactions: Synthesis of ab initio calculations and combination rules. International Journal of Quantum Chemistry 1983, 23, 821. (32) Stone, A. J.; Tong, C. S. Anisotropy of atom–atom repulsions. Journal of Computational Chemistry 1994, 15, 1377. (33) Brdarski, S.; Karlström, G. Modeling of the Exchange Repulsion Energy. J.Phys.Chem.A 1998, 102, 8182. (34) Aquilanti, V.; Bartolomei, M.; Cappelletti, D.; Carmona-Novillo, E.; Pirani, F. The N2–N2 system: An experimental potential energy surface and calculated rotovibrational levels of the molecular nitrogen dimer. The Journal of Chemical Physics 2002, 117, 615. (35) 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.et al.; Gaussian, Inc.: Wallingford CT, 2009. (36) AIMAll Todd A. Keith, TK Gristmill Software, Overland Park KS, USA, (aim.tkgristmill.com). 2016. (37) Maxwell, P.; Martin Pendas, A.; Popelier, P. L. A. Extension of the interacting quantum atoms (IQA) approach to B3LYP level density functional theory. PhysChemChemPhys 2016, 18, 20986. (38) Bondi, A. van der Waals Volumes and Radii. The Journal of Physical Chemistry 1964, 68, 441. (39) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. Journal of Computational Chemistry 2004, 25, 1157. (40) Lorentz, H. A. Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Annalen der Physik 1881, 248, 127. (41) Berthelot Sur le Mélange des Gaz. Comptes rendus hebdomadaires des séances de l'Académie des Sciences 1898, 126, 1703. (42) Raymo, F. M.; Bartberger, M. D.; Houk, K. N.; Stoddart, J. F. The Magnitude of [C−H···O] Hydrogen Bonding in Molecular and Supramolecular Assemblies. Journal of the American Chemical Society 2001, 123, 9264. (43) Popelier, P. L. A. New insights in atom-atom interactions for future drug design. Current Topics in Med.Chem. 2012, 12, 1924. (44) Delhommelle, J.; MilliÉ, P. Inadequacy of the Lorentz-Berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation. Molecular Physics 2001, 99, 619. (45)

Waldman, M.; Hagler, A. T.

New combining rules for rare gas van der waals parameters.

J.Comput.Chem. 1993, 14, 1077. (46) Kong, C. L. Combining rules for intermolecular potential parameters. II. Rules for the Lennard-Jones (12–6) potential and the Morse potential. J.Chem.Phys. 1973, 59, 2464. (47) Gilbert, T. L. Soft-Sphere Model for Closed-Shell Atoms and Ions. The Journal of Chemical Physics 1968, 49, 2640. (48) Gilbert, T. L.; Simpson, O. C.; Williamson, M. A. Relation between charge and force parameters of closed-shell atoms and ions. J.Chem.Phys. 1975, 63, 4061. 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(49) Boehm, H.-J.; Ahlrichs, R. A study of short-range repulsions. J.Chem.Phys. 1982, 77, 2028.


Page 30 of 31

Page 31 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The intra-atomic energies of topological atoms behave exponentially upon molecular compression. 88x66mm (96 x 96 DPI)

ACS Paragon Plus Environment

Exponential Relationships Capturing Atomistic ... - ACS Publications

Recommend Documents