Analytical harmonic vibrational frequencies with VV10-containing density functionals: Theory, efficient implementation, and benchmark assessments.

VV10 is a powerful nonlocal density functional for long-range correlation that is used to include dispersion effects in many modern density functionals, such as the meta-generalized gradient approximation (mGGA), B97M-V, the hybrid GGA, ωB97X-V, and the hybrid mGGA, ωB97M-V. While energies and analytical gradients for VV10 are already widely available, this study reports the first derivation and efficient implementation of the analytical second derivatives of the VV10 energy. The additional compute cost of the VV10 contributions to analytical frequencies is shown to be small in all but the smallest basis sets for recommended grid sizes. This study also reports the assessment of VV10-containing functionals for predicting harmonic frequencies using the analytical second derivative code. The contribution of VV10 to simulating harmonic frequencies is shown to be small for small molecules but important for systems where weak interactions are important, such as water clusters. In the latter cases, B97M-V, ωB97M-V, and ωB97X-V perform very well. The convergence of frequencies with respect to the grid size and atomic orbital basis set size is studied, and recommendations are reported. Finally, scaling factors to allow comparison of scaled harmonic frequencies with experimental fundamental frequencies and to predict zero-point vibrational energy are presented for some recently developed functionals (including r2SCAN, B97M-V, ωB97X-V, M06-SX, and ωB97M-V).

Approaching the basis set limit in Gaussian-orbital-based periodic calculations with transferability: Performance of pure density functionals for simple semiconductors.

Simulating solids with quantum chemistry methods and Gaussian-type orbitals (GTOs) has been gaining popularity. Nonetheless, there are few systematic studies that assess the basis set incompleteness error (BSIE) in these GTO-based simulations over a variety of solids. In this work, we report a GTO-based implementation for solids and apply it to address the basis set convergence issue. We employ a simple strategy to generate large uncontracted (unc) GTO basis sets that we call the unc-def2-GTH sets. These basis sets exhibit systematic improvement toward the basis set limit as well as good transferability based on application to a total of 43 simple semiconductors. Most notably, we found the BSIE of unc-def2-QZVP-GTH to be smaller than 0.7 mEh per atom in total energies and 20 meV in bandgaps for all systems considered here. Using unc-def2-QZVP-GTH, we report bandgap benchmarks of a combinatorially designed meta-generalized gradient approximation (mGGA) functional, B97M-rV, and show that B97M-rV performs similarly (a root-mean-square-deviation of 1.18 eV) to other modern mGGA functionals, M06-L (1.26 eV), MN15-L (1.29 eV), and Strongly Constrained and Appropriately Normed (SCAN) (1.20 eV). This represents a clear improvement over older pure functionals such as local density approximation (1.71 eV) and Perdew-Burke-Ernzerhof (PBE) (1.49 eV), although all these mGGAs are still far from being quantitatively accurate. We also provide several cautionary notes on the use of our uncontracted bases and on future research on GTO basis set development for solids.

Transition states, reaction paths, and thermochemistry using the nuclear-electronic orbital analytic Hessian.

The nuclear-electronic orbital (NEO) method is a multicomponent quantum chemistry theory that describes electronic and nuclear quantum effects simultaneously while avoiding the Born-Oppenheimer approximation for certain nuclei. Typically specified hydrogen nuclei are treated quantum mechanically at the same level as the electrons, and the NEO potential energy surface depends on the classical nuclear coordinates. This approach includes nuclear quantum effects such as zero-point energy and nuclear delocalization directly into the potential energy surface. An extended NEO potential energy surface depending on the expectation values of the quantum nuclei incorporates coupling between the quantum and classical nuclei. Herein, theoretical methodology is developed to optimize and characterize stationary points on the standard or extended NEO potential energy surface, to generate the NEO minimum energy path from a transition state down to the corresponding reactant and product, and to compute thermochemical properties. For this purpose, the analytic coordinate Hessian is developed and implemented at the NEO Hartree-Fock level of theory. These NEO Hessians are used to study the SN2 reaction of ClCH3Cl- and the hydride transfer of C4H9 +. For each system, analysis of the single imaginary mode at the transition state and the intrinsic reaction coordinate along the minimum energy path identifies the dominant nuclear motions driving the chemical reaction. Visualization of the electronic and protonic orbitals along the minimum energy path illustrates the coupled electronic and protonic motions beyond the Born-Oppenheimer approximation. This work provides the foundation for applying the NEO approach at various correlated levels of theory to a wide range of chemical reactions.

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

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

