Density-functional theory-symmetry-adapted intermolecular perturbation theory with density fitting: A new efficient method to study intermolecular interaction energies

Energy component analysis of π interactions

Fundamental features of biomolecules, such as their structure, solvation, and crystal packing and even the docking of drugs, rely on noncovalent interactions. Theory can help elucidate the nature of these interactions, and energy component analysis reveals the contributions from the various intermolecular forces: electrostatics, London dispersion terms, induction (polarization), and short-range exchange-repulsion. Symmetry-adapted perturbation theory (SAPT) provides one method for this type of analysis. In this Account, we show several examples of how SAPT provides insight into the nature of noncovalent π-interactions. In cation-π interactions, the cation strongly polarizes electrons in π-orbitals, leading to substantially attractive induction terms. This polarization is so important that a cation and a benzene attract each other when placed in the same plane, even though a consideration of the electrostatic interactions alone would suggest otherwise. SAPT analysis can also support an understanding of substituent effects in π-π interactions. Trends in face-to-face sandwich benzene dimers cannot be understood solely in terms of electrostatic effects, especially for multiply substituted dimers, but SAPT analysis demonstrates the importance of London dispersion forces. Moreover, detailed SAPT studies also reveal the critical importance of charge penetration effects in π-stacking interactions. These effects arise in cases with substantial orbital overlap, such as in π-stacking in DNA or in crystal structures of π-conjugated materials. These charge penetration effects lead to attractive electrostatic terms where a simpler analysis based on atom-centered charges, electrostatic potential plots, or even distributed multipole analysis would incorrectly predict repulsive electrostatics. SAPT analysis of sandwich benzene, benzene-pyridine, and pyridine dimers indicates that dipole/induced-dipole terms present in benzene-pyridine but not in benzene dimer are relatively unimportant. In general, a nitrogen heteroatom contracts the electron density, reducing the magnitude of both the London dispersion and the exchange-repulsion terms, but with an overall net increase in attraction. Finally, using recent advances in SAPT algorithms, researchers can now perform SAPT computations on systems with 200 atoms or more. We discuss a recent study of the intercalation complex of proflavine with a trinucleotide duplex of DNA. Here, London dispersion forces are the strongest contributors to binding, as is typical for π-π interactions. However, the electrostatic terms are larger than usual on a fractional basis, which likely results from the positive charge on the intercalator and its location between two electron-rich base pairs. These cation-π interactions also increase the induction term beyond those of typical noncovalent π-interactions.

Dispersion-corrected density functional theory for aromatic interactions in complex systems

Aromatic interactions play a key role in many chemical and biological systems. However, even if very simple models are chosen, the systems of interest are often too large to be handled with standard wave function theory (WFT). Although density functional theory (DFT) can easily treat systems of more than 200 atoms, standard semilocal (hybrid) density functional approximations fail to describe the London dispersion energy, a factor that is essential for accurate predictions of inter- and intramolecular noncovalent interactions. Therefore dispersion-corrected DFT provides a unique tool for the investigation and analysis of a wide range of complex aromatic systems. In this Account, we start with an analysis of the noncovalent interactions in simple model dimers of hexafluorobenzene (HFB) and benzene, with a focus on electrostatic and dispersion interactions. The minima for the parallel-displaced dimers of HFB/HFB and HFB/benzene can only be explained when taking into account all contributions to the interaction energy and not by electrostatics alone. By comparison of saturated and aromatic model complexes, we show that increased dispersion coefficients for sp(2)-hybridized carbon atoms play a major role in aromatic stacking. Modern dispersion-corrected DFT yields accurate results (about 5-10% error for the dimerization energy) for the relatively large porphyrin and coronene dimers, systems for which WFT can provide accurate reference data only with huge computational effort. In this example, it is also demonstrated that new nonlocal, density-dependent dispersion corrections and atom pairwise schemes mutually agree with each other. The dispersion energy is also important for the complex inter- and intramolecular interactions that arise in the molecular crystals of aromatic molecules. In studies of hexahelicene, dispersion-corrected DFT yields "the right answer for the right reason". By comparison, standard DFT calculations reproduce intramolecular distances quite accurately in single-molecule calculations while inter- and intramolecular distances become too large when dispersion-uncorrected solid-state calculations are carried out. Dispersion-corrected DFT can fix this problem, and these results are in excellent agreement with experimental structure and energetic (sublimation) data. Uncorrected treatments do not even yield a bound crystal state. Finally, we present calculations for the formation of a cationic, quadruply charged dimer of a porphyrin derivative, a case where dispersion is required in order to overcome strong electrostatic repulsion. A combination of dispersion-corrected DFT with an adequate continuum solvation model can accurately reproduce experimental free association enthalpies in solution. As in the previous examples, consideration of the electrostatic interactions alone does not provide a qualitatively or quantitatively correct picture of the interactions of this complex.

