Benchmarking Basis Sets for Density Functional Theory Thermochemistry Calculations: Why Unpolarized Basis Sets and the Polarized 6-311G Family Should Be Avoided.

Basis sets are a crucial but often largely overlooked choice in setting up quantum chemistry calculations. The choice of the basis set can be critical in determining the accuracy and calculation time of your quantum chemistry calculations. Clear recommendations based on thorough benchmarking are essential but not readily available currently. This study investigates the relative quality of basis sets for general properties by benchmarking basis set performance for a diverse set of 139 reactions (from the diet-150-GMTKN55 data set). In our analysis, we find the distributions of errors are often significantly non-Gaussian, meaning that the joint consideration of median errors, mean absolute errors, and outlier statistics is helpful to provide a holistic understanding of basis set performance. Our direct comparison of performance between most modern basis sets provides quantitative evidence for basis set recommendations that broadly align with the established understanding of basis set experts and is evident in the design of modern basis sets. For example, while zeta is a good measure of quality, it is not the only determining factor for an accurate calculation with unpolarized double- and triple-ζ basis sets (like 6-31G and 6-311G) having very poor performance. Appropriate use of polarization functions (e.g., 6-31G*) is essential to obtain the accuracy offered by double- or triple-ζ basis sets. In our study, the best performances for double- and triple-ζ basis sets are 6-31++G** and pcseg-2, respectively. However, the performances of singly polarized double-ζ and doubly polarized triple-ζ basis sets are quite similar with one key exception: the polarized 6-311G basis set family has poor parametrization, which means its performance is more like a double-ζ than a triple-ζ basis set. All versions of the 6-311G basis set family should be avoided entirely for valence chemistry calculations moving forward.

On the specialization of Gaussian basis sets for core-dependent properties.

Despite the fact that most quantum chemistry basis sets are designed for accurately modeling valence chemistry, these general-purpose basis sets continue to be widely used to model core-dependent properties. Core-specialized basis sets are designed with specific features to accurately represent the behavior of the core region. This design typically incorporates Gaussian primitives with higher exponents to capture core behavior effectively, as well as some decontraction of basis functions to provide flexibility in describing the core electronic wave function. The highest Gaussian exponent and the degree of contraction for both s- and p-basis functions effectively characterize these design aspects. In this study, we compare the design and performance of general-purpose basis sets against several literature-based basis sets specifically designed for three core-dependent properties: J coupling constants, hyperfine coupling constants, and magnetic shielding constants (used for calculating chemical shifts). Our findings consistently demonstrate a significant reduction in error when employing core-specialized basis sets, often at a marginal increase in computational cost compared to the popular 6-31G** basis set. Notably, for expedient calculations of J coupling, hyperfine coupling, and magnetic shielding constants, we recommend the use of the pcJ-1, EPR-II, and pcSseg-1 basis sets, respectively. For higher accuracy, the pcJ-2, EPR-III, and pcSseg-2 basis sets are recommended.

Model Chemistry Recommendations for Scaled Harmonic Frequency Calculations: A Benchmark Study.

