ASSESSMENT OF EQUATION-OF-MOTION COUPLED-CLUSTER METHODS WITH APPROXIMATE TREATMENTS OF HIGHER-ORDER EXCITATIONS AND DEVELOPMENT OF NOVEL SCHEMES FOR ACCURATE CALCULATIONS OF DIRADICAL ELECTRONIC SPECTRA AND BOND BREAKING By Adeayo Olayinka Ajala A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Chemistry - Doctor of Philosophy 2017 ABSTRACT ASSESSMENT OF EQUATION-OF-MOTION COUPLED-CLUSTER METHODS WITH APPROXIMATE TREATMENTS OF HIGHER-ORDER EXCITATIONS AND DEVELOPMENT OF NOVEL SCHEMES FOR ACCURATE CALCULATIONS OF DIRADICAL ELECTRONIC SPECTRA AND BOND BREAKING By Adeayo Olayinka Ajala The development and implementation of electronic structure methods based on the exponential wave function ansatz of the single-reference coupled-cluster (CC) theory and its extensions to excited states exploiting the equation-of-motion (EOM) and linear response frameworks have witnessed great success in a wide range of applications, but there are areas of chemistry, especially studies of chemical reaction pathways and photochemistry, where further improvements in the existing CC and EOMCC methodologies are needed. In order to make progress in this area, it is important to evaluate the quality of the results that the existing CC/EOMCC methods provide, particularly in applications involving the interpretation and prediction of photochemical phenomena and electronic excitations spectra involving closed- and open-shell molecules. Thus, in the first part of this PhD project we use a database set of 28 organic molecules ranging from linear polyenes, unsaturated cyclic hydrocarbons, aromatic hydrocarbons, and heterocycles to aldehydes, ketones, amides, and nucleobases to examine the performance of the completely renormalized (CR) EOMCC approaches for excited electronic states, in which the relatively inexpensive non-iterative corrections due to triple excitations are added to the energies obtained with the standard EOMCC approach with singles and doubles, abbreviated as EOMCCSD. We focus on two variants of the approximately size-intensive CR-EOMCC methodology with singles, doubles, and noniterative triples, abbreviated as δ-CR-EOMCCSD(T), and the analogous two variants of the newer, rigorously size-intensive, left-eigenstate δ-CR-EOMCC(2,3) approach based on the biorthogonal formulation of the method of moments of CC equations. In the second part of this dissertation, we focus on the development of new EOMCC methods that are particularly well-suited for accurate calculations of diradical electronic spectra and single bond breaking. They are the cost-effective variants of the doubly electronattached (DEA) EOMCC methodologies with up to 3-particle–1-hole (3p-1h) or 4-particle– 2-hole (4p-2h) excitations, abbreviated as DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC (3p-1h,4p-2h){Nu }, respectively, which utilize the idea of applying a linear electron-attaching operator to the correlated CC ground state of an (N − 2)-electron closed-shell reference system in order to generate the ground and excited states of the N -electron open-shell species of interest, while using Nu active unoccupied orbitals to select the dominant 3p-1h and 4p-2h terms. We demonstrate that the relatively inexpensive DEA-EOMCC(3p-1h,4p-2h){Nu } method greatly reduces the computational costs of the parent active-space DEA-EOMCC (4p-2h){Nu } and full DEA-EOMCC(4p-2h) approaches, needed to obtain highly accurate results for open-shell systems having two electrons outside the closed-shell cores, with virtually no loss in accuracy of the resulting excitation and dissociation energies. We also show that the active-space DEA-EOMCC(3p-1h){Nu } method accurately reproduces the results of the parent DEA-EOMCC(3p-1h) calculations at the small fraction of the cost. In addition to a series of benchmark examples that illustrate the performance of the DEA-EOMCC(3p-1h){Nu }, DEA-EOMCC(3p-1h,4p-2h){Nu }, and other DEA-EOMCC approaches with 3p-1h and 4p-2h excitations, including singlet–triplet gaps in methylene, trimethylenemethane, and several antiaromatic diradicals and bond breaking in the fluorine molecule, we provide the most essential details of DEA-EOMCC equations with an active-space treatment of 3p-1h and 4p-2h terms, as implemented in our codes and interfaced with the GAMESS package. Copyright by ADEAYO OLAYINKA AJALA 2017 This is dedicated to my lovely wife Danielle and my daughter Dara, my supportive parents Oluwatoyin and Abayomi (deceased), and my siblings Babatunde (deceased) and Olanrewaju. v ACKNOWLEDGMENTS I owe a great deal of gratitude to my PhD advisor, Professor Piotr Piecuch, for his patience and thorough guidance throughout my graduate career. Since my first day in his group he has driven me to be a better scientist in paying attention to details and discovering my potential. I am also indebted to him for his continuous support throughout my academic journey at Michigan State University. I would like to thank the other members of my committee, Professor Benjamin G. Levine, Professor Robert I. Cukier, and Professor Rémi Beaulac for all their support, advice, and patience in overseeing my graduate study. It is also important to mention the valuable effort of Dr. Jun Shen from our group. As a result of his guidance and several discussions, I progressed much more quickly than I would have without him, especially in the development work presented in this dissertation. I would also like to thank Dr. Jun Shen for checking some of the equations presented in this work. My sincere appreciation also goes to Dr. Jared Hansen for his great assistance in one of the projects presented in this dissertation. Lastly, I would like to thank other members of the Piecuch group, both past and present, for the help and friendship they provided, including Dr. Nicholas P. Bauman, Mr. Jorge Emiliano Deustua, Mr. Ilias Magoulas, and Mr. Stephen H. Yuwono. vi TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Chapter 2 Project Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Chapter 3 3.1 3.2 3.3 3.4 3.5 3.6 Benchmarking Completely Renormalized Equation-of-Motion Coupled-Cluster Methods with Triple Excitations . . . . . . . . 10 Background Information and Motivation . . . . . . . . . . . . . . . . . . . . 10 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Benchmark Molecules and Computational Details . . . . . . . . . . . . . . . 18 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4.1 Geometries and their Effect on the Calculated Vertical Excitation Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4.2 Remarks on the Calculated Vertical Excitation Spectra for Molecules in the Database of Ref. [17] . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.2.1 Linear Polyenes: Ethene, E-butadiene, all-E-hexatriene, and all-E-octatetraene . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4.2.2 Unsaturated Cyclic Hydrocarbons: Cyclopropene, Cyclopentadiene, and Norbornadiene . . . . . . . . . . . . . . . . . . 50 3.4.2.3 Aromatic Hydrocarbons: Benzene and Naphthalene . . . . . 52 3.4.2.4 Heterocycles: Furan, Pyrrole, Imidazole, Pyridine, Pyrazine, Pyrimidine, Pyridazine, s-Triazine, and s-Tetrazine . . . . . 53 3.4.2.5 Aldehydes and Ketones: Formaldehyde, Acetone, and p-Benzoquinone . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.2.6 Amides: Formamide, Acetamide, and Propanamide . . . . . 61 3.4.2.7 Nucleobases: Cytosine, Uracil, Thymine, and Adenine . . . 64 Statistical Error Analysis for the Entire Benchmark Set . . . . . . . . . . . . 65 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Chapter 4 4.1 4.2 Economical Doubly Electron-Attached Equation-of-Motion Coupled-Cluster Methods with an Active-Space Treatment of 3-particle–1-hole and 4-particle–2-hole Excitations . . . . . . . Background Information and Motivation . . . . . . . . . . . . . . . . . . . . Theory and Algorithmic Details . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Basic Elements of the DEA-EOMCC Formalism and the Previously Developed DEA-EOMCC(3p-1h), DEAEOMCC(4p-2h), and DEA-EOMCC(4p-2h){Nu } Approximations . . vii 107 108 114 114 The Active-Space DEA-EOMCC(3p-1h){Nu } and DEAEOMCC(3p-1h,4p-2h){Nu } Approaches Developed in this Thesis Project119 4.2.3 Key Details of the Efficient Computer Implementation of the ActiveSpace DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC (3p-1h,4p-2h){Nu } Methods . . . . . . . . . . . . . . . . . . . . . . . 122 4.2.4 The Remaining Algorithmic Details and Illustrative Timings of DEAEOMCC(3p-1h){Nu } and DEA-EOMCC (3p-1h, 4p-2h){Nu } Calculations . . . . . . . . . . . . . . . . . . . . . 137 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 4.3.1 Adiabatic Energy Gaps Involving Low-Lying Singlet and Triplet States of Methylene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 4.3.2 The Complete Basis Set Limit Investigation of Singlet and Triplet States of Methylene . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 4.3.3 Adiabatic Singlet–Triplet Gap in Trimethylenemethane . . . . . . . . 154 4.3.4 Vertical Singlet–Triplet Gaps in Antiaromatic Diradicals . . . . . . . 163 4.3.5 F–F Bond Dissociation in the Fluorine Molecule . . . . . . . . . . . . 171 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 4.2.2 4.3 4.4 Chapter 5 Summary and Future Outlook . . . . . . . . . . . . . . . . . . . . 176 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 APPENDIX A Derivation of the Non-Iterative Energy Correction Defining the Biorthogonal MMCC Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 APPENDIX B Factorized Form of the DEA-EOMCC(4p-2h) Equations Based on CCSD Reference Wave Functions . . . . . . . . . . . . . . . . . . . . . . . . 186 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 viii LIST OF TABLES Table 3.1: Symmetry unique Cartesian coordinates of the ground-state geometries for the 28 molecules comprising the benchmark set resulting from the MP2/631G∗ optimizations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Table 3.2: Symmetry unique Cartesian coordinates of the ground-state geometries of the 28 molecules comprising the benchmark set resulting from the CRCC(2,3),D/TZVP optimizations carried out in this work. . . . . . . . . . 29 Table 3.3: Vertical excitation energies (in eV) for singlet states of all molecules in the test set using the geometries optimized at the MP2/6-31G∗ level.a ,b . . . 36 Table 3.4: Vertical excitation energies (in eV) for singlet states of all molecules in the test set using, in the case of the EOMCCSD and δ-CR-EOMCC results, the geometries optimized at the CR-CC(2,3),D/TZVP level.a . . . . . . 41 Table 3.5: Statistical error analyses of the various CC/EOMCC data, including EOMCCSD, δ-CR-EOMCCSD(T),IA (δ-CR(T),IA), δ-CR-EOMCCSD(T),ID (δ-CR(T),ID), δ-CR-EOMCC(2,3),A (δ-CR(2,3),A), δ-CR-EOMCC(2,3),D (δ-CR(2,3),D), CC3, and EOMCCSDT-3, from their CASPT2 counterparts. 68 Table 3.6: Statistical error analyses of the various CC/EOMCC data, including EOMCCSD, δ-CR-EOMCCSD(T),IA (δ-CR(T),IA), δ-CR-EOMCCSD(T),ID (δ-CR(T),ID), δ-CR-EOMCC(2,3),A (δ-CR(2,3),A), δ-CR-EOMCC(2,3),D (δ-CR(2,3),D), CC3, and EOMCCSDT-3, from their TBE-2 counterparts. 68 Table 3.7: Statistical error analyses of the EOMCCSD, δ-CR-EOMCCSD(T),IA (δ-CR(T),IA), δ-CR-EOMCCSD(T),ID (δ-CR(T),ID), δ-CR-EOMCC(2,3),A (δ-CR(2,3),A), and δ-CR-EOMCC(2,3),D (δ-CR(2,3),D) data from their CC3 counterparts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Table 3.8: Statistical error analyses of the EOMCCSD, δ-CR-EOMCCSD(T),IA (δ-CR(T),IA), δ-CR-EOMCCSD(T),ID (δ-CR(T),ID), δ-CR-EOMCC(2,3),A (δ-CR(2,3),A), and δ-CR-EOMCC(2,3),D (δ-CR(2,3),D) data from their EOMCCSDT-3 counterparts. . . . . . . . . . . . . . . . . . . . . . . . . . 69 Table 4.1: Explicit algebraic expressions for the one- and two-body matrix elements (CCSD) of H̄N,open . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 ix Table 4.2: The various classes of restricted projections that must be considered when generating the computationally efficient form of the equations defining the DEA-EOMCC(3p-1h, 4p-2h){Nu } eigenvalue problem. . . . . . . . . . . . 128 Table 4.3: A comparison of CPU times required by the various DEA-EOMCC calculations characterizing the X 3 A2 state of TMM, as described by the cc-pVDZ and, in parentheses, cc-pVTZ basis sets, along with the formal scalings of the most expensive steps in the diagonalization of H̄N,open with no , nu , and Nu .a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Table 4.4: A comparison of the full CI and various DEA-EOMCC adiabatic excitation energies, along with the corresponding MaxUE and NPE values relative to full CI, characterizing the low-lying states of methylene, as described by the TZ2P basis set.a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 Table 4.5: Comparison of the adiabatic excitation energies (in eV) for the low-lying states of CH2 , as obtained with various DEA-EOMCC approaches using the ACxZ (x=T, Q, and 5) basis sets and (N − 2)-electron RHF orbitals, and extrapolating to the CBS limit, with the corresponding QMC results and experiment.a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 Table 4.6: The selected DEA-EOMCC results for the adiabatic singlet–triplet separation, ∆ES−T = E(B 1 A1 ) − E(X 3 A02 ), in kcal/mol, in trimethylenemethane (TMM), as described by the cc-pVDZ and, in parentheses, ccpVTZ basis sets, calculated using the SF-DFT/6-31G(d) geometries.a . . 158 Table 4.7: The DEA-EOMCC results for the vertical singlet–triplet gaps, ∆ES−T = E(S) − E(T), in kcal/mol, in systems 1–7, as described by the VDZ and mATZ basis sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 Table 4.8: A comparison of the ∆ES−T values (in kcal/mol) characterizing systems 1–7 obtained with the DEA-EOMCC[4p-2h] extrapolation in the DEAEOMCC(4p-2h){Nu }/VDZ calculations with the UCCSD(T), UBD(T), ROHF-MkCCSD, and CASSCF-MkCCSD results reported by Saito et al. 171 Table 4.9: A comparison of the full CI and various DEA-EOMCC ground-state energies of F2 at the equilibrium (Re = 2.66816 bohr) and a few other F–F distances, along with the corresponding MUE and NPE values relative to full CI, obtained with the DZ basis set.a . . . . . . . . . . . . . . . . . . 172 x LIST OF FIGURES Figure 3.1: Benchmark set of molecules considered in this work. . . . . . . . . . . . 19 Figure 3.2: Correlation plot of the δ-CR-EOMCC(2,3),D/TZVP excitation energies (in eV) computed using the MP2/6-31G∗ and CR-CC(2,3),D/TZVP geometries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 3.3: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . . . . . . . . . . . . . . . . . 78 Figure 3.4: Correlation plot for the calculated singlet excited states: δ-CR- EOMCCSD(T),IA/TZVP (δ-CR-(T),IA) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . 79 Figure 3.5: Correlation plot for the calculated singlet excited states: δ-CR- EOMCCSD(T),ID/TZVP (δ-CR-(T),ID) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . 80 Figure 3.6: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . . . . . . 81 Figure 3.7: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . . . . . . 82 Figure 3.8: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . . . . . . . . . . . . . . . . . . . . . . 83 Figure 3.9: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR-(T),IA) vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . xi 84 Figure 3.10: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR-(T),ID) vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . 85 Figure 3.11: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR-(2,3),A) vs TBE-2 data. The absence of the table on the right hand side indicates that there are no outliers. . . . . . . . . . . . . 86 Figure 3.12: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR-(2,3),D) vs TBE-2 data. The absence of the table on the right hand side indicates that there are no outliers. . . . . . . . . . . . . 87 Figure 3.13: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . . . . . . . . . . . . . . . . . . . . . . 88 Figure 3.14: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP (δ-CR-(T),IA) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . 89 Figure 3.15: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP (δ-CR-(T),ID) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . 90 Figure 3.16: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . . . . . . . . . 91 Figure 3.17: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . . . 92 Figure 3.18: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/631G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. . . . . . . . . . . . . . . . . 93 xii Figure 3.19: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP (δ-CR-(T),IA) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . 94 Figure 3.20: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP (δ-CR-(T),ID) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . 95 Figure 3.21: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. . . . . 96 Figure 3.22: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 97 Figure 3.23: Normal distribution curves for the deviation of the δ-CR-EOMCCSD(T),IA excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 3.24: Normal distribution curves for the deviation of the δ-CR-EOMCCSD(T),ID excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Figure 3.25: Normal distribution curves for the deviation of the δ-CR-EOMCC(2,3),A excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Figure 3.26: Normal distribution curves for the deviation of δ-CR-EOMCC(2,3),D excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Figure 4.1: The key elements of the algorithm used to compute (CCSD) (+2) hΦAbc )C |Φi in the efficient implementation of the DEAk |(H̄N,open Rµ EOMCC(3p-1h,4p-2h){Nu } method. . . . . . . . . . . . . . . . . . . . . 136 xiii Figure 4.2: Diradical systems considered in this work. 1: C4 H4 , 2: [C5 H5 ]+ , 3: C4 H3 NH2 , 4: C4 H3 CHO, 5: C4 H2 NH2 CHO, 6: C4 H2 -1, 2-(CH2 )2 , and 7: C4 H2 -1, 3-(CH2 )2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 xiv Chapter 1 Introduction Quantum chemistry has become an indispensable part of modern molecular science, providing us with concepts and ideas which play a central role in our understanding of chemical systems and processes. The development of quantum chemistry has also provided us with versatile computational tools, whose application in many areas of physical and biological sciences offer quantitative data and useful insights, helping us to predict, verify, and understand experimental observations and measurements. Given the plethora of electronic structure methods, accompanied by the emergence and development of program suites, it is essential to critically evaluate the performance and accuracy of existing quantum chemistry approximations and to formulate and implement new and improved ideas. Both of these aspects are reflected in this dissertation. One of the most effective approaches to the critical assessment of the existing methods is by testing these methods on databases that contain larger numbers of well-characterized molecular species. A few examples of the databases that exist for benchmarking a broad range of ground-state properties of ab initio and density functional theory (DFT) methodologies are the Gaussian Gn test sets [1–6] for atomization energies, ionization potentials, electron affinities, proton affinities, and enthalpies of formation, the S22 [7], S26 [8], and S66 [9] training sets for non-covalent interactions, the ATcT [10–12] benchmark set for thermochemical data, the DBH24 [13] test set for barrier heights, the W4 [14] training set for 1 total atomization energies, the 3dBE70 [15] benchmark set for average bond energies of 3d transition-metal-containing compounds, and the “mindless” benchmark set [16], which is a diversity-oriented collection of randomly generated molecules for main group thermochemical properties. Unfortunately, for the corresponding excited-state properties, where theoretical and computational quantum chemistry methods have great potential to make significant contributions to interpretation and prediction of photochemical phenomena and electronic molecular spectra, only one such comprehensive test set has been proposed so far [17]. As a result, fewer methods for excited states have been systematically tested on large datasets. The test set developed in Ref. [17] consists of 28 organic molecules, 149 singlet vertical excitation energies, and 72 triplet vertical excitation energies, with the majority of excited states being dominated by one-electron transitions and some having more substantial two-electron excitation components. In Refs. [17–25], several quantum chemistry approaches, including the complete-active-space self-consistent field [26, 27] based second-order perturbation theories, such as CASPT2 [28, 29] and NEVPT2 [30–33], a variety of coupled-cluster (CC) [34–37] linear-response [38–47] and equation-of-motion (EOM) [48–54] CC methods, time-dependent DFT [55, 56], and the DFT-based multi-reference configuration interaction (CI) [57] approximation, have been tested using this set [17–25], providing useful information. We should also mention the analogous work from Pal and co-workers [58], which highlights a benchmark investigation of the ionized (IP) EOMCC approach with singles and doubles (IP-EOMCCSD) [59–62] using the geometries and selected spectroscopic properties of a variety of doublet radicals, comparing its performance with experiment and other CC approaches. All of the above examples illustrate that, benchmarking using larger datasets has become an important procedure in the evaluation and assessment of the quality and validation of 2 ab initio electronic structure methods, including the higher-level methods of CC theory and its extensions to excited and open-shell states. Following on this theme, the first part of this dissertation will focus on using the comprehensive dataset of Ref. [17] to examine the performance of new generations of EOMCC methods developed by our group based on the idea of correcting the excitation energies obtained with EOMCCSD for the effects of triple excitations, abbreviated as δ-CR-EOMCCSD(T) [63] (see, also, Refs. [64–66]) and δ-CREOMCC(2,3) [67–69]. In doing this, we will rely on our exhaustive study published in Ref. [70] and the accompanying high-level EOMCC and linear response CC data reported in Refs. [17, 19, 20, 22, 24, 25]. Our effort to benchmark the non-iterative δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) triples corrections to EOMCCSD will show that these approaches are capable of greatly improving the results of EOMCCSD calculations without the need to go all the way to prohibitively expensive EOMCC levels, such as EOMCCSDT [71–73], where triple excitations in the cluster and EOM excitation operators are treated fully, but there are situations, especially the excitation spectra of open-shell species, such as radicals and diradicals, where it is worth considering alternative approaches. The δ-CR-EOMCCSD(T) and δ-CREOMCC(2,3) methods are cost efficient and capable of handling excited states dominated by one- as well as two-electron transitions, but like all single-reference particle-conserving CC/EOMCC schemes implemented using the spin-integrated spin-orbital equations they break the spin-symmetry of the non-relativistic Hamiltonian in open-shell systems, even when the reference determinant is of the restricted (e.g, restricted open-shell Hartree–Fock or ROHF) type. The δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) methods, being highly correlated, eliminate much of this problem numerically, offering accurate results for radical [66, 68] and diradical [74] electronic spectra, even when the excited states of interest have 3 a multi-reference character, but the problem of spin contamination in open-shell systems remains. Within the EOMCC framework that interests us in this thesis, the cleanest approach to electronic spectra of open-shell systems, resulting in a rigorously spin-adapted description, is to consider the electron-attached (EA) [75–81] and IP [59–62, 77, 79–85] EOMCC theories and their extensions to multiple electron-attachment and multiple ionization cases, such as the case of doubly electron-attached (DEA) EOMCC framework and its doubly ionized (DIP) counterpart [53, 86–94]. To appreciate these kinds of mehods, we begin our discussion with the EA- and IP-EOMCC methodologies. In the EA- and IP-EOMCC theories, the ground and excited states of an (N ± 1)-electron system are generated by applying a linear electronattaching or ionizing operator to the correlated CC ground state of the related N -electron closed-shell core. The basic EA-EOMCCSD [75, 76] and IP-EOMCCSD [59–62] approximations, in which the electron-attaching and ionizing operators of EOMCC are truncated at the 2-particle–1-hole (2p-1h) and 2-hole–1-particle (2h-1p) terms, respectively, have difficulties with accurately describing excitation spectra of radicals [75, 79, 80, 95, 96, 96, 97], but, in analogy to particle-conserving EOMCC schemes, where one needs to go beyond doubles, one can resolve this inadequacy through the inclusion of 3p-2h/3h-2p [79, 80] and 4p-3h/4h-3p [85, 96] components of the electron-attaching and ionizing operators of EOMCC. The resulting EA-EOMCCSDT [78], IP-EOMCCSDT [82, 83], EA-EOMCCSD(3p-2h) [79, 80], IP-EOMCCSD(3h-2p) [79, 80], EA-EOMCCSD(4p-3h) [96], IP-EOMCCSD(4h-3p) [96], EAEOMCCSDTQ [85], and IP-EOMCCSDTQ [85] schemes greatly improve the EA-EOMCCSD and IP-EOMCCSD results even when the excited states of radicals of interest gain a significant multi-determinantal character, but computer costs associated with such high-level EAand IP-EOMCC methods are quite high. For example, the EA-EOMCCSDT calculations, 4 where the cluster operator is truncated at the three-body terms and the EOM electronattaching operator at 3p-2h excitations, are characterized by the expensive CPU iterative steps which scale as n3o n5u in the underlying CCSDT computations for the closed-shell core and n2o n5u in the steps related to the diagonalization of the similarity-transformed Hamiltonian (no and nu are the numbers of occupied and unoccupied orbitals used in the post-SCF calculations). These kinds of iterative steps are not computationally affordable for routine chemical applications for larger molecular problems. In order to address the issue of large costs of high-order EA- and IP-EOMCC calculations, the active-space CC [98–109] and EOMCC [71, 72, 110–114] approaches (see Ref. [115] for a review) have been extended to the EA-EOMCC and IP-EOMCC [79–81, 96] methodologies. The examples of the active-space EA- and IP-EOMCC methods are EA-EOMCCSDt and IP-EOMCCSDt [79–81]. The EA-EOMCCSDt and IP-EOMCCSDt methods are obtained by reducing the numbers of 3p-2h and 3h-2p amplitudes in the parent EA-EOMCCSD(3p-2h) and IP-EOMCCSD(3h-2p) calculations with the help of active orbitals. As shown in Refs. [79–81, 96, 116], this offers major savings in the computational effort compared to the EAEOMCCSD(3p-2h) and IP-EOMCCSD(3h-2p) approaches and their EA-EOMCCSDT and IP-EOMCCSDT counterparts, which treat 3p-2h and 3h-2p terms fully, with virtually no loss of accuracy in the calculated radical excitation spectra. This is telling us that one should be able to use the idea of active orbitals to select higher-order terms in the DEAand DIP-EOMCC frameworks, which we discuss next. In the DEA- and DIP-EOMCC approaches [53, 86–94], which are particularly well-suited to describe systems having two electrons outside the N -electron closed-shell core, such as diradicals, the ground and excited states of (N + 2)- or (N − 2)-electron species are obtained by applying suitably defined operators that attach two electrons to (the DEA case) 5 or remove two electrons from (the DIP case) the N -electron reference system, while relaxing the remaining electrons. The levels of theory that leads to the desired accuracies in describing electronic spectra of diradicals (errors in the energy gaps of 1 kcal/mol or smaller) are DEA-EOMCC(4p-2h) in the DEA case and DIP-EOMCC(4h-2p) in the DIP case. In the DEA-EOMCC(4p-2h) approach, we truncate the electron-attaching operator of EOMCC at 4p-2h terms. In the DIP-EOMCC(4h-2p) scheme, the corresponding ionizing operator is truncated at 4h-2p component. As demonstrated in Refs. [92–94], the incorporation of 4p-2h and 4h-2p excitations in the DEA- and DIP-EOMCC frameworks greatly improving the accuracy compared to the basic DEA-EOMCC(3p-1h) and DIP-EOMCC(3h-1p) models, where the EOM operators attaching or removing two electrons from the related closed-shell cores are truncated at 3p-1h and 3h-1p levels. Unfortunately, as in all CC/EOMCC methods with higher-than-double excitations, the full implementation of the DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p) approaches that give these high accuracies comes at a high price, resulting in schemes that have very expensive iterative n2o n6u (the DEA case) and n4o n4u (the DIP case) steps. One way to incorporate 4p-2h and 4h-2p excitations without running into the prohibitive costs of the full DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p) computations is to employ the previously discussed active-space ideas to select the dominant higher-rank excitation amplitudes [92, 93]. This is particularly true in the DEA-EOMCC(4p-2h) case, where the use of active orbitals to select the dominant 4p-2h components reduces the iterative n2o n6u steps of full DEA-EOMCC(4p-2h) to a much more acceptable Nu2 n2o n4u level, where Nu ( nu ) is the number of active orbitals unoccupied in the N -electron closed-shell reference system. Similar reduction in the computational effort is observed when one selects the dominant 4h-2p terms in the DIP-EOMCC(4h-2p) approach using No (< no ) active occupied orbitals. This allows us to replace the original n4o n4u scaling of DIP-EOMCC(4h-2p) by much 6 less expensive No2 n2o n4u . The results reported in Refs. [92, 93] show that the DEA- and DIPEOMCC methods with an active-space treatment of 4p-2h and 4h-2p excitations accurately reproduce the parent, nearly exact DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p) data at the fraction of the computational cost. The active-space DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p) approaches of Refs. [92, 93] have been very successful in accurately describing diradical electronic spectra and single bond breaking, but one issue that prevents such methods from becoming more popular is the high cost of handling the lower-order 3p-1h terms within the active-space DEAEOMCC(4p-2h) framework. Indeed, the cost of computing 3p-1h terms fully scales as iterative no n5u , making the full treatment of 3p-1h contributions very expensive when basis sets used in the calculations are larger (so that nu is large). This issue is addressed in this dissertation by developing the new form of the DEA-EOMCC approach, in which both 3p-1h and 4p-2h terms are treated with the help of active orbitals, allowing us to reduce the scalings of the CPU steps associated with the full treatment of these terms from the original and prohibitively expensive no n5u and n2o n6u levels to Nu no n4u and Nu2 n2o n4u , respectively [94]. As shown in this thesis research and as demonstrated in our recently published studies [94, 117], the DEA-EOMCC calculations with an active-space treatment of 3p-1h and 4p-2h contributions provide highly accurate results for electronic spectra of diradicals and single bond breaking, which can compete with those obtained with the parent full DEA-EOMCC(4p-2h) approach, at the small fraction of the computational cost of the latter method. The byproduct of this work is the implementation of the active-space DEA-EOMCC(3p-1h) scheme, which is not as accurate as its higher-level DEA-EOMCC(4p-2h) counterpart, but still quite useful in diradical applications. The active-space DEA-EOMCC(3p-1h) approach, where 4p-2h correlations are neglected, uses active orbitals to select the dominant 3p-1h contri7 butions, reducing the expensive no n5u steps of full DEA-EOMCC(3p-1h) to a considerably more practical Nu no n4u level. As shown in this dissertation, and as demonstrated in our recent work [94, 117], the active-space DEA-EOMCC(3p-1h) method allows us to accurately reproduce the parent DEA-EOMCC(3p-1h) data, where 3p-1h terms are treated fully at the small fraction of the computational costs of full DEA-EOMCC(3p-1h) calculations. 8 Chapter 2 Project Objectives The main objectives of this dissertation work are A. Benchmark the δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) approaches for vertical excitation energies using a database of 28 small to medium sized organic molecules. B. Develop and apply reduced-cost DEA-EOMCC(4p-2h) method with an active-space treatment of 3p-1h as well as 4p-2h excitations and its lower-level counterpart truncated at 3p-1h terms to study the electronic spectra of diradicals and single bond dissociations. C. Describe the details of our implementations of the spin- and symmetry-adapted DEAEOMCC methods based on the CCSD reference wave function with an active-space treatment of 3p-1h and 4p-2h contributions and how the DEA-EOMCC equations were derived, factorized, and translated into FORTRAN code. 9 Chapter 3 Benchmarking Completely Renormalized Equation-of-Motion Coupled-Cluster Methods with Triple Excitations 3.1 Background Information and Motivation As mentioned in the Introduction, the first part of this dissertation focuses on using the test set introduced in Ref. [17], with additional information and updates provided in Refs. [19, 21, 24, 25], to examine the performance of the newer generations of the non-iterative triples corrections to the vertical excitation energies obtained in the EOMCCSD [48–50, 118–121] calculations, abbreviated as δ-CR-EOMCCSD(T) [63] and δ-CR-EOMCC(2,3) [69]. These corrections, resulting from the excited-state extensions [67, 122–124] of the method of moments of CC equations (MMCC) [65, 68, 122, 125–129], are the approximately sizeintensive [42, 130] (δ-CR-EOMCCSD(T)) or rigorously size-intensive (δ-CR-EOMCC(2,3)) modifications of the CR-EOMCC approaches called CR-EOMCCSD(T) [63, 65, 66] and CREOMCC(2,3) [67, 68]. Herein lies our motivation for benchmarking the δ-CR-EOMCCSD(T) 10 and δ-CR-EOMCC(2,3) methods using the test set of Ref. [17]. It is generally thought that the basic EOMCCSD [48–50] and linear-response CCSD [41, 42] approximations, which construct the desired excited-state information on top of the conventional CCSD [131–134] ground state, and which are characterized by the relatively inexpensive iterative CPU steps that scale as n2o n4u or N 6 , where N is a measure of the system size, provide an accurate description of excited states dominated by one-electron transitions. However, the EOMCCSD method is often not accurate enough to obtain a quantitative description of such states, especially when larger polyatomic species are examined (cf., e.g., Refs. [69, 135–137]; for a thorough evaluation of a number of EOMCC methods, including EOMCCSD, illustrating the same, see Refs. [17, 19, 20, 22, 24, 25]). It also fails to describe states with significant two-electron excitation contributions [63–68, 70–72, 74, 122]. One can address these shortcomings by including the effects of connected triple excitations, as is done in full EOMCCSDT [71–73]. While the full treatment of triple excitations substantially improves the description of excited electronic states, often providing virtually exact results (see, e.g., Refs. [64, 65, 71–73, 115, 138–141]), it is also accompanied by a steep increase in the CPU times characterizing the EOMCCSDT computations (the iterative n3o n5u or N 8 steps), limiting its applicability to systems with up to a dozen or so correlated electrons and smaller basis sets. Thus, if one is to make use of the EOMCC methodologies in accurate calculations of molecular electronic spectra in medium size and larger systems, including vertical excitation energies that interest us in this work, EOMCC schemes which can account for the effects of triples in an approximate, cost effective, and yet reliable manner need to be employed. There are several ways of incorporating triple excitations in the EOMCC and linearresponse CC formalisms without running into the prohibitive computational costs of full 11 EOMCCSDT. For example, one can select the dominant triply excited components of the cluster operator T that defines the underlying ground-state CC wave function and threebody components of the linear excitation operator Rµ in the EOMCC wave function ansatz through the use of active orbitals, as is done in the active-space EOMCCSDt method [71, 72, 110–112] (see Refs. [115, 129] for reviews; cf., also, Refs. [98–109] for the closely related activespace CC methods for the ground electronic states). While this allows for full-EOMCCSDTquality results at the cost of EOMCCSD times a prefactor proportional to the numbers of active occupied and active unoccupied orbitals used to select the triples, the approach is no longer strictly speaking black-box as one has to select the active orbitals. One can also contemplate approaches for identifying the most important triples contributions through the many-body perturbation theory or the aforementioned MMCC analysis. Some examples of these types of non-iterative triples methods are the EOMCC(2)PT(2) approach [142] and its size-intensive EOMCCSD(2)T modification [143], the linear-response CCSDR(3) method [46, 47], the EOMCCSD(T) [51], EOMCCSD(T̃) [52], and EOMCCSD(T0 ) [52] hierarchy obtained from the perturbative analysis of the EOMCCSDT equations, the CR-EOMCC family, such as the original CR-EOMCCSD(T) schemes [63, 65, 66], the newer CR-EOMCC(2,3) approaches [67, 68], and their approximately and rigorously size-intensive δ-CR-EOMCC counterparts [63, 69], as well as the related N-EOMCCSD(T) scheme [144], the spin-flip [145–147] extensions of variants A and D of CR-EOMCC(2,3) implemented in Ref. [148] and abbreviated as EOMCCSD(fT) and EOMCCSD(dT), respectively, and the iterative EOMCCSDT-n (n = 1, 2, 3) [51, 52] and CC3 [44–47] methodologies. These various approaches account for the leading triples effects, while replacing the iterative n3o n5u (N 8 ) CPU steps of full EOMCCSDT by the iterative (EOMCCSDT-n and CC3) or even less expensive non-iterative (EOMCC(2)PT(2), EOMCCSD(2)T , CCSDR(3), EOMCCSD(T), 12 0 EOMCCSD(T̃), EOMCCSD(T ), CR-EOMCCSD(T), N-EOMCCSD(T), CR-EOMCC(2,3)) n3o n4u (N 7 ) CPU costs. This, combined with their black-box character, makes the above approximate treatment of triple excitations within the EOMCC and linear-response CC frameworks attractive candidates for routine use in photochemistry and other areas where accurate excited-state information is called for. Another promising approach in this category is the similarity-transformed EOMCC (STEOMCC) methodology [24, 53, 54], which incorporates higher-than-double excitations into the EOMCC framework through the suitable transformation of the similarity-transformed Hamiltonian of EOMCCSD. The performance of many of the above methods is already well established. In particular, the EOMCCSDT-3 and CC3 approaches and their non-iterative EOMCCSD(T), EOMCCSD(T̃), and CCSDR(3) counterparts have been thoroughly examined in Refs. [17, 19–22, 24, 25] using the database developed in Ref. [17]. The same applies to the STEOMCC approaches, which have been tested against the EOMCCSDT-3 and CC3 data that are generally recognized as accurate, using the singlet and triplet excited states of the 28 molecules constituting the database of Ref. [17]. Although a number of successful applications of the CR-EOMCCSD(T), CR-EOMCC(2,3), δ-CR-EOMCCSD(T), and δ-CR-EOMCC(2,3) methods have been published (see, e.g., Refs. [63–69, 74, 111, 115, 135, 136, 149–186]), showing considerable promise in applications involving singly and more multi-reference doubly excited states, none of the CR-EOMCC approaches have been subjected to a comprehensive statistical evaluation of the type used in Refs. [17–25]. This present work addresses this issue by testing the approximately size-intensive δ-CR-EOMCCSD(T) method and its biorthogonal and strictly size-intensive δ-CR-EOMCC(2,3) counterpart against the previously published EOMCCSDT-3, CC3, CASPT2, and theoretical best estimate (TBE) data using the database of excited states developed in Ref. [17] and the subsequent studies [19–22, 24, 25]. Com13 parisons with the EOMCCSDT-3 and CC3 results reported in Refs. [17, 19, 21, 22, 24, 25] are particularly useful, since the δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) approaches replace the iterative n3o n4u (N 7 ) CPU steps and ∼ N 6 storage requirements by the much less expensive iterative n2o n4u (N 6 ) and non-iterative n3o n4u (N 7 ) steps and ∼ N 4 storage requirements. Comparisons with the CASPT2 and TBE data are useful too, since δ-CR-EOMCC methods are computational black boxes and the multi-reference CASPT2 approach is not. Furthermore, the existing TBE data can presently be regarded as some of the best estimates of the excitation energies for the molecules comprising the database of Ref. [17], which one would like to reproduce with a reasonable accuracy. 3.2 Theory As already alluded to in Section 3.1, in the δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) methods of Refs. [63] and [69], we correct the vertical excitation energies obtained in the EOMCCSD calculations for the leading triples effects extracted from the MMCC considerations. Thus, if ωµ = Eµ − E0 represents the vertical excitation energy from the ground state |Ψ0 i to the excited state |Ψµ i, the δ-CR-EOMCCSD(T) and δ-CR-EOMCC(2,3) values of ωµ can be given the following general form: (CCSD) ωµ = ωµ (CCSD) where ωµ + δµ , (3.1) is the EOMCCSD excitation energy obtained by diagonalizing the similarity- transformed Hamiltonian of CCSD, i.e., H̄ (CCSD) = e−T1 −T2 HeT1 +T2 , 14 (3.2) in the space of singly and doubly excited determinants, |Φai i and |Φab ij i, respectively, and (CCSD) δµ is the triples correction to ωµ . As usual, T1 and T2 are the singly and doubly excited components of the cluster operator T obtained by solving the ground-state CCSD equations and i, j, . . . (a, b, . . . ) are the spin-orbitals which are occupied (unoccupied) in the reference determinant |Φi (in the case of the δ-CR-EOMCC calculations reported in this dissertation, the restricted Hartree–Fock (RHF) determinant). Both δ-CR-EOMCC approaches use the same general expression for δµ , namely, δµ = X ijk `abc µ,ijk Mµ,abc , (3.3) i 0.75 eV CASPT2 EOMCCSD REL 6.62 7.42 1.22 2 1A ( g ! *) 2 1A ( g ! *) 5.42 6.61 1.23 Octatetraene 2 1A ( g ! *) 4.64 5.98 1.21 Benzene 2 1E2g( ! *) 8.18 9.21 1.17 Naphthalene 3 1Ag( ! *) 6.67 7.77 1.15 Naphthalene 3 1B3u( ! *) 7.74 9.03 1.18 Pyridine 3 1B2( ! *) 8.60 9.64 1.17 Pyrazine 1 1B3g( ! *) 8.47 9.74 1.18 Butadiene Hexatriene 12 R2 = 0.9783 EOMCCSD (eV) 10 8 2 1A ( g ! *) 8.61 9.54 1.16 3 1A ( 1 ! *) 7.21 7.97 1.09 1 1B ( 2 ! *) 6.31 7.09 1.06 1 1E (!" !*) 7.50 8.28 1.08 2 1E (!" !*) 8.99 10.28 1.16 2 1B (n 1g !*) 6.45 7.25 1.11 Tetrazine 3 1B (n 1g !*) 6.73 8.36 1.16 Tetrazine 2 1B3g( ! *) 8.34 9.43 1.17 5.64 6.55 1.12 Pyrazine Pyrimidine 6 Pyridazine Triazine 4 Triazine MaxE 1.63 2 MSE 0.49 Tetrazine MUE 0.50 0 0 2 4 6 8 10 12 CASPT2 (eV) Benzoquinone 1 1B3u(n !*) Formamide 3 1A ( ! *) 10.54 11.40 1.10 Thymine 3 1A (n! "*) 6.85 7.68 1.12 Adenine 6 1 A (!"!*) 6.87 7.72 1.10 7.56 8.48 1.12 Adenine 7 1A (!" !*) Figure 3.3: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 78 Outliers > 0.75 eV 12 R2 = 0.9830 -CR-(T),IA (eV) 10 8 State 6 -CR-(T),IA REL Octatetraene 2 1Ag( ! *) 4.64 5.41 1.21 Naphthalene 3 1B3u( ! *) 7.74 8.61 1.18 8.47 9.27 1.18 8.99 9.85 1.16 6.73 7.89 1.16 Pyrazine Triazine 4 CASPT2 Tetrazine 1 1B ( 3g 2 1E (!" 3 1B (n 1g ! *) !*) !*) MaxE MSE MUE 1.16 0.20 0.24 2 0 0 2 4 6 8 10 12 CASPT2 (eV) Figure 3.4: Correlation plot for the calculated singlet excited states: δ-CREOMCCSD(T),IA/TZVP (δ-CR-(T),IA) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 79 Outliers > 0.75 eV 12 R2 = 0.9842 -CR-(T),ID (eV) 10 8 State 6 CASPT2 -CR-(T),ID REL Naphthalene 3 1B3u( ! *) 7.74 8.50 1.18 Tetrazine 3 1B1g(n 6.73 7.68 1.16 !*) 4 MaxE MSE MUE 0.95 0.12 0.21 2 0 0 2 4 6 8 10 12 CASPT2 (eV) Figure 3.5: Correlation plot for the calculated singlet excited states: δ-CREOMCCSD(T),ID/TZVP (δ-CR-(T),ID) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 80 Outliers > 0.75 eV 12 R2 = 0.9863 -CR-(2,3),A (eV) 10 8 6 4 MaxE MSE MUE 0.68 0.00 0.17 2 0 0 2 4 6 8 10 12 CASPT2 (eV) Figure 3.6: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 81 Outliers > 0.75 eV 12 R2 = 0.9869 -CR-(2,3),D (eV) 10 8 6 4 MaxE MSE MUE 0.53 -0.06 0.18 2 0 0 2 4 6 8 10 12 CASPT2 (eV) Figure 3.7: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs CASPT2/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 82 Outliers > 0.75 eV 12 R2 = 0.9721 State EOMCCSD (eV) 10 8 6.55 7.42 1.22 Hexatriene 2 1Ag( ! *) 5.09 6.61 1.23 2 1A ( g ! *) 4.47 5.98 1.21 Cyclopentadiene 2 1A ( 1 ! *) 6.28 7.05 1.15 Benzene 2 1E ( 2g ! *) 8.15 9.21 1.17 1 1B ( 1g ! *) 5.75 6.53 1.11 3 1A ( g ! *) 6.49 7.77 1.15 3 1A ( 1 ! *) 7.18 7.94 1.07 Benzoquinone 1 1B (n 3u 5.55 6.55 1.12 Propanamide 2 1A ( ! *) 7.09 7.87 1.09 Uracil 3 1A (n! "*) 6.56 7.69 1.12 Naphthalene Pyridine MaxE MSE MUE 1.52 0.45 0.45 2 0 0 2 4 6 8 10 REL 2 1Ag( ! *) Naphthalene 4 EOMCCSD Butadiene Octatetraene 6 TBE-2 !*) 12 TBE-2 (eV) Figure 3.8: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 83 Outliers > 0.75 eV 12 R2 = 0.9764 -CR-(T),IA (eV) 10 8 State 6 -CR-(T),IA REL Hexatriene 2 1Ag( ! *) 5.09 6.03 1.23 Octatetraene 2 1Ag( ! *) 4.47 5.41 1.21 6.49 7.35 1.15 !*) 5.55 6.34 1.12 (n! "*) 6.56 7.40 1.12 Naphthalene Benzoquinone 4 TBE-2 Uracil 3 1A ( g 1 1B (n 3u 3 1A ! *) MaxE MSE MUE 0.94 0.18 0.22 2 0 0 2 4 6 8 10 12 TBE-2 (eV) Figure 3.9: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR(T),IA) vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 84 Outliers > 0.75 eV 12 R2 = 0.9769 -CR-(T),ID (eV) 10 8 State 6 Hexatriene Octatetraene Naphthalene 4 2 1Ag( ! *) TBE-2 -CR-(T),ID REL 5.09 5.91 1.23 2 1A ( g ! *) 4.47 5.32 1.21 3 1A ( g ! *) 6.49 7.25 1.15 MaxE MSE MUE 0.85 0.11 0.19 2 0 0 2 4 6 8 10 12 TBE-2 (eV) Figure 3.10: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR(T),ID) vs TBE-2 data. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 85 Outliers > 0.75 eV 12 R2 = 0.9794 -CR-(2,3),A (eV) 10 8 6 4 MaxE MSE MUE 0.70 0.00 0.19 2 0 0 2 4 6 8 10 12 TBE-2 (eV) Figure 3.11: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR(2,3),A) vs TBE-2 data. The absence of the table on the right hand side indicates that there are no outliers. 86 Outliers > 0.75 eV 12 R2 = 0.9797 -CR-(2,3),D (eV) 10 8 6 4 MaxE MSE MUE 0.58 -0.07 0.19 2 0 0 2 4 6 8 10 12 TBE-2 (eV) Figure 3.12: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries (δ-CR(2,3),D) vs TBE-2 data. The absence of the table on the right hand side indicates that there are no outliers. 87 Outliers > 0.50 eV State 12 R2 = 0.9835 EOMCCSD (eV) 10 8 6.77 7.42 1.22 Hexatriene 2 1Ag( ! *) 5.72 6.61 1.23 Octatetraene 2 1Ag( ! *) 4.97 5.98 1.21 Benzene 2 1E2g( ! *) 8.43 9.21 1.17 Naphthalene 3 1Ag( ! *) 6.90 7.77 1.15 Naphthalene 3 1B3u( ! *) 8.12 9.03 1.18 4 1A ( 1 ! *) 8.68 9.44 1.14 3 1B ( 2 ! *) 8.77 9.64 1.17 1 1B ( 3g ! *) 8.77 9.74 1.18 2 1A ( g ! *) 8.69 9.54 1.16 2 1E (!" 9.44 10.28 1.16 2 1B (n 2g !*) 6.23 6.77 1.12 Tetrazine 3 1B (n 1g !*) 7.08 8.36 1.16 Tetrazine 2 1B3g( ! *) 8.47 9.43 1.17 Benzoquinone 1 1B3u(n !*) 5.82 6.55 1.12 Benzoquinone 2 1B1u ( ! *) 7.82 8.47 1.10 Uracil 3 1A (n! "*) 6.87 7.69 1.12 Uracil 4 1A (n! "*) 7.12 7.74 1.11 Pyrazine Pyrazine 4 Triazine Tetrazine MaxE MSE MUE 1.28 0.30 0.30 0 0 2 4 6 8 10 12 CC3 (eV) REL 2 1Ag( ! *) Pyridine 2 EOMCCSD Butadiene Pyridine 6 CC3 !*) Figure 3.13: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 88 Outliers > 0.50 eV 12 R2 = 0.9866 -CR-(T),IA (eV) 10 8 State 6 -CR-(T),IA REL 3 1B1g(n !*) 7.08 7.89 1.16 Benzoquinone 1 1B (n 3u !*) 5.82 6.34 1.12 Uracil 1A (n! "*) 6.87 7.40 1.12 Tetrazine 4 CC3 3 MaxE MSE MUE 0.81 0.01 0.13 2 0 0 2 4 6 8 10 12 CC3 (eV) Figure 3.14: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP (δ-CR-(T),IA) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 89 Outliers > 0.50 eV 12 R2 = 0.9879 -CR-(T),ID (eV) 10 8 6 State Tetrazine 3 1B (n 1g !*) CC3 -CR-(T),ID REL 7.08 7.68 1.16 4 MaxE MSE MUE 0.60 -0.06 0.15 2 0 0 2 4 6 8 10 12 CC3 (eV) Figure 3.15: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP (δ-CR-(T),ID) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 90 Outliers > 0.50 eV 12 R2 = 0.9902 -CR-(2,3),A (eV) 10 8 6 4 MaxE MSE MUE 0.44 -0.17 0.19 2 0 0 2 4 6 8 10 12 CC3 (eV) Figure 3.16: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 91 Outliers > 0.50 eV 12 R2 = 0.9906 -CR-(2,3),D (eV) 10 8 6 State Tetrazine 1 1B ( 2u ! *) CC3 -CR-(2,3),D REL 5.12 4.58 1.11 4 MaxE MSE MUE 0.54 -0.24 0.25 2 0 0 2 4 6 8 10 12 CC3 (eV) Figure 3.17: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs CC3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 92 Outliers > 0.50 eV 12 State R2 = 0.9860 EOMCCSD (eV) 10 8 2 1Ag( ! *) 6.89 7.42 1.22 Hexatriene 2 1Ag( ! *) 5.88 6.61 1.23 Octatetraene 2 1Ag( ! *) 5.17 5.98 1.21 8.60 9.21 1.17 Naphthalene Naphthalene Pyridine Pyridine 4 Pyrazine MaxE MSE MUE 0.93 0.22 0.22 2 0 0 2 4 6 8 10 12 REL Butadiene Benzene 6 EOMCCSDT-3 EOMCCSD 2 1E ( 2g 3 1A ( g ! *) 7.14 7.77 1.15 3 1B ( 3u ! *) 8.33 9.03 1.18 4 1A ( 1 ! *) 8.86 9.44 1.14 3 1B ( 2 ! *) 8.97 9.64 1.17 1 1B ( 3g ! *) 9.00 9.74 1.18 1A ( g ! *) Pyrazine 2 ! *) 8.90 9.54 1.16 Triazine 2 1E (!" !*) 9.64 10.28 1.16 Tetrazine 3 1B1g(n !*) 7.43 8.36 1.16 Tetrazine 2 1B3g( ! *) 8.72 9.43 1.17 EOMCCSDT-3 (eV) Figure 3.18: Correlation plot for the calculated singlet excited states: EOMCCSD/TZVP vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 93 Outliers > 0.50 eV 12 R2 = 0.9885 -CR-(T),IA (eV) 10 8 6 4 MaxE MSE MUE 0.46 -0.08 0.14 2 0 0 2 4 6 8 10 12 EOMCCSDT-3 (eV) Figure 3.19: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),IA/TZVP (δ-CR-(T),IA) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 94 Outliers > 0.50 eV 12 R2 = 0.9893 -CR-(T),ID (eV) 10 8 6 4 MaxE MSE MUE 0.37 -0.15 0.18 2 0 0 2 4 6 8 10 12 EOMCCSDT-3 (eV) Figure 3.20: Correlation plot for the calculated singlet excited states: δ-CR-EOM CCSD(T),ID/TZVP (δ-CR-(T),ID) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 95 Outliers > 0.50 eV 12 R2 = 0.9904 -CR-(2,3),A (eV) 10 8 6 4 MaxE MSE MUE 0.48 -0.26 0.26 2 0 0 2 4 6 8 10 12 EOMCCSDT-3 (eV) Figure 3.21: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),A/TZVP (δ-CR-(2,3),A) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The absence of the table on the right hand side indicates that there are no outliers. 96 Outliers > 0.50 eV 12 R2 = 0.9904 -CR-(2,3),D (eV) 10 8 State 6 Tetrazine 1 1B2u( ! *) EOMCCSDT-3 -CR-(2,3),D REL 5.16 4.58 1.11 4 MaxE MSE MUE 0.58 -0.33 0.33 2 0 0 2 4 6 8 10 12 EOMCCSDT-3 (eV) Figure 3.22: Correlation plot for the calculated singlet excited states: δ-CR-EOM CC(2,3),D/TZVP (δ-CR-(2,3),D) vs EOMCCSDT-3/TZVP vertical excitation energies (in eV) at MP2/6-31G∗ geometries. The table on the right hand side shows the list of outliers, marked with open circles in the plot. 97 variant A of δ-CR-EOMCC(2,3) offering a similar performance as long as the excited states of interest are not dominated by two-electron transitions, which is the case here. As already alluded to above, using the 1 1 B3g (n2 → π ∗2 ) state of s-tetrazine as an example, and as shown in our past work [67, 68], variants D of the CR-EOMCC(2,3) and δ-CR-EOMCC(2,3) methodologies are generally more robust than other CR-EOMCC/δ-CR-EOMCC corrections when the doubly excited states are considered. Most of the above remarks related to statistical comparisons of the δ-CR-EOMCC and CC3 excitation energies apply to the analogous comparisons of the δ-CR-EOMCC and EOMCCSDT-3 data. With the cutoff threshold of 0.50 eV used in the examination of the EOMCCSD and δ-CR-EOMCC vs CC3 and EOMCCSDT-3 results, we cannot say as much about the relative performance of the various non-iterative triples δ-CR-EOMCC corrections relative to EOMCCSDT-3 as in the case of the analogous comparisons with CC3, since all four corrections work equally well, producing no outliers, but we continue to observe a steady increase in the R2 correlation factor when going from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC(2,3),A, to δCR-EOMCC(2,3),D, from 0.9860 when the EOMCCSD and EOMCCSDT-3 excitation energies are compared to 0.9885, 0.9893, 0.9904, and 0.9904 when we compare the EOMCCSDT3 data with their δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC(2,3),A, and δ-CR-EOMCC(2,3),D counterparts, respectively. This is reflected in the error distribution curves shown in Figs. 3.23–3.26, where we can see that the error distribution relative to the EOMCCSDT-3 data becomes increasingly narrower as we go from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC(2,3),A, to δ-CR-EOMCC(2,3),D. Once again, based on the MUE and MSE values collected in Table 3.8, one might argue that δ-CR-EOMCCSD(T),IA and δ-CR-EOMCCSD(T),ID approaches 98 are more effective in reproducing the EOMCCSDT-3 results than their biorthogonal δ-CREOMCC(2,3),A and δ-CR-EOMCC(2,3),D counterparts, but knowing that the error distribution relative to the excitation energies obtained in the EOMCCSDT-3 calculations is narrowest in the δ-CR-EOMCC(2,3),D case and keeping in mind that the EOMCCSDT-3 method is not necessarily more accurate than the δ-CR-EOMCC(2,3),D approach, particularly when the excited states in question have some contributions from double excitations, reinforces our belief that the overall best approach among the four types of triples corrections to EOMCCSD excitation energies investigated in this work is δ-CR-EOMCC(2,3),D, with δ-CR-EOMCC(2,3),A offering similar accuracies as long as the excited states of interest are not dominated by two-electron transitions. 99 70 60 Population 50 40 30 20 10 0 -0.5 0 0.5 1 1.5 Excitation Energy Deviation (eV) Figure 3.23: Normal distribution curves for the deviation of the δ-CR-EOMCCSD(T),IA excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. 100 70 60 Population 50 40 30 20 10 0 -0.5 0 0.5 1 1.5 Excitation Energy Deviation (eV) Figure 3.24: Normal distribution curves for the deviation of the δ-CR-EOMCCSD(T),ID excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. 101 70 60 Population 50 40 30 20 10 0 -0.5 0 0.5 1 1.5 Excitation Energy Deviation (eV) Figure 3.25: Normal distribution curves for the deviation of the δ-CR-EOMCC(2,3),A excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. 102 70 60 Population 50 40 30 20 10 0 -0.5 0 0.5 1 1.5 Excitation Energy Deviation (eV) Figure 3.26: Normal distribution curves for the deviation of δ-CR-EOMCC(2,3),D excitation energies at MP2/6-31G∗ geometries from the CASPT2/TZVP (green), CC3/TZVP (black), EOMCCSDT-3/TZVP (red) and TBE-2 (blue) results. 103 3.6 Conclusions We have used a comprehensive test set of 148 singlet excited states of 28 medium-size organic molecules taken from Ref. [17] to benchmark two variants of the approximately size-intensive CR-EOMCC method with singles, doubles, and non-iterative triples, abbreviated as δ-CREOMCCSD(T),IA and δ-CR-EOMCCSD(T),ID [63], derived from the MMCC formalism [65, 68, 122, 125–129], and the analogous two variants of the rigorously size-intensive δ-CREOMCC(2,3) approach, designated as δ-CR-EOMCC(2,3),A and δ-CR-EOMCC(2,3),D, respectively [69], based on the generalization of the biorthogonal MMCC formalism [127–129] to excited states [67, 122–124], against the previously published CASPT2 [17, 21], TBE [17, 21] (especially, TBE-2 [21]), CC3 [17, 25] and EOMCCSDT-3 [22, 24] results. The list of 148 singlet excited states used to initiate the various statistical error analyses reported in this study has been obtained by excluding the doubly excited 1 1 B3g (n2 → π ∗2 ) state of s-tetrazine, for which the CC3 and EOMCCSDT-3 reference data are unavailable and the existing CASPT2 and TBE excitation energies reported in Refs. [17, 21] unreliable, from the 149 singlet excitations collected in Table I of Ref. [17]. We have, however, determined the δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC(2,3),A, and δ-CR-EOMCC(2,3),D vertical excitation energies, as well as their EOMCCSD counterparts, for all 149 excited states listed in Table I of Ref. [17], along with the 54 additional excitations. All of the δ-CR-EOMCC calculations performed in this work and the underlying EOMCCSD computations were performed using the TZVP basis set used in the previous studies [17– 20, 22–25] and two sets of the ground-state nuclear geometries, including the MP2/6-31G∗ geometries taken from Ref. [17] and the CR-CC(2,3),D geometries obtained in this work, were used to determine the EOMCCSD, δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, 104 δ-CR-EOMCC(2,3),A, and δ-CR-EOMCC(2,3),D vertical excitation spectra. By comparing our δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR- EOMCC(2,3),A, and δ-CR-EOMCC(2,3),D excitation energies for the 148 excited states of 28 molecules taken from Table I of Ref. [17], as described above, or their appropriate subsets for which the relevant reference CASPT2, TBE-2, CC3, and EOMCCSDT-3 data are available, we have shown that the non-iterative triples corrections to the EOMCCSD excitation energies defining the relatively inexpensive, single-reference, black-box δ-CR-EOMCC approaches provide significant improvements in the EOMCCSD data, while closely matching the results of the iterative and considerably more expensive CC3 and EOMCCSDT-3 calculations and their CASPT2 and TBE counterparts, typically to within ∼ 0.1 − 0.2 eV, i.e., to within intrinsic errors of the CC3, EOMCCSDT-3, CASPT2, and TBE estimates. We have also demonstrated that the δ-CR-EOMCC methods, especially the most robust δ-CR-EOMCC(2,3),D approach that works well for singly as well as doubly excited states, are capable of bringing the results of the CC3 and EOMCCSDT-3 calculations to a closer agreement with the CASPT2 and TBE data, demonstrating the utility of the cost effective δ-CR-EOMCC methods in applications involving molecular electronic spectra. This has allowed us to conclude that the overall best balanced approach among the four types of triples corrections to EOMCCSD excitation energies investigated in this work is δ-CREOMCC(2,3),D, with δ-CR-EOMCC(2,3),A offering similar accuracies as long as the excited states of interest are not dominated by two-electron transitions. We have reached these conclusions by performing a variety of full and partial statistical error analyses and examining the suitably designed correlation and error distribution plots. We have also used the four δ-CR-EOMCC approaches considered in this study to identify and accurately characterize 54 additional singlet excited states in the energy range covered by 105 Table I of Ref. [17], including five states that can be found in the Supporting Information to Ref. [17] and 49 states that have not been considered in the earlier benchmark work [17–25], which can be used in future benchmark studies. The aforementioned 1 1 B3g (n2 → π ∗2 ) state of s-tetrazine, listed in Table I of Ref. [17], which could not be used in our overall statistical error analyses due to the absence of the reliable benchmark data to judge our δ-CR-EOMCC results, and six other states among the 54 states outside the set of 149 states listed in Table I of Ref. [17] are almost pure two-electron transitions, which many quantum chemistry methods have problems with, but we have provided arguments, based on the successful track record involving various CR-EOMCC or δ-CR-EOMCC calculations, including quasi-degenerate excited states dominated by double excitations [63–68, 74, 111, 149, 152, 155, 162, 179, 182] and the comparison of our best δ-CR-EOMCC(2,3),D excitation energies for the 1 1 B3g (n2 → π ∗2 ) state of s-tetrazine with the recently published NEVPT2 data [23], that our δ-CREOMCC calculations for the doubly excited states found in this work and other additional states that have not been considered in the prior work [17–25] are accurate to within ∼ 0.2 − 0.3 eV. We have suggested full EOMCCSDT, active-space EOMCCSDt, or accurate multi-reference CI calculations for all of the additional excited states found in our calculations to verify if our assessment of the accuracy of the δ-CR-EOMCC calculations for these extra states is correct. In summary, we have identified the δ-CR-EOMCC(2,3) methodology, especially its δCR-EOMCC(2,3),D variant, as a useful and, at the same time, rigorously size-intensive approach for the routine and highly accurate calculations of molecular electronic spectra, even when the excited states of interest have more substantial two-electron contributions, with the δ-CR-EOMCC(2,3),A approximation offering an equally good description as long as the excited states of interest are dominated by one-electron transitions. 106 Chapter 4 Economical Doubly Electron-Attached Equation-of-Motion Coupled-Cluster Methods with an Active-Space Treatment of 3-particle–1-hole and 4-particle–2-hole Excitations The second part of this dissertation is concerned with the development and application of economical DEA-EOMCC approximations that have emerged following our group’s initial implementation of the full DEA-EOMCC(3p-1h) and DEA-EOMCC(4p-2h) methods and the active-space DEA-EOMCC(4p-2h) approach, in which 4p-2h terms are treated using active orbitals [92, 93]. We begin by reviewing fundamental elements of the DEAEOMCC theory defining the existing full DEA-EOMCC(3p-1h) and full and active-space DEA-EOMCC(4p-2h) approaches introduced prior to my method development work [92, 93]. Then, we discuss the new generation of DEA-EOMCC approaches truncated at either 4p-2h 107 or 3p-1h excitations, where both 3p-1h and 4p-2h terms are treated using active orbitals, resulting in major savings in the computational effort compared to their all-orbital counterparts. In addition to the discussion of the key equations, details of computer codes, and examples of CPU timings, we present the results of the various DEA-EOMCC calculations, including methods developed in this thesis project, focusing on the determination of electronic spectra of diradicals, especially their singlet–triplet gaps, and one example of single bond breaking where the DEA-EOMCC approaches can be useful too. Much of our discussion is tied up to our original work published in Refs. [94] and [117]. 4.1 Background Information and Motivation Quantum chemistry methods based on the exponential wave function ansatz [198, 199] of single-reference CC theory [34–37, 200, 201] and their extensions to excited states and properties other than energy exploiting the EOM [48–50, 118, 119] and linear response [38–42, 120, 121, 202] frameworks have witnessed considerable success in a wide range of molecular applications (cf., e.g., Refs. [203] and [204] for selected reviews). As pointed out in the Introduction, this includes extensions of the EOMCC formalism to open-shell systems obtained by adding electron(s) to or removing electron(s) from the corresponding closed-shell cores via the EA [75–81] or IP [59–62, 77, 79–85] methodologies, their linear response [205] and SAC-CI [96, 206–208] counterparts, and their multiply attached/ionized generalizations, such as the DEA- and DIP-EOMCC schemes [53, 86–94] or the EOMCC approach to triple electron attachment [209]. There is growing interest in the EA/IP-EOMCC, DEA/DIP-EOMCC, and similar approaches, as a way to handle ground and excited states of open-shell species around closed shells, such as radicals and diradicals, since they offer 108 several advantages over the conventional particle-conserving CC/EOMCC treatments that rely on the spin-integrated, but not spin-adapted, spin-orbital formulation employing the unrestricted or restricted open-shell references. The EA- and IP-EOMCC methods and their multiply electron-attached and multiply ionized extensions, in which one diagonalizes the similarity-transformed Hamiltonian of the CC theory obtained in the calculations for the reference closed-shell system in the appropriate sector of the Fock space corresponding to the open-shell species of interest, provide a rigorously spin-adapted description, while offering a potential of being very accurate when suitable approximations, including those developed in this work are applied. Recently, our group has developed high-level variants of the DEA- and DIP-EOMCC approaches with up to 4p-2h and 4h-2p excitations, abbreviated as DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p), respectively, and their less expensive active-space counterparts, designated as DEA-EOMCC(4p-2h){Nu } and DIP-EOMCC(4h-2p){No }, where Nu and No indicate the numbers of active unoccupied and active occupied orbitals used to select the corresponding 4p-2h and 4h-2p contributions, as promising new ways to describe multireference systems having two electrons outside the corresponding closed-shell cores [92, 93]. The active-space DEA-EOMCC(4p-2h){Nu } and DIP-EOMCC(4h-2p){No } approaches of Refs. [92] and [93] have been shown to be highly successful in challenging test cases involving single bond breaking in closed-shell molecules leading to doublet dissociation fragments and electronic spectra of diradicals, producing excellent results when compared to full CI or experiment, independent of the type of molecular orbitals (MOs) employed in the calculations, and almost perfectly reproducing the parent full DEA-EOMCC(4p-2h) and DIP-EOMCC(4h-2p) data at the fraction of the computer cost. Unfortunately, even with the help of active orbitals to select the dominant 4p-2h excitations, calculations at the DEA-EOMCC(4p-2h) 109 level remain quite expensive, especially when larger basis sets have to be employed. This has prompted our interest in investigating new types of approximations to the existing activespace DEA-EOMCC(4p-2h){Nu } and full DEA-EOMCC(4p-2h) methods [92, 93], with the goal of developing more economical and cost-effective approaches [94] that can reduce the computer effort further without sacrificing high accuracy the DEA-EOMCC(4p-2h)-level theories offer. This work reports my contributions to this area. To appreciate the new DEA-EOMCC methods proposed and tested in this work, and the advantages they offer compared to the existing schemes in this category, we first examine the computer costs of the high-level DEA-EOMCC(4p-2h) calculations. The most expensive CPU steps of the DEA-EOMCC(4p-2h) computations with a full treatment of 3p-1h and 4p-2h contributions scale as n2o n6u or N 8 (see Refs. [92] and [93]). As a result, the DEA-EOMCC(4p-2h) approach, although highly accurate in applications involving single bond breaking in closed-shell molecules and diradical electronic spectra [92, 93], is usually prohibitively expensive (for the examples of timings, see Section II.B of Ref. [92]). The active-space analog of the DEA-EOMCC(4p-2h) method, designated as DEAEOMCC(4p-2h){Nu }, which uses a subset of Nu unoccupied orbitals to select the dominant 4p-2h excitations, developed in Refs. [92] and [93] reduces the most expensive n2o n6u steps of the full DEA-EOMCC(4p-2h) approach to a considerably more manageable Nu2 n2o n4u or ∼ N 6 level, which is equivalent to costs of the standard EOMCCSD calculations or costs of the ground-state CCSD computations times a relatively small prefactor if Nu  nu , but does not solve the problem in its entirety. Indeed, although the resulting active- space DEA-EOMCC(4p-2h){Nu } approach offers substantial savings in the computer effort compared to its full DEA-EOMCC(4p-2h) parent without loss of accuracy [92, 93], the DEA-EOMCC(4p-2h){Nu } calculations remain expensive when larger basis sets are em110 ployed. This is due to the fact that in the existing implementation of the active-space DEA-EOMCC(4p-2h){Nu } method described in Refs. [92] and [93] the lower-rank 3p-1h components are still treated fully using all orbitals in the basis set, and this becomes a serious problem in applications using larger bases, since computer costs associated with 3p-1h contributions can be high too. Indeed, the full treatment of 3p-1h excitations within the DEA-EOMCC framework requires CPU steps that scale as no n5u , which can be as demanding as or, in some cases, more time consuming than the Nu2 n2o n4u steps of the DEAEOMCC(4p-2h){Nu } method associated with 4p-2h terms, especially when nu becomes larger, since typical values of Nu and no are much smaller than nu . Clearly, the same analysis applies to the lower-level DEA-EOMCC(3p-1h) approach [53, 86, 88, 89, 92, 93], in which the electron-attaching operator of the DEA-EOMCC formalism is truncated at 3p-1h component. The DEA-EOMCC(3p-1h) calculations, in which 4p-2h contributions are neglected, can be quite expensive too due to the no n5u steps resulting from a full treatment of 3p-1h excitations. There clearly is a need to address the above concerns if we want the DEA-EOMCC methodology, especially its presently highest DEA-EOMCC(4p-2h) level, to be more widely exploited in molecular applications. It is, therefore, essential that the expensive CPU steps of the no n5u and n2o n6u types, which originate from the presence of 3p-1h and 4p-2h components in the DEA-EOMCC wave function expansions, are replaced by steps that are considerably more manageable. The previous DEA-EOMCC(4p-2h){Nu } method described in Refs. [92] and [93] addresses this issue, but only partly, since it focuses on 4p-2h excitations without doing anything about their 3p-1h counterparts, which lead to high costs too if treated fully. In this dissertation, we present a solution to the above problems. Thus, we propose and test a new class of computationally affordable variants of the DEA-EOMCC approach, 111 abbreviated as DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC(3p-1h,4p-2h){Nu }. In analogy to the DEA-EOMCC(4p-2h){Nu } and DIP-EOMCC(4h-2p){No } methods of Refs. [92] and [93] and their active-space CC [98–109], EOMCC [71, 72, 110–112, 114], and EA/IPEOMCC [79–81, 96] predecessors (see Ref. [115] for a review), the DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC(3p-1h,4p-2h){Nu } schemes developed in this work are based on the idea of employing active orbitals to capture dominant excitation (in this case, electron attaching) amplitudes. Thus, the DEA-EOMCC(3p-1h){Nu } approach uses Nu active unoccupied orbitals to select a small subset of the dominant 3p-1h amplitudes within the standard DEA-EOMCC(3p-1h) framework, in which the operator attaching two electrons to the corresponding closed-shell core is truncated at 3p-1h component. Similarly, the DEAEOMCC(3p-1h,4p-2h){Nu } method uses Nu active unoccupied orbitals to select the dominant 3p-1h and 4p-2h amplitudes within the higher-level DEA-EOMCC(4p-2h) scheme. In other words, the DEA-EOMCC(3p-1h){Nu } approach is a natural approximation to its DEA-EOMCC(3p-1h) parent, which becomes full DEA-EOMCC(3p-1h) when all unoccupied orbitals in the basis set are active (i.e., Nu = nu ). Similarly, the DEA-EOMCC (3p-1h,4p-2h){Nu } approach is a natural approximation to its DEA-EOMCC(4p-2h) parent or to the previously developed DEA-EOMCC(4p-2h){Nu } method of Refs. [92] and [93], becoming full DEA-EOMCC(4p-2h) when Nu = nu . The DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC(3p-1h,4p-2h){Nu } approaches proposed in this work offer significant savings in the computer effort compared to their full DEA-EOMCC(3p-1h) and DEAEOMCC(4p-2h) counterparts by reducing the expensive N 6 -like no n5u steps associated with 3p-1h excitations to a N 5 -like Nu no n4u level, which for larger systems is less expensive than costs of the underlying CCSD calculations. As in the case of the previously proposed DEAEOMCC(4p-2h){Nu } approximation [92, 93], the DEA-EOMCC(3p-1h,4p-2h){Nu } method 112 replaces the N 8 -like n2o n6u steps associated with 4p-2h contributions by the much more affordable, CCSD-type, Nu2 n2o n4u steps, in addition to using the relatively inexpensive, N 5 -like, Nu no n4u steps in handling the lower-order 3p-1h terms. In order to test the performance of the DEA-EOMCC(3p-1h){Nu } and DEA-EOMCC (3p-1h,4p-2h){Nu } approaches developed in this work, especially when compared to the previously examined [92, 93] DEA-EOMCC(3p-1h), DEA-EOMCC(4p-2h){Nu }, and full DEA-EOMCC(4p-2h) methods, we investigate adiabatic excitation energies characterizing low-lying states of methylene, singlet–triplet gaps in trimethylenemethane (TMM), a series of cyclobutadiene and its derivatives and cyclopentadienyl cation, and bond breaking in the F2 molecule. We show that the new DEA-EOMCC(3p-1h,4p-2h){Nu } approach with the activespace treatment of 3p-1h and 4p-2h excitations and its lower-level DEA-EOMCC(3p-1h){Nu } counterpart ignoring 4p-2h contributions, while using active orbitals to select the dominant 3p-1h components, accurately reproduce the results obtained with the considerably more expensive parent DEA-EOMCC methods with a full treatment of 3p-1h and full or activespace treatment of 4p-2h excitations at the small fraction of the computer effort. This is particularly valuable when the higher-level DEA-EOMCC(3p-1h,4p-2h){Nu } approach is employed, since the explicit inclusion of 4p-2h excitations in the DEA-EOMCC wave function expansions leads to a robust and highly accurate description of the electronic structure and spectra of diradicals and single bond breaking in closed shells leading to doublet dissociation fragments, independent of the MO basis exploited in such calculations [92–94]. 113 4.2 4.2.1 Theory and Algorithmic Details Basic Elements of the DEA-EOMCC Formalism and the Previously Developed DEA-EOMCC(3p-1h), DEAEOMCC(4p-2h), and DEA-EOMCC(4p-2h){Nu } Approximations In the DEA-EOMCC formalism exploited in this work, one represents the ground (µ = 0) (N ) and excited (µ > 0) states |Ψµ i of the N -electron system obtained by adding two electrons to the closed-shell core using the following wave function ansatz [53, 86, 88, 89, 92–94]: (N ) (+2) |Ψµ i = Rµ (N −2) |Ψ0 i, (4.1) where (N −2) |Ψ0 i = eT |Φ(N −2) i (4.2) is the CC ground state of the (N −2)-electron closed-shell species, with |Φ(N −2) i designating the corresponding reference determinant that serves as the Fermi vacuum. The operator T entering Eq. (4.2) is the usual particle-conserving cluster operator, obtained in the groundstate CC calculations for the (N − 2)-electron reference system, and (+2) Rµ MR = X Rµ,np-(n−2)h , (4.3) n=2 where MR = N in the exact case and MR < N in approximate schemes, is the EOM operator attaching two electrons to the corresponding (N − 2)-electron closed-shell core, 114 while allowing excitations of the remaining electrons via the Rµ,np-(n−2)h components with n > 2. Once we decide on specific truncations in the many-body expansions representing (+2) the T and Rµ operators (assuming that the highest many-body rank in T (abbreviated as MT ) is, as explained in Refs. [92] and [93], at least (MR − 2)) and once the relevant components of T are determined by solving the ground-state CC equations for the (N − 2)(+2) electron reference system, we obtain the Rµ,np-(n−2)h components of the Rµ operator and the corresponding vertical electron-attachment energies (N ) ωµ (N ) where Eµ (N ) = Eµ (N −2) − E0 , (N ) (4.4) (N −2) is the energy of the N -electron state |Ψµ i and E0 is the ground-state CC energy of the (N − 2)-electron reference system, by solving the following non-Hermitian eigenvalue problem [92–94]: (+2) (H̄N,open Rµ (N ) (+2) )C |Φ(N −2) i = ωµ Rµ |Φ(N −2) i (4.5) in the space of N -electron determinants corresponding to the Rµ,np-(n−2)h components in(+2) cluded in Rµ . Here, H̄N,open is the open part of the similarity-transformed form of the Hamiltonian H, written in the normal-ordered representation HN = H−hΦ(N −2) |H|Φ(N −2) i, i.e., the open part of the H̄N = e−T HN eT = (HN eT )C (4.6) operator obtained in the underlying CC calculations for the (N −2)-electron reference system, and subscript C designates the connected operator product. Thus, H̄N,open is this part of H̄N , Eq. (4.6), which corresponds to diagrams of (HN eT )C that have external fermion lines. 115 (N −2) It is easy to show that H̄N,open is equivalent to H̄ − E0 1, where H̄ = e−T HeT is the similarity-transformed form of H and 1 is the unit operator. the aforementioned condition MR − 2 6 MT has to be satisfied to obtain the connected form of the eigenvalue problem represented by Eq. (4.5), which is, in turn, key to retaining the desired property of size (N ) intensivity of the resulting electron-attachment energies ωµ . (+2) Different truncations in the cluster and electron attaching operators, T and Rµ , respec- tively, which enter the above equations and which satisfy the condition MR − 2 6 MT , lead to various DEA-EOMCC schemes. In particular, in the full DEA-EOMCC(4p-2h) method developed in Refs. [92] and [93], which presently is the highest implemented level of the DEA-EOMCC theory, we truncate the cluster operator T at double excitations, i.e., we use the standard CCSD approach to determine the (N − 2)-electron reference ground state (N −2) |Ψ0 i, and set MR in Eq. (4.3) at 4, obtaining (+2) Rµ = Rµ,2p + Rµ,3p-1h + Rµ,4p-2h , (4.7) where Rµ,2p = X rab (µ) aa ab , (4.8) al,al,AF e>F e>F 134 and ef Iabce = h̄ab rcf . (4.64) The antisymmetrizer Apq = A pq , which enters the above equations is defined as Apq ≡ A pq = 1 − (pq), (4.65) with (pq) representing the transposition of indices p and q. Next, we consider some key components of the efficient computer implementation of the DEA-EOMCC(3p-1h,4p-2h){Nu } eigenvalue equations, namely, Eqs. (4.29) - (4.64). Figure (4.1) gives the important details of the algorithm that is used to compute the projection of the DEA-EOMCC eigenvalue problem on the selected 3p-1h determinants. The algorithm for calculating the remaining 2p and 4p-2h projections is similar but we do not discuss it in this dissertation. One important element of our algorithm is that the explicit loops that are used to construct the DEA-EOMCC(3p-1h,4p-2h){Nu } equations projected on |ΦAbck i, Eq. 4.37, range over active indices represented by the use of bold, uppercase letters for variables in the loop. Within these loops, we have utilized a high degree of vectorization while exploiting efficient matrix multiplication routines from the BLAS library to perform essential computations. In doing so, we have maintained the Einstein summation convention throughout and also imposed the traditional summation symbol P free index belongs to the active set of orbitals. 135 in instances where at least one unoccupied Calculate I˜am for all values of e, f Eq. (4.56) Set Iam = I˜am for all values of a, m LOOP OVER D ef n Calculate IDm = IDm + vmn rDef for all values of m, Eq. (4.54) X X Df n vmn raDf for all values of a, m, Eq. (4.55) Calculate Iam = Iam + 2 n f (>D) END OF LOOP OVER D LOOP OVER A Calculate χAbck and ∆Abck for all values b, c, k, Eqs. (4.41) and (4.42) Set ∆Abck = χAbck for all values of b, c, k LOOP OVER E kE m k k k + h̄E Calculate ∆Abc = ∆Abc A rEbc + h̄Am rEbc + X E m . Furthermore, amplitudes defining the Lµ deexcitation operator, `µ,i A 1 ...in comparison of Eq. (A.21) with Eq. (3.4) reveals the presence of the generalized moments of i ...i a ...a (A) 1 n (m ) = hΦ 1 n |H̄ (A) R the CC/EOMCC equations Mµ,a µ |Φi. It should be noted A i1 ...in 1 ...an i ...i 1 n (m ) with n > m that for a given CC/EOMCC method A, not all of the moments Mµ,a A A 1 ...an are non-zero. Indeed, for a given approximation, there is generally a value of n above which 184 i ...i 1 n (m ) are zero. For instance, in the case of CCSD (the m = 2 all moments Mµ,a A A 1 ...an case), only the triply, quadruply, pentuply, and hextuply excited moments, i.e. moments with n = 3 − 6, are nonzero when the Hamiltonian contains pairwise interactions only, and thus N0,A = 6 (the CCSD equations are solved by zeroing the singly and doubly excited ij moments, Mi0,a (2) and M0,ab (2), respectively, hence the triply excited moments are the first to be nonzero). Similarly, in the EOMCCSD case with µ > 0, Nµ,A = 8. Taking advantage of these observations, along with the fact that the moments are zero for n > Nµ,A , Eq. (A.21) can be rewritten as (A) δµ ≡ (A) Eµ − Eµ = Nµ,A X X n=mA +1 i1 <···