Implementation of analytic gradients for CCSD and EOM-CCSD using Cholesky decomposition of the electron-repulsion integrals and their derivatives: Theory and benchmarks.

We present a general formulation of analytic nuclear gradients for the coupled-cluster with single and double substitution (CCSD) and equation-of-motion (EOM) CCSD energies computed using Cholesky decomposition (CD) representations of the electron repulsion integrals. By rewriting the correlated energy and response equations such that the storage of the largest four-index intermediates is eliminated, CD leads to a significant reduction in disk storage requirements, reduced I/O penalties, and an improved parallel performance. CD thus extends the scope of the systems that can be treated by (EOM-)CCSD methods, although analytic gradients in the framework of CD are needed to extend the applicability of (EOM-)CCSD methods in the context of geometry optimizations. This paper presents a formulation of analytic (EOM-)CCSD gradient within the CD framework and reports on the salient details of the corresponding implementation. The accuracy and the capabilities of analytic CD-based (EOM-)CCSD gradients are illustrated by benchmark calculations and several illustrative examples.

Linker-Dependent Singlet Fission in Tetracene Dimers.

Separation of triplet excitons produced by singlet fission is crucial for efficient application of singlet fission materials. While earlier works explored the first step of singlet fission, the formation of the correlated triplet pair state, the focus of recent studies has been on understanding the second step of singlet fission, the formation of independent triplets from the correlated pair state. We present the synthesis and excited-state dynamics of meta- and para-bis(ethynyltetracenyl)benzene dimers that are analogues to the ortho-bis(ethynyltetracenyl)benzene dimer reported by our groups previously. A comparison of the excited-state properties of these dimers allows us to investigate the effects of electronic conjugation and coupling on singlet fission between the ethynyltetracene units within a dimer. In the para isomer, in which the two chromophores are conjugated, the singlet exciton yields the correlated triplet pair state, from which the triplet excitons can decouple via molecular rotations. In contrast, the meta isomer in which the two chromophores are cross-coupled predominantly relaxes via radiative decay. We also report the synthesis and excited-state dynamics of two para dimers with different bridging units joining the ethynyltetracenes. The rate of singlet fission is found to be faster in the dimer with the bridging unit that has orbitals closer in energy to that of the ethynyltetracene chromophores.

Singlet Fission in a Covalently Linked Cofacial Alkynyltetracene Dimer.

Singlet fission is a process in which a singlet exciton converts into two triplet excitons. To investigate this phenomenon, we synthesized two covalently linked 5-ethynyl-tetracene (ET) dimers with differing degrees of intertetracene overlap: BET-X, with large, cofacial overlap of tetracene π-orbitals, and BET-B, with twisted arrangement between tetracenes exhibits less overlap between the tetracene π-orbitals. The two compounds were crystallographically characterized and studied by absorption and emission spectroscopy in solution, in PMMA and neat thin films. The results show that singlet fission occurs within 1 ps in an amorphous thin film of BET-B with high efficiency (triplet yield: 154%). In solution and the PMMA matrix the S1 of BET-B relaxes to a correlated triplet pair (1)(T1T1) on a time scale of 2 ps, which decays to the ground state without forming separated triplets, suggesting that triplet energy transfer from (1)(T1T1) to a nearby chromophore is essential for producing free triplets. In support of this hypothesis, selective excitation of BET-B doped into a thin film of diphenyltetracene (DPT) leads to formation of the (1)(T1T1) state of BET-B, followed by generation of both DPT and BET-B triplets. For the structurally cofacial BET-X, an intermediate forms in <180 fs and returns to the ground state more rapidly than BET-B. First-principles calculations predict a 2 orders of magnitude faster rate of singlet fission to the (1)(T1T1) state in BET-B relative to that of crystalline tetracene, attributing the rate increase to greater coupling between the S1 and (1)(T1T1) states and favorable energetics for formation of the separated triplets.

On couplings and excimers: lessons from studies of singlet fission in covalently linked tetracene dimers.