Despite the widespread popularity of scaled harmonic frequency calculations to predict experimental fundamental frequencies in chemistry, sparse benchmarking is available to guide users on the appropriate level of theory and basis set choices (model chemistry) or deep understanding of expected errors. An updated assessment of the best approach for scaling to minimize errors is also overdue. Here, we assess the performance of over 600 popular, contemporary, and robust model chemistries in the calculation of scaled harmonic frequencies, evaluating different scaling factor types and their implications in the scaled harmonic frequencies and model chemistry performance. We can summarize our results into three main findings: (1) Using model-chemistry-specific scaling factors optimized for three different frequency regions (low (<1,000 cm-1), mid (1,000-2,000 cm-1), and high (>2,000 cm-1)) results in substantial improvements in the agreement between the scaled harmonic and experimental frequencies compared to other choices. (2) Larger basis sets and more robust levels of theory generally lead to superior performance; however, the particular model chemistry choice matters and poor choices lead to significantly reduced accuracies. (3) Outliers are expected in routine calculations regardless of the model chemistry choice. Our benchmarking results here do not consider the intensity of vibrational transitions; however, we draw upon previous benchmarking results for dipole moments that highlight the importance of diffuse functions (i.e., augmented basis sets) in high-quality intensity predictions. In terms of specific recommendations, overall, the highest accuracy model chemistries are double-hybrid density functional approximations with a non-Pople augmented triple-ζ basis set, which can produce median frequency errors down to 7.6 cm-1 (DSD-PBEP86/def2-TZVPD), which is very close to the error in the harmonic approximation, i.e., the anharmonicity error. Double-ζ basis sets should not be used with double-hybrid functionals as there is no improvement compared to hybrid functional results (unlike for double-hybrid triple-ζ model chemistries). Note that 6-311G* and 6-311+G* basis sets perform like a double-ζ basis set for vibrational frequencies. After scaling, all studied hybrid functionals with non-Pople triple-ζ basis sets will produce median errors of less than 15 cm-1, with the best result of 9.9 cm-1 with B97-1/def2-TZVPD. Appropriate matching of double-ζ basis sets with hybrid functionals can produce high-quality results, but the precise choice of functional and basis set is more important. The B97-1, TPSS0-D3(BJ), or ωB97X-D hybrid density functionals with 6-31G*, pc-1, or pcseg-1 are recommended for fast routine calculations, all delivering median errors of 11-12 cm-1. Note that dispersion corrections are not easily available for B97-1; given its strong performance here, we recommend these be added to major programs in coming updates.

The performance of CCSD(T) for the calculation of dipole moments in diatomics.

This work analyzes the accuracy of the coupled cluster with single, double, and perturbative triple excitation [CCSD(T)] method for predicting dipole moments. In particular, we benchmark CCSD(T) predictions for the equilibrium bond length, vibrational frequency, and dipole moment versus accurate experimental data. As a result, we find that CCSD(T) leads to accurate dipole moments. However, in some cases, it disagrees with the experimental values, and the disagreement can not be satisfactorily explained via relativistic or multi-reference effects. Therefore, our results indicate that benchmark studies for energy and geometry properties do not accurately describe other electron density magnitudes.

VIBFREQ1295: A New Database for Vibrational Frequency Calculations.

High-throughput approaches for producing approximate vibrational spectral data for molecules of astrochemistry interest rely on harmonic frequency calculations using computational quantum chemistry. However, model chemistry recommendations (i.e., a level of theory and basis set pair) for these calculations are not yet available and, thus, thorough benchmarking against comprehensive benchmark databases is needed. Here, we present a new database for vibrational frequency calculations (VIBFREQ1295) storing 1295 experimental fundamental frequencies and CCSD(T)(F12*)/cc-pVDZ-F12 ab initio harmonic frequencies from 141 molecules. VIBFREQ1295's experimental data was complied through a comprehensive review of contemporary experimental data, while the ab initio data was computed here. The chemical space spanned by the molecules chosen is considered in-depth and is shown to have good representation of common organic functional groups and vibrational modes. Scaling factors are routinely used to approximate the effect of anharmonicity and convert computed harmonic frequencies to predicted fundamental frequencies. With our experimental and high-level ab initio data, we find that a single global uniform scaling factor of 0.9617(3) results in median differences of 15.9(5) cm-1. A far superior performance with a median difference of 7.5(5) cm-1 can be obtained, however, by using separate scaling factors (SFs) for three regions: frequencies less than 1000 cm-1 (SF = 0.987(1)), between 1000 and 2000 cm-1 (SF = 0.9727(6)), and above 2000 cm-1 (SF = 0.9564(4)). This sets a lower bound for the performance that could be reliably obtained using scaling of harmonic frequency calculations to predict experimental fundamental frequencies. VIBFREQ1295's most important purpose is to provide a robust database for benchmarking the performance of any vibrational frequency calculations. VIBFREQ1295 data could also be used to train machine-learning models for the prediction of vibrational spectra and as a reference and data starting point for more detailed spectroscopic modeling of particular molecules. The database can be found as part of the Supporting Information for this paper or in the Harvard DataVerse at https://doi.org/10.7910/DVN/VLVNU7.

Photodissociation of dicarbon: How nature breaks an unusual multiple bond.

