Toward Accurate Quantum Mechanical Thermochemistry: (1) Extensible Implementation and Comparison of Bond Additivity Corrections and Isodesmic Reactions.

Obtaining accurate enthalpies of formation of chemical species, ΔHf, often requires empirical corrections that connect the results of quantum mechanical (QM) calculations with the experimental enthalpies of elements in their standard state. One approach is to use atomization energy corrections followed by bond additivity corrections (BACs), such as those defined by Petersson et al. or Anantharaman and Melius. Another approach is to utilize isodesmic reactions (IDRs) as shown by Buerger et al. We implement both approaches in Arkane, an open-source software that can calculate species thermochemistry using results from various QM software packages. In this work, we collect 421 reference species from the literature to derive ΔHf corrections and fit atomization energy corrections and BACs for 15 commonly used model chemistries. We find that both types of BACs yield similar accuracy, although Anantharaman- and Melius-type BACs appear to generalize better. Furthermore, BACs tend to achieve better accuracy than IDRs for commonly used model chemistries, and IDRs can be less robust because of the sensitivity to the chosen reference species and reactions. Overall, Anantharaman- and Melius-type BACs are our recommended approach for achieving accurate QM corrections for enthalpies.

Kinetic Modeling of API Oxidation: (2) Imipramine Stress Testing.

Gauging the chemical stability of active pharmaceutical ingredients (APIs) is critical at various stages of pharmaceutical development to identify potential risks from drug degradation and ensure the quality and safety of the drug product. Stress testing has been the major experimental method to study API stability, but this analytical approach is time-consuming, resource-intensive, and limited by API availability, especially during the early stages of drug development. Novel computational chemistry methods may assist in screening for API chemical stability prior to synthesis and augment contemporary API stress testing studies, with the potential to significantly accelerate drug development and reduce costs. In this work, we leverage quantum chemical calculations and automated reaction mechanism generation to provide new insights into API degradation studies. In the continuation of part one in this series of studies [Grinberg Dana et al., Mol. Pharm. 2021 18 (8), 3037-3049], we have generated the first ab initio predictive chemical kinetic model of free-radical oxidative degradation for API stress testing. We focused on imipramine oxidation in an azobis(isobutyronitrile) (AIBN)/H2O/CH3OH solution and compared the model's predictions with concurrent experimental observations. We analytically determined iminodibenzyl and desimipramine as imipramine's two major degradation products under industry-standard AIBN stress testing conditions, and our ab initio kinetic model successfully identified both of them in its prediction for the top three degradation products. This work shows the potential and utility of predictive chemical kinetic modeling and quantum chemical computations to elucidate API chemical stability issues. Further, we envision an automated digital workflow that integrates first-principle models with data-driven methods that, when actively and iteratively combined with high-throughput experiments, can substantially accelerate and transform future API chemical stability studies.

Kinetic Modeling of API Oxidation: (1) The AIBN/H2O/CH3OH Radical "Soup".

Stress testing of active pharmaceutical ingredients (API) is an important tool used to gauge chemical stability and identify potential degradation products. While different flavors of API stress testing systems have been used in experimental investigations for decades, the detailed kinetics of such systems as well as the chemical composition of prominent reactive species, specifically reactive oxygen species, are unknown. As a first step toward understanding and modeling API oxidation in stress testing, we investigated a typical radical "soup" solution an API is subject to during stress testing. Here we applied ab initio electronic structure calculations to automatically generate and refine a detailed chemical kinetics model, taking a fresh look at API oxidation. We generated a detailed kinetic model for a representative azobis(isobutyronitrile) (AIBN)/H2O/CH3OH stress-testing system with a varied cosolvent ratio (50%/50%-99.5%/0.5% vol water/methanol) for 5.0 mM AIBN and representative pH values of 4-10 at 40 °C that was stirred and open to the atmosphere. At acidic conditions, hydroxymethyl alkoxyl is the dominant alkoxyl radical, and at basic conditions, for most studied initial methanol concentrations, cyanoisopropyl alkoxyl becomes the dominant alkoxyl radical, albeit at an overall lower concentration. At acidic conditions, the levels of cyanoisopropyl peroxyl, hydroxymethyl peroxyl, and hydroperoxyl radicals are relatively high and comparable, while, at both neutral and basic pH conditions, superoxide becomes the prominent radical in the system. The present work reveals the prominent species in a common model API stress testing system at various cosolvent and pH conditions, sets the stage for an in-depth quantitative API kinetic study, and demonstrates the usage of novel software tools for automated chemical kinetic model generation and ab initio refinement.