Electronic factors controlling singlet fission (SF) rates are investigated in covalently linked dimers of tetracene. Using covalent linkers, relative orientation of the individual chromophores can be controlled, maximizing the rates of SF. Structures with coplanar and staggered arrangements of tetracene moieties are considered. The electronic structure calculations and three-state kinetic model for SF rates provide explanations for experimentally observed low SF yields in coplanar dimers and efficient SF in staggered dimers. The calculations illuminate the role of the excimer formation in SF process. The structural relaxation in the S1 state leads to the increased rate of the multi-exciton (ME) state formation, but impedes the second step, separation of the ME state into independent triplets. The slower second step reduces SF yield by allowing other processes, such as radiationless relaxation, to compete with triplet generation. The calculations of electronic couplings also suggest an increased rate of radiationless relaxation at the excimer geometries. Thus, the excimer serves as a trap of the ME state. The effect of covalent linkers on the electronic factors and SF rates is investigated. In all considered structures, the presence of the linker leads to larger couplings, however, the effect on the overall rate is less straightforward, since the linkers generally result in less favorable energetics. This complex behavior once again illustrates the importance of integrative approaches that evaluate the overall rate, rather than focusing on specific electronic factors such as energies or couplings.

Quantifying charge resonance and multiexciton character in coupled chromophores by charge and spin cumulant analysis.

We extend excited-state structural analysis to quantify the charge-resonance and multi-exciton character in wave functions of weakly interacting chromophores such as molecular dimers. The approach employs charge and spin cumulants which describe inter-fragment electronic correlations in molecular complexes. We introduce indexes corresponding to the weights of local, charge resonance, and biexciton (with different spin structure) configurations that can be computed for general wave functions thus allowing one to quantify the character of doubly excited states. The utility of the approach is illustrated by applications to several small dimers, e.g., He-H2, (H2)2, and (C2H4)2, using full and restricted configuration interaction schemes. In addition, we present calculations for several systems relevant to singlet fission, such as tetracene, 1,6-diphenyl-1,3,5-hexatriene, and 1,3-diphenylisobenzofuran dimers.

What we can learn from the norms of one-particle density matrices, and what we can't: some results for interstate properties in model singlet fission systems.

The utility of the norms of one-particle density matrices, ∥γ∥, for understanding the trends in electronic properties is discussed. Using several model systems that are relevant in the context of singlet fission (butadiene, octatetraene, and ethylene dimer), the dependence of interstate properties (such as transition dipole moments and nonadiabatic couplings, NACs) on molecular geometries is investigated. ∥γ∥ contains the principal information about the changes in electronic states involved, such as varying degree of one-electron character of the transition; thus, it captures leading trends in one-electron interstate properties (i.e., when ∥γ∥ is small, the respective interstate matrix elements are also small). However, finer variations in properties that arise due to the dependence of the matrix elements of the respective operators may not be reproduced. Analysis of NACs in ethylene dimer reveals that intermolecular components of NACs follow the trends in ∥γ∥ well, as they are determined primarily by the characters of the two wave functions; however, intramolecular components depend on the relative orientation of the two moieties via the dependence in the derivative of the electron-nuclear Coulomb operator. Therefore, intramolecular NACs may exhibit large variations even when the changes in ∥γ∥ are small. We observe large NACs at perfectly stacked geometry; however, larger values (by a factor of 1.6) are observed at slip-stacked (along the long axis) geometries. Larger values of NACs at slip-stacked configurations are due to the breaking of symmetry of the local environment of the heavy atoms and not due to the wave function composition. We found that the variations in ∥γ∥ for ethylene dimer are due to a varying admixture of the charge-resonance configurations in the S1 state, whereas the (1)ME state retains its pure multiexciton character.

General implementation of the resolution-of-the-identity and Cholesky representations of electron repulsion integrals within coupled-cluster and equation-of-motion methods: theory and benchmarks.