The dicarbon molecule (C2) is found in flames, comets, stars, and the diffuse interstellar medium. In comets, it is responsible for the green color of the coma, but it is not found in the tail. It has long been held to photodissociate in sunlight with a lifetime precluding observation in the tail, but the mechanism was not known. Here we directly observe photodissociation of C2 From the speed of the recoiling carbon atoms, a bond dissociation energy of 602.804(29) kJ·mol[Formula: see text] is determined, with an uncertainty comparable to its more experimentally accessible N2 and O2 counterparts. The value is within 0.03 kJ·mol-1 of high-level quantum theory. This work shows that, to break the quadruple bond of C2 using sunlight, the molecule must absorb two photons and undergo two "forbidden" transitions.

Full spectroscopic model and trihybrid experimental-perturbative-variational line list for CN.

Accurate line lists are important for the description of the spectroscopic nature of small molecules. While a line list for CN (an important molecule for chemistry and astrophysics) exists, no underlying energy spectroscopic model has been published, which is required to consider the sensitivity of transitions to a variation of the proton-to-electron mass ratio. Here we have developed a Duo energy spectroscopic model as well as a novel hybrid style line list for CN and its isotopologues, combining energy levels that are derived experimentally (Marvel), using the traditional/perturbative approach (Mollist), and the variational approach (from a Duo spectroscopic model using standard ExoMol methodology). The final Trihybrid ExoMol-style line list for 12C14N consists of 28 004 energy levels (6864 experimental, 1574 perturbative, the rest variational) and 2285 103 transitions up to 60 000 cm-1 between the three lowest electronic states (X 2Σ+, A 2Π, and B 2Σ+). The spectroscopic model created is used to evaluate CN as a molecular probe to constrain the variation of the proton-to-electron mass ratio; no overly promising sensitive transitions for extragalactic study were identified.

Experimental energy levels of 12C14N through marvel analysis.

The cyano radical (CN) is a key molecule across many different factions of astronomy and chemistry. Accurate, empirical rovibronic energy levels with uncertainties are determined for eight doublet states of CN using the marvel (Measured Active Rotational-Vibrational Energy Levels) algorithm. 40 333 transitions were validated from 22 different published sources to generate 8083 spin-rovibronic energy levels. The empirical energy levels obtained from the marvel analysis are compared to current energy levels from the mollist line list. The mollist transition frequencies are updated with marvel energy level data which brings the frequencies obtained through experimental data up to 77.3 per cent from the original 11.3 per cent, with 92.6 per cent of the transitions with intensities over 10-23 cm molecule-1 at 1000 K now known from experimental data. At 2000 K, 100.0 per cent of the partition function is recovered using only marvel energy levels, while 98.2 per cent is still recovered at 5000 K.

Computation of Dipole Moments: A Recommendation on the Choice of the Basis Set and the Level of Theory.

Together with experimental data, theoretically predicted dipole moments represent a valuable tool for different branches in the chemical and physical sciences. With the diversity of levels of theory and basis sets available, a reliable combination must be carefully chosen in order to achieve accurate predictions. In a recent publication (J. Chem. Theory Comput. 2018, 14 (4), 1969-1981), Hait and Head-Gordon took a first step in this regard by providing recommendations on the best density functionals suitable for these purposes. However, no extensive study has been performed to provide recommendations on the basis set choice. Here, we shed some light into this matter by evaluating the performance of 38 general-purpose basis sets of single- up to triple-ζ-quality, when coupled with nine different levels of theory, in the computation of dipole moments. The calculations were performed on a data set with 114 small molecules containing second- and third-row elements. We based our analysis in regularized root-mean-square errors (regularized RMSE), in which the difference between the calculated µcalc and benchmark µbmk dipole moment values is derived as (µcalc[D] - µbmk[D])/(max(µbmk[D],1[D])). This procedure ensures relative errors for ionic species and absolute errors for species with small dipole moment values. Our results indicate that the best compromise between accuracy and computational efficiency is achieved by performing the computations with an augmented double-ζ-quality basis set (i.e., aug-pc-1, aug-pcseg-1, aug-cc-pVDZ) together with a hybrid functional (e.g., ωB97X-V, SOGGA11-X). Augmented triple-ζ basis sets could enhance the accuracy of the computations, but the computational cost of introducing such a basis set is substantial compared with the small improvement provided. These findings also highlight the crucial role that augmentation of the basis set with diffuse functions on both hydrogen and non-hydrogen atoms plays in the computation of dipole moments.

