Mapping spin contamination-free potential energy surfaces using restricted open-shell methods with Grassmannians.

The Lagrange-based Grassmann interpolation (G-Int) method has been extended for open-shell systems using restricted open-shell (RO) methods. The performance of this method was assessed in constructing potential energy surfaces (PESs) for vanadium(II) oxide, benzyl radical, and methanesulfenyl chloride radical cation. The density matrices generated by G-Int when used as initial guesses for self-consistent field (SCF) calculations, exhibit superior performance compared to other traditional SCF initial guess schemes, such as SADMO, GWH, and CORE. Additionally, the energy obtained from the G-Int scheme satisfies the variational principle and outperforms the direct energy-based Lagrange interpolation approach. In the case of methanesulfenyl chloride radical cation, a unique example with a flat PES at the end region along the H-C-S-Cl dihedral angle, the use of an equally-spaced grid sampling leads to significant oscillations near the end of the interval due to the effects of Runge's phenomenon. Introducing an unequally-spaced grid sampling based on a scaled Gauss-Chebyshev quadrature effectively mitigated the Runge's phenomenon, making it suitable for combining with G-Int in constructing PESs for general applications. Thus, G-Int provides an efficient and robust strategy for building spin contamination-free PESs with consistent accuracy.

Analysis of two overlapping fragmentation approaches in density matrix construction: GMBE-DM vs. ADMA.

In this study, we conduct a comparative analysis of two density matrix construction methods: the generalized many-body expansion for building density matrices (GMBE-DM) based on the set-theoretical principle of inclusion/exclusion and the adjustable density matrix assembler (ADMA) based on the Mulliken-Mezey ansatz. We apply these methods to various noncovalent clusters, including water clusters, ion-water clusters, and ion-pair clusters, using both small 6-31G(d) and large def2-TZVPPD basis sets. Our findings reveal that the GMBE-DM method, particularly when combined with the purification scheme and truncation at the one-body level [GMBE(1)-DM-P], exhibits superior performance across all test systems and basis sets. In contrast, all ADMA set of methods show reasonable results only with small and compact basis sets. For example, GMBE(1)-DM-P outperforms the best ADMA method by at least 4 and 16 times with small and large basis sets, respectively, in the case of (H2O)N=6-55. This highlights the significance of the basis set choice for ADMA, which is even more critical than the fragmentation scheme, such as the size of subsystems, while GMBE-DM consistently produces accurate results irrespective of the chosen basis set. Consequently, the efficient and robust GMBE(1)-DM-P approach is recommended as a fragmentation method for generating accurate absolute and relative energies across different binding patterns and basis sets for noncovalent clusters.

Photodissociation Dynamics of the Highly Stable ortho-Nitroaniline Cation.

0rtho-Nitroaniline (ONA) is a model for the insensitive high explosive 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) that shares strong hydrogen bonding character between adjacent nitro and amino groups. This work reports femtosecond time-resolved mass spectrometry (FTRMS) measurements and theoretical calculations that explain the high stability of the ONA cation compared with related nitroaromatic molecules. Ab initio calculations found that the lowest-lying electronic excited state of the ONA cation, D1, lies more than 2 eV above the ground state, and the energetic barriers to rearrangement and dissociation reactions exceed this D1 energy. These theoretical results were confirmed by FTRMS pump-probe measurements showing that (1) fragment ions represented less than 30% of the total ion yield when a 1014 W cm-2, 1300 nm, 20 fs pump pulse was used to ionize ONA; and (2) 3.1 eV (400 nm) photons were required to induce dissociation of the ONA cation. Stronger coupling between the ground D0 and excited D4 states of the ONA cation at the geometry of neutral ONA resulted in a transient enhancement of fragment ion yields at <300 fs pump-probe delay times, prior to relaxation of the ONA cation to its optimal geometry.

Improving second-order Møller-Plesset perturbation theory for noncovalent interactions with the machine learning-corrected ab initio dispersion potential.

