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

Benchmarking Swaths of Intermolecular Interaction Components with Symmetry-Adapted Perturbation Theory

A benchmark database for interaction energy components of various noncovalent interactions (NCIs) along their dissociation curve is one of the essential needs in theoretical chemistry, especially for the development of force fields and machine-learning methods. We utilize DFT-SAPT or SAPT(DFT) as one of the most accurate methods to generate an extensive stock of the energy components, including dispersion energies extrapolated to the complete basis set limit (CBS). Precise analyses of the created data, and benchmarking the total interaction energies against the best available CCSD(T)/CBS values, reveal different aspects of the methodology and the nature of NCIs. For example, error cancellation effects between the S2 approximation and nonexact xc-potentials occur, and large charge transfer energies in some systems, including heavy atoms, can explain the lower accuracy of DFT-SAPT. This method is perfect for neutral complexes containing light nonmetals, while other systems with heavier atoms should be treated carefully. In the last part, a representative data set for all NCIs is extracted from the original data.

PhyNEO: A Neural-Network-Enhanced Physics-Driven Force Field Development Workflow for Bulk Organic Molecule and Polymer Simulations

An accurate, generalizable, and transferable force field plays a crucial role in the molecular dynamics simulations of organic polymers and biomolecules. Conventional empirical force fields often fail to capture precise intermolecular interactions due to their negligence of important physics, such as polarization, charge penetration, many-body dispersion, etc. Moreover, the parameterization of these force fields relies heavily on top-down fittings, limiting their transferabilities to new systems where the experimental data are often unavailable. To address these challenges, we introduce a general and fully ab initio force field construction strategy, named PhyNEO. It features a hybrid approach that combines both the physics-driven and the data-driven methods and is able to generate a bulk potential with chemical accuracy using only quantum chemistry data of very small clusters. Careful separations of long-/short-range interactions and nonbonding/bonding interactions are the key to the success of PhyNEO. By such a strategy, we mitigate the limitations of pure data-driven methods in long-range interactions, thus largely increasing the data efficiency and the scalability of machine learning models. The new approach is thoroughly tested on poly(ethylene oxide) and polyethylene glycol systems, giving superior accuracies in both microscopic and bulk properties compared to conventional force fields. This work thus offers a promising framework for the development of advanced force fields in a wide range of organic molecular systems.

Reliable Dimerization Energies for Modeling of Supramolecular Junctions

Accurate estimates of intermolecular interaction energy, ΔE, are crucial for modeling the properties of organic electronic materials and many other systems. For a diverse set of 50 dimers comprising up to 50 atoms (Set50-50, with 7 of its members being models of single-stacking junctions), benchmark ΔE data were compiled. They were obtained by the focal-point strategy, which involves computations using the canonical variant of the coupled cluster theory with singles, doubles, and perturbative triples [CCSD(T)] performed while applying a large basis set, along with extrapolations of the respective energy components to the complete basis set (CBS) limit. The resulting ΔE data were used to gauge the performance for the Set50-50 of several density-functional theory (DFT)-based approaches, and of one of the localized variants of the CCSD(T) method. This evaluation revealed that (1) the proposed "silver standard" approach, which employs the localized CCSD(T) method and CBS extrapolations, can be expected to provide accuracy better than two kJ/mol for absolute values of ΔE, and (2) from among the DFT techniques, computationally by far the cheapest approach (termed "ωB97X-3c/vDZP" by its authors) performed remarkably well. These findings are directly applicable in cost-effective yet reliable searches of the potential energy surfaces of noncovalent complexes.

Using Noncovalent Interactions to Test the Precision of Projector-Augmented Wave Data Sets

The projector-augmented wave (PAW) method is one of the approaches that are widely used to approximately treat core electrons and thus to speed up plane-wave basis set electronic structure calculations. However, PAW involves approximations, and it is thus important to understand how they affect the results. Tests of the precision of PAW data sets often use the properties of isolated atoms or atomic solids. While this is sufficient to identify problematic PAW data sets, little information has been gained to understand the origins of the errors and suggest ways to correct them. Here, we show that the interaction energies of molecular dimers are very useful not only to identify problematic PAW data sets but also to uncover the origin of the errors. Using dimers from the S22 and S66 test sets and other dimers, we find that the error in the interaction energy is composed of a short-range component with an exponential decay and a long-range electrostatic part caused by an error in the total charge density. We propose and evaluate a simple improvable scheme to correct the long-range error and find that even in its simple and readily usable form, it is able to reduce the interaction energy errors to less than half on average for hydrogen-bonded dimers.

H2O-HCN complex: A new potential energy surface and intermolecular rovibrational states from rigorous quantum calculations

In this work the H2O-HCN complex is quantitatively characterized in two ways. First, we report a new rigid-monomer 5D intermolecular potential energy surface (PES) for this complex, calculated using the symmetry-adapted perturbation theory based on density functional theory method. The PES is based on 2833 ab initio points computed employing the aug-cc-pVQZ basis set, utilizing the autoPES code, which provides a site-site analytical fit with the long-range region given by perturbation theory. Next, we present the results of the quantum 5D calculations of the fully coupled intermolecular rovibrational states of the H2O-HCN complex for the total angular momentum J values of 0, 1, and 2, performed on the new PES. These calculations rely on the quantum bound-state methodology developed by us recently and applied to a variety of noncovalently bound binary molecular complexes. The vibrationally averaged ground-state geometry of H2O-HCN determined from the quantum 5D calculations agrees very well with that from the microwave spectroscopic measurements. In addition, the computed ground-state rotational transition frequencies, as well as the B and C rotational constants calculated for the ground state of the complex, are in excellent agreement with the experimental values. The assignment of the calculated intermolecular vibrational states of the H2O-HCN complex is surprisingly challenging. It turns out that only the excitations of the intermolecular stretch mode can be assigned with confidence. The coupling among the angular degrees of freedom (DOFs) of the complex is unusually strong, and as a result most of the excited intermolecular states are unassigned. On the other hand, the coupling of the radial, intermolecular stretch mode and the angular DOFs is weak, allowing straightforward assignment of the excitation of the former.

Accurate three-body noncovalent interactions: the insights from energy decomposition