General mathematical formulation of scattering processes in atom-diatomic collisions in the RmatReact methodology.

Accurately modelling cold and ultracold reactive collisions occurring over deep potential wells, such as [Formula: see text], requires the development of new theoretical and computational methodologies. One potentially useful framework is the R-matrix method adopted widely for electron-molecule collisions which has more recently been applied to non-reactive heavy-particle collisions such as Ar-Ar. The existing treatment of non-reactive elastic and inelastic scattering needs to be substantially extended to enable modelling of reactive collisions: this is the subject of this paper. Herein, we develop the general mathematical formulation for non-reactive elastic and inelastic scattering, photoassociation, photodissociation, charge exchange and reactive scattering using the R-matrix method. Of particular note is that the inner region, of central importance to calculable R-matrix methodologies, must be finite in all scattering coordinates rather than a single scattering coordinate as for non-reactive scattering. This article is part of a discussion meeting issue 'Advances in hydrogen molecular ions: H3+, H5+ and beyond'.

Low-temperature chemistry using the R-matrix method.

Techniques for producing cold and ultracold molecules are enabling the study of chemical reactions and scattering at the quantum scattering limit, with only a few partial waves contributing to the incident channel, leading to the observation and even full control of state-to-state collisions in this regime. A new R-matrix formalism is presented for tackling problems involving low- and ultra-low energy collisions. This general formalism is particularly appropriate for slow collisions occurring on potential energy surfaces with deep wells. The many resonance states make such systems hard to treat theoretically but offer the best prospects for novel physics: resonances are already being widely used to control diatomic systems and should provide the route to steering ultracold reactions. Our R-matrix-based formalism builds on the progress made in variational calculations of molecular spectra by using these methods to provide wavefunctions for the whole system at short internuclear distances, (a regime known as the inner region). These wavefunctions are used to construct collision energy-dependent R-matrices which can then be propagated to give cross sections at each collision energy. The method is formulated for ultracold collision systems with differing numbers of atoms.

Accurate Electron Densities at Nuclei Using Small Ramp-Gaussian Basis Sets.

Electron densities at nuclei are difficult to calculate accurately with all-Gaussian basis sets because they lack an electron-nuclear cusp. The newly developed mixed ramp-Gaussian basis sets, such as R-31G, possess electron-nuclear cusps due to the presence of ramp functions in the basis. The R-31G basis set is a general-purpose mixed ramp-Gaussian basis set modeled on the 6-31G basis set. The prediction of electron densities at nuclei using R-31G basis sets for Li-F outperforms Dunning, Pople, and Jensen general purpose all-Gaussian basis sets of triple-ζ quality or lower and the cc-pVQZ basis set. It is of similar quality to the specialized pcJ-0 basis set which was developed with partial decontraction of core functions and extra high exponent s-Gaussians to predict electron density at the nucleus. These results show significant advantages in the properties of mixed ramp-Gaussian basis sets compared to all-Gaussian basis sets.

Electron-vibration entanglement in the Born-Oppenheimer description of chemical reactions and spectroscopy.

Entanglement is sometimes regarded as the quintessential measure of the quantum nature of a system and its significance for the understanding of coupled electronic and vibrational motions in molecules has been conjectured. Previously, we considered the entanglement developed in a spatially localized diabatic basis representation of the electronic states, considering design rules for qubits in a low-temperature chemical quantum computer. We extend this to consider the entanglement developed during high-energy processes. We also consider the entanglement developed using adiabatic electronic basis, providing a novel way for interpreting effects of the breakdown of the Born-Oppenheimer (BO) approximation. We consider: (i) BO entanglement in the ground-state wavefunction relevant to equilibrium thermodynamics, (ii) BO entanglement associated with low-energy wavefunctions relevant to infrared and tunneling spectroscopies, (iii) BO entanglement in high-energy eigenfunctions relevant to chemical reaction processes, and (iv) BO entanglement developed during reactive wavepacket dynamics. A two-state single-mode diabatic model descriptive of a wide range of chemical phenomena is used for this purpose. The entanglement developed by BO breakdown correlates simply with the diameter of the cusp introduced by the BO approximation, and a hierarchy appears between the various BO-breakdown correction terms, with the first-derivative correction being more important than the second-derivative correction which is more important than the diagonal correction. This simplicity is in contrast to the complexity of BO-breakdown effects on thermodynamic, spectroscopic, and kinetic properties. Further, processes poorly treated at the BO level that appear adequately treated using the Born-Huang adiabatic approximation are found to have properties that can only be described using a non-adiabatic description. For the entanglement developed between diabatic electronic states and the nuclear motion, qualitatively differently behavior is found compared to traditional properties of the density matrix and hence entanglement provides new information about system properties. For chemical reactions, this type of entanglement simply builds up as the transition-state region is crossed. It is robust to small changes in parameter values and is therefore more attractive for making quantum qubits than is the related fragile ground-state entanglement, provided that coherent motion at the transition state can be sustained.