Symmetry-adapted perturbation theory based on unrestricted Kohn-Sham orbitals for high-spin open-shell van der Waals complexes

Two open-shell formulations of the symmetry-adapted perturbation theory are presented. They are based on the spin-unrestricted Kohn-Sham (SAPT(UKS)) and unrestricted Hartree-Fock (SAPT(UHF)) descriptions of the monomers, respectively. The key reason behind development of SAPT(UKS) is that it is more compatible with density functional theory (DFT) compared to the previous formulation of open-shell SAPT based on spin-restricted Kohn-Sham method of Żuchowski et al. [J. Chem. Phys. 129, 084101 (2008)]. The performance of SAPT(UKS) and SAPT(UHF) is tested for the following open-shell van der Waals complexes: He···NH, H(2)O···HO(2), He···OH, Ar···OH, Ar···NO. The results show an excellent agreement between SAPT(UKS) and SAPT(ROKS). Furthermore, for the first time SAPT based on DFT is shown to be suitable for the treatment of interactions involving Π-state radicals (He···OH, Ar···OH, Ar···NO). In the interactions of transition metal dimers ((3)Σ(u)(+))Au(2) and ((13)Σ(g)(+))Cr(2) we show that SAPT is incompatible with the use of effective core potentials. The interaction energies of both systems expressed instead as supermolecular UHF interaction plus dispersion from SAPT(UKS) result in reasonably accurate potential curves.

Physically-motivated force fields from symmetry-adapted perturbation theory

We present a general methodology for generating accurate and transferable ab initio force fields, employing the framework of symmetry adapted perturbation theory (SAPT). The resulting force fields are "physically-motivated" in that they contain separate, explicit terms to account for the various fundamental intermolecular interactions, such as exchange, electrostatics, induction, and dispersion, with each term parametrized to a corresponding term in the SAPT energy decomposition. Crucially, the resulting force fields are largely compatible with existing, standard simulation packages, requiring only minimal modifications. We present several novel parametrization techniques that yield robust, physically meaningful atomic parameters that are transferable between molecular environments. We demonstrate the accuracy and generality of our method by validating against experimental second virial coefficients for a variety of small molecules. We then show that the resulting atomic parameters can be combined using physically motivated ansatzes to accurately predict arbitrary heteromolecular interaction energies, with example applications including prediction of gas adsorption in functionalized metal-organic framework materials.

Quantum-mechanical analysis of the energetic contributions to π stacking in nucleic acids versus rise, twist, and slide

Symmetry-adapted perturbation theory (SAPT) is applied to pairs of hydrogen-bonded nucleobases to obtain the energetic components of base stacking (electrostatic, exchange-repulsion, induction/polarization, and London dispersion interactions) and how they vary as a function of the helical parameters Rise, Twist, and Slide. Computed average values of Rise and Twist agree well with experimental data for B-form DNA from the Nucleic Acids Database, even though the model computations omitted the backbone atoms (suggesting that the backbone in B-form DNA is compatible with having the bases adopt their ideal stacking geometries). London dispersion forces are the most important attractive component in base stacking, followed by electrostatic interactions. At values of Rise typical of those in DNA (3.36 Å), the electrostatic contribution is nearly always attractive, providing further evidence for the importance of charge-penetration effects in π-π interactions (a term neglected in classical force fields). Comparison of the computed stacking energies with those from model complexes made of the "parent" nucleobases purine and 2-pyrimidone indicates that chemical substituents in DNA and RNA account for 20-40% of the base-stacking energy. A lack of correspondence between the SAPT results and experiment for Slide in RNA base-pair steps suggests that the backbone plays a larger role in determining stacking geometries in RNA than in B-form DNA. In comparisons of base-pair steps with thymine versus uracil, the thymine methyl group tends to enhance the strength of the stacking interaction through a combination of dispersion and electrosatic interactions.

Dispersion corrected hartree-fock and density functional theory for organic crystal structure prediction

We present and evaluate dispersion corrected Hartree-Fock (HF) and Density Functional Theory (DFT) based quantum chemical methods for organic crystal structure prediction. The necessity of correcting for missing long-range electron correlation, also known as van der Waals (vdW) interaction, is pointed out and some methodological issues such as inclusion of three-body dispersion terms are discussed. One of the most efficient and widely used methods is the semi-classical dispersion correction D3. Its applicability for the calculation of sublimation energies is investigated for the benchmark set X23 consisting of 23 small organic crystals. For PBE-D3 the mean absolute deviation (MAD) is below the estimated experimental uncertainty of 1.3 kcal/mol. For two larger π-systems, the equilibrium crystal geometry is investigated and very good agreement with experimental data is found. Since these calculations are carried out with huge plane-wave basis sets they are rather time consuming and routinely applicable only to systems with less than about 200 atoms in the unit cell. Aiming at crystal structure prediction, which involves screening of many structures, a pre-sorting with faster methods is mandatory. Small, atom-centered basis sets can speed up the computation significantly but they suffer greatly from basis set errors. We present the recently developed geometrical counterpoise correction gCP. It is a fast semi-empirical method which corrects for most of the inter- and intramolecular basis set superposition error. For HF calculations with nearly minimal basis sets, we additionally correct for short-range basis incompleteness. We combine all three terms in the HF-3c denoted scheme which performs very well for the X23 sublimation energies with an MAD of only 1.5 kcal/mol, which is close to the huge basis set DFT-D3 result.

Ab initio calculations of the interaction between CO2 and the acetate ion

A series of ab initio calculations designed to investigate the interaction of CO(2) with acetate are presented. The lowest energy structure, AC-CO(2)-η(2), is predicted by CCSD(T)/aVTZ to be bound by -10.6 kcal/mol. Six of the bound complexes have binding energies on the order of -8 kcal/mol, but analysis shows that the η(1)-CT complex is fundamentally different from the others. The η(1)-CT complex is characterized by geometric distortion, large polarization and induction effects and charge transfer whereas the other five complexes have little geometric distortion and negligible charge transfer. The amount of charge that is transferred from the anion to the CO(2) in the η(1)-CT complex is estimated to be about half an electron by NPA, DMA, CHELPG, and Mulliken analyses, whereas the EDA-ALMO-CTA (B3LYP) approach predicts a charge transfer of 75 me(-). However, the transfer of this small amount of charge leads to an energy lowering of -56 kcal/mol, without which the complex would not be bound. The RI-MP2 geometries closely approximate those resulting from the CCSD optimizations, and the optimized second-order opposite spin (O2) method performs well for all the complexes except for the η(1)-CT complex. DFT methods do not reproduce all the ab initio geometries, binding energies and/or energy ordering of these complexes although the range-separated hybrid meta-GGA (M11) and nonlocal (VV10 and vdwDF10) functionals are shown to yield results significantly better than other functionals considered for this system. The fact that there is such variation among DFT methods has implications for DFT-based ab initio molecular dynamics simulations and for the parametrization of classical force fields based on DFT calculations.

Factors that distort the heme structure in Heme-Nitric Oxide/OXygen-Binding (H-NOX) protein domains. A theoretical study

DFT and dispersion-corrected DFT calculations were carried out to probe the factors that distort the heme structure in Heme-Nitric oxide/OXygen-binding (H-NOX) protein domains. Various model systems that include heme, heme+surrounding residues, and heme+surrounding residues+additional protein environment were examined; the latter system was calculated with a quantum mechanics/molecular mechanics (QM/MM) method. The computations were extended to a myoglobin (Mb) protein, in which the heme structure is quite planar, in contrast to that in H-NOX. The natural tendency of the heme is to be planar. The strong structural distortion in H-NOX is mainly brought about by the intermolecular interactions between the whole heme molecule (heme ring plus its peripheral substituents) and the surrounding residues, among which the polar residues (Tyr140, Pro115, Mse98) play major roles in distorting the heme structure. The two peripheral propionate substituents that are oriented on the same side of the heme plane can also make the molecule distort, but the distortion caused by this factor is not significant. In Mb, the surrounding residues considered are all nonpolar and do not cause a structural distortion. The different structural features of the heme macrocycle in the different proteins (H-NOX and Mb) are reproduced by the calculations. The dispersion correction is necessary, since it improves the calculated structures. The effects of the distortion on the binding affinity of the axial ligand to the heme were also examined.

Is it possible to derive quantitative information on polarization of electron density from the multipolar model?

The accuracy of electrostatic properties estimated from the Hansen-Coppens multipolar model was verified. Tests were carried out to determine whether the multipolar model is accurate enough to study changes of electrostatic properties under the influence of a crystal field. Perturbed and unperturbed electron densities of individual molecules of amino acids and dipeptides were obtained from cluster and perturbation theory calculations. This enabled the changes in electrostatic properties values caused by polarization of the electron density to be characterized. Multipolar models were then fitted to the subsequent theoretical electron densities. The study revealed that electrostatic properties obtained from the multipolar models are significantly different from those obtained directly from the theoretical densities. The electrostatic properties of isolated molecules are reproduced better by multipolar models than the electrostatic properties of molecules in a crystal. Changes of electrostatic properties caused by perturbation of electron density due to the crystal environment are barely described by the multipolar model. As a consequence, the electrostatic properties obtained from multipolar models fitted to the perturbed theoretical densities derived either from cluster or periodic calculations do not differ much from those estimated from multipolar models fitted to densities of isolated molecules. The main reason for this seems to be related to an inadequate description of electron-density polarization in the vicinity of the nuclei by the multipolar model.

Large-scale symmetry-adapted perturbation theory computations via density fitting and Laplace transformation techniques: investigating the fundamental forces of DNA-intercalator interactions

Symmetry-adapted perturbation theory (SAPT) provides a means of probing the fundamental nature of intermolecular interactions. Low-orders of SAPT (here, SAPT0) are especially attractive since they provide qualitative (sometimes quantitative) results while remaining tractable for large systems. The application of density fitting and Laplace transformation techniques to SAPT0 can significantly reduce the expense associated with these computations and make even larger systems accessible. We present new factorizations of the SAPT0 equations with density-fitted two-electron integrals and the first application of Laplace transformations of energy denominators to SAPT. The improved scalability of the DF-SAPT0 implementation allows it to be applied to systems with more than 200 atoms and 2800 basis functions. The Laplace-transformed energy denominators are compared to analogous partial Cholesky decompositions of the energy denominator tensor. Application of our new DF-SAPT0 program to the intercalation of DNA by proflavine has allowed us to determine the nature of the proflavine-DNA interaction. Overall, the proflavine-DNA interaction contains important contributions from both electrostatics and dispersion. The energetics of the intercalator interaction are are dominated by the stacking interactions (two-thirds of the total), but contain important contributions from the intercalator-backbone interactions. It is hypothesized that the geometry of the complex will be determined by the interactions of the intercalator with the backbone, because by shifting toward one side of the backbone, the intercalator can form two long hydrogen-bonding type interactions. The long-range interactions between the intercalator and the next-nearest base pairs appear to be negligible, justifying the use of truncated DNA models in computational studies of intercalation interaction energies.

An automated approach for the parameterization of accurate intermolecular force-fields: pyridine as a case study

An automated protocol is proposed and validated, which integrates accurate quantum mechanical calculations with classical numerical simulations. Intermolecular force fields, (FF) suitable for molecular dynamics (MD) and Monte Carlo simulations, are parameterized through a novel iterative approach, fully based on quantum mechanical data, which has been automated and coded into the PICKY software, here presented. The whole procedure is tested and validated for pyridine, whose bulk phase, described through MD simulations performed with the specifically parameterized FF, is characterized by computing several of its thermodynamic, structural, and transport properties, comparing them with their experimental counterparts.