SAPT codes for calculations of intermolecular interaction energies.

Symmetry-adapted perturbation theory (SAPT) is a method for calculations of intermolecular (noncovalent) interaction energies. The set of SAPT codes that is described here, the current version named SAPT2020, includes virtually all variants of SAPT developed so far, among them two-body SAPT based on perturbative, coupled cluster, and density functional theory descriptions of monomers, three-body SAPT, and two-body SAPT for some classes of open-shell monomers. The properties of systems governed by noncovalent interactions can be predicted only if potential energy surfaces (force fields) are available. SAPT is the preferred approach for generating such surfaces since it is seamlessly connected to the asymptotic expansion of interaction energy. SAPT2020 includes codes for automatic development of such surfaces, enabling generation of complete dimer surfaces with a rigid monomer approximation for dimers containing about one hundred atoms. These codes can also be used to obtain surfaces including internal degrees of freedom of monomers.

Correcting long-range electrostatics in DFTB.

We demonstrate that the atom-based charge model implemented in the current versions of the density functional tight binding (DFTB) method fails to reproduce the correct charge distribution of a range of systems, including homonuclear molecules, graphene, and nanotubes, resulting in serious distortions in the electrostatic interactions for such systems caused by the missing quadrupole moments. In particular, this failure seriously impacts the long- and medium-range interaction energies of the DFTB plus dispersion (DFTB-D) model, leading to incorrect predictions of translational or rotational barriers in such systems. We show explicitly on examples of H2 and N2 that correct quadrupole moments-and consequently correct electrostatic interactions-can be restored in such systems by adding additional bond (ghost) sites to the homonuclear molecules. Attempts to determine the point charges associated with the additional sites using the usual Mulliken population analysis lead to unphysical results. Instead, these charges can be determined using the actual DFTB densities used in the parameterization process. For homonuclear molecules, we propose an extension to the DFTB-D model by adding charges that reproduce the physically correct quadrupolar charge distribution. The resulting DFTB-D-Q model greatly improves the rotational barriers for interactions of molecular hydrogen and nitrogen with benzene.

Dispersion-Corrected DFT Struggles with Predicting Three-Body Interaction Energies.

We demonstrate that the dispersion-corrected density functional theory (DFT-D) schemes fall short of predicting reliable three-body interaction energies. This concerns also a popular variant of DFT-D called the "many-body dispersion" (MBD) method, which might seem surprising in the light of the fact that its name contains the very phrase "many-body". The main reason for the inaccuracy of the three-body interaction energies in the DFT-D schemes can be attributed to internal deficiencies of the standard DFT functionals that the existing "-D" methods are incapable of correcting since the main problems emerge from the terms not related to the dispersion component. At present, it seems that none of the a posteriori dispersion techniques are able to predict accurately the total interaction energy for a supermolecular system together with its simultaneous decomposition into the many-body components. On the other hand, if one is interested only in the three-body interaction energies, we propose an adjustment to the MBD approach that achieves good accuracy in conjunction with the supermolecular MP2.

Blind test of density-functional-based methods on intermolecular interaction energies.

In the past decade, a number of approaches have been developed to fix the failure of (semi)local density-functional theory (DFT) in describing intermolecular interactions. The performance of several such approaches with respect to highly accurate benchmarks is compared here on a set of separation-dependent interaction energies for ten dimers. Since the benchmarks were unknown before the DFT-based results were collected, this comparison constitutes a blind test of these methods.

Predictions for water clusters from a first-principles two- and three-body force field.

A new rigid-monomer three-body potential has been developed for water by fitting it to more than 70 thousand trimer interaction energies computed ab initio using coupled-cluster methods and augmented triple-zeta-quality basis sets. This potential was used together with a modified form of a previously developed two-body potential and with a polarization model of four- and higher-body interactions to predict the energetics of the water trimer, hexamer, and 24-mer. Despite using the rigid-monomer approximation, these predictions agree better with flexible-monomer benchmarks than published results obtained with flexible-monomer force fields. An unexpected finding of our work is that simple polarization models predict four-body interactions to within a few percent, whereas for three-body interactions these models are known to have errors on the order of 50%.

Communication: resolving the three-body contribution to the lattice energy of crystalline benzene: benchmark results from coupled-cluster theory.

Coupled-cluster theory including single, double, and perturbative triple excitations [CCSD(T)] has been applied to trimers that appear in crystalline benzene in order to resolve discrepancies in the literature about the magnitude of non-additive three-body contributions to the lattice energy. The present results indicate a non-additive three-body contribution of 0.89 kcal mol(-1), or 7.2% of the revised lattice energy of -12.3 kcal mol(-1). For the trimers for which we were able to compute CCSD(T) energies, we obtain a sizeable difference of 0.63 kcal mol(-1) between the CCSD(T) and MP2 three-body contributions to the lattice energy, confirming that three-body dispersion dominates over three-body induction. Taking this difference as an estimate of three-body dispersion for the closer trimers, and adding an Axilrod-Teller-Muto estimate of 0.13 kcal mol(-1) for long-range contributions yields an overall value of 0.76 kcal mol(-1) for three-body dispersion, a significantly smaller value than in several recent studies.

Localized overlap algorithm for unexpanded dispersion energies.

First-principles-based, linearly scaling algorithm has been developed for calculations of dispersion energies from frequency-dependent density susceptibility (FDDS) functions with account of charge-overlap effects. The transition densities in FDDSs are fitted by a set of auxiliary atom-centered functions. The terms in the dispersion energy expression involving products of such functions are computed using either the unexpanded (exact) formula or from inexpensive asymptotic expansions, depending on the location of these functions relative to the dimer configuration. This approach leads to significant savings of computational resources. In particular, for a dimer consisting of two elongated monomers with 81 atoms each in a head-to-head configuration, the most favorable case for our algorithm, a 43-fold speedup has been achieved while the approximate dispersion energy differs by less than 1% from that computed using the standard unexpanded approach. In contrast, the dispersion energy computed from the distributed asymptotic expansion differs by dozens of percent in the van der Waals minimum region. A further increase of the size of each monomer would result in only small increased costs since all the additional terms would be computed from the asymptotic expansion.

Communication: Density functional theory overcomes the failure of predicting intermolecular interaction energies.

Density-functional theory (DFT) revolutionized the ability of computational quantum mechanics to describe properties of matter and is by far the most often used method. However, all the standard variants of DFT fail to predict intermolecular interaction energies. In recent years, a number of ways to go around this problem has been proposed. We show that some of these approaches can reproduce interaction energies with median errors of only about 5% in the complete range of intermolecular configurations. Such errors are comparable to typical uncertainties of wave-function-based methods in practical applications. Thus, these DFT methods are expected to find broad applications in modelling of condensed phases and of biomolecules.

Efficient Calculations of Dispersion Energies for Nanoscale Systems from Coupled Density Response Functions.

Dispersion energies computed from coupled Kohn-Sham (CKS) dynamic density-density response functions are known to be highly accurate. At the same time, the computational algorithm is of only modest complexity compared to other accurate methods of dispersion energy calculation. We present a new implementation of this algorithm that removes several computational barriers present in current implementations and enables calculations of dispersion energies for systems with more than 200 atoms using more than 5000 basis functions. The improvements were mainly achieved by reorganizing the algorithm to minimize memory and disk usage. We present applications to two systems: the buckycatcher complex with fullerene and the vancomycin complex with a diacetyl-Lys-d-Ala-d-Ala bacterial wall precursor, both calculations performed with triple-ζ-quality basis sets. Our implementation makes it possible to use ab initio computed dispersion energies in popular "density functional theory plus dispersion" approaches.

Interaction energies of large clusters from many-body expansion.

In the canonical supermolecular approach, calculations of interaction energies for molecular clusters involve a calculation of the whole cluster, which becomes expensive as the cluster size increases. We propose a novel approach to this task by demonstrating that interaction energies of such clusters can be constructed from those of small subclusters with a much lower computational cost by applying progressively lower-level methods for subsequent terms in the many-body expansion. The efficiency of such "stratified approximation" many-body approach (SAMBA) is due to the rapid convergence of the many-body expansion for typical molecular clusters. The method has been applied to water clusters (H(2)O)(n), n = 6, 16, 24. For the hexamer, the best results that can be obtained with current computational resources in the canonical supermolecular method were reproduced to within about one tenth of the uncertainty of the canonical approach while using 24 times less computer time in the many-body expansion calculations. For (H(2)O)(24), SAMBA is particularly beneficial and we report interaction energies with accuracy that is currently impossible to obtain with the canonical supermolecular approach. Moreover, our results were computed using two orders of magnitude smaller computer resources than used in the previous best calculations for this system. We also show that the basis-set superposition errors should be removed in calculations for large clusters.

A molecular dynamics study of 1,1-diamino-2,2-dinitroethylene (FOX-7) crystal using a symmetry adapted perturbation theory-based intermolecular force field.

A dimer potential energy function for 1,1-diamino-2,2-dinitroethylene (FOX-7) has been developed using symmetry adapted perturbation theory based on a Kohn-Sham density functional theory description of the monomers [SAPT(DFT)]. Interaction energies of 1008 dimer configurations were computed in an augmented double zeta basis set and fitted to an atom-atom intermolecular potential energy function of Coulomb plus Buckingham exp-6 form. The potential was used in isothermal-isostress molecular dynamics simulations to study the structure and thermal/pressure response of FOX-7 crystal. The simulated structure is in very good agreement with experiment and the computed thermal/pressure response of the crystal shows significant anisotropy with respect to crystallographic direction, in-line with experimental observations. It is concluded that SAPT(DFT) is an excellent method for development of intermolecular potentials for energetic molecular crystals.

Vibration-rotation-tunneling states of the benzene dimer: an ab initio study.

An improved intermolecular potential surface for the benzene dimer is constructed from interaction energies computed by symmetry-adapted perturbation theory, SAPT(DFT), with the inclusion of third-order contributions. Twelve characteristic points on the surface have been investigated also using the coupled-cluster method with single, double, and perturbative triple excitations, CCSD(T), and triple-zeta quality basis sets with midbond functions. The SAPT and CCSD(T) results are in close agreement and provide the best representation of these points to date. The potential was used in calculations of vibration-rotation-tunneling (VRT) levels of the dimer by a method appropriate for large amplitude intermolecular motions and tunneling between multiple equivalent minima in the potential. The resulting VRT levels were analyzed with the use of the permutation-inversion full cluster tunneling (FCT) group G(576) and a chain of subgroups that starts from the molecular symmetry group C(s)(M) of the rigid dimer at its equilibrium C(s) geometry and leads to G(576) if all possible intermolecular tunneling mechanisms are feasible. Further information was extracted from the calculated wave functions. It was found, in agreement with the experimental data, that for all of the 54 G(576) symmetry species (with different nuclear spin statistical weights) the lower VRT states have a tilted T-shape (TT) structure; states with the parallel-displaced structure are higher in energy than the ground state of A symmetry by at least 30 cm(-1). The dissociation energy D(0) equals 870 cm(-1), while the depth D(e) of the TT minimum in the potential is 975 cm(-1). Hindered rotation of the cap in the TT structure and tilt tunneling lead to level splittings on the order of 1 cm(-1). Also intermolecular vibrations with excitation energies starting at a few cm(-1) were identified. A further small, but probably significant, level splitting was assigned to cap turnover, although in scans of the potential surface we could not find a plausible 'reaction path' for this process. Rotational constants were extracted from energy levels calculated for total angular momentum J = 0 and 1, and from expectation values of the inertia tensor. Although the end-over-end rotational constant B + C agrees well with the measured microwave spectra, there is disagreement with the measurements concerning the (a)symmetric rotor character of the benzene dimer. It is concluded from calculations for the 54 nuclear spin species that the microwave spectrum should show overlapping contributions from many different species. Another interesting conclusion regards the role of the quantum number K, for a prolate near-symmetric rotor the projection of the total angular momentum on the prolate axis. For the benzene dimer, K has a substantial effect on the energy levels associated with the intermolecular motions of the complex.

Improved interaction energy benchmarks for dimers of biological relevance.

The set of interaction energies for 22 systems of biological importance, developed by Jurecka et al. [Phys. Chem. Chem. Phys., 2006, 8, 1985] and called S22, became an often used benchmark for evaluating the performance of various computational methods. As the quality of such methods improves, the uncertainties of the S22 energies are becoming too large to enable meaningful comparisons. We therefore improved the benchmarks by performing calculations in larger basis sets than used by Jurecka et al. The basis extensions included additions of higher angular momentum, diffuse, and midbond functions. The percentage deviations of the original S22 interaction energies from our values are up to 15.4%. We have estimated that the average (unsigned) uncertainty of our results should be about 1.0%. This estimate includes contributions from basis set truncation, frozen-core approximation, and neglected electron excitations.

Interactions of graphene sheets deduced from properties of polycyclic aromatic hydrocarbons.

Intermolecular interactions of coronene dimer were studied with symmetry-adapted perturbation theory based on the density functional theory description of the monomers [SAPT(DFT)]. The most stable stacked structure was found to have the interaction energy of -17.45 kcal/mol, slightly lower than the structure analogous to graphite (-17.36 kcal/mol). The latter energy was extrapolated to the interaction energy of two graphene sheets. The effects of interactions of multiple layers were also estimated leading to the exfoliation energy of graphite equal to 45.3 meV per carbon atom. The SAPT(DFT)-based decomposition into physical quantities of the interaction energies shows the dominant effect of the dispersion interactions with a weaker electrostatic contribution due to penetration effects. The extrapolated physical picture of the graphene-graphene interaction is very similar to that of smaller stacked polycyclic aromatic hydrocarbons.

Ab initio potential energy surfaces for NH((3)sigma(-))-NH((3)sigma(-)) with analytical long range.

We present four-dimensional ab initio potential energy surfaces for the three different spin states of the NH((3)Sigma(-))-NH((3)Sigma(-)) complex. The potentials are partially based on the work of Dhont et al. [J. Chem. Phys. 123, 184302 (2005)]. The surface for the quintet state is obtained at the RCCSD(T)/augmented correlation-consistent polarized valence triple-zeta (aug-cc-pVTZ) level of theory and the energy differences with the singlet and triplet states are calculated at the complete active space with nth-order perturbation theory/aug-cc-pVTZ (n=2,3) level of theory. The ab initio potentials are fitted to coupled spherical harmonics in the angular coordinates, and the long range is further expanded as a power series in 1/R. The RCCSD(T) potential is corrected for a size-consistency error of about 0.5x10(-6) E(h) prior to fitting. The long-range coefficients obtained from the fit are found to be in good agreement with first and second-order perturbation theory calculations.

Crystal structure prediction for cyclotrimethylene trinitramine (RDX) from first principles.

Crystal structure prediction and molecular dynamics methods were applied to the cyclotrimethylene trinitramine (RDX) crystal to explore the stability rankings of various polymorphs using a recently developed nonempirical potential energy function that describes the RDX dimer interactions. The energies of 500 high-density structures resulting from molecular packing were minimized and the 14 lowest-energy structures were subjected to isothermal-isostress molecular dynamics (NsT-MD) simulations. For both crystal structure prediction methods and molecular dynamics simulations, the lowest-energy polymorph corresponded to the experimental structure; furthermore, the lattice energy of this polymorph was lower than that of the other polymorphs by at least 1.1 kcal mol(-1). Crystal parameters and densities of the low-energy crystal produced by the NsT-MD simulations matched those of the experimental crystal to within 1% of density and cell edge lengths and 0.01 degrees of the cell angle. The arrangement of the molecules within the time-averaged unit cell were in equally outstanding agreement with experiment, with the largest deviation of the location of the molecular mass centers being less than 0.07 A and the largest deviation in molecular orientation being less than 2.8 degrees . NsT-MD simulations were also used to calculate crystallographic parameters as functions of temperature and pressure and the results were in a reasonable agreement with experiment.

Dispersionless density functional theory.

A new density functional (DF) method is proposed for calculations of intermolecular interaction energies. The exchange-correlation functional was optimized in such a way that the method recovers the interaction energies with the dispersion (including exchange-dispersion) component subtracted and therefore our approach is named the dispersionless DF (dlDF) method. The dlDF method is shown to predict very well the dispersionless part of the interaction energy for all types of intermolecular interactions. Thus, if combined with a dispersion component, computed ab initio or from a simple function fitted to ab initio values, it provides accurate and physically justified interaction energies in the whole range of intermolecular separations. Our dispersion function is significantly more accurate than the published ones.

Symmetry-adapted perturbation theory utilizing density functional description of monomers for high-spin open-shell complexes.

We present an implementation of symmetry-adapted perturbation theory (SAPT) to interactions of high-spin open-shell monomers forming high-spin dimers. The monomer spin-orbitals used in the expressions for the electrostatic and exchange contributions to the interaction energy are obtained from density functional theory using a spin-restricted formulation of the open-shell Kohn-Sham (ROKS) method. The dispersion and induction energies are expressed through the density-density response functions predicted by the time-dependent ROKS theory. The method was applied to several systems: NH...He, CN...Ne, H2O...HO2, and NH...NH. It provides accuracy comparable to that of the best previously available methods such as the open-shell coupled-cluster method with single, double, and noniterative triple excitations, RCCSD(T), with a significantly reduced computational cost.