We present a general implementation of the resolution-of-the-identity (RI) and Cholesky decomposition (CD) representations of electron repulsion integrals within the coupled-cluster with single and double substitutions (CCSD) and equation-of-motion (EOM) family of methods. The CCSD and EOM-CCSD equations are rewritten to eliminate the storage of the largest four-index intermediates leading to a significant reduction in disk storage requirements, reduced I/O penalties, and, as a result, improved parallel performance. In CCSD, the number of rate-determining contractions is also reduced; however, in EOM the number of operations is increased because the transformed integrals, which are computed once in the canonical implementation, need to be reassembled at each Davidson iteration. Nevertheless, for large jobs the effect of the increased number of rate-determining contractions is surpassed by the significantly reduced memory and disk usage leading to a considerable speed-up. Overall, for medium-size examples, RI/CD CCSD calculations are approximately 40% faster compared with the canonical implementation, whereas timings of EOM calculations are reduced by a factor of two. More significant speed-ups are obtained in larger bases, i.e., more than a two-fold speed-up for CCSD and almost five-fold speed-up for EOM-EE-CCSD in cc-pVTZ. Even more considerable speedups (6-7-fold) are achieved by combining RI/CD with the frozen natural orbitals approach. The numeric accuracy of RI/CD approaches is benchmarked with an emphasis on energy differences. Errors in EOM excitation, ionization, and electron-attachment energies are less than 0.001 eV with typical RI bases and with a 10(-4) threshold in CD. Errors with 10(-2) and 10(-3) thresholds, which afford more significant computational savings, are less than 0.04 and 0.008 eV, respectively.

Analytic Evaluation of Nonadiabatic Couplings within the Complex Absorbing Potential Equation-of-Motion Coupled-Cluster Method.

We present the theory for the evaluation of nonadiabatic couplings (NACs) involving resonance states within the complex absorbing potential equation-of-motion coupled-cluster (CAP-EOM-CC) framework implemented within the singles and doubles approximation. Resonance states are embedded in the continuum and undergo rapid decay through autodetachment. In addition, nuclear motion can facilitate transitions between different resonances and between resonances and bound states. These nonadiabatic transitions affect the chemical fate of resonances and have distinct spectroscopic signatures. The NAC vector is a central quantity needed to model such effects. In the CAP-EOM-CC framework, resonance states are treated on the same footing as bound states. Using the example of fumaronitrile, which supports a bound radical anion and several anionic resonances, we analyze the NAC between bound states and pseudocontinuum states, between bound states and resonances, and between two resonances. We find that the NAC between a bound state and a resonance is nearly independent of the CAP strength and thus straightforward to evaluate, whereas the NAC between two resonance states or between a bound state and a pseudocontinuum state is more difficult to evaluate.

Revisiting the Performance of Time-Dependent Density Functional Theory for Electronic Excitations: Assessment of 43 Popular and Recently Developed Functionals from Rungs One to Four.

In this paper, the performance of more than 40 popular or recently developed density functionals is assessed for the calculation of 463 vertical excitation energies against the large and accurate QuestDB benchmark set. For this purpose, the Tamm-Dancoff approximation offers a good balance between computational efficiency and accuracy. The functionals ωB97X-D and BMK are found to offer the best performance overall with a root-mean square error (RMSE) of around 0.27 eV, better than the computationally more demanding CIS(D) wave function method with a RMSE of 0.36 eV. The results also suggest that Jacob's ladder still holds for time-dependent density functional theory excitation energies, though hybrid meta generalized-gradient approximations (meta-GGAs) are not generally better than hybrid GGAs. Effects of basis set convergence, gauge invariance correction to meta-GGAs, and nonlocal correlation (VV10) are also studied, and practical basis set recommendations are provided.

Faster Exact Exchange for Solids via occ-RI-K: Application to Combinatorially Optimized Range-Separated Hybrid Functionals for Simple Solids with Pseudopotentials Near the Basis Set Limit.

In this work, we developed and showcased the occ-RI-K algorithm to compute the exact exchange contribution in density functional calculations of solids near the basis set limit. Within the Gaussian planewave (GPW) density fitting, our algorithm achieves a 1-2 orders of magnitude speedup compared to conventional GPW algorithms. Since our algorithm is well suited for simulations with large basis sets, we applied it to 12 hybrid density functionals with pseudopotentials and a large uncontracted basis set to assess their performance on band gaps of 25 simple solids near the basis set limit. The largest calculation performed in this work involves 16 electrons and 350 basis functions in the unit cell utilizing a 6 × 6 × 6 k-mesh. With 20-27% exact exchange, global hybrid functionals (B3LYP, PBE0, revPBE0, B97-3, SCAN0) perform similarly with a root-mean-square deviation (RMSD) of 0.61-0.77 eV, while other global hybrid functionals such as M06-2X (2.02 eV) and MN15 (1.05 eV) show higher RMSD due to their increased fraction of exact exchange. A short-range hybrid functional, HSE achieves a similar RMSD (0.76 eV) but shows a notable underestimation of band gaps due to the complete lack of long-range exchange. We found that two combinatorially optimized range-separated hybrid functionals, ωB97X-rV (3.94 eV) and ωB97M-rV (3.40 eV), and the two other range-separated hybrid functionals, CAM-B3LYP (2.41 eV) and CAM-QTP01 (4.16 eV), significantly overestimate the band gap because of their high fraction of long-range exact exchange. Given the failure of ωB97X-rV and ωB97M-rV, we have yet to find a density functional that offers consistent performance for both molecules and solids. Our algorithm development and density functional assessment will serve as a stepping stone toward developing more accurate hybrid functionals and applying them to practical applications.