A unified diabatic description for electron transfer reactions, isomerization reactions, proton transfer reactions, and aromaticity.

While diabatic approaches are ubiquitous for the understanding of electron-transfer reactions and have been mooted as being of general relevance, alternate applications have not been able to unify the same wide range of observed spectroscopic and kinetic properties. The cause of this is identified as the fundamentally different orbital configurations involved: charge-transfer phenomena involve typically either 1 or 3 electrons in two orbitals whereas most reactions are typically closed shell. As a result, two vibrationally coupled electronic states depict charge-transfer scenarios whereas three coupled states arise for closed-shell reactions of non-degenerate molecules and seven states for the reactions implicated in the aromaticity of benzene. Previous diabatic treatments of closed-shell processes have considered only two arbitrarily chosen states as being critical, mapping these states to those for electron transfer. We show that such effective two-state diabatic models are feasible but involve renormalized electronic coupling and vibrational coupling parameters, with this renormalization being property dependent. With this caveat, diabatic models are shown to provide excellent descriptions of the spectroscopy and kinetics of the ammonia inversion reaction, proton transfer in N2H7(+), and aromaticity in benzene. This allows for the development of a single simple theory that can semi-quantitatively describe all of these chemical phenomena, as well as of course electron-transfer reactions. It forms a basis for understanding many technologically relevant aspects of chemical reactions, condensed-matter physics, chemical quantum entanglement, nanotechnology, and natural or artificial solar energy capture and conversion.

Non-adiabatic effects in thermochemistry, spectroscopy and kinetics: the general importance of all three Born-Oppenheimer breakdown corrections.

Using a simple model Hamiltonian, the three correction terms for Born-Oppenheimer (BO) breakdown, the adiabatic diagonal correction (DC), the first-derivative momentum non-adiabatic correction (FD), and the second-derivative kinetic-energy non-adiabatic correction (SD), are shown to all contribute to thermodynamic and spectroscopic properties as well as to thermal non-diabatic chemical reaction rates. While DC often accounts for >80% of thermodynamic and spectroscopic property changes, the commonly used practice of including only the FD correction in kinetics calculations is rarely found to be adequate. For electron-transfer reactions not in the inverted region, the common physical picture that diabatic processes occur because of surface hopping at the transition state is proven inadequate as the DC acts first to block access, increasing the transition state energy by (ℏω)(2)λ/16J(2) (where λ is the reorganization energy, J the electronic coupling and ω the vibration frequency). However, the rate constant in the weakly-coupled Golden-Rule limit is identified as being only inversely proportional to this change rather than exponentially damped, owing to the effects of tunneling and surface hopping. Such weakly-coupled long-range electron-transfer processes should therefore not be described as "non-adiabatic" processes as they are easily described by Born-Huang ground-state adiabatic surfaces made by adding the DC to the BO surfaces; instead, they should be called just "non-Born-Oppenheimer" processes. The model system studied consists of two diabatic harmonic potential-energy surfaces coupled linearly through a single vibration, the "two-site Holstein model". Analytical expressions are derived for the BO breakdown terms, and the model is solved over a large parameter space focusing on both the lowest-energy spectroscopic transitions and the quantum dynamics of coherent-state wavepackets. BO breakdown is investigated pertinent to: ammonia inversion, aromaticity in benzene, the Creutz-Taube ion, the bacterial photosynthetic reaction centre, BNB, the molecular conductor Alq3, and inverted-region charge recombination in a ferrocene-porphyrin-fullerene triad photosynthetic model compound. Throughout, the fundamental nature of BO breakdown is linked to the properties of the cusp catastrophe: the cusp diameter is shown to determine the magnitudes of all couplings, numerical basis-set and trajectory-integration requirements, and to determine the transmission coefficient κ used to understand deviations from transition-state theory.

Bond angle variations in XH3 [X = N, P, As, Sb, Bi]: the critical role of Rydberg orbitals exposed using a diabatic state model.

Ammonia adopts sp(3) hybridization (HNH bond angle 108°) whereas the other members of the XH3 series PH3, AsH3, SbH3, and BiH3 instead prefer octahedral bond angles of 90-93°. We use a recently developed general diabatic description for closed-shell chemical reactions, expanded to include Rydberg states, to understand the geometry, spectroscopy and inversion reaction profile of these molecules, fitting its parameters to results from Equation of Motion Coupled-Cluster Singles and Doubles (EOM-CCSD) calculations using large basis sets. Bands observed in the one-photon absorption spectrum of NH3 at 18.3 eV, 30 eV, and 33 eV are reassigned from Rydberg (formally forbidden) double excitations to valence single-excitation resonances. Critical to the analysis is the inclusion of all three electronic states in which two electrons are placed in the lone-pair orbital n and/or the symmetric valence σ* antibonding orbital. An illustrative effective two-state diabatic model is also developed containing just three parameters: the resonance energy driving the high-symmetry planar structure, the reorganization energy opposing it, and HXH bond angle in the absence of resonance. The diabatic orbitals are identified as sp hybrids on X; for the radical cations XH3(+) for which only 2 electronic states and one conical intersection are involved, the principle of orbital following dictates that the bond angle in the absence of resonance is acos(-1/5) = 101.5°. The multiple states and associated multiple conical intersection seams controlling the ground-state structure of XH3 renormalize this to acos[3 sin(2)(2(1/2)atan(1/2))/2 - 1/2] = 86.7°. Depending on the ratio of the resonance energy to the reorganization energy, equilibrium angles can vary from these limiting values up to 120°, and the anomalously large bond angle in NH3 arises because the resonance energy is unexpectedly large. This occurs as the ordering of the lowest Rydberg orbital and the σ* orbital swap, allowing Rydbergization to compresses σ* to significantly increase the resonance energy. Failure of both the traditional and revised versions of the valence-shell electron-pair repulsion (VSEPR) theory to explain the ground-state structures in simple terms is attributed to exclusion of this key physical interaction.

Efficient calculation of integrals in mixed ramp-Gaussian basis sets.

Algorithms for the efficient calculation of two-electron integrals in the newly developed mixed ramp-Gaussian basis sets are presented, alongside a Fortran90 implementation of these algorithms, RampItUp. These new basis sets have significant potential to (1) give some speed-up (estimated at up to 20% for large molecules in fully optimised code) to general-purpose Hartree-Fock (HF) and density functional theory quantum chemistry calculations, replacing all-Gaussian basis sets, and (2) give very large speed-ups for calculations of core-dependent properties, such as electron density at the nucleus, NMR parameters, relativistic corrections, and total energies, replacing the current use of Slater basis functions or very large specialised all-Gaussian basis sets for these purposes. This initial implementation already demonstrates roughly 10% speed-ups in HF/R-31G calculations compared to HF/6-31G calculations for large linear molecules, demonstrating the promise of this methodology, particularly for the second application. As well as the reduction in the total primitive number in R-31G compared to 6-31G, this timing advantage can be attributed to the significant reduction in the number of mathematically complex intermediate integrals after modelling each ramp-Gaussian basis-function-pair as a sum of ramps on a single atomic centre.

Mixed Ramp-Gaussian Basis Sets.

We discuss molecular orbital basis sets that contain both Gaussian and polynomial (ramp) functions. We show that, by modeling ramp-Gaussian products as sums of ramps, all of the required one- and two-electron integrals can be computed quickly and accurately. To illustrate our approach, we construct R-31+G, a mixed ramp-Gaussian basis in which the core basis functions of the 6-31+G basis are replaced by ramps. By performing self-consistent Hartree-Fock calculations, we show that the thermochemical predictions of R-31+G and 6-31+G are similar but the former has the potential to be significantly faster.