Non-covalent interactions from a Quantum Chemical Topology perspective

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter's role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Using the Relative Energy Gradient Method with Interacting Quantum Atoms to Determine the Reaction Mechanism and Catalytic Effects in the Peptide Hydrolysis in HIV-1 Protease

The reaction mechanism in an active site is of the utmost importance when trying to understand the role that an enzyme plays in biological processes. In a recently published paper [Theor. Chem. Acc. 2017, 136, 86], we formalised the Relative Energy Gradient (REG) method for automating an Interacting Quantum Atoms (IQA) analysis. Here, the REG method is utilised to determine the mechanism of peptide hydrolysis in the aspartic active site of the enzyme HIV-1 Protease. Using the REG method along with the IQA approach we determine the mechanism of peptide hydrolysis without employing any arbitrary parameters and with remarkable ease (albeit at large computational cost: the system contains 133 atoms, which means that there are 17 689 individual IQA terms to be calculated). When REG and IQA work together it is possible to determine a reaction mechanism at atomistic resolution from data directly derived from quantum calculations, without arbitrary parameters. Moreover, the mechanism determined by this novel method gives concrete insight into how the active site residues catalyse peptide hydrolysis.

Fluorine Gauche Effect Explained by Electrostatic Polarization Instead of Hyperconjugation: An Interacting Quantum Atoms (IQA) and Relative Energy Gradient (REG) Study

We present an interacting quantum atoms (IQA) study of the gauche effect by comparing 1,2-difluoroethane, 1,2-dichloroethane, and three conformers of 1,2,3,4,5,6-hexafluorocyclohexane. In the 1,2-difluoroethane, the gauche effect is observed in that the gauche conformation is more stable than the anti, whereas in 1,2-dichloroethane the opposite is true. The analysis performed here is exhaustive and unbiased thanks to using the recently introduced relative energy gradient (REG) method [ Thacker , J. C. R. ; Popelier , P. L. A. Theor. Chem. Acc . 2017 , 136 , 86 ], as implemented in the in-house program ANANKE. We challenge the common explanation that hyperconjugation is responsible for the gauche stability in 1,2-difluoroethane and instead present electrostatics as the cause of gauche stability. Our explanation of the gauche effect is also is seen in other molecules displaying local gauche conformations, such as the recently synthesized "all-cis" hexafluorocyclohexane and its conformers where all the fluorine atoms are in the equatorial positions. Using our extension of the traditional IQA methodology that allows for the partitioning of electrostatic terms into polarization and charge transfer, we propose that the cause of gauche stability is 1,3 C···F electrostatic polarization interactions. In other words, if a number of fluorine atoms are aligned, then the stability due to polarization of nearby carbon atoms is increased.

Machine Learning of Dynamic Electron Correlation Energies from Topological Atoms

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 H₂···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.

Realistic sampling of amino acid geometries for a multipolar polarizable force field

The Quantum Chemical Topological Force Field (QCTFF) uses the machine learning method kriging to map atomic multipole moments to the coordinates of all atoms in the molecular system. It is important that kriging operates on relevant and realistic training sets of molecular geometries. Therefore, we sampled single amino acid geometries directly from protein crystal structures stored in the Protein Databank (PDB). This sampling enhances the conformational realism (in terms of dihedral angles) of the training geometries. However, these geometries can be fraught with inaccurate bond lengths and valence angles due to artefacts of the refinement process of the X-ray diffraction patterns, combined with experimentally invisible hydrogen atoms. This is why we developed a hybrid PDB/nonstationary normal modes (NM) sampling approach called PDB/NM. This method is superior over standard NM sampling, which captures only geometries optimized from the stationary points of single amino acids in the gas phase. Indeed, PDB/NM combines the sampling of relevant dihedral angles with chemically correct local geometries. Geometries sampled using PDB/NM were used to build kriging models for alanine and lysine, and their prediction accuracy was compared to models built from geometries sampled from three other sampling approaches. Bond length variation, as opposed to variation in dihedral angles, puts pressure on prediction accuracy, potentially lowering it. Hence, the larger coverage of dihedral angles of the PDB/NM method does not deteriorate the predictive accuracy of kriging models, compared to the NM sampling around local energetic minima used so far in the development of QCTFF.

Accurate prediction of polarised high order electrostatic interactions for hydrogen bonded complexes using the machine learning method kriging

As intermolecular interactions such as the hydrogen bond are electrostatic in origin, rigorous treatment of this term within force field methodologies should be mandatory. We present a method able of accurately reproducing such interactions for seven van der Waals complexes. It uses atomic multipole moments up to hexadecupole moment mapped to the positions of the nuclear coordinates by the machine learning method kriging. Models were built at three levels of theory: HF/6-31G(**), B3LYP/aug-cc-pVDZ and M06-2X/aug-cc-pVDZ. The quality of the kriging models was measured by their ability to predict the electrostatic interaction energy between atoms in external test examples for which the true energies are known. At all levels of theory, >90% of test cases for small van der Waals complexes were predicted within 1 kJ mol(-1), decreasing to 60-70% of test cases for larger base pair complexes. Models built on moments obtained at B3LYP and M06-2X level generally outperformed those at HF level. For all systems the individual interactions were predicted with a mean unsigned error of less than 1 kJ mol(-1).