In this work, we utilize our recently developed machine learning (ML)-corrected ab initio dispersion (aiD) potential, known as D3-ML, which is based on the comprehensive SAPT10K dataset and relies solely on Cartesian coordinates as input, to address the dispersion deficiencies in second-order Møller-Plesset perturbation theory (MP2) by replacing its problematic dispersion and exchange-dispersion terms with D3-ML. This leads to the development of a new dispersion-corrected MP2 method, MP2+aiD(CCD), which outperforms other spin-component-scaled and dispersion-corrected MP2 methods as well as popular ML models for predicting noncovalent interactions across various datasets, including S66 × 8, NAP6 (containing 6 naphthalene dimers), L7, S12L, DNA-ellipticine, the C60 dimer, and C60[6]CPPA. In addition, MP2+aiD(CCD) exhibits comparable or even superior performance compared to the contemporary ωB97M-V functional. The limited performance of pure ML models for systems outside the training set or larger than those in the training set highlights their instability and unpredictability. Conversely, the outstanding performance and transferability of the hybrid MP2+aiD(CCD) method can be attributed to the fusion of the physical electronic structure method and a data-driven ML model, combining the strengths of both sides. This investigation firmly establishes MP2+aiD(CCD) as one of the most accurate and reliable fifth-order scaling correlated wave function methods currently available for modeling noncovalent interactions, even for large complexes. MP2+aiD(CCD) is expected to be reliably applicable in investigating real-life complexes at the hundred-atom scale.

Ab initio dispersion potentials based on physics-based functional forms with machine learning.

In this study, we introduce SAPT10K, a comprehensive dataset comprising 9982 noncovalent interaction energies and their binding energy components (electrostatics, exchange, induction, and dispersion) for diverse intermolecular complexes of 944 unique dimers. These complexes cover significant portions of the intermolecular potential energy surface and were computed using higher-order symmetry-adapted perturbation theory, SAPT2+(3)(CCD), with a large aug-cc-pVTZ basis set. The dispersion energy values in SAPT10K serve as crucial inputs for refining the ab initio dispersion potentials based on Grimme's D3 and many-body dispersion (MBD) models. Additionally, Δ machine learning (ML) models based on newly developed intermolecular features, which are derived from intermolecular histograms of distances for element/substructure pairs to simultaneously account for local environments as well as long-range correlations, are also developed to address deficiencies of the D3/MBD models, including the inflexibility of their functional forms, the absence of MBD contributions in D3, and the standard Hirshfeld partitioning scheme used in MBD. The developed dispersion models can be applied to complexes involving a wide range of elements and charged monomers, surpassing other popular ML models, which are limited to systems with only neutral monomers and specific elements. The efficient D3-ML model, with Cartesian coordinates as the sole input, demonstrates promising results on a testing set comprising 6714 dimers, outperforming another popular ML model, component-based machine-learned intermolecular force field (CLIFF), by 1.5 times. These refined D3/MBD-ML models have the capability to replace the time-consuming dispersion components in symmetry-adapted perturbation theory-based calculations and can promptly illustrate the dispersion contribution in noncovalent complexes for supramolecular assembly and chemical reactions.

Generating accurate density matrices on the tangent space of a Grassmann manifold.

Interpolating a density matrix from a set of known density matrices is not a trivial task. This is because a linear combination of density matrices does not necessarily correspond to another density matrix. In this Communication, density matrices are examined as objects of a Grassmann manifold. Although this manifold is not a vector space, its tangent space is a vector space. As a result, one can map the density matrices on this manifold to their corresponding vectors in the tangent space and then perform interpolations on that tangent space. The resulting interpolated vector can be mapped back to the Grassmann manifold, which can then be utilized (1) as an optimal initial guess for a self-consistent field (SCF) calculation or (2) to derive energy directly without time-consuming SCF iterations. Such a promising approach is denoted as Grassmann interpolation (G-Int). The hydrogen molecule has been used to illustrate that the described interpolated method in this work preserves the essential attributes of a density matrix. For phosphorus mononitride and ferrocene, it was demonstrated numerically that reference points for the definition of the corresponding tangent spaces can be chosen arbitrarily. In addition, the interpolated density matrices provide a superior and essentially converged initial guess for an SCF calculation to make the SCF procedure itself unnecessary. Finally, this accurate, efficient, robust, and systematically improved G-Int strategy has been used for the first time to generate highly accurate potential energy surfaces with fine details for the difficult case, ferrocene.