Why iron? A spin-polarized conceptual density functional theory study on metal-binding specificity of porphyrin.

Heme is a key cofactor of hemoproteins in which porphyrin is often found to be preferentially metalated by the iron cation. In our previous work [Feng, X. T.; Yu, J. G.; Lei, M.; Fang, W. H.; Liu, S. B. J. Phys. Chem. B 2009, 113, 13381], conceptual density functional theory (CDFT) descriptors have been applied to understand the metal-binding specificity of porphyrin. We found that the iron-porphyrin complex significantly differs in many aspects from porphyrin complexes with other metal cations except Ru, for which similar behaviors for the reactivity descriptors were discovered. In this study, we employ the spin-polarized version of CDFT to investigate the reactivity for a series of (pyridine)(n)-M(ll)-porphyrin complexes-where M = Mg, Ca, Cr, Mn, Co, Ni, Cu, Zn, Ru, and Cd, and n = 0, 1, and 2-to further appreciate the metal-binding specificity of porphyrin. Both global and local descriptors were examined within this framework. We found that, within the spin resolution, not only chemical reactivity descriptors from CDFT of the iron complex are markedly different from that of other metal complexes, but we also discovered substantial differences in reactivity descriptors between Fe and Ru complexes. These results confirm that spin properties play a highly important role in physiological functions of hemoproteins. Quantitative reactivity relationships have been revealed between global and local spin-polarized reactivity descriptors. These results contribute to our better understanding of the metal binding specificity and reactivity for heme-containing enzymes and other metalloproteins alike.

Toward understanding metal-binding specificity of porphyrin: a conceptual density functional theory study.

Porphyrin is a key cofactor of hemoproteins. The complexes it forms with divalent metal cations such as Fe, Mg, and Mn compose an important category of compounds in biological systems, serving as a reaction center for a number of essential life processes. Employing density functional theory (DFT) and conceptual DFT approaches, the structural properties and reactivity of (pyridine)(n)-M-porphyrin complexes were systematically studied for the following selection of divalent metal cations: Mg, Ca, Cr, Mn, Co, Ni, Cu, Zn, Ru, and Cd with n varying from 0, 1, to 2. Metal selectivity and porphyrin specificity were investigated from the perspective of both structural and reactivity properties. Quantitative structural and reactivity relationships have been discovered between bonding interactions, charge distributions, and DFT chemical reactivity descriptors. These results are beneficial to our understanding of the chemical reactivity and metal cation specificity for heme-containing enzymes and other metalloproteins alike.

New and Efficient Equation-of-Motion Coupled-Cluster Framework for Core-Excited and Core-Ionized States.

We present a fully analytical implementation of the core-valence separation (CVS) scheme for the equation-of-motion (EOM) coupled-cluster singles and doubles (CCSD) method for calculations of core-level states. Inspired by the CVS idea as originally formulated by Cederbaum, Domcke, and Schirmer, pure valence excitations are excluded from the EOM target space and the frozen-core approximation is imposed on the reference-state amplitudes and multipliers. This yields an efficient, robust, practical, and numerically balanced EOM-CCSD framework for calculations of excitation and ionization energies as well as state and transition properties (e.g., spectral intensities, natural transition, and Dyson orbitals) from both the ground and excited states. The errors in absolute excitation/ionization energies relative to the experimental reference data are on the order of 0.2-3.0 eV, depending on the K-edge considered and on the basis set used, and the shifts are systematic for each edge. Compared to a previously proposed CVS scheme where CVS was applied as a posteriori projection only during the solution of the EOM eigenvalue equations, the new scheme is computationally cheaper. It also achieves better cancellation of errors, yielding similar spectral profiles but with absolute core excitation and ionization energies that are systematically closer to the corresponding experimental data. Among the presented results are calculations of transient-state X-ray absorption spectra, relevant for interpretation of UV-pump/X-ray probe experiments.