Vertical valence ionization potential benchmarks from equation-of-motion coupled cluster theory and QTP functionals.

The ionization potential (IP) of a molecule quantifies the energy required to remove an electron from the system. As such, it is a fundamental quantity in the context of redox chemistry, charge transfer, and molecular electronics. The accurate theoretical prediction of this property is therefore highly desirable for virtual materials design. Furthermore, vertical IPs are of interest in the development of many-body Green's function methods like the GW formalism, as well as density functionals and semiempirical methods. In this contribution, we report over 1468 vertical valence IPs calculated with the IP variant of equation-of-motion coupled cluster theory with singles and doubles (IP-EOM-CCSD) covering 155 molecules. The purpose of this is two-fold: First, the quality of the predicted IPs is compared with respect to experiments and higher-order coupled cluster theory. This confirms the overall high accuracy and robustness of this method, with some outliers which are discussed in detail. Second, a large set of consistent theoretical reference values for vertical valence IPs are generated. This addresses a lack of reliable reference data for lower-lying valence IPs, where experimental data are often unavailable or of dubious quality. The benchmark set is then used to assess the quality of the eigenvalues predicted by different density functional approximations (via Bartlett's IP-eigenvalue theorem) and the extended Koopmans' theorem approach. The QTP family of functionals are found to be remarkably accurate, low-cost alternatives to IP-EOM-CCSD.

Valence and charge-transfer optical properties for some SinCm (m, n ≤ 12) clusters: Comparing TD-DFT, complete-basis-limit EOMCC, and benchmarks from spectroscopy.

Accurate optical characterization of the closo-Si12C12 molecule is important to guide experimental efforts toward the synthesis of nano-wires, cyclic nano-arrays, and related array structures, which are anticipated to be robust and efficient exciton materials for opto-electronic devices. Working toward calibrated methods for the description of closo-Si12C12 oligomers, various electronic structure approaches are evaluated for their ability to reproduce measured optical transitions of the SiC2, Si2Cn (n = 1-3), and Si3Cn (n = 1, 2) clusters reported earlier by Steglich and Maier [Astrophys. J. 801, 119 (2015)]. Complete-basis-limit equation-of-motion coupled-cluster (EOMCC) results are presented and a comparison is made between perturbative and renormalized non-iterative triples corrections. The effect of adding a renormalized correction for quadruples is also tested. Benchmark test sets derived from both measurement and high-level EOMCC calculations are then used to evaluate the performance of a variety of density functionals within the time-dependent density functional theory (TD-DFT) framework. The best-performing functionals are subsequently applied to predict valence TD-DFT excitation energies for the lowest-energy isomers of SinC and Sin-1C7-n (n = 4-6). TD-DFT approaches are then applied to the SinCn (n = 4-12) clusters and unique spectroscopic signatures of closo-Si12C12 are discussed. Finally, various long-range corrected density functionals, including those from the CAM-QTP family, are applied to a charge-transfer excitation in a cyclic (Si4C4)4 oligomer. Approaches for gauging the extent of charge-transfer character are also tested and EOMCC results are used to benchmark functionals and make recommendations.

Vibrational Characterization of Radical Ion Adducts between Imidazole and CO2.