The Grassmann interpolation method for spin-unrestricted open-shell systems.

The recently reported Grassmann interpolation (G-Int) method [J. A. Tan and K. U. Lao, J. Chem. Phys. 158, 051101 (2023)] has been extended to spin-unrestricted open-shell systems. In contrast to closed-shell systems, where G-Int has to be performed only once since the α and ß density matrices are the same, spin-unrestricted open-shell systems require G-Int to be performed twice-one for the α spin and another for the ß spin density matrix. In this work, we tested the performance of G-Int to the carbon monoxide radical cation CO●+ and nickelocene complex, which have the doublet and triple ground states, respectively. We found that the Frobenius norm errors associated with the interpolations for the α and ß spin density matrices are comparable for a given molecular geometry. These G-Int density matrices, when used as an initial guess for a self-consistent field (SCF) calculation, outperform the conventional SCF guess schemes, such as the superposition of atomic densities, purified superposition of atomic densities, core Hamiltonian, and generalized Wolfsberg-Helmholtz approximation. Depending on the desired accuracy, these G-Int density matrices can be used to directly evaluate the SCF energy without performing SCF iterations. In addition, the spin-unrestricted G-Int density matrices have been used for the first time to directly calculate the atomic charges using the Mulliken and ChElPG population analysis.

Electronic structure theory on modeling short-range noncovalent interactions between amino acids.

While short-range noncovalent interactions (NCIs) are proving to be of importance in many chemical and biological systems, these atypical bindings happen within the so-called van der Waals envelope and pose an enormous challenge for current computational methods. We introduce SNCIAA, a database of 723 benchmark interaction energies of short-range noncovalent interactions between neutral/charged amino acids originated from protein x-ray crystal structures at the "gold standard" coupled-cluster with singles, doubles, and perturbative triples/complete basis set [CCSD(T)/CBS] level of theory with a mean absolute binding uncertainty less than 0.1 kcal/mol. Subsequently, a systematic assessment of commonly used computational methods, such as the second-order Møller-Plesset theory (MP2), density functional theory (DFT), symmetry-adapted perturbation theory (SAPT), composite electronic-structure methods, semiempirical approaches, and the physical-based potentials with machine learning (IPML) on SNCIAA is carried out. It is shown that the inclusion of dispersion corrections is essential even though these dimers are dominated by electrostatics, such as hydrogen bondings and salt bridges. Overall, MP2, ωB97M-V, and B3LYP+D4 turned out to be the most reliable methods for the description of short-range NCIs even in strongly attractive/repulsive complexes. SAPT is also recommended in describing short-range NCIs only if the δMP2 correction has been included. The good performance of IPML for dimers at close-equilibrium and long-range conditions is not transferable to the short-range. We expect that SNCIAA will assist the development/improvement/validation of computational methods, such as DFT, force-fields, and ML models, in describing NCIs across entire potential energy surfaces (short-, intermediate-, and long-range NCIs) on the same footing.

Accurate prediction of global-density-dependent range-separation parameters based on machine learning.