Application of multipolar charge models and molecular dynamics simulations to study stark shifts in inhomogeneous electric fields

Atomic multipole moments are used to investigate vibrational frequency shifts of CO and H(2) in uniform and inhomogeneous electric fields using ab initio calculations and Molecular Dynamics (MD) simulations. The importance of using atomic multipole moments that can accurately represent both molecular electrostatics and the vibrational response of the molecule to changes in the local electric field is highlighted. The vibrational response of CO to applied uniform and inhomogeneous electric fields is examined using Density Functional Theory calculations for a range of test fields, and the results are used to assess the performance of different atomic multipole models. In uniform fields, the calculated Stark tuning rates of Deltamu = 0.52 cm(-1)/(MV/cm) (DFT), Deltamu = 0.55 cm(-1)/(MV/cm) (fluctuating three-point charge model), and Deltamu = 0.64 cm(-1)/(MV/cm) (Multipole model up to octupole), compare favorably with the experimentally measured value of 0.67 cm(-1)/(MV/cm). For H(2), which has no permanent dipole moment, CCSD(T) calculations demonstrate the importance of bond-weakening effects in force fields in response to the applied inhomogeneous electric field. Finally, CO in hexagonal ice is considered as a test system to highlight the performance of selected multipolar models in MD simulations. The approach discussed here can be applied to calibrate a range of multipolar charge models for diatomic probes, with applications to interpret Stark spectroscopy measurements in protein active sites.

Prediction of interaction energies of substituted hydrogen-bonded Watson-Crick cytosine:guanine(8X) base pairs

We investigated the variation in the interaction energy between the Watson-Crick hydrogen-bonded DNA base pairs guanine and cytosine (G(8X):C), where guanine is substituted in the C8 position by 37 different functional groups. Base pairs were optimized at the B3LYP/6-311+G(2d,p) level. A base pair complex containing a more strongly electron-withdrawing group remarkably forms a more stable base pair with C. Multivariate linear regression provided a quantitative relationship between the interaction energies and descriptors generated by the quantum chemical topology (QCT) approach. The descriptors were sampled from the monomers only, not the supermolecular base pair complexes. A model with r2 = 0.96 and a root-mean-square (rms) value of 0.6 kJ/mol was obtained for a training set of 28 base pair complexes. The model was tested by an external test set of 9 complexes, yielding r2 = 0.99 and an rms value of 0.2 kJ/mol. The results indicated that the bonds C6=O6 and N2-H2 at the hydrogen-bonded frontier of the guanine derivatives play an important role in transmitting the substituent effects. A linear correlation between substitution energies and Hammett constants (sigma(m)) was also obtained for all 37 substituents, yielding r2 = 0.82 and an rms value of 1.2 kJ/mol. The model based on QCT descriptors can therefore be used for the prediction of the interaction energy of the base pair G(8x):C, strictly based on data for the G(8x) monomers only.

Interaction energy contributions of H-bonded and stacked structures of the AT and GC DNA base pairs from the combined density functional theory and intermolecular perturbation theory approach

Stacked and Watson-Crick structures of DNA base pairs are investigated with the DFT-SAPT variant of intermolecular perturbation theory, yielding a rigorous decomposition of the interaction energy into electrostatic, induction, dispersion, and exchange contributions. Their interplay in the various structures is analyzed. Total interaction energies extrapolated to the complete basis set limit are compared with corresponding second-order Møller-Plesset and estimated coupled-cluster theory results.

An analysis of the electrostatic interaction between nucleic acid bases

Results from several commonly used approximate methods of evaluating electrostatic interactions have been compared to the rigorous, nonexpanded electrostatic energies at both uncorrelated and correlated levels of theory. We examined a number of energy profiles for both hydrogen bonded and stacked configurations of the nucleic acid base pairs. We found that the penetration effects play an extremely important role and the expanded electrostatic energies are significantly underestimated with respect to the ab initio values. Apart from the inability to reproduce the magnitudes of the ab initio electrostatic energy, there are other problems with the available approximate electrostatic models. For example, the distributed multipole analysis, one of the most advanced methods, is extremely sensitive to the basis set and level of theory used to evaluate the multipole moments. Detailed ab initio results are provided that other researchers could use to test their approximate models.

Long range behavior of high-rank topological multipole moments