We address the molecular level origins of the dramatic difference in the catalytic mechanisms of CO2 activation by the seemingly similar molecules pyridine (Py) and imidazole (Im). This is accomplished by comparing the fundamental interactions of CO2 radical anions with Py and Im in the isolated, gas phase PyCO2- and ImCO2- complexes. These species are prepared by condensation of the neutral compounds onto a (CO2) n- cluster ion beam by entrainment in a supersonic jet ion source. The structures of the anionic complexes are determined by theoretical analysis of their vibrational spectra, obtained by IR photodissociation of weakly bound CO2 molecules in a photofragmentation mass spectrometer. Although the radical PyCO2- system adopts a carbamate-like configuration corresponding to formation of an N-C covalent bond, the ImCO2- species is revealed to be best described as an ion-molecule complex in which an oxygen atom in the CO2- radical anion is H-bonded to the NH group. Species that feature a covalent N-C interaction in ImCO2- are calculated to be locally stable structures, but are much higher in energy than the largely electrostatically bound ion-molecule complex. These results support the suggestion from solution phase electrochemical studies (Bocarsly et al. ACS Catal. 2012, 2, 1684-1692) that the N atoms are not directly involved in the catalytic activation of CO2 by Im.

A note on the accuracy of KS-DFT densities.

The accuracy of the density of wave function methods and Kohn-Sham (KS) density functionals is studied using moments of the density, ⟨rn⟩=∫ ρ rrndτ=∫0∞4πr2ρ(r)rndr,where n=-1,-2,0,1,2,and 3 provides information about the short- and long-range behavior of the density. Coupled cluster (CC) singles, doubles, and perturbative triples (CCSD(T)) is considered as the reference density. Three test sets are considered: boron through neon neutral atoms, two and four electron cations, and 3d transition metals. The total density and valence only density are distinguished by dropping appropriate core orbitals. Among density functionals tested, CAMQTP00 and ωB97x show the least deviation for boron through neon neutral atoms. They also show accurate eigenvalues for the HOMO indicating that they should have a more correct long-range behavior for the density. For transition metals, some density functional approximations outperform some wave function methods, suggesting that the KS determinant could be a better starting point for some kinds of correlated calculations. By using generalized many-body perturbation theory (MBPT), the convergence of second-, third-, and fourth-order KS-MBPT for the density is addressed as it converges to the infinite-order coupled cluster result. For the transition metal test set, the deviations in the KS density functional theory methods depend on the amount of exact exchange the functional uses. Functionals with exact exchange close to 25% show smaller deviations from the CCSD(T) density.

Automatic generation of reaction energy databases from highly accurate atomization energy benchmark sets.

In this contribution, we discuss how reaction energy benchmark sets can automatically be created from arbitrary atomization energy databases. As an example, over 11 000 reaction energies derived from the W4-11 database, as well as some relevant subsets are reported. Importantly, there is only very modest computational overhead involved in computing >11 000 reaction energies compared to 140 atomization energies, since the rate-determining step for either benchmark is performing the same 140 quantum chemical calculations. The performance of commonly used electronic structure methods for the new database is analyzed. This allows investigating the relationship between the performances for atomization and reaction energy benchmarks based on an identical set of molecules. The atomization energy is found to be a weak predictor for the overall usefulness of a method. The performance of density functional approximations in light of the number of empirically optimized parameters used in their design is also discussed.

Does the ionization potential condition employed in QTP functionals mitigate the self-interaction error?