An impressive collection of accurate two-body interaction energies for small complexes has been assembled into benchmark databases and used to improve the performance of multiple density functional, semiempirical, and machine learning methods. Similar benchmark data on nonadditive three-body energies in molecular trimers are comparatively scarce, and the existing ones are practically limited to homotrimers. In this work, we present a benchmark dataset of 20 equilibrium noncovalent interaction energies for a small but diverse selection of 10 heteromolecular trimers. The new 3BHET dataset presents complexes that combine different interactions including π-π, anion-π, cation-π, and various motifs of hydrogen and halogen bonding in each trimer. A detailed symmetry-adapted perturbation theory (SAPT)-based energy decomposition of the two- and three-body interaction energies shows that 3BHET consists of electrostatics- and dispersion-dominated complexes. The nonadditive three-body contribution is dominated by induction, but its influence on the overall bonding type in the complex (as exemplified by its position on the ternary diagram) is quite small. We also tested the extended SAPT (XSAPT) approach which is capable of including some nonadditive interactions in clusters of any size. The resulting three-body dispersion term (obtained from the many-body dispersion formalism) is mostly in good agreement with the supermolecular CCSD(T)-MP2 values and the nonadditive induction term is similar to the three-body SAPT(DFT) data, but the overall three-body XSAPT energies are not very accurate as they are missing the first-order exchange terms.

Efficient Calculation of the Dispersion Energy for Multireference Systems with Cholesky Decomposition: Application to Excited-State Interactions

Accurate and efficient prediction of dispersion interactions in excited-state complexes poses a challenge due to the complex nature of electron correlation effects that need to be simultaneously considered. We propose an algorithm for computing the dispersion energy in nondegenerate ground- or excited-state complexes with arbitrary spin. The algorithm scales with the fifth power of the system size due to employing Cholesky decomposition of Coulomb integrals and a recently developed recursive formula for density response functions of the monomers. As a numerical illustration, we apply the new algorithm in the framework of multiconfigurational symmetry adapted perturbation theory, SAPT(MC), to study interactions in dimers with localized excitons. The SAPT(MC) analysis reveals that the dispersion energy may be the main force stabilizing excited-state dimers.

Graphene as Nanocarrier for Gold(I)-Monocarbene Complexes: Strength and Nature of Physisorption

Gold(I) metal complexes are finding increasing applications as therapeutic agents against a variety of diseases. As their potential use as effective metallodrugs is continuously confirmed, the issue of their administration, distribution and delivery to desired biological targets emerges. Graphene and its derivatives possess attractive properties in terms of high affinity and low toxicity, suggesting that they can efficaciously be used as drug nanocarriers. In the present study, we computationally address the adsorption of a gold(I) N-heterocyclic monocarbene, namely, IMeAuCl (where IMe = 1,3-dimethylimidazol-2-ylidene), on graphene. The Au(I) N-heterocyclic carbene family has indeed shown promising anticancer activity and the N-heterocyclic ring could easily interact with planar graphene nanostructures. By means of high-level electronic structure approaches, we investigated the strength and nature of the involved interaction using small graphene prototypes, which allow us to benchmark the best-performing DFT functionals as well as assess the role of the different contributions to total interaction energies. Moreover, realistic adsorption enthalpies and free energy values are obtained by exploiting the optimal DFT method to describe the drug adsorption on larger graphene models. Such values (ΔH_ads = -18.4 kcal/mol and ΔG_ads= -7.20 kcal/mol for the largest C₁₅₀H₃₀ model) indicate a very favorable adsorption, mainly arising from the dispersion component of the interaction, with the electrostatic attraction also playing a non-negligible role.

High-Order Quantum-Mechanical Analysis of Hydrogen Bonding in Hachimoji and Natural DNA Base Pairs

High-order quantum chemistry is applied to hydrogen-bonded natural DNA nucleobase pairs [adenine:thymine (A:T) and guanine:cytosine (G:C)] and non-natural Hachimoji nucleobase pairs [isoguanine:1-methylcytosine (B:S) and 2-aminoimidazo[1,2a][1,3,5]triazin-4(1H)-one:6-amino-5-nitropyridin-2-one (P:Z)] to see how the intermolecular interaction energies and their energetic components (electrostatics, exchange-repulsion, induction/polarization, and London dispersion interactions) vary among the base pairs. We examined the Hoogsteen (HG) geometries in addition to the traditional Watson-Crick (WC) geometries. Coupled-cluster theory through perturbative triples [CCSD(T)] extrapolated to the complete basis set (CBS) limit and high-order symmetry-adapted perturbation theory (SAPT) at the SAPT2+(3)(CCD)δMP2/aug-cc-pVTZ level are used to estimate highly accurate noncovalent interaction energies. Electrostatic interactions are the most attractive component of the interaction energies, but the sum of induction/polarization and London dispersion is nearly as large, for all base pairs and geometries considered. Interestingly, the non-natural Hachimoji base pairs interact more strongly than the corresponding natural base pairs, by -21.8 (B:S) and -0.3 (P:Z) kcal mol^-1 in the WC geometries, according to CCSD(T)/CBS. This is consistent with the H-bond distances being generally shorter in the non-natural base pairs. The natural base pairs are energetically more stabilized in their Hoogsteen geometries than in their WC geometries. The Hoogsteen geometry makes the A:T base pair slightly more stable, by -0.8 kcal mol^-1, and it greatly stabilizes the G:C⁺ base pair, by -15.3 kcal mol^-1. The G:C⁺ stabilization is mainly due to the fact that C has typically added a proton when found in Hoogsteen geometries. By contrast, Hoogsteen geometries are substantially less favorable than WC geometries for non-natural Hachimoji base pairs, by 17.3 (B:S) and 13.8 (P:Z) kcal mol^-1.

Long-range interactions of aromatic molecules with alkali-metal and alkaline-earth-metal atoms