The construction of a high-rank multipolar force field (for peptides) is a complex task, leading to several intermediate questions in need of a clear answer. Here we focus on the convergence of the (electrostatic) multipolar expansion at medium and long range. Using molecular electron densities, quantum chemical topology (QCT) defines the atoms as finite volumes, each endowed with multipole moments. The terms in the multipole expansion are grouped according to powers of the internuclear distance, R(-L). Given two atom types at a given distance, we determine which rank (L) is necessary for the electrostatic energy to converge to the exact interaction energy within a certain error. With this information, the rank of the expansion for each interaction can be adapted to the required accuracy and the available computing power.

Transferability of Quantum Topological Atoms in Terms of Electrostatic Interaction Energy

Understanding atomic transferability is important to guide the design of a force field. Atoms in molecules are defined and computed according to the theory of quantum chemical topology (QCT). The electron density associated with such topological atoms is conveniently described by high-rank multipole moments. Here, we assess the transferability of atoms by means of their electrostatic interaction energy, using a convergent multipole expansion. The test systems are (H2O)3 and serine...(H2O)5. The effect of a varying electron density cutoff (i.e., truncating the atoms) is discussed and the effect of polarization is quantified.

The electrostatic potential generated by topological atoms. II. Inverse multipole moments

Quantum chemical topology defines finite atoms, whose bounded electron density generates a well-defined electrostatic potential. A multipole expansion based on spherical tensors provides a potential that is formally convergent outside the divergence sphere. Part I of this series [P. L. A. Popelier and M. Rafat, Chem. Phys. Lett.376, 148 (2003)] showed that a continuous multipole expansion expands the convergence region, thereby allowing the electrostatic potential to be evaluated at short range. Here, we propose a different method, based on "inverse" multipole moments, enabling an expansion that converges everywhere. These moments are defined by inverse (i.e., negative) powers of the magnitude of the position vector describing the electron density inside the atom. We illustrate this technique on nitrogen in N(2), oxygen in H(2)O, and oxygen in the phenolic group of the amino acid tyrosine. The proposed method constitutes a considerable advance over the method presented in Part I.

A TDDFT Study of the Optical Response of DNA Bases, Base Pairs, and Their Tautomers in the Gas Phase

We present calculations of the optical response of the DNA bases and base pairs both in their normal and tautomeric forms in the gas phase, using time-dependent density functional theory (TDDFT). These calculations are performed in real time within the adiabatic approximation with a basis of local orbitals. Our results for the individual bases are in good agreement with experiment and computationally more demanding calculations of chemical accuracy. The optical response of base pairs indicates that the differences between normal and tautomeric forms in certain cases are significant enough to provide a means of identification.

Ab Initio Study of the Complexes of Halogen-Containing Molecules RX (X=Cl, Br, and I) and NH3: Towards Understanding the Nature of Halogen Bonding and the Electron-Accepting Propensities of Covalently Bonded Halogen Atoms

Ab initio calculations have been performed on a series of complexes formed between halogen-containing molecules and ammonia to gain a deeper insight into the nature of halogen bonding. It appears that the dihalogen molecules form the strongest halogen-bonded complexes with ammonia, followed by HOX; the charge-transfer-type contribution has been demonstrated to dominate the halogen bonding in these complexes. For the complexes involving carbon-bound halogen molecules, our calculations clearly indicate that electrostatic interactions are mainly responsible for their binding energies. Whereas the halogen-bond strength is significantly enhanced by progressive fluorine substitution, the substitution of a hydrogen atom by a methyl group in the CH(3)X...NH(3) complex weakened the halogen bonding. Moreover, remote substituent effects have also been noted in the complexes of halobenzenes with different para substituents. The influence of the hybridization state of the carbon atom bonded to the halogen atom has also been examined and the results reveal that halogen-bond strengths decrease in the order HC triple bond CX > H(2)C=CHX approximately O=CHX approximately C(6)H(5)X > CH(3)X. In addition, several excellent linear correlations have been established between the interaction energies and both the amount of charge transfer and the electrostatic potentials corresponding to an electron density of 0.002 au along the R-X axis; these correlations provide good models with which to evaluate the electron-accepting abilities of the covalently bonded halogen atoms. Finally, some positively charged halogen-bonded systems have been investigated and the effect of the charge has been discussed.

Hydrogen Bonding and Covalent Effects in Binding of Cisplatin to Purine Bases: Ab Initio and Atoms in Molecules Studies