Though contrary to conventional wisdom, the interpretation of all occupied Kohn-Sham eigenvalues as vertical ionization potentials is justified by several formal and numerical arguments. Similarly, the performance of density functional approximations (DFAs) for fractionally charged systems has been extensively studied as a measure of one- and many-electron self-interaction errors (MSIEs). These complementary perspectives (initially recognized in ab initio dft) are shown to lead to the unifying concept that satisfying Bartlett's IP theorem in DFA's mitigates self-interaction errors. In this contribution, we show that the IP-optimized QTP functionals (reparameterization of CAM-B3LYP where all eigenvalues are approximately equal to vertical IPs) display reduced self-interaction errors in a variety of tests including the He2+ potential curve. Conversely, the MSIE-optimized rCAM-B3LYP functional also displays accurate orbital eigenvalues. It is shown that the CAM-QTP and rCAM-B3LYP functionals show improved dissociation limits, fundamental gaps and thermochemical accuracy compared to their parent functional CAM-B3LYP.

Predictive coupled-cluster isomer orderings for some SinCm (m, n ≤ 12) clusters: A pragmatic comparison between DFT and complete basis limit coupled-cluster benchmarks.

The accurate determination of the preferred Si12C12 isomer is important to guide experimental efforts directed towards synthesizing SiC nano-wires and related polymer structures which are anticipated to be highly efficient exciton materials for the opto-electronic devices. In order to definitively identify preferred isomeric structures for silicon carbon nano-clusters, highly accurate geometries, energies, and harmonic zero point energies have been computed using coupled-cluster theory with systematic extrapolation to the complete basis limit for set of silicon carbon clusters ranging in size from SiC3 to Si12C12. It is found that post-MBPT(2) correlation energy plays a significant role in obtaining converged relative isomer energies, suggesting that predictions using low rung density functional methods will not have adequate accuracy. Utilizing the best composite coupled-cluster energy that is still computationally feasible, entailing a 3-4 SCF and coupled-cluster theory with singles and doubles extrapolation with triple-ζ (T) correlation, the closo Si12C12 isomer is identified to be the preferred isomer in the support of previous calculations [X. F. Duan and L. W. Burggraf, J. Chem. Phys. 142, 034303 (2015)]. Additionally we have investigated more pragmatic approaches to obtaining accurate silicon carbide isomer energies, including the use of frozen natural orbital coupled-cluster theory and several rungs of standard and double-hybrid density functional theory. Frozen natural orbitals as a way to compute post-MBPT(2) correlation energy are found to be an excellent balance between efficiency and accuracy.

Core-core and core-valence correlation energy atomic and molecular benchmarks for Li through Ar.

We have established benchmark core-core, core-valence, and valence-valence absolute coupled-cluster single double (triple) correlation energies (±0.1%) for 210 species covering the first- and second-rows of the periodic table. These species provide 194 energy differences (±0.03 mEh) including ionization potentials, electron affinities, and total atomization energies. These results can be used for calibration of less expensive methodologies for practical routine determination of core-core and core-valence correlation energies.

A density functional for core-valence correlation energy.

A density functional, εCV-DFT(ρc, ρv), describing the core-valence correlation energy has been constructed as a linear combination of εLY P (corr)(ρc), εV WN5 (corr)(ρc, ρv), εPBE (corr)(ρc, ρv), εSlater (ex)(ρc, ρv), εHCTH (ex)(ρc, ρv), εHF (ex)(ρc, ρv), and FCV-DFTNi,Zi, a function of the nuclear charges. This functional, with 6 adjustable parameters, reproduces (±0.27 kcal/mol rms error) a benchmark set of 194 chemical energy changes including 9 electron affinities, 18 ionization potentials, and 167 total atomization energies covering the first- and second-rows of the periodic table. This is almost twice the rms error (±0.16 kcal/mol) obtained with CCSD(T)/MTsmall calculations, but less than half the rms error (±0.65 kcal/mol) obtained with MP2/GTlargeXP calculations, and somewhat smaller than the rms error (±0.39 kcal/mol) obtained with CCSD/MTsmall calculations. The largest positive and negative errors from εCV-DFT(ρc, ρv) were 0.88 and -0.75 kcal/mol with the set of 194 core-valence energy changes ranging from +3.76 kcal/mol for the total atomization energy of propyne to -9.05 kcal/mol for the double ionization of Mg. Evaluation of the εCV-DFT(ρc, ρv) functional requires less time than a single SCF iteration, and the accuracy is adequate for any model chemistry based on the CCSD(T) level of theory.