In this work, we develop an accurate and efficient XGBoost machine learning model for predicting the global-density-dependent range-separation parameter, ωGDD, for long-range corrected functional (LRC)-ωPBE. This ωGDDML model has been built using a wide range of systems (11 466 complexes, ten different elements, and up to 139 heavy atoms) with fingerprints for the local atomic environment and histograms of distances for the long-range atomic correlation for mapping the quantum mechanical range-separation values. The promising performance on the testing set with 7046 complexes shows a mean absolute error of 0.001 117 a0-1 and only five systems (0.07%) with an absolute error larger than 0.01 a0-1, which indicates the good transferability of our ωGDDML model. In addition, the only required input to obtain ωGDDML is the Cartesian coordinates without electronic structure calculations, thereby enabling rapid predictions. LRC-ωPBE(ωGDDML) is used to predict polarizabilities for a series of oligomers, where polarizabilities are sensitive to the asymptotic density decay and are crucial in a variety of applications, including the calculations of dispersion corrections and refractive index, and surpasses the performance of all other popular density functionals except for the non-tuned LRC-ωPBE. Finally, LRC-ωPBE (ωGDDML) combined with (extended) symmetry-adapted perturbation theory is used in calculating noncovalent interactions to further show that the traditional ab initio system-specific tuning procedure can be bypassed. The present study not only provides an accurate and efficient way to determine the range-separation parameter for LRC-ωPBE but also shows the synergistic benefits of fusing the power of physically inspired density functional LRC-ωPBE and the data-driven ωGDDML model.

An accurate and efficient fragmentation approach via the generalized many-body expansion for density matrices.

With relevant chemical space growing larger and larger by the day, the ability to extend computational tractability over that larger space is of paramount importance in virtually all fields of science. The solution we aim to provide here for this issue is in the form of the generalized many-body expansion for building density matrices (GMBE-DM) based on the set-theoretical derivation with overlapping fragments, through which the energy can be obtained by a single Fock build. In combination with the purification scheme and the truncation at the one-body level, the DM-based GMBE(1)-DM-P approach shows both highly accurate absolute and relative energies for medium-to-large size water clusters with about an order of magnitude better than the corresponding energy-based GMBE(1) scheme. Simultaneously, GMBE(1)-DM-P is about an order of magnitude faster than the previously proposed MBE-DM scheme [F. Ballesteros and K. U. Lao, J. Chem. Theory Comput. 18, 179 (2022)] and is even faster than a supersystem calculation without significant parallelization to rescue the fragmentation method. For even more challenging systems including ion-water and ion-pair clusters, GMBE(1)-DM-P also performs about 3 and 30 times better than the energy-based GMBE(1) approach, respectively. In addition, this work provides the first overlapping fragmentation algorithm with a robust and effective binning scheme implemented internally in a popular quantum chemistry software package. Thus, GMBE(1)-DM-P opens a new door to accurately and efficiently describe noncovalent clusters using quantum mechanics.

Predicting and Understanding Non-Covalent Interactions Using Novel Forms of Symmetry-Adapted Perturbation Theory.

Although sometimes derided as "weak" interactions, non-covalent forces play a critical role in ligand binding and crystal packing and in determining the conformational landscape of flexible molecules. Symmetry-adapted perturbation theory (SAPT) provides a framework for accurate ab initio calculation of intermolecular interactions and furnishes a natural decomposition of the interaction energy into physically meaningful components: semiclassical electrostatics (rigorously obtained from monomer charge densities), Pauli or steric repulsion, induction (including both polarization and charge transfer), and dispersion. This decomposition helps to foster deeper understanding of non-covalent interactions and can be used to construct transferable, physics-based force fields. Separability of the SAPT interaction energy also provides the flexibility to construct composite methods, a feature that we exploit to improve the description of dispersion interactions. These are challenging to describe accurately because they arise from nonlocal electron correlation effects that appear for the first time at second order in perturbation theory but are not quantitatively described at that level.As with all quantum-chemical methods, a major limitation of SAPT is nonlinear scaling of the computational cost with respect to system size. This cost can be significantly mitigated using "SAPT0(KS)", which incorporates monomer electron correlation by means of Kohn-Sham (KS) molecular orbitals from density functional theory (DFT), as well as by an "extended" theory called XSAPT, developed by the authors. XSAPT generalizes traditional dimer SAPT to many-body systems, so that a ligand-protein interaction (for example) can be separated into contributions from individual amino acids, reducing the cost of the calculation below that of even supramolecular DFT while retaining the accuracy of high-level ab initio quantum chemistry.This Account provides an overview of the SAPT0(KS) approach and the XSAPT family of methods. Several low-cost variants are described that provide accuracy approaching that of the best ab initio benchmarks yet are affordable enough to tackle ligand-protein binding and sizable host-guest complexes. These variants include SAPT+aiD, which uses ab initio atom-atom dispersion potentials ("+aiD") in place of second-order SAPT dispersion, and also SAPT+MBD, which incorporates many-body dispersion (MBD) effects that are important in the description of nanoscale materials. Applications to drug binding highlight the size-extensive nature of dispersion, which is not a weak interaction in large systems. Other applications highlight how a physics-based analysis can sometimes upend conventional wisdom regarding intermolecular forces. In particular, careful reconsideration of π-π interactions makes clear that the quadrupolar electrostatics (or "Hunter-Sanders") model of π-π stacking should be replaced by a "van der Waals model" in which conformational preferences arise from a competition between dispersion and Pauli repulsion. Our analysis also suggests that molecular shape, rather than aromaticity per se, is the key factor driving strong stacking interactions. Looking forward, we anticipate that XSAPT-based methods can play a role in screening of drug candidates and in materials design.