The isotropic and anisotropic coefficients C_n ^l,m of the long-range spherical expansion ∼1/Rⁿ (R-the intermolecular distance) of the dispersion and induction intermolecular energies are calculated using the first principles for the complexes containing an aromatic molecule (benzene, pyridine, furan, and pyrrole) and alkali-metal (Li, Na, K, Rb, and Cs) or alkaline-earth-metal (Be, Mg, Ca, Sr, and Ba) atoms in their electronic ground states. The values of the first- and second-order properties of the aromatic molecules are calculated using the response theory with the asymptotically corrected LPBE0 functional. The second-order properties of the closed-shell alkaline-earth-metal atoms are obtained using the expectation-value coupled cluster theory and of the open-shell alkali-metal atoms using analytical wavefunctions. These properties are used for the calculation of the dispersion C_n,disp ^l,m and induction C_n,ind ^l,m coefficients (C_n ^l,m=C_n,disp ^l,m+C_n,ind ^l,m) with n up to 12 using the available implemented analytical formulas. It is shown that the inclusion of the coefficients with n > 6 is important for reproducing the interaction energy in the van der Waals region at R ≈ 6 Å. The reported long-range potentials should be useful for constructing the analytical potentials valid for the whole intermolecular interaction range, which are needed for spectroscopic and scattering studies.

Assessment of three-body dispersion models against coupled-cluster benchmarks for crystalline benzene, carbon dioxide, and triazine

To study the contribution of three-body dispersion to crystal lattice energies, we compute the three-body contributions to the lattice energies for crystalline benzene, carbon dioxide, and triazine using various computational methods. We show that these contributions converge quickly as the intermolecular distances between the monomers grow. In particular, the smallest value among the three pairwise intermonomer closest-contact distances, R_min, shows a strong correlation with the three-body contribution to the lattice energy, and, here, the largest of the closest-contact distances, R_max, serves as a cutoff criterion to limit the number of trimers to be considered. We considered all trimers up to R_max=15Å. The trimers with R_min<4Å contribute 90.4%, 90.6%, and 93.9% of the total three-body contributions for crystalline benzene, carbon dioxide, and triazine, respectively, for the coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] method. For trimers with R_min>4Å, the second-order Møller-Plesset perturbation theory (MP2) supplemented with the Axilrod-Teller-Muto (ATM) three-body dispersion correction reproduces the CCSD(T) values for the cumulative three-body contributions with errors of less than 0.1 kJ mol^-1. Moreover, three-body contributions are converged within 0.15 kJ mol^-1 by R_max=10Å. From these results, it appears that in molecular crystals where dispersion dominates the three-body contribution to the lattice energy, the trimers with R_min>4Å can be computed with the MP2+ATM method to reduce the computational cost, and those with R_max>10Å appear to be basically negligible.

Understanding the Role of Intramolecular Ion-Pair Interactions in Conformational Stability Using an Ab Initio Thermodynamic Cycle

Intramolecular ion-pair interactions yield shape and functionality to many molecules. With proper orientation, these interactions overcome steric factors and are responsible for the compact structures of several peptides. In this study, we present a thermodynamic cycle based on isoelectronic and alchemical mutation to estimate the intramolecular ion-pair interaction energy. We determine these energies for 26 benchmark molecules with common ion-pair combinations and compare them with results obtained using intramolecular symmetry-adapted perturbation theory. For systems with long linkers, the ion-pair energies evaluated using both approaches deviate by less than 2.5% in the vacuum phase. The thermodynamic cycle based on density functional theory facilitates calculations of salt-bridge interactions in model tripeptides with continuum/microsolvation modeling and four large peptides: 1EJG (crambin), 1BDK (bradykinin), 1L2Y (a mini-protein with a tryptophan cage), and 1SCO (a toxin from the scorpion venom).

Overcoming Artificial Multipoles in Intramolecular Symmetry-Adapted Perturbation Theory

Intramolecular symmetry-adapted perturbation theory (ISAPT) is a method to compute and decompose the noncovalent interaction energy between two molecular fragments A and B covalently connected via a linker C. However, the existing ISAPT algorithm displays several issues for many fragmentation patterns (that is, specific assignments of atoms to the A/B/C subsystems), including an artificially repulsive electrostatic energy (even when the fragments are hydrogen-bonded) and very large and mutually cancelling induction and exchange-induction terms. We attribute those issues to the presence of artificial dipole moments at the interfragment boundary, as the atoms of A and B directly connected to C are missing electrons on one of their hybrid orbitals. Therefore, we propose several new partitioning algorithms which reassign one electron, on a singly occupied link hybrid orbital, from C to each of A/B. Once the contributions from these link orbitals are added to fragment density matrices, the computation of ISAPT electrostatic, induction, and dispersion energies proceeds exactly as normal, and the exchange energy expressions need only minor modifications. Among the link partitioning algorithms introduced, the so-called ISAPT(SIAO1) approach (in which the link orbital is obtained by a projection onto the intrinsic atomic orbitals (IAOs) of a given fragment followed by orthogonalization to this fragment's occupied space) leads to reasonable values of all ISAPT corrections for all fragmentation patterns, and exhibits a fast and systematic basis set convergence. This improvement is made possible by a significant reduction in magnitude (even though not a complete elimination) of the unphysical dipole moments at the interfragment boundaries. We demonstrate the utility of the improved ISAPT partitioning by examining intramolecular interactions in several pentanediol isomers, examples of linear and branched alkanes, and the open and closed conformations of a family of N-arylimide molecular torsion balances.

Physical mechanisms of intermolecular interactions from symmetry-adapted perturbation theory

Symmetry-adapted perturbation theory (SAPT) is a method for computational studies of noncovalent interactions between molecules. This method will be discussed here from the perspective of establishing the paradigm for understanding mechanisms of intermolecular interactions. SAPT interaction energies are obtained as sums of several contributions. Each contribution possesses a clear physical interpretation as it results from some specific physical process. It also exhibits a specific dependence on the intermolecular separation R. The four major contributions are the electrostatic, induction, dispersion, and exchange energies, each due to a different mechanism, valid at any R. In addition, at large R, SAPT interaction energies are seamlessly connected with the corresponding terms in the asymptotic multipole expansion of interaction energy in inverse powers of R. Since such expansion explicitly depends on monomers' multipole moments and polarizabilities, this connection provides additional insights by rigorously relating interaction energies to monomers' properties.

Potential energy surfaces of a stacked dimer of benzene and its radical cation: what remains and what appears

A stepwise qualitative consideration of the reduction in symmetry for the highly symmetrical "right sandwich" and "twisted sandwich" structures of the stacked benzene dimer caused by the pseudo-Jahn-Teller effect made it possible to construct a scheme of the potential energy surface (PES) of this dimer. Thirty-six equivalent structures of minimum energy are ordered on this extremely flat surface, transforming into each other in an almost barrier-free manner. There are two kinds of these transformations, both of which are pseudorotation. The transformation pathways inherent in this neutral dimer are also characteristic of its radical cation (RC). Structural transformations of the RC are inextricably linked with changes in its electronic state since its stationary structures relate to different electronic states. (C₆H₆)₂⁺˙ exists in the form of two orbital isomers, each of which "pseudorotates" on its own area of the PES. During the pseudorotation, the SOMO distribution on the fragments of (C₆H₆)₂⁺˙ changes, which is identical to what occurs during the pseudorotation of the Jahn-Teller benzene RC. These areas are connected by pairwise interconversions of their minimum energy structures. The interconversion of the orbital isomers is a bypassing of conical intersections between corresponding electronic states, which occurs through a synchronized pseudorotation of the fragments. The conclusions of the qualitative consideration and the results of quantum chemical calculations of various levels performed for (C₆H₆)₂⁺˙ are in full agreement with each other. The revealed features of the structure of the PES of the two reference systems in studying the intermolecular and ion-molecule interactions are the basis for considering their more complex analogs, primarily their less symmetrical ones.

A general tight-binding based energy decomposition analysis scheme for intermolecular interactions in large molecules

In this work, a general tight-binding based energy decomposition analysis (EDA) scheme for intermolecular interactions is proposed. Different from the earlier version [Xu et al., J. Chem. Phys. 154, 194106 (2021)], the current tight-binding based density functional theory (DFTB)-EDA is capable of performing interaction analysis with all the self-consistent charge (SCC) type DFTB methods, including SCC-DFTB2/3 and GFN1/2-xTB, despite their different formulas and parameterization schemes. In DFTB-EDA, the total interaction energy is divided into frozen, polarization, and dispersion terms. The performance of DFTB-EDA with SCC-DFTB2/3 and GFN1/2-xTB for various interaction systems is discussed and assessed.

Implementation of Symmetry-Adapted Perturbation Theory based on density functional theory and using hybrid exchange-correlation kernels for dispersion terms

We report the implementation of a symmetry-adapted perturbation theory algorithm based on a density functional theory description of the monomers [SAPT(DFT)]. The implementation adopts a density-fitting treatment of hybrid exchange-correlation kernels to enable the description of monomers with hybrid functionals, as in the algorithm by Bukowski, Podeszwa, and Szalewicz [Chem. Phys. Lett. 414, 111(2005)]. We have improved the algorithm by increasing numerical stability with QR factorization, and optimized the computation of the exchange-correlation kernel with its 2-index density-fitted representation. The algorithm scales as O(N5) formally and is usable for systems with up to ∼3000 basis functions, as demonstrated forthe C60-buckycatcher complex with the aug-cc-pVDZ basis set. The hybrid-kernel-based SAPT(DFT) algorithm is shown to be as accurate as SAPT(DFT) implementations based on local effective exact exchange potentials obtained from the local Hartree-Fock (LHF) method, while avoiding the lower-scaling [O(N4)] but iterative and sometimes hard-to-converge LHF process. The hybrid-kernel algorithm outperforms Hartree-Fock-based SAPT (SAPT0) for the S66 test set, and its accuracy is comparable to the many-body perturbation theory based SAPT2+ approach, which scales as O(N7), although SAPT2+ exhibits a more narrow distribution of errors.

Reliable crystal structure predictions from first principles

An inexpensive and reliable method for molecular crystal structure predictions (CSPs) has been developed. The new CSP protocol starts from a two-dimensional graph of crystal’s monomer(s) and utilizes no experimental information. Using results of quantum mechanical calculations for molecular dimers, an accurate two-body, rigid-monomer ab initio-based force field (aiFF) for the crystal is developed. Since CSPs with aiFFs are essentially as expensive as with empirical FFs, tens of thousands of plausible polymorphs generated by the crystal packing procedures can be optimized. Here we show the robustness of this protocol which found the experimental crystal within the 20 most stable predicted polymorphs for each of the 15 investigated molecules. The ranking was further refined by performing periodic density-functional theory (DFT) plus dispersion correction (pDFT+D) calculations for these 20 top-ranked polymorphs, resulting in the experimental crystal ranked as number one for all the systems studied (and the second polymorph, if known, ranked in the top few). Alternatively, the polymorphs generated can be used to improve aiFFs, which also leads to rank one predictions. The proposed CSP protocol should result in aiFFs replacing empirical FFs in CSP research.

Developing theoretical frameworks to predict new polymorphs is highly desirable. Here the authors present an ab initio based force-field approach for crystal structure prediction offering a dramatic computational speed-up over fully ab initio schemes.

A computational inspection of the dissociation energy of mid-sized organic dimers

The gas-phase value of the dissociation energy ( D0) is a key parameter employed in both experimental and theoretical descriptions of noncovalent complexes. The D0 data were obtained for a set of mid-sized organic dimers in their global minima which was located using geometry optimizations that applied ample basis sets together with either the conventional second-order Møller–Plesset (MP2) method or several dispersion-corrected density-functional theory (DFT-D) schemes. The harmonic vibrational zero-point (VZP) and deformation energies from the MP2 calculations were combined with electronic energies from the coupled cluster theory with singles, doubles, and iterative triples [CCSD(T)] extrapolated to the complete basis set (CBS) limit to estimate D0 with the aim of inspecting values that were most recently measured, and an analogous comparison was performed using the DFT-D data. In at least one case (namely, for the aniline⋯methane cluster), the D0 estimate that employed the CCSD(T)/CBS energies differed from experiment in the way that could not be explained by a possible deficiency in the VZP contribution. Curiously, one of the DFT-D schemes (namely, the B3LYP-D3/def2-QZVPPD) was able to reproduce all measured D0 values to within 1.0 kJ/mol from experimental error bars. These findings show the need for further measurements and computations of some of the complexes. In order to facilitate such studies, the physical nature of intermolecular interactions in the investigated dimers was analyzed by means of the DFT-based symmetry-adapted perturbation theory.

Dispersion Interactions in Exciton-Localized States. Theory and Applications to π-π* and n-π* Excited States

We address the problem of intermolecular interaction energy calculations in molecular complexes with localized excitons. Our focus is on the correct representation of the dispersion energy. We derive an extended Casimir-Polder formula for direct computation of this contribution through second order in the intermolecular interaction operator V̂. An alternative formula, accurate to infinite order in V̂, is derived within the framework of the adiabatic connection (AC) theory. We also propose a new parametrization of the VV10 nonlocal correlation density functional, so that it corrects the CASSCF energy for the dispersion contribution and can be applied to excited-state complexes. A numerical investigation is carried out for benzene, pyridine, and peptide complexes with the local exciton corresponding to the lowest π–π* or n– π* states. The extended Casimir-Polder formula is implemented in the framework of multiconfigurational symmetry-adapted perturbation theory, SAPT(MC). A SAPT(MC) analysis shows that the creation of a localized exciton affects mostly the electrostatic component of the interaction energy of investigated complexes. Nevertheless, the changes in Pauli repulsion and dispersion energies cannot be neglected. We verify the performance of several perturbation- and AC-based methods. Best results are obtained with a range-separated variant of an approximate AC approach employing extended random phase approximation and CASSCF wave functions.

Factors influencing hydrogen peroxide versus water inclusion in molecular crystals

Hydrate formation is often unavoidable during crystallization, leading to performance degradation of pharmaceuticals and energetics. In some cases, water molecules trapped within crystal lattices can be substituted for hydrogen peroxide, improving the solubility of drugs and detonation performance of explosives. The present work compares hydrates and hydrogen peroxide solvates in two ways: (1) analyzing structural motifs present in crystal structures accessed from the Cambridge Structural Database and (2) developing potential energy surfaces for water and hydrogen peroxide interacting with functional groups of interest at geometries relevant to the solid state. By elucidating fundamental differences in local interactions that can be formed with molecules of hydrogen peroxide and/or water, the analyses presented here provide a foundation for the design and selection of candidate molecules for the formation of hydrogen peroxide solvates.

First-principles modelling of the new generation of subnanometric metal clusters: Recent case studies

The very recent development of highly selective techniques making possible the synthesis and experimental characterization of subnanometric (subnanometer-sized) metal clusters (even single atoms) is pushing our understanding far beyond the present knowledge in materials science, driving these clusters as a new generation of quantum materials at the lower bounds of nanotechnology. When the size of the metal cluster is reduced to a small number of atoms, the d-band of the metal splits into a subnanometric d-type molecular orbitals network in which all metal atoms are inter-connected, with the inter-connections having the length of a chemical bond (1-2 Å). These molecular characteristics are at the very core of the high stability and novel properties of the smallest metal clusters, with their integration into colloidal materials interacting with the environment having the potential to further boost their performance in applications such as luminescence, sensing, bioimaging, theranostics, energy conversion, catalysis, and photocatalysis. Through the presentation of very recent case studies, this Feature Article is aimed to illustrate how first-principles modelling, including methods beyond the state-of-the-art and an interplay with cutting-edge experiments, is helping to understand the special properties of these clusters at the most fundamental level. Moreover, it will be discussed how superfluid helium droplets can act both as nano-reactors and carriers to achieve the synthesis and surface deposition of metal clusters. This concept will be illustrated with the quantum simulation of the helium droplet-assisted soft-landing of a single Au atom onto a titanium dioxide (TiO₂) surface. Next, it will be shown how the application of first-principles methods have disclosed the fundamental reasons why subnanometric Cu₅ clusters are resistant to irreversible oxidation, and capable of increasing and extending into the visible region the solar absorption of TiO₂, of augmenting its efficiency for photo-catalysis beyond a factor of four, also considering the decomposition and photo-activation of CO₂ as a prototypical (photo-) catalytic reaction. Finally, I will discuss how the modification of the same material with subnanometric Ag₅ clusters has converted it into a "reporter" of a surface polaron property as well as a novel two-dimensional polaronic material.

Open-Shell Variant of the London Dispersion-Corrected Hartree-Fock Method (HFLD) for the Quantification and Analysis of Noncovalent Interaction Energies

The London dispersion (LD)-corrected Hartree–Fock (HF) method (HFLD) is an ab initio approach for the quantification and analysis of noncovalent interactions (NCIs) in large systems that is based on the domain-based local pair natural orbital coupled-cluster (DLPNO-CC) theory. In the original HFLD paper, we discussed the implementation, accuracy, and efficiency of its closed-shell variant. Herein, an extension of this method to open-shell molecular systems is presented. Its accuracy is tested on challenging benchmark sets for NCIs, using CCSD(T) energies at the estimated complete basis set limit as reference. The HFLD scheme was found to be as accurate as the best-performing dispersion-corrected exchange-correlation functionals, while being nonempirical and equally efficient. In addition, it can be combined with the well-established local energy decomposition (LED) for the analysis of NCIs, thus yielding additional physical insights.

Comprehensive Basis-Set Testing of Extended Symmetry-Adapted Perturbation Theory and Assessment of Mixed-Basis Combinations to Reduce Cost

Hybrid or "extended" symmetry-adapted perturbation theory (XSAPT) replaces traditional SAPT's treatment of dispersion with better performing alternatives while at the same time extending two-body (dimer) SAPT to a many-body treatment of polarization using a self-consistent charge embedding procedure. The present work presents a systematic study of how XSAPT interaction energies and energy components converge with respect to the choice of Gaussian basis set. Errors can be reduced in a systematic way using correlation-consistent basis sets, with aug-cc-pVTZ results converged within <0.1 kcal/mol. Similar (if slightly less systematic) behavior is obtained using Karlsruhe basis sets at much lower cost, and we introduce new versions with limited augmentation that are even more efficient. Pople-style basis sets, which are more efficient still, often afford good results if a large number of polarization functions are included. The dispersion models used in XSAPT afford much faster basis-set convergence as compared to the perturbative description of dispersion in conventional SAPT, meaning that "compromise" basis sets (such as jun-cc-pVDZ) are no longer required and benchmark-quality results can be obtained using triple-ζ basis sets. The use of diffuse functions proves to be essential, especially for the description of hydrogen bonds. The "δ(Hartree-Fock)" correction for high-order induction can be performed in double-ζ basis sets without significant loss of accuracy, leading to a mixed-basis approach that offers 4× speedup over the existing (cubic scaling) XSAPT approach.

Mini Review: Quantum Confinement of Atomic and Molecular Clusters in Carbon Nanotubes

We overview our recent developments on a computational approach addressing quantum confinement of light atomic and molecular clusters (made of atomic helium and molecular hydrogen) in carbon nanotubes. We outline a multi-scale first-principles approach, based on density functional theory (DFT)-based symmetry-adapted perturbation theory, allowing an accurate characterization of the dispersion-dominated particle–nanotube interaction. Next, we describe a wave-function-based method, allowing rigorous fully coupled quantum calculations of the pseudo-nuclear bound states. The approach is illustrated by showing the transition from molecular aggregation to quasi-one-dimensional condensed matter systems of molecular deuterium and hydrogen as well as atomic ⁴He, as case studies. Finally, we present a perspective on future-oriented mixed approaches combining, e.g., orbital-free helium density functional theory (He-DFT), machine-learning parameterizations, with wave-function-based descriptions.

Symmetry-Adapted Perturbation Theory Based on Multiconfigurational Wave Function Description of Monomers

We present a formulation of the multiconfigurational (MC) wave function symmetry-adapted perturbation theory (SAPT). The method is applicable to noncovalent interactions between monomers which require a multiconfigurational description, in particular when the interacting system is strongly correlated or in an electronically excited state. SAPT(MC) is based on one- and two-particle reduced density matrices of the monomers and assumes the single-exchange approximation for the exchange energy contributions. Second-order terms are expressed through response properties from extended random phase approximation (ERPA). The dispersion components of SAPT(MC) have been introduced in our previous works [HapkaM.J. Chem. Theory Comput.2019, 15, 1016−102730525591; HapkaM.J. Chem. Theory Comput.2019, 15, 6712–672331670950]. SAPT(MC) is applied either with generalized valence bond perfect pairing (GVB) or with complete active space self-consistent field (CASSCF) treatment of the monomers. We discuss two model multireference systems: the H₂ ··· H₂ dimer in out-of-equilibrium geometries and interaction between the argon atom and excited state of ethylene. Using the C₂H₄* ··· Ar complex as an example, we examine second-order terms arising from negative transitions in the linear response function of an excited monomer. We demonstrate that the negative-transition terms must be accounted for to ensure qualitative prediction of induction and dispersion energies and develop a procedure allowing for their computation. Factors limiting the accuracy of SAPT(MC) are discussed in comparison with other second-order SAPT schemes on a data set of small single-reference dimers.

Towards Quantum-Chemical Modeling of the Activity of Anesthetic Compounds

The modeling of the activity of anesthetics is a real challenge because of their unique electronic and structural characteristics. Microscopic approaches relevant to the typical features of these systems have been developed based on the advancements in the theory of intermolecular interactions. By stressing the quantum chemical point of view, here, we review the advances in the field highlighting differences and similarities among the chemicals within this group. The binding of the anesthetics to their partners has been analyzed by Symmetry-Adapted Perturbation Theory to provide insight into the nature of the interaction and the modeling of the adducts/complexes allows us to rationalize their anesthetic properties. A new approach in the frame of microtubule concept and the importance of lipid rafts and channels in membranes is also discussed.

A Scalable Graph Neural Network Method for Developing an Accurate Force Field of Large Flexible Organic Molecules

An accurate force field is the key to the success of all molecular mechanics simulations on organic polymers and biomolecules. Accurate correlated wave function (CW) methods scale poorly with system size, so this poses a great challenge to the development of an extendible ab initio force field for large flexible organic molecules at the CW level of accuracy. In this work, we combine the physics-driven nonbonding potential with a data-driven subgraph neural network bonding model (named sGNN). Tests on polyethylene glycol, polyethene, and their block polymers show that our strategy is highly accurate and robust for molecules of different sizes and chemical compositions. Therefore, one can develop a parameter library of small molecular fragments (with sizes easily accessible to CW methods) and assemble them to predict the energy of large polymers, thus opening a new path to next-generation organic force fields.

Simplified tuning of long-range corrected density functionals for use in symmetry-adapted perturbation theory

Long considered a failure, second-order symmetry-adapted perturbation theory (SAPT) based on Kohn-Sham orbitals, or SAPT0(KS), can be resurrected for semiquantitative purposes using long-range corrected density functionals whose asymptotic behavior is adjusted separately for each monomer. As in other contexts, correct asymptotic behavior can be enforced via "optimal tuning" based on the ionization energy theorem of density functional theory, but the tuning procedure is tedious, expensive for large systems, and comes with a troubling dependence on system size. Here, we show that essentially identical results are obtained using a fast, convenient, and automated tuning procedure based on the size of the exchange hole. In conjunction with "extended" (X)SAPT methods that improve the description of dispersion, this procedure achieves benchmark-quality interaction energies, along with the usual SAPT energy decomposition, without the hassle of system-specific tuning.

Optimized damping parameters for empirical dispersion corrections to symmetry-adapted perturbation theory

Symmetry-adapted perturbation theory (SAPT) has become an invaluable tool for studying the fundamental nature of non-covalent interactions by directly computing the electrostatics, exchange (steric) repulsion, induction (polarization), and London dispersion contributions to the interaction energy using quantum mechanics. Further application of SAPT is primarily limited by its computational expense, where even its most affordable variant (SAPT0) scales as the fifth power of system size [O(N⁵)] due to the dispersion terms. The algorithmic scaling of SAPT0 is reduced from O(N⁵)→O(N⁴) by replacing these terms with the empirical D3 dispersion correction of Grimme and co-workers, forming a method that may be termed SAPT0-D3. Here, we optimize the damping parameters for the -D3 terms in SAPT0-D3 using a much larger training set than has previously been considered, namely, 8299 interaction energies computed at the complete-basis-set limit of coupled cluster through perturbative triples [CCSD(T)/CBS]. Perhaps surprisingly, with only three fitted parameters, SAPT0-D3 improves on the accuracy of SAPT0, reducing mean absolute errors from 0.61 to 0.49 kcal mol^-1 over the full set of complexes. Additionally, SAPT0-D3 exhibits a nearly 2.5× speedup over conventional SAPT0 for systems with ∼300 atoms and is applied here to systems with up to 459 atoms. Finally, we have also implemented a functional group partitioning of the approach (F-SAPT0-D3) and applied it to determine important contacts in the binding of salbutamol to G-protein coupled β₁-adrenergic receptor in both active and inactive forms. SAPT0-D3 capabilities have been added to the open-source Psi4 software.

Efficient Density-Fitted Explicitly Correlated Dispersion and Exchange Dispersion Energies

The leading-order dispersion and exchange-dispersion terms in symmetry-adapted perturbation theory (SAPT), E_disp⁽²⁰⁾ and E_exch-disp⁽²⁰⁾, suffer from slow convergence to the complete basis set limit. To alleviate this problem, explicitly correlated variants of these corrections, E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12, have been proposed recently. However, the original formalism (M., Kodrycka , J. Chem. Theory Comput. 2019, 15, 5965-5986), while highly successful in terms of improving convergence, was not competitive to conventional orbital-based SAPT in terms of computational efficiency due to the need to manipulate several kinds of two-electron integrals. In this work, we eliminate this need by decomposing all types of two-electron integrals using robust density fitting. We demonstrate that the error of the density fitting approximation is negligible when standard auxiliary bases such as aug-cc-pVXZ/MP2FIT are employed. The new implementation allowed us to study all complexes in the A24 database in basis sets up to aug-cc-pV5Z, and the E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 values exhibit vastly improved basis set convergence over their conventional counterparts. The well-converged E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 numbers can be substituted for conventional E_disp⁽²⁰⁾ and E_exch-disp⁽²⁰⁾ ones in a calculation of the total SAPT interaction energy at any level (SAPT0, SAPT2+3, ...). We show that the addition of F12 terms does not improve the accuracy of low-level SAPT treatments. However, when the theory errors are minimized in high-level SAPT approaches such as SAPT2+3(CCD)δMP2, the reduction of basis set incompleteness errors thanks to the F12 treatment substantially improves the accuracy of small-basis calculations.

A non-empirical intermolecular force-field for trinitrobenzene and its application in crystal structure prediction

An anisotropic atom-atom distributed intermolecular force-field (DIFF) for rigid trinitrobenzene (TNB) is developed using distributed multipole moments, dipolar polarizabilities, and dispersion coefficients derived from the charge density of the isolated molecule. The short-range parameters of the force-field are fitted to first- and second-order symmetry-adapted perturbation theory dimer interaction energy calculations using the distributed density-overlap model to guide the parameterization of the short-range anisotropy. The second-order calculations are used for fitting the damping coefficients of the long-range dispersion and polarization and also for relaxing the isotropic short-range coefficients in the final model, DIFF-srL2(rel). We assess the accuracy of the unrelaxed model, DIFF-srL2(norel), and its equivalent without short-range anisotropy, DIFF-srL0(norel), as these models are easier to derive. The model potentials are contrasted with empirical models for the repulsion-dispersion fitted to organic crystal structures with multipoles of iterated stockholder atoms (ISAs), FIT(ISA,L4), and with Gaussian Distributed Analysis (GDMA) multipoles, FIT(GDMA,L4), commonly used in modeling organic crystals. The potentials are tested for their ability to model the solid state of TNB. The non-empirical models provide more reasonable relative lattice energies of the three polymorphs of TNB and propose more sensible hypothetical structures than the empirical force-field (FIT). The DIFF-srL2(rel) model successfully has the most stable structure as one of the many structures that match the coordination sphere of form III. The neglect of the conformational flexibility of the nitro-groups is a significant approximation. This methodology provides a step toward force-fields capable of representing all phases of a molecule in molecular dynamics simulations.

Ab initio effective one-electron potential operators: Applications for charge-transfer energy in effective fragment potentials

The concept of effective one-electron potentials (EOPs) has proven to be extremely useful in efficient description of electronic structure of chemical systems, especially extended molecular aggregates such as interacting molecules in condensed phases. Here, a general method for EOP-based elimination of electron repulsion integrals is presented, that is tuned toward the fragment-based calculation methodologies such as the second generation of the effective fragment potentials (EFP2) method. Two general types of the EOP operator matrix elements are distinguished and treated either via the distributed multipole expansion or the extended density fitting (DF) schemes developed in this work. The EOP technique is then applied to reduce the high computational costs of the effective fragment charge-transfer (CT) terms being the bottleneck of EFP2 potentials. The alternative EOP-based CT energy model is proposed, derived within the framework of intermolecular perturbation theory with Hartree-Fock noninteracting reference wavefunctions, compatible with the original EFP2 formulation. It is found that the computational cost of the EFP2 total interaction energy calculation can be reduced by up to 38 times when using the EOP-based formulation of CT energy, as compared to the original EFP2 scheme, without compromising the accuracy for a wide range of weakly interacting neutral and ionic molecular fragments. The proposed model can thus be used routinely within the EFP2 framework.

Extension of an Atom-Atom Dispersion Function to Halogen Bonds and Its Use for Rational Design of Drugs and Biocatalysts

A dispersion function Das in the form of a damped atom-atom asymptotic expansion fitted to ab initio dispersion energies from symmetry-adapted perturbation theory was improved and extended to systems containing heavier halogen atoms. To illustrate its performance, the revised Das function was implemented in the multipole first-order electrostatic and second-order dispersion (MED) scoring model. The extension has allowed applications to a much larger set of biocomplexes than it was possible with the original Das. A reasonable correlation between MED and experimentally determined inhibitory activities was achieved in a number of test cases, including structures featuring nonphysically shortened intermonomer distances, which constitute a particular challenge for binding strength predictions. Since the MED model is also computationally efficient, it can be used for reliable and rapid assessment of the ligand affinity or multidimensional scanning of amino acid side-chain conformations in the process of rational design of novel drugs or biocatalysts.

Nonapproximated third-order exchange induction energy in symmetry-adapted perturbation theory

The exchange terms in symmetry-adapted perturbation theory (SAPT) are normally calculated within the so-called S² or single exchange approximation, which approximates the all-electron antisymmetrizer by interchanges of at most one electron pair between the interacting molecules. This approximation is typically very accurate at the van der Waals minimum separation and at larger intermolecular distances but begins to deteriorate at short range. Nonapproximated expressions for the second-order SAPT exchange corrections have been derived some time ago by Schäffer and Jansen [Mol. Phys. 111, 2570 (2013)]. In this work, we extend Schäffer and Jansen's formalism to derive and implement a nonapproximated expression for the third-order exchange-induction correction. Numerical tests on several representative noncovalent databases show that the S² approximation underestimates the exchange-induction contributions in both second and third orders. This underestimation is very similar in relative terms, but the larger absolute values of the third-order exchange-induction effects, and their near complete cancellation with the corresponding induction energies, make the third-order errors more severe. In the worst-case scenario of interactions involving ions, the breakdown of the S² approximation can result in a qualitatively wrong, attractive character of SAPT total energies at short range {as first observed by Lao and Herbert [J. Phys. Chem. A 116, 3042 (2012)]}. As expected, the inclusion of the full third-order exchange-induction energy in place of its S²-approximated counterpart restores the correct, repulsive short-range behavior of the SAPT potential energy curves computed through the third order.

A nuclear spin and spatial symmetry-adapted full quantum method for light particles inside carbon nanotubes: clusters of 3He, 4He, and para-H2

A new nuclear spin and spatial symmetry-adapted full quantum method for light fermionic and bosonic particles under cylindrical carbon nanotube confinement.

Many-Body Dispersion

A broad range of approaches to many-body dispersion are discussed, including empirical approaches with multiple fitted parameters, augmented density functional-based approaches, symmetry adapted perturbation theory, and a supermolecule approach based on coupled cluster theory. Differing definitions of "body" are considered, specifically atom-based vs molecule-based approaches.

The Nature of Strong Chalcogen Bonds Involving Chalcogen-Containing Heterocycles

Chalcogen bonds are σ hole interactions and have been used in recent years as an alternative to hydrogen bonds. In general, the electrostatic potential at the chalcogen atom and orbital delocalization effects are made responsible for the orientation of the chalcogen bond. Here, we were able to show by means of SAPT calculations that neither the induction (orbital delocalization effects) nor the electrostatic term is causing the spatial orientation of strong chalcogen bonds in tellurium-containing aromatics. Instead, steric interactions (Pauli repulsion) are responsible for the orientation. Against chemical intuition the dispersion energies of the examined tellurium-containing aromatics are far less important for the net attractive forces compared to the energies in the corresponding sulfur and selenium compounds. Our results underline the importance of often overlooked steric interactions (Pauli repulsion) in conformational control of σ hole interactions.

Parametrizing the Spatial Dependence of 1H NMR Chemical Shifts in π-Stacked Molecular Fragments

Most recently a renewed interest in several areas has arisen in factors governing the 1H NMR chemical shift (1H CS) of protons in aromatic systems. Therefore, it is important to describe how 1H CS values are affected by π-stacking intermolecular interactions. The parametrization of radial and angular dependences of the 1H CS is proposed, which is based on conventional gauge-independent atomic orbital (GIAO) calculations of explicit molecular fragments. Such a parametrization is exemplified for a benzene dimer with intermonomer vertical and horizontal distances which are in the range of values often found in crystals of organic compounds. Results obtained by the GIAO calculations combined with B3LYP and MP2 methods were compared, and revealed qualitatively the same trends in the 1H CS data. The parametrization was found to be quantitatively correct for the T-shaped benzene dimers, and its limitations were discussed. Parametrized 1H CS surfaces should become useful for providing additional restraints in the search of site-specific information through an analysis of structurally induced 1H CS changes.

How Much Dispersion Energy Is Included in the Multiconfigurational Interaction Energy?

We demonstrate how to quantify the amount of dispersion interaction recovered by supermolecular calculations with the multiconfigurational self-consistent field (MCSCF) wave functions. For this purpose, we present a rigorous derivation which connects the portion of dispersion interaction captured by the assumed wave function model-the residual dispersion interaction-with the size of the active space. Based on the obtained expression for the residual dispersion contribution, we propose a dispersion correction for the MCSCF that avoids correlation double counting. Numerical demonstration for model four-electron dimers in both ground and excited states described with the complete active space self-consistent field (CASSCF) reference serves as a proof-of-concept for the method. Accurate results, largely independent of the size of the active space, are obtained. For many-electron systems, routine CASSCF interaction energy calculations recover a tiny fraction of the full second-order dispersion energy. We found that the residual dispersion is non-negligible only for purely dispersion-bound complexes.

Description of Chiral Complexes within Functional-Group Symmetry-Adapted Perturbation Theory-The Case of (S/R)-Carvone with Derivatives of (-)-Menthol

Symmetry-adapted perturbation theory (SAPT) and functional-group SAPT (F-SAPT) are applied to examine differences in interaction energies of diastereoisomeric complexes of two chiral molecules of natural origin: (S/R)-carvone with (-)-menthol. The study is extended by including derivatives of menthol with its hydroxy group exchanged by another functional group, thus examining the substituent effect of the interaction and the interaction differences between diastereoisomers. The partitioning of the interaction energy into functional-group components allows one to explain this phenomenon by the mutual cancellation of attractive and repulsive interactions between functional groups. In some cases, one can identify dominant chiral interactions between groups of atoms of carvone and menthol derivatives, while in many other instances, no major interaction can be distinguished and the net chiral difference results from subtle near cancellation of several smaller terms. Our results indicate that the F-SAPT method can be faithfully utilized for such analyses.

From Molecular Aggregation to a One-Dimensional Quantum Crystal of Deuterium Inside a Carbon Nanotube of 1 nm Diameter

The quantum motion of clusters of up to four deuterium molecules under confinement in a single-wall (1 nm diameter) carbon nanotube is investigated by applying a highly accurate full quantum treatment of the most relevant nuclear degrees of freedom and an ab initio-derived potential model of the underlying dispersion-dominated intermolecular interactions. The wave functions and energies are calculated using an ad hoc-developed discrete variable representation (DVR) numerical approach in internal coordinates, with the space grid approaching a few billion grid points. We unambiguously demonstrate the formation of a solid-like pyramidal one-dimensional chain structure of molecules under the cylindrical nanotube confinement. The onset of solid-like packing is explained by analyzing the potential minima landscape. The stabilization of collective rotational motion through "rigid rotations" of four deuterium molecules provides conclusive evidence for the onset of a quantum solid-like behavior resembling that of quantum rings featuring persistent current (charged particles) or persistent flow (neutral particles).