An efficient extrapolation to the (T)/CBS limit.

We extrapolate to the perturbative triples (T)/complete basis set (CBS) limit using double ζ basis sets without polarization functions (Wesleyan-1-Triples-2ζ or "Wes1T-2Z") and triple ζ basis sets with a single level of polarization functions (Wesleyan-1-Triples-3ζ or "Wes1T-3Z"). These basis sets were optimized for 102 species representing the first two rows of the Periodic Table. The species include the entire set of neutral atoms, positive and negative atomic ions, as well as several homonuclear diatomic molecules, hydrides, rare gas dimers, polar molecules, such as oxides and fluorides, and a few transition states. The extrapolated Wes1T-(2,3)Z triples energies agree with (T)/CBS benchmarks to within ±0.65 mEh, while the rms deviations of comparable model chemistries W1, CBS-APNO, and CBS-QB3 for the same test set are ±0.23 mEh, ±2.37 mEh, and ±5.80 mEh, respectively. The Wes1T-(2,3)Z triples calculation time for the largest hydrocarbon in the G2/97 test set, C6H5Me(+), is reduced by a factor of 25 when compared to W1. The cost-effectiveness of the Wes1T-(2,3)Z extrapolation validates the usefulness of the Wes1T-2Z and Wes1T-3Z basis sets which are now available for a more efficient extrapolation of the (T) component of any composite model chemistry.

Evaluation of the heats of formation of corannulene and C60 by means of inexpensive theoretical procedures.

Inexpensive ab initio procedures that employ homologous sequences of isodesmic reactions for the calculation of enthalpies of formation of moderate-sized organic molecules were tested with benzene, naphthalene, phenanthrene, and triphenylene. Two size-consistent adjustable parameters were found to bring the calculated values within the uncertainty of the experimental values. These procedures were then applied to C20H10 (corannulene) and C60 (buckminsterfullerene). The results, specifically, Δ(f)H(298)(0)(C20H10) = 484 ± 4 kJ mol(-1) and Δ(f)H(298)(0)(C60) = 2531 ± 15 kJ mol(-1), are in excellent agreement with both the recent definitive W1h calculations of Karton et al. for corannulene [Δ(f)H(298)(0)(C20H10) = 485.2 ± 7.9 kJ mol(-1)] and their estimated value for buckminsterfullerene [Δ(f)H(298)(0)(C60) = 2521.6 ± 13.6 kJ mol(-1)] ( J. Phys. Chem. A 2013, 117, 1834-1842). We support their conclusion that the experimental values should be reexamined.

CCSD(T)/CBS atomic and molecular benchmarks for H through Ar.

We extrapolate to the coupled cluster single and double excitation and the perturbative triples (CCSD(T))/complete basis set (CBS) limit with a sequence of optimized n-tuple-ζ augmented polarization augmented (nZaPa) basis sets (n = 4, 5, 6, and 7) for 115 species representing the first two rows of the Periodic Table. The species include the entire set of atoms, positive and negative atomic ions, homonuclear diatomic molecules, and hydrides. The benchmark set also includes the rare gas dimers, polar molecules such as oxides and fluorides, and a few transition states for chemical reactions. The CCSD correlation energies agree with available CCSD-F12b/3C(FIX) values to within ±0.18 mEh root-mean-square (rms) deviation. The (T) components agree to within ±0.10 mEh and the total CCSD(T) correlation energies to within ±0.26 mEh or 0.1% rms deviation, which is probably the better measure, since the largest deviation is 0.43 mEh or 0.13%. These CBS limits can now be used as benchmarks to calibrate more approximate calculations using smaller basis sets. The sequence of basis sets provides data on convergence patterns for each component of the correlation energy.