Buffer concentration dramatically affects the stability of S-nitrosothiols in aqueous solutions.

S-nitrosothiols (RSNOs) are an important group of nitric oxide (NO)-donating compounds with low toxicity and wide biomedical applications. In this paper, we, for the first time, demonstrate that the concentration of buffer remarkably affects the stability of RSNOs including naturally occurring S-nitrosoglutathione (GSNO) and synthetic S-nitroso-N-acetylpenicillamine (SNAP). For a solution with a high concentration of GSNO (e.g., 50 mM) and an initial near-neutral pH, the optimal buffer concentration is close to the GSNO concentration under our experimental conditions. A lower buffer concentration does not have adequate buffer capacity to resist the pH drop caused by GSNO decomposition. The decreased solution pH further accelerates GSNO decomposition because GSNO is most stable at near-neutral pH according to our density functional theory (DFT) calculations. A higher-than-optimal buffer concentration also reduces the GSNO stability because buffer ingredients including phosphate, Tris base, and HEPES consume NO/N2O3. In contrast to GSNO, the highest SNAP stability is obtained when the starting solution at a neutral pH does not contain buffer species, and the stability decreases as the buffer concentration increases. This is because SNAP is more stable at mildly acidic pH and the SNAP decomposition-induced pH drop stabilizes the donor. When the RSNO concentration is low (e.g., 1 mM), the buffer concentration also matters because any excess buffer accelerates the donor decomposition. Since the effect of buffer concentration was previously overlooked and suboptimal buffer concentrations were often used, this paper will aid in the formulation of RSNO solutions to obtain the maximum stability for prolonged storage and sustained NO release.

Coupled Cluster Benchmarking of Large Noncovalent Complexes in L7 and S12L as Well as the C60 Dimer, DNA-Ellipticine, and HIV-Indinavir.

In this work, we report the benchmark binding energies of the seven complexes within the L7 data set, six host-guest complexes from the S12L data set, a C60 dimer, the DNA-ellipticine intercalation complex, and the largest system of the study, the HIV-indinavir system, which contained 343 atoms or 139 heavy atoms. The high-quality values reported were obtained via a focal point method that relies on the canonical form of second-order Møller-Plesset theory and the domain-based local pair natural orbital scheme for the coupled cluster with single double and perturbative triple excitations [DLPNO-CCSD(T)] extrapolated to the complete basis set (CBS) limit. The results in this work not only corroborate but also improve upon some previous benchmark values for large noncovalent complexes albeit at a relatively steep cost. Although local CCSD(T) and the largely successful fixed-node diffusion Monte Carlo (FN-DMC) have been shown to generally agree for small- to medium-size systems, a discrepancy in their reported binding energy values arises for large complexes, where the magnitude of the disagreement is a definite cause for concern. For example, the largest deviation in the L7 data set was 2.8 kcal/mol (∼10%) on the low end in C3GC. Such a deviation only grows worse in the S12L set, which showed a difference of up to 10.4 kcal/mol (∼25%) by a conservative estimation in buckycatcher-C60. The DNA-ellipticine complex also generated a disagreement of 4.4 kcal/mol (∼10%) between both state-of-the-art methods. The disagreement between local CCSD(T) and FN-DMC in large noncovalent complexes shows that it is urgently needed to have the canonical CCSD(T), the Monte Carlo CCSD(T), or the full configuration interaction quantum Monte Carlo approaches available to large systems on the hundred-atom scale to solve this dilemma. In addition, the performances of cheaper popular computational methods were assessed for the studied complexes with respect to DLPNO-CCSD(T)/CBS. r2SCAN-3c, B97M-V, and PBE0+D4 work well in large noncovalent complexes in this work, and GFN2-xTB performs well in π-π stacking complexes. B97M-V is the most reliable computationally efficient approach to predicting noncovalent interactions for large complexes, being the only one to have binding errors within the so-called 1 kcal/mol "chemical accuracy". The benchmark interaction energies of these host-guest complexes, molecular materials, and biological systems with electronic and medicinal implications provide crucial reference data for the improvement of current and future lower-cost methods.

Ultrafast Dynamics of Nitro-Nitrite Rearrangement and Dissociation in Nitromethane Cation.

We report new insights into the ultrafast rearrangement and dissociation dynamics of nitromethane cation (NM+) using pump-probe measurements, electronic structure calculations, and ab initio molecular dynamics simulations. The "roaming" nitro-nitrite rearrangement (NNR) pathway involving large-amplitude atomic motion, which has been previously described for neutral nitromethane, is demonstrated for NM+. Excess energy resulting from initial population of the electronically excited D2 state of NM+ upon strong-field ionization provides the necessary energy to initiate NNR and subsequent dissociation into NO+. Both pump-probe measurements and molecular dynamics simulations are consistent with the completion of NNR within 500 fs of ionization with dissociation into NO+ and OCH3 occurring ∼30 fs later. Pump-probe measurements indicate that NO+ formation is in competition with the direct dissociation of NM+ to CH3+ and NO2. Electronic structure calculations indicate that a strong D0 → D1 transition can be excited at 650 nm when the C-N bond is stretched from its equilibrium value (1.48 Å) to 1.88 Å. On the other hand, relaxation of the NM+ cation after ionization into D0 occurs in less than 50 fs and results in observation of intact NM+. Direct dissociation of the equilibrium NM+ to produce NO2+ and CH3 can be induced with 650 nm excitation via a weakly allowed D0 → D2 transition.

Accurate molecular polarizabilities with coupled cluster theory and machine learning.

The molecular dipole polarizability describes the tendency of a molecule to change its dipole moment in response to an applied electric field. This quantity governs key intra- and intermolecular interactions, such as induction and dispersion; plays a vital role in determining the spectroscopic signatures of molecules; and is an essential ingredient in polarizable force fields. Compared with other ground-state properties, an accurate prediction of the molecular polarizability is considerably more difficult, as this response quantity is quite sensitive to the underlying electronic structure description. In this work, we present highly accurate quantum mechanical calculations of the static dipole polarizability tensors of 7,211 small organic molecules computed using linear response coupled cluster singles and doubles theory (LR-CCSD). Using a symmetry-adapted machine-learning approach, we demonstrate that it is possible to predict the LR-CCSD molecular polarizabilities of these small molecules with an error that is an order of magnitude smaller than that of hybrid density functional theory (DFT) at a negligible computational cost. The resultant model is robust and transferable, yielding molecular polarizabilities for a diverse set of 52 larger molecules (including challenging conjugated systems, carbohydrates, small drugs, amino acids, nucleobases, and hydrocarbon isomers) at an accuracy that exceeds that of hybrid DFT. The atom-centered decomposition implicit in our machine-learning approach offers some insight into the shortcomings of DFT in the prediction of this fundamental quantity of interest.

Electron confinement meet electron delocalization: non-additivity and finite-size effects in the polarizabilities and dispersion coefficients of the fullerenes.

In this work, we used finite-field derivative techniques and density functional theory (DFT) to compute the static isotropic polarizability series (αl with l = 1, 2, 3) for the C60-C84 fullerenes and quantitatively assess the intrinsic non-additivity in these fundamental response properties. By comparing against classical models of the fullerenes as conducting spherical shells (or solid spheres) of uniform electron density, a detailed critical analysis of the derived effective scaling laws (α1 ∼ N1.2, α2 ∼ N2.0, α3 ∼ N2.7) demonstrates that the electronic structure of finite-sized fullerenes-a unique dichotomy of electron confinement and delocalization effects due to their quasi-spherical cage-like structures and encapsulated void spaces-simultaneously limits and enhances their quantum mechanical response to electric field perturbations. Corresponding frequency-dependent polarizabilities were obtained by inputting the αl series into the hollow sphere model (within the modified single frequency approximation), and used to compute the molecular dispersion coefficients (Cn with n = 6, 8, 9, 10) needed to describe the non-trivial van der Waals (vdW) interactions in fullerene-based systems. Using first-order perturbation theory in conjunction with >140 000 DFT calculations, we also computed the non-negligible zero-point vibrational contributions to α1 in C60 and C70, thereby enabling a more accurate and direct comparison between theory and experiment for these quintessential nanostructures.

Coupled cluster benchmarks of large noncovalent complexes: The L7 dataset as well as DNA-ellipticine and buckycatcher-fullerene.

In this work, benchmark binding energies for dispersion-bound complexes in the L7 dataset, the DNA-ellipticine intercalation complex, and the buckycatcher-C60 complex with 120 heavy atoms using a focal-point method based on the canonical form of second-order Møller-Plesset theory (MP2) and the domain based local pair natural orbital scheme for the coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] extrapolated to the complete basis set (CBS) limit are reported. This work allows for increased confidence given the agreement with respect to values recently obtained using the local natural orbital CCSD(T) for L7 and the canonical CCSD(T)/CBS result for the coronene dimer (C2C2PD). Therefore, these results can be considered pushing the CCSD(T)/CBS binding benchmark to the hundred-atom scale. The disagreements between the two state-of-the-art methods, CCSD(T) and fixed-node diffusion Monte Carlo, are substantial with at least 2.0 (∼10%), 1.9 (∼5%), and 10.3 kcal/mol (∼25%) differences for C2C2PD in L7, DNA-ellipticine, and buckycatcher-C60, respectively. Such sizable discrepancy above "chemical accuracy" for large noncovalent complexes indicates how challenging it is to obtain benchmark binding interactions for systems beyond small molecules, although the three up-to-date density functionals, PBE0+D4, ωB97M-V, and B97M-V, agree better with CCSD(T) for these large systems. In addition to reporting these values, different basis sets and various CBS extrapolation parameters for Hartree-Fock and MP2 correlation energies were tested for the first time in large noncovalent complexes with the goal of providing some indications toward optimal cost effective routes to approach the CBS limit without substantial loss in quality.

NENCI-2021. I. A large benchmark database of non-equilibrium non-covalent interactions emphasizing close intermolecular contacts.

In this work, we present NENCI-2021, a benchmark database of ∼8000 Non-Equilibirum Non-Covalent Interaction energies for a large and diverse selection of intermolecular complexes of biological and chemical relevance. To meet the growing demand for large and high-quality quantum mechanical data in the chemical sciences, NENCI-2021 starts with the 101 molecular dimers in the widely used S66 and S101 databases and extends the scope of these works by (i) including 40 cation-π and anion-π complexes, a fundamentally important class of non-covalent interactions that are found throughout nature and pose a substantial challenge to theory, and (ii) systematically sampling all 141 intermolecular potential energy surfaces (PESs) by simultaneously varying the intermolecular distance and intermolecular angle in each dimer. Designed with an emphasis on close contacts, the complexes in NENCI-2021 were generated by sampling seven intermolecular distances along each PES (ranging from 0.7× to 1.1× the equilibrium separation) and nine intermolecular angles per distance (five for each ion-π complex), yielding an extensive database of 7763 benchmark intermolecular interaction energies (Eint) obtained at the coupled-cluster with singles, doubles, and perturbative triples/complete basis set [CCSD(T)/CBS] level of theory. The Eint values in NENCI-2021 span a total of 225.3 kcal/mol, ranging from -38.5 to +186.8 kcal/mol, with a mean (median) Eint value of -1.06 kcal/mol (-2.39 kcal/mol). In addition, a wide range of intermolecular atom-pair distances are also present in NENCI-2021, where close intermolecular contacts involving atoms that are located within the so-called van der Waals envelope are prevalent-these interactions, in particular, pose an enormous challenge for molecular modeling and are observed in many important chemical and biological systems. A detailed symmetry-adapted perturbation theory (SAPT)-based energy decomposition analysis also confirms the diverse and comprehensive nature of the intermolecular binding motifs present in NENCI-2021, which now includes a significant number of primarily induction-bound dimers (e.g., cation-π complexes). NENCI-2021 thus spans all regions of the SAPT ternary diagram, thereby warranting a new four-category classification scheme that includes complexes primarily bound by electrostatics (3499), induction (700), dispersion (1372), or mixtures thereof (2192). A critical error analysis performed on a representative set of intermolecular complexes in NENCI-2021 demonstrates that the Eint values provided herein have an average error of ±0.1 kcal/mol, even for complexes with strongly repulsive Eint values, and maximum errors of ±0.2-0.3 kcal/mol (i.e., ∼±1.0 kJ/mol) for the most challenging cases. For these reasons, we expect that NENCI-2021 will play an important role in the testing, training, and development of next-generation classical and polarizable force fields, density functional theory approximations, wavefunction theory methods, and machine learning based intra- and inter-molecular potentials.

Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package.

This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange-correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear-electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an "open teamware" model and an increasingly modular design.

Attracting Opposites: Promiscuous Ion-π Binding in the Nucleobases.

Ion-π interactions between the face of a molecular π-system and a cation or anion are among the strongest noncovalent interactions known, with applications throughout biochemistry and structural biology, molecular recognition and host-guest chemistry, as well as enzyme kinetics and organocatalysis. In this work, we examine the competing notions of selectivity and flexibility in this class of noncovalent interactions by investigating how certain π-systems can be promiscuous ion-π binders with the versatility to form favorable cation- and anion-π complexes. We focus our efforts on a detailed theoretical case study of the DNA/RNA nucleobases by first demonstrating that these π-systems are promiscuous ion-π binders with the biologically relevant Li+/Na+ cations and F-/Cl- anions via benchmark-quality quantum-mechanical binding energy curves computed at the CCSD(T)/CBS level of theory. Using a symmetry-adapted perturbation theory (SAPT)-based energy decomposition analysis, we explore the different physicochemical driving forces underlying the formation of cation- and anion-π complexes, as well as the crucial role played by charge penetration effects in determining the nontrivial (and often counterintuitive) electrostatics in anion-π systems. In doing so, a unified view of these rather distinct noncovalent binding motifs emerges with the finding that both cation- and anion-π complexes are strongly stabilized by an essentially ring-independent potential that can only be overcome by substantially unfavorable electrostatics. This work furnishes a more comprehensive explanation for decades of observed correlations between the equilibrium binding energy and the electrostatic potential above the ring and provides new insight into the nature of selectivity and flexibility in this important class of noncovalent interactions. Quite interestingly, the analysis presented herein demonstrates that π-systems have an inherent propensity to bind both cations and anions, thereby implying that promiscuous ion-π binding is not an exotic property of the nucleobases and should be common in nature.