Ab initio and density functional calculations are employed to investigate the role of hydrogen bonding in the binding of cisplatin to the purine bases guanine and adenine. Through the use of the theory of atoms in molecules (AIM), it is shown that hydrogen bonds are ubiquitous in such systems, with N-H...N and N-H...Cl interactions present in addition to the expected N-H...O. This in turn means that the known stability of cisplatin-guanine complexes cannot be ascribed solely to hydrogen bonding and allows decomposition of total binding energy into contributions from covalent and hydrogen bonds. To do so, a new method for predicting hydrogen bond energies from bond critical point properties is proposed, employing partial least-squares analysis to remove the family dependence of simple models. Still more hydrogen bond motifs are found in bifunctional complexes of the general type purine-[Pt(NH(3))(2)](2+)-purine, including purine...purine contacts, though again the energetics of these are insufficient to explain the observed trends in stability. Finally, the effect of platination on the pairing of guanine with cytosine is studied in a similar manner, revealing large redistributions of hydrogen bonding but surprisingly small overall changes in pairing energy.

Accurate interaction energies of hydrogen-bonded nucleic acid base pairs

Hydrogen-bonded nucleic acids base pairs substantially contribute to the structure and stability of nucleic acids. The study presents reference ab initio structures and interaction energies of selected base pairs with binding energies ranging from -5 to -47 kcal/mol. The molecular structures are obtained using the RI-MP2 (resolution of identity MP2) method with extended cc-pVTZ basis set of atomic orbitals. The RI-MP2 method provides results essentially identical with the standard MP2 method. The interaction energies are calculated using the Complete Basis Set (CBS) extrapolation at the RI-MP2 level. For some base pairs, Coupled-Cluster corrections with inclusion of noniterative triple contributions (CCSD(T)) are given. The calculations are compared with selected medium quality methods. The PW91 DFT functional with the 6-31G basis set matches well the RI-MP2/CBS absolute interaction energies and reproduces the relative values of base pairing energies with a maximum relative error of 2.6 kcal/mol when applied with Becke3LYP-optimized geometries. The Becke3LYP DFT functional underestimates the interaction energies by few kcal/mol with relative error of 2.2 kcal/mol. Very good performance of nonpolarizable Cornell et al. force field is confirmed and this indirectly supports the view that H-bonded base pairs are primarily stabilized by electrostatic interactions.

Calculation of electrostatic interaction energies in molecular dimers from atomic multipole moments obtained by different methods of electron density partitioning

Accurate and fast evaluation of electrostatic interactions in molecular systems is still one of the most challenging tasks in the rapidly advancing field of macromolecular chemistry, including molecular recognition, protein modeling and drug design. One of the most convenient and accurate approaches is based on a Buckingham-type approximation that uses the multipole moment expansion of molecular/atomic charge distributions. In the mid-1980s it was shown that the pseudoatom model commonly used in experimental X-ray charge density studies can be easily combined with the Buckingham-type approach for calculation of electrostatic interactions, plus atom-atom potentials for evaluation of the total interaction energies in molecular systems. While many such studies have been reported, little attention has been paid to the accuracy of evaluation of the purely electrostatic interactions as errors may be absorbed in the semiempirical atom-atom potentials that have to be used to account for exchange repulsion and dispersion forces. This study is aimed at the evaluation of the accuracy of the calculation of electrostatic interaction energies with the Buckingham approach. To eliminate experimental uncertainties, the atomic moments are based on theoretical single-molecule electron densities calculated at various levels of theory. The electrostatic interaction energies for a total of 11 dimers of alpha-glycine, N-acetylglycine and L-(+)-lactic acid structures calculated according to Buckingham with pseudoatom, stockholder and atoms-in-molecules moments are compared with those evaluated with the Morokuma-Ziegler energy decomposition scheme. For alpha-glycine a comparison with direct "pixel-by-pixel" integration method, recently developed Gavezzotti, is also made. It is found that the theoretical pseudoatom moments combined with the Buckingham model do predict the correct relative electrostatic interactions energies, although the absolute interaction energies are underestimated in some cases. The good agreement between electrostatic interaction energies computed with Morokuma-Ziegler partitioning, Gavezzotti's method, and the Buckingham approach with atoms-in-molecules moments demonstrates that reliable and accurate evaluation of electrostatic interactions in molecular systems of considerable complexity is now feasible.

Atomic properties of selected biomolecules: quantum topological atom types of carbon occurring in natural amino acids and derived molecules

We seek to recover rigorous atom types from amino acid wave functions. The atom types emerge from a cluster analysis operating on a set of seven atomic properties, including kinetic energy, volume, population, and dipole, quadrupole, octupole, and hexadecapole moments. These properties are acquired by partitioning the molecular electron density into quantum topological atoms. Wave functions are generated at the B3LYP/6-311+G(2d,p)//HF/6-31G(d) level for a sensible conformation of each of the 20 naturally occurring amino acids and smaller derived molecules, which together constitute a data set of 57 molecules. From this set 213 unique quantum topological carbons are obtained, which are linked according to the similarity of their properties. After introducing a statistical separation criterion, our cluster analysis proposes two representations: a cruder one with 5 atom types and a finer one with 21 atom types. The immediate coordination of the central carbon plays a major role in labeling the atom types.