SQMBox: Interfacing a semiempirical integral library to modular ab initio electronic structure enables new semiempirical methods

Numerous methodological developments as well as the increased availability of high-performance computing resources have contributed to computational chemistry becoming an integral discipline among the modern chemical sciences.^1–3 A large number of different ab initio as well as semiempirical electronic structure codes are now in widespread use. A combination of these two electronic structure realms typically takes place in multi-level simulation workflows, e.g., to explore the chemical^4,5 and conformer space^6,7 or in free energy calculations.^8,9 While some quantum chemistry programs are “general purpose”, many software packages have a certain application-specific focus with respect to the implemented methods, and a semiempirical equivalent to an ab initio method is generally not available in the same software package. Although many quantum chemistry software packages include functionality for both ab initio and semiempirical Hamiltonians,^10–18 these packages typically provide support for semiempirical methods using separate code paths with restricted functionality. Hence, if an electronic structure method is intended to be used in combination with a different Hamiltonian, such as a semiempirical Hamiltonian, a reimplementation of the respective electronic structure routine is typically required. A unified implementation would have numerous advantages: (a) it reduces the implementation effort to develop semiempirical variants of ab initio methods; (b) multi-level electronic structure workflows become naturally accessible since the identical code path guarantees the same functionality for ab initio and semiempirical Hamiltonians; and (c) due to their lower cost, semiempirical Hamiltonians can be used to efficiently test new electronic structure implementations (e.g., when analytic and numerical gradients are compared).

In this work, we present the implementation of a novel semiempirical integral library into the GPU-accelerated electronic structure code TeraChem. We will explain in detail how the modular setup of TeraChem and the new SQMBox library enable the straightforward extension to semiempirical Hamiltonians. The benefits of this approach are demonstrated by combining TeraChem’s spin-restricted ensemble-referenced Kohn–Sham (REKS) method with the semiempirical tight-binding Hamiltonian GFN1-xTB. For complete compatibility with the ab initio code, we add Mulliken-approximated¹⁹ Fock exchange to this semiempirical library. The inclusion of this term was previously reported to increase the computational cost of semiempirical methods by almost an order of magnitude.²⁰ We show how to efficiently compute semiempirical exchange by utilizing an earlier presented refactorization,²¹ reuse of stored intermediates, and GPU acceleration. We also generalize the exchange algorithm to encompass non-identical and non-symmetric density matrices and present formulas for the efficient evaluation of these terms in the supplementary material.

We first outline the key implementation aspects of electronic structure methods in TeraChem, followed by a recap of the extended tight-binding methods GFN1-xTB and GFN2-xTB. We then describe in more detail how we develop the Mulliken approximation-based semiempirical integrals by decomposing the individual xTB Hamiltonian elements according to their order with respect to the density matrix. Finally, we demonstrate the performance of our GPU-accelerated Mulliken-approximated Fock exchange implementation and show exploratory results for the GFN1-xTB Hamiltonian in combination with REKS and complete active space (CAS) electronic structure methods, demonstrating that gradients and nonadiabatic coupling elements are readily available through the newly established setup.

In a seminal paper by Almlöf and co-workers, the superiority of an integral-direct scheme for the electron repulsion integrals (ERIs) in self-consistent field (SCF) procedures was demonstrated.²² Here, the fourth-order tensors (μν|κλ) are recomputed whenever required and contracted immediately with the density matrix. This is in contrast to “conventional” integral handling, where the integrals are precomputed and stored. Due to its lower memory requirements, the integral-direct scheme is applicable to large molecular systems. Thus, it has become the accepted standard for the computation of ERIs and is implemented in the majority of ab initio electronic structure codes. However, most programs, specifically in post-SCF calculations, further benefit from efficient integral approximations that are performed in a somewhat semi-direct scheme. Popular variants include the resolution-of-the-identity (RI)^23,24 and Cholesky decomposition (CD)²⁵ methods. Both schemes factorize the four-index tensor as follows:

Here, μ, ν, κ, and λ denote the atomic orbitals, while P is an auxiliary index corresponding to the auxiliary or Cholesky basis. The necessary contractions with the atomic orbital (AO) density matrix (or orbital coefficients^24,26) are then carried out on these three-index tensors. This typically reduces the formal scaling of the calculation. For optimum performance, a lower dimensional intermediate (i.e., a second- or third-order tensor) can be precomputed and stored. Such integral handling schemes have found particular use in post-SCF calculations, where they reduce the formal scaling of ERI calculations on the molecular orbital (MO) basis. Consequently, the integral handling in many quantum chemistry codes depends, at least in part, on the electronic structure calculation to be performed.

TeraChem, a GPU-accelerated ab initio electronic structure code, approaches this problem differently.^27,28 In TeraChem, ERIs are always evaluated in a fully AO integral-direct manner in both SCF and post-SCF methods. This approach is rooted in the availability of a very efficient AO integral library (called “IntBox,” see Fig. 1). Its efficiency stems from a highly optimized GPU code that exploits primitive pair screening as well as the architecture specifics of common consumer-grade GPUs. As a consequence, every electronic structure method in TeraChem is formulated in terms of the basic quantities (i.e., Fock matrix intermediates such as Coulomb and exchange operators) that are provided by the IntBox library. Implemented electronic structure methods directly “inherit” excellent performance by leveraging this integral library. Effectively, this approach creates a clean separation of the electronic Schrödinger equation. While the wavefunction ansatz is encoded in the electronic structure routines inside the main TeraChem program, the Hamiltonian represented in the Gaussian-type orbital (GTO) basis is evaluated in the IntBox library.

In general, the energy of a single-determinant electronic wavefunction expanded in an orbital basis can be formulated in terms of energy contributions that depend on the one-particle density matrix D,

Here, σ, τ = α, β denote the electron spin and μ, ν, κ, and λ denote the contracted GTOs. The prefactors c_K and c_XC are method-specific scaling factors of the Fock exchange and the density-functional exchange-correlation energy. For the ab initio methods implemented in TeraChem, the IntBox library provides the Fock matrix terms given in brackets in the last line, as well as the nuclear gradients for the corresponding energy terms given in the first line (the XC term is handled separately in the TeraChem code and is not part of the IntBox library due to the grid dependence and nonlinear potential in the density). The same intermediates found in the last line of Eq. (2) can also be used to compute multi-determinant wavefunctions. The most naive way would be to sum up the respective terms for individual determinant contributions, but in practice, it is often convenient to generate the two-electron integral tensor in the molecular orbital (MO) basis from transformed J and K intermediates.²⁹ 

Originally, the electronic structure code made direct calls to the ab initio integral library. To achieve a more modular workflow in which the electronic structure procedure is truly agnostic to the explicit form of the underlying Hamiltonian, we have now set-up a wrapper routine (“IntWrapper”), which communicates with the electronic structure code and the respective integral library at the backend. Schematically, this workflow is shown in Fig. 1. The electronic structure code passes a transition or particle density matrix to IntWrapper and obtains a basic quantity (e.g., a J matrix) in return. Depending on the run mode, the received basic quantity has been computed at the backend by the ab initio (IntBox) or the semiempirical (SQMBox) integral library. In this way, unmodified electronic structure routines can carry out ab initio and semiempirical electronic structure methods using the same code path. Once the respective quantities [cf. the last line in Eq. (2) and gradients] for a semiempirical Hamiltonian are available and coded up, the entire functionality of the existing electronic structure code becomes available. This is the approach we follow in this work, and the details are given in Sec. II B.

We developed a semiempirical integral library (SQMBox), written in C/C++ and CUDA, and interfaced with the quantum chemistry software TeraChem. It is interfaced to the main code via a C interface, i.e., only standard C variable types and pointers are passed between the main program and the library. These design decisions ensure that SQMBox can be easily interfaced to other quantum chemistry codes. We take the semiempirical tight-binding Hamiltonians GFN1-xTB (originally GFN-xTB) and GFN2-xTB as blueprints to develop the library. In contrast to other popular semiempirical methods (see Ref. 30 for an overview), these tight-binding methods are formulated in terms of GTOs, allowing us to leverage existing routines (the overlap matrix and its gradient) from TeraChem. The parameters for these Hamiltonians are currently hard-coded in the SQMBox library, but user-specified parameters could be read from an external parameter file. The starting point for SQMBox development is the extended tight-binding energy expressions previously given,^31,32 which we briefly review here. The energy expression for either of these Hamiltonians can be formulated in a concise way as follows:

Here, n = 1, 2 refer to the GFN1-xTB and GFN2-xTB Hamiltonians, respectively. The first term in Eq. (3) refers to a purely classical short-range repulsion function. The dispersion (E_disp) is evaluated by means of the D3 model^33,34 in GFN1-xTB and by means of a self-consistent D4 model^35,36 in GFN2-xTB. The covalent bond formation is mainly driven by the third term, which is the extended Hückel-type (EHT) energy. In GFN1- and GFN2-xTB, E_EHT carries most of the element-specific parameters. The last term describes the electrostatic and exchange-correlation effects in both xTB methods. Most of the terms in Eq. (3) are similar in form for GFN1-xTB and GFN2-xTB. However, major differences are found in E_ES+XC. In GFN2-xTB, this takes the following form:

Here, the shell l of an atom denotes the summed quantity (here the electron population) of the atomic orbitals κ of a specific principal quantum number and angular momentum. The Δq_l are the shell charges and are given as the differences between the actual and reference electron population of the shell,

The first term in Eq. (4) refers to the isotropic electrostatic interaction energy, often denoted as E_γ, which results from the second-order tight-binding expansion in terms of density fluctuations. The sum goes over all shells in the system, which are distinguished by their atom index as well as their principal and angular momentum quantum numbers. The second term in Eq. (4) originates from the third-order term in the TB expansion. In GFN2-XTB, this is restricted to same-shell contributions. In GFN1-xTB, this term is evaluated in terms of atomic partitioning. The last five terms are absent in GFN1-xTB and refer to the multipole electrostatic terms and XC terms, which are considered up to second order in GFN2-xTB. All these terms are based on cumulative atomic multipole moments,³⁷ which are obtained analogously to Eq. (5) using atom-centered multipole moment integrals instead of overlap integrals. We refer the interested reader to the original publication,³¹ and for the purpose of this paper, we simply mention that the dependence on the density matrix is the same in the isotropic and anisotropic terms. Hence, we restrict our discussion to the isotropic Mulliken-approximated integral terms.

To obtain a semiempirical integral library, the individual energy terms in the xTB methods are grouped according to their dependence on the electron density matrix D. Next, these terms are juxtaposed to the respective ab initio terms (see Fig. 2). The energy terms with no dependence on D already in the tight-binding formulation are the repulsion and D3 dispersion. These are entirely classical and do not depend on the electronic density matrix. The energy terms with a linear dependence on D contribute to the electronic energy, but do not require a self-consistent solution. In the tight-binding energy expression, this directly applies to the extended Hückel-type energy, which can be handled equivalently to the kinetic energy in ab initio theories (see the top entry in Fig. 2).

The second- and third-order tight-binding terms, however, can be decomposed into terms that are of the same or lower order in the density matrix. We will first look at the second-order term: by inserting Eq. (5) into the first term on the right-hand side (rhs) of Eq. (4) and collecting all terms with no dependence on D, we obtain the reference density Coulomb interaction, which is purely classical and solely depending on q_l,0. This plain q_l,0 electrostatic term is treated like the nuclear repulsion energy in ab initio theories and denoted as $Ẽnn$ in the following (see ρ₀–ρ₀ interaction in Fig. 2). Furthermore, we obtain a term where the density matrix-dependent electronic charges q_l are interacting with the reference charges q_l,0. This term can be treated equivalently to the nuclear-electron interaction integrals in ab initio methods (denoted as $Ẽne$in the following). In case of GFN2-xTB, a multipole component (arising from $EESqμ$ and $EESqΘ$ due to q_l,0-dipole/quadrupole interactions) contributes to this term as well.

The electrostatic interactions between q_l correspond to a Mulliken-approximated interelectronic Coulomb interaction (J),

In our implementation, the Mataga–Nishimoto–Ohno–Klopman-type^38–40 Coulomb interaction terms are generalized to allow for range-separation as used in the context of exchange (see below),

R_AB is the interatomic distance and η_ll’ is the average of the effective chemical hardness of the two shells involved. Note that GFN1- and GFN2-xTB both use a shell-wise isotropic electrostatic treatment but use different averaging schemes for the hardness (harmonic and arithmetic, respectively). In Eq. (7), the c_FR and c_LR are the full- and long-range scaling coefficients and ω is a range-separation parameter analogous to range-separated hybrid density functional approximations. To the best of our knowledge, this is the first suggestion to introduce range separation in the context of empirical xTB treatment of the Coulomb interaction. The second term in Eq. (7) is motivated by the fact that the ω parameter is of inverse distance (bohr⁻¹) units. Hence, the second term in Eq. (7) can be regarded as a modified harmonic weighting and may be a better approximation than the error function scaling proposed in Ref. 21, as it still gives the maximum magnitude for the Coulomb interaction at zero distance.

The q_l,0-independent multipole electrostatic and the XC contributions, which occur in GFN2-xTB, are quadratic in D and also contribute to the interelectronic electrostatic energy $ẼJ$ in our implementation. Decomposing the third-order tight binding terms leads to energy terms that range from zeroth to third order in the electron density. The decomposed terms read as

Equation (8) is of zeroth order in the electron density and just an atomic (shell-wise), geometry-invariant constant. We add this contribution to the classical energy $Ẽnn$ to be in agreement with the energies obtained with the original xTB implementation.⁴¹ The first order term in Eq. (9) is added to the semiempirical equivalent of nuclear–electron interaction (⁠$Ẽne$⁠), and likewise, the second order term in Eq. (10) is added to $ẼJ$⁠. The last term, which is the third order in the electron density has no equivalent in ab initio wavefunction theory. We handle this term like the XC energy (with corresponding XC potential and kernel) in first principles Kohn–Sham DFT codes.

We mention at this point that Nishimoto has previously applied a similar decomposition procedure when deriving the TD-DFTB3 equations and also drew similar connections between the tight-binding and ab initio Hamiltonians.⁴² 

It should be emphasized that the semiempirical integral library presented here has its origin in semiempirical tight-binding methods. Hence, integrals are calculated in the nonorthogonal atomic orbital basis set, just like in the ab initio library. That means that the zero-differential overlap (ZDO) approximation⁴³ is not employed, and all Fock matrix quantities are used in the generalized eigenvalue problem,

All quantities and results presented in this paper refer to this non-ZDO treatment. However, the semiempirical library could be used in combination with the ZDO approximation by setting the overlap matrix to a unit matrix.

Meanwhile, the previously presented xTB Hamiltonians did not include any Fock exchange, the semiempirical integral library presented here does include an implementation of Mulliken-approximated Fock exchange. The matrix to be computed is as follows:

Our library is set up to use the same hardness parameters in the Mataga–Nishimoto–Ohno–Klopman-type term for Coulomb (J) and exchange (K) [see Eq. (7) for the definition of γ]. In principle, this provides a self-interaction error free treatment within the monopole approximation if Coulomb and exchange enter the Hamiltonian with identical scaling and range-separation parameters. As written in Eq. (13), the evaluation of the Mulliken-approximated Fock exchange would scale as O(N⁴). Initial implementations of long-range Fock exchange corrected DFTB were found to be considerably slower (by about a factor of 10).²⁰ This computational complexity is one of the reasons why this term is absent in most DFTB methods. However, it is often added when excited states are of interest.^21,44–46 Employing the well-known neglect of diatomic differential overlap (NDDO) approximation would bring down the scaling of this term to N², but that needs to be accompanied with the ZDO of the overlap matrix in Eq. (12), which is incompatible with most DFTB and GFNn-xTB Hamiltonians. Fortunately, the computational complexity of the Mulliken-approximated exchange can be significantly lowered, leading to a form with O(N³) scaling,

Here, the entire build for exchange contribution to the Fock matrix can be formulated in terms of cubically-scaling matrix–matrix multiplications and quadratically-scaling element-wise products (denoted by _°). This allows efficient evaluation of computing architectures like graphics processing units (GPUs). This rearrangement was pointed out earlier by Humeniuk and Mitrić,²¹ but no GPU implementation has been reported up to this point. We implemented the expression in Eq. (14) efficiently on the CPU as well as on the GPU (single and double precision), making use of existing BLAS and cuBLAS libraries. By storing intermediates, a total of six matrix–matrix multiplications need to be performed. If D is symmetric, i.e., in ground state SCF calculations, Eq. (14) can be further reduced to only four matrix–matrix multiplications, greatly reducing the computational effort. The computational efficiency of evaluating this term on the GPU will be demonstrated in Sec. IV.

The gradient expressions corresponding to the energy/potential terms listed in Fig. 2 and generalized to non-identical and non-symmetric density matrices X, Y, and Z are given as

Here, α_C is the Cartesian direction of nucleus C. The latter equation for the derivative of the density functional can be expressed as a derivative of an exchange-correlation kernel, which is evaluated for the reference density (Y in this case) and contracted with two other (transition) density matrices X and Z. The same routine can then be used for both ground and excited state cases (for the ground state, X = Y = Z = D). Essentially all electronic terms in the semiempirical Hamiltonian make use of the overlap matrix. For optimum efficiency in the gradient calculation, we want to avoid multiple evaluations of overlap gradients. Therefore, we split each of the gradient terms into two components: an overlap or basis function component and a Hamiltonian or Hellmann–Feynman-type component,

Here, we use the chain rule to differentiate the overlap matrix elements first and then factor out the nuclear derivative of the overlap matrix element. The second term on the right-hand side is the Hellmann–Feynman-type term and includes all partial derivatives with constant overlap. For the first term, this leads to a derivative contribution of the form

For each energy contribution to the electronic energy, we form the corresponding intermediate and store it in a matrix (denoted W^S here). In our implementation, these are added together and then added to the energy-weighted density matrix, which appears in the derivative expression of the orbital orthonormality constraint in both ab initio and semiempirical tight-binding calculations. This ensures that the derivative of the overlap matrix is computed only once per gradient evaluation.

From the above-mentioned description, it should be clear that the semiempirical integral library never computes any overlap integrals or their derivatives but rather obtains overlap integrals on input and returns energy/potential terms and corresponding W^S matrices or gradients as output. Given that the evaluation of overlap integrals and their gradients is standard in electronic structure codes, these are then evaluated by the integral library of the electronic structure program, in our case IntBox/TeraChem. This has a practical implication: up to this point, our semiempirical integral library setup is agnostic to the type of basis function employed, as long as they are centered on the individual atoms. In TeraChem, we employ Gaussian-type orbitals (GTOs), but also Slater-type orbitals (STOs) or less standard functions like polynomial “ramp” functions⁴⁷ could, in principle, be employed.

The GFN2-xTB Hamiltonian contains energy terms that depend on the electric dipole moment and quadrupole moment integrals. For variational ground state solutions, the computation of nuclear derivatives has been described in the Supplementary Material of Ref. 31. In our general semiempirical integral library, different and non-symmetric density matrices (X and Y) may enter the corresponding energy expression. Just like in the overlap case, we collect the corresponding intermediates for the dipole and quadrupole moment integrals,

Here, M^β and M^βγ denote the Cartesian dipole and quadrupole moment integrals, with components along the β and γ axes, respectively. The extra index A indicates that these are evaluated with the origin at atom A (with μ ∈ A; see Ref. 31 for details). For this reason, the individual W matrices for multipole integrals are non-symmetric with respect to exchanging the indices μν. The evaluation of the multipole gradients in this fashion is certainly non-standard in electronic structure packages, and a GTO-based routine has been added to the TeraChem integral library.

To demonstrate the advantages of our approach, we combine the GFN1-xTB Hamiltonian with the two-state interaction two-state-averaged spin-restricted ensemble-referenced Kohn–Sham method with two active electrons and orbitals, denoted SI2-SA2-REKS(2,2) in the following.^48–50 Here, the ground and first excited states are obtained as solutions from the state interaction Hamiltonian,

Here, the E_REKS(2,2) state is the energy of the perfectly spin-paired configuration, and E_OSS is the energy of the open-shell singlet configuration. The coupling parameter ∆_rs is related to the Lagrange multiplier for the active orbital pair r and s as well as to their fractional occupation number difference. The orbitals used in SI2-SA2-REKS(2,2) are obtained by minimizing the SA2-REKS(2,2) energy expression,

Here, the sum goes over all possible Slater determinants that can be generated in the (2,2) active space, and c_L are their weights in the ensemble. We note that GFN1- and GFN2-xTB have been specifically parameterized for the purpose of obtaining molecular geometries and describing noncovalent interactions, hence we do not expect these parameters to be ideal for excited state properties. Nevertheless, we refrain from modifying the semiempirical parameters since our goal here is to demonstrate how easily the semiempirical Hamiltonians can be used with ab initio electronic structure methods. In the future, it will be important to optimize the semiempirical parameters.

Adding Fock exchange is crucial to improve the description of excited states in both density functional and density functional tight-binding (DFTB) Hamiltonians. In the latter, it is typically added in a range-separated form approaching 100% Fock exchange at long distances.^20,21 We augment the GFN1-xTB′ Hamiltonian (i.e., GFN1-xTB without third-order term) with long-range Mulliken-approximated Fock exchange such that the energy becomes

Here, the damped and range-separated Coulomb terms are from Eq. (7). We will use ω = 0.4 a.u. (as well as c_FR = 0 and c_LR = 1) in the following (see Sec. III). Following Ref. 21 (but at variance with Ref. 20), we add the exchange directly to the Hamiltonian (as opposed to using the difference density). This is motivated by the fact that, though the GFN1-xTB Hamiltonian has been fitted empirically, it is formally based on the density fluctuation expansion of a functional without Fock exchange. It is noted for completeness that implementing a difference density formulation is possible by decomposing the exchange contributions into reference and fluctuation contributions, similar to what we did for the second- and third-order tight binding terms (see above). The corresponding reference contributions would then also be added to $Ẽnn$ and $Ẽne$⁠, respectively.

We noticed that the long-range correction alone is not sufficient for a reasonable description of the electronic states near intersections (see Sec. IV). This is in agreement with findings for SI2-SA2-REKS(2,2) in combination with long-range corrected DFTB (LC-DFTB).⁵¹ In that study, the authors resorted to spin-polarized DFTB Hamiltonians with corresponding parameters. Therein, an atomic on-site correction proportional to the spin density matrix is added, which leads to an improved description of the branching plane near minimum energy conical intersections (MECIs). In this work, we suggest a parameter-economic scheme that is inspired by spin-polarized DFTB but in line with our setup. We use the ab initio integral library in TeraChem to compute the atomic on-site ab initio K contributions and add them to the Hamiltonian. The corresponding energy then becomes

For the purpose of this paper, the scaling factor c_K = 0.05 is chosen based on benchmarks for vertical excitation energies (see Sec. III). But in future efforts, which should involve full reparameterization, this scaling factor can be treated it as an element-specific parameter as well. A schematic of the Fock exchange contributions added to the GFN1-xTB Hamiltonian is given in Fig. 3.

We combine this Hamiltonian with the SI2-SA2-REKS(2,2) approach, which is described in detail elsewhere.^49,50 Using this method, we compute the ground and first excited electronic states, their nuclear gradients, and nonadiabatic coupling vectors using the existing implementation in TeraChem⁵⁰ but with SQMBox semiempirical integrals (cf. Fig. 1). Finally, we also show that we can use the GFN1-xTB-type Hamiltonians in combination with other excited state wavefunction approaches, such as CASCI, CASSCF, RPA, TDA, and hh-TDA, demonstrating the strength of our approach while revealing limitations of the current GFN1-xTB parametrization for excited states.

We added the semiempirical integral library presented in this work to the current development version of the GPU-accelerated electronic structure program TeraChem.^27,28 All calculations presented in this work were carried out with this development version. We implemented the GFN1-xTB³² and GFN2-xTB³¹ Hamiltonian expressions into our semiempirical integral library, following the decomposition strategy described earlier. For now, we have implemented GFN2-xTB without the self-consistent D4 dispersion model^35,36 but rather in combination with the D3(BJ)^33,34 London dispersion correction and corresponding Axilrod–Teller–Muto correction, which had already been used in the original parametrization process (damping parameters: a₁ = 0.52, a₂ = 5.00, s₈ = 2.70, s₉ = 5.00, all other parameters remain unchanged).³¹ We will add support for self-consistent D4 dispersion at a later stage and note here that, due to the exponential dependence on the density matrix, this term will lead to contributions to the XC energy/potential equivalent in the semiempirical integral library presented here.

We tested our implementation first by comparing energies and gradients computed for ground state SCF solutions with the GFN1-xTB and GFN2-xTB Hamiltonians against the standalone xtb code.⁴¹ The implementation of Mulliken-approximated Fock exchange to SQMBox and the compatibility of excited state workflows were tested by comparison of analytic and numerical gradients. Among the considered electronic structure routines were configuration interaction singles (CIS), random-phase approximation (RPA), complete active space configuration interaction (CASCI) and self-consistent field (CASSCF), and spin-restricted ensemble-referenced Kohn–Sham (REKS) theory in combination with GFN1-xTB and GFN2-xTB Hamiltonians, with the restriction that third-order terms [see Eqs. (8)–(11)] of the tight binding expansion have been neglected up to this point and will be added for post-SCF/response-type schemes in the future. Detailed derivations and expressions for the individual Fock matrix quantities and gradients are given in the supplementary material.

In the computational performance assessment of the Mulliken-approximated Fock exchange, we performed GFN1-xTB SCF calculations augmented with 100% of the semiempirical Mulliken-approximated full-range Fock exchange (c_FR = 1 and c_LR = 0). Fullerene systems were taken from Ref. 52, and larger systems were constructed as assemblies of those. The considered geometries are given in the supplementary material and range from C₁₀₀ to a dimer of C₂₄₀@C₇₂₀. The reported timings in Fig. 4 are the timings for a single K build, averaged over the different SCF cycles. An incremental Fock build in terms of different density matrices was used in the SCF.⁵³ This allowed us to design an efficient mixed precision GPU-based scheme that is detailed in Sec. IV A, while the numerical precision assessment is shown in the supplementary material (Table S1 and Fig. S2). We note that the third-order term of the GFN1-xTB Hamiltonian has a nonlinear potential in the density matrix, hence, its mean-field potential is built using the full density matrix, similar to the XC potential in Kohn–Sham DFT calculations.⁵⁴ 

To demonstrate the strength of the TeraChem/SQMBox setup in enabling new method combinations, we combine the SI2-SA2-REKS(2,2) approach with the LC-GFN1-xTB′(K_at) Hamiltonian described in Sec. II E. The individual Fock exchange components described in Sec. II are added to the GFN1-xTB Hamiltonian, and we scanned for suitable range-separation and atomic exchange scaling factors ω and c_K, respectively. We used the lowest vertical excitation benchmark sets collected in Ref. 55 (systems originally taken from Refs. 56–58) and tested possible parameter combinations ω = 0.1, 0.2, 0.3, 0.4, and 0.5 and c_K = 0.005, 0.01, and 0.05, respectively. We selected ω = 0.4 and c_K = 0.05 in the following. We then optimized minimum energy conical intersections (MECIs) with this method using a well-established benchmark set.⁵⁹ It is emphasized that the standard TeraChem architecture (run conical) has been used, which employs DL-FIND (using hybrid delocalized internal coordinates).⁶⁰ In other words, no further modifications of the code were made after setting up and interfacing with the semiempirical integral library. All structures were converged to an S₁/S₀ gap of at most 1.5 meV, except for the ethylene (ethylidene) and methyliminium (methylamine) structures, for which the gap could not be reduced below about 0.2 eV (see Fig. 5 for structures and supplementary material for explicit energy values).

Figure 6 shows the S₁/S₀ potential energy surfaces around the twisted-pyramidal conical intersection geometry of ethylene with LC-GFN1-xTB′(K_at)/SI2-SA2-REKS(2,2), demonstrating that the topography is properly conical around the MECI. Here, potential energy scans are performed close to the MECI by sampling ±0.002 Å displacements along the branching plane vectors g and h. These are determined as the normalized difference gradient of the two states (g) and as the normalized nonadiabatic coupling vector between the two states (h) after Gram–Schmidt-type orthogonalization against g. To complement this, the GFN1-xTB′ Hamiltonian with 100% Mulliken-approximated Fock exchange [termed GFN1-xTB′(100) in the following] is combined with the floating-occupation molecular orbital complete-active space configuration interaction and the two state-averaged complete active space self-consistent field method, with both methods using a (2,2) active space [termed FOMO-CASCI(2,2) and SA2-CASSCF(2,2), respectively, in the following]. For FOMO-CASCI(2,2), a constant fractional occupation of 0.5 for the active spin-MOs is used.

A potential energy surface plot of ethylene using a rigid scan along the torsion and pyramidalization angles is shown to demonstrate the importance of adding Fock exchange (Fig. 7). Here, we used the SI2-SA2-REKS(2,2) approach with the LC-GFN1-xTB′(K_at) and GFN1-xTB′ Hamiltonians.

To identify some characteristics of the GFN1-xTB′ Hamiltonians in the context of excited state calculations, we additionally computed the lowest four vertical excitations of ethylene with the time-dependent density functional theory (TD-DFT) or random-phase approximation (RPA),⁶¹ its Tamm–Dancoff approximated (TDA) variant,⁶² as well as the recently presented variant of the hole–hole Tamm–Dancoff approximation (hh-TDA).⁵⁵ These were compared to the excited state obtained from SI2-SA2-REKS(2,2). We employ both the Fock exchange-inclusive LC-GFN1-xTB′(K_at) and the exchange-free GFN1-xTB′ Hamiltonians in this context (see Table I).

We first demonstrate the performance for the evaluation of the semiempirical Mulliken-approximated Fock exchange in our implementation. In this context, we show how single-precision arithmetic on GPUs can efficiently be exploited while retaining sufficient numerical precision in the total energies. This is followed by presenting results for calculations that were enabled by the newly presented modular setup: we combine the GFN1-xTB Hamiltonian with the multiconfigurational electronic structure method REKS and showcase this for the optimization of minimum energy crossing point geometries on an established benchmark set. The straightforward combination with CAS-type methods is then shown for the description of the conical intersection geometry of ethylene.

For full generalizability to ab initio codes, a semiempirical Mulliken-approximated Fock exchange has been added to the semiempirical integral library. Here, we assess the computational performance of our implementation in the computation of semiempirical K matrices. In initial implementations of long-range Fock exchange corrected DFTB, this became one of the computational bottlenecks, making LC-DFTB calculations slower by about a factor of 10 compared to DFTB.²⁰ In Fig. 4, we plot the average computational time for the formation of the semiempirical Mulliken-approximated Fock exchange matrix as a function of basis set size. For this purpose, we use the otherwise unmodified GFN1-xTB Hamiltonian with its corresponding basis set but include the Mulliken-approximated Fock exchange (100%), as described in Eq. (14). Fullerenes and assemblies of these ranging from 100 to 1920 atoms are used (see supplementary material for details). We report the time required to form the semiempirical K matrix, averaged over all SCF iterations. Even though we use a highly efficient implementation, which only requires four matrix–matrix multiplications [i.e., O(N³) steps], the computational expense becomes considerable for several thousand basis functions on the CPU (almost 80 s for close to 8000 basis functions using 4 CPU cores). Evaluation on the GPU in double precision significantly reduces the computational cost by at least a factor of two. We mention that we used a GeForce GTX 970 card, i.e., a consumer grade GPU from 2014. Such cards are known to perform dramatically better if they work in single-precision. This is confirmed here, and we find that the single-precision implementation allows us to form the semiempirical K matrix in not even 4 s for the largest system considered with almost 8000 basis functions. We use incremental Fock builds here, so a difference density matrix is used in Eq. (14) for all but the first SCF cycle.⁵³ In this way, only the very first formation of the semiempirical K matrix needs to be performed in double precision, while lower precision is sufficient for the remaining SCF iterations. We use lower precision with a hybrid CPU/mixed precision GPU scheme based on the following decomposition of the overlap matrix:

Here, we decompose the overlap matrix into a diagonal part (unit matrix I) and an off-diagonal part $S̄$⁠. In our mixed precision scheme, terms that depend on $S̄$ are evaluated in single precision. Terms that are proportional to the diagonal part only are equivalent to a quadratically scaling ZDO-type evaluation of K. These are always evaluated on the CPU in double precision [see Eq. (64) in the supplementary material for details]. The latter does not add any significant computational cost to the overall computation time. Based on our numerical findings, we furthermore suggest to evaluate K entirely in double precision on the CPU if n_AO < 300 + 100 × n_CPU, where the latter is the number of available CPU cores. For larger systems, we switch between double and mixed precision on the GPU, depending on the type of density matrix passed in. If the density matrix is a full particle density matrix, the code runs in double precision, while the mixed precision is used if a difference or transition density matrix is used. This is tested by forming the trace, i.e., the semiempirical K matrix is formed in double precision, if

where we keep the threshold parameter ϵ_thr close to zero. This identifies if a full or a difference/transition density is passed in. The extra computational cost of forming the trace, which scales with n_AO × n_AO, is negligible compared to the overall cubic formation of the Mulliken-approximated K matrix. This dynamic switching in combination with the mixed precision evaluation of K, allows us to exploit the extremely fast evaluation of Mulliken-approximated Fock exchange with consumer-grade GPUs. In combination with an incremental Fock matrix SCF procedure,⁵³ we note only minor numerical differences in the final single point energies due to the formation of the Mulliken-approximated exchange matrix in single precision, which becomes reduced by an order of magnitude by each using double precision in the first iteration and by resorting to the mixed precision scheme (see Table S1 in the supplementary material). With this implementation, the rate-determining computational step in exchange-augmented xTB-type calculations remains the diagonalization of the Hamiltonian matrix (see Fig. 4). With the new algorithm, there is no longer any reason to avoid Mulliken-approximated Fock exchange.

The availability of semiempirical Hamiltonian elements via the semiempirical integral library enables new method combinations, as demonstrated in the following sections. The LC-GFN1-xTB′(K_at) Hamiltonian is combined with the REKS method and applied in the optimization of minimum energy conical intersections. We also show that the GFN1-xTB′ Hamiltonian with a 100% Mulliken-approximated Fock exchange can be straightforwardly combined with the CASCI and CASSCF methods. These are then used, along with the aforementioned REKS scheme, in the description of the minimum energy conical intersection region. Not only energies and gradients, but also the analytic derivative coupling implementations in TeraChem are now accessible with semiempirical Hamiltonians through the established framework. In the context of ethylene, we point out peculiarities related to the purpose-specific parametrization of the GFN1-xTB Hamiltonian, using additional electronic structure methods for further reference in the computation of the ππ* excited state energy. Finally, we will point out the importance of Fock exchange (on-site and long-range) in the description of the excited state potential energy surface with REKS.

We combine the two-state interaction two-state-averaged spin-restricted ensemble-averaged Kohn–Sham DFT method [SI2-SA2-REKS(2,2)]⁴⁸ with the semiempirical LC-GFN1-xTB′(K_at) Hamiltonian (see Sec. II E). In Fig. 5, S₁/S₀ MECI geometries that are optimized at this level of theory are shown (green) and overlayed with the reference geometries computed at the multi-reference configuration interaction singles and doubles (MRCISD) level of theory (taken from Ref. 59). We find that the optimized structures agree well with the reference structures, suggesting that this method combination could have potential in the explorations and preoptimizations of MECIs in photochemical studies. We note, however, that the relative energies for the MECIs are underestimated for all systems but ketene. This is likely due to the parametrization of the GFN1-xTB Hamiltonian, which has been parameterized (without Fock exchange) for a different purpose, i.e., minimum geometries, nuclear Hessians, and noncovalent interactions.³² At this point, it is emphasized that the MECI optimization procedure not only makes use of adiabatic nuclear gradients but also of nonadiabatic coupling vectors between both states. Due to the established setup (see Fig. 1), gradients and nonadiabatic coupling vectors for this new method combination become available simply by providing the aforementioned semiempirical integral components (Fock matrix equivalents and gradients) in place of the ab initio components. This way, elaborate re-derivation and implementation becomes unnecessary, and the existing implementation in TeraChem can be used.^49,50

To highlight this further, the S₀ and S₁ potential energy surfaces around the twisted-pyramidal MECI of ethylene, which were computed with LC-GFN1-xTB′(K_at)/SI2-SA2-REKS(2,2), are shown in Fig. 8. Here, we used the optimized MECI geometry and computed the S₁ and S₀ surfaces on a grid along the branching plane vectors g and h, which are obtained from the gradient difference and nonadiabatic coupling vectors between both states, respectively. We observe a peaked conical intersection along these vectors. The branching plane vectors are visualized in the insets and indicate contributions from the pyramidalization, torsion, and also from bond C–C bond stretching modes, which are in good agreement with the branching plane vectors from previously reported MRCISD and SI2-SA2-REKS(2,2)/ωPBEh calculations.^51,59 Consequently, the presented framework provided us with easy access to gradients and nonadiabatic coupling, enabling photochemical studies in the near future.

In Fig. 8, we also plot the corresponding potential energy surfaces computed with FOMO-CASCI(2,2) and SA2-CASSCF(2,2) based on the GFN1-xTB′(100) Hamiltonian (see Sec. III). We find again that we retrieve the full functionality of the available implementations in TeraChem, which initially had been formulated only in terms of ab initio theories.^63–65 Thanks to the cleanly separated Hamiltonian/wavefunction setup and the SQMBox library, these novel method combinations become accessible with semiempirical Hamiltonians, like the GFN1-xTB Hamiltonian used here.

We have already noted before that we are using a purpose-specifically fitted Hamiltonian (original purpose: ground state geometries, noncovalent interactions, and nuclear Hessians) in an off-purpose setting simply to demonstrate our implementation. Unsurprisingly, we find that the GFN1-xTB parametrization is not optimum for excited state purposes. Using the linear-response time-dependent density functional theory (TD-DFT) implementation,⁶¹ its Tamm–Dancoff approximated variant⁶² and the recently added hole–hole Tamm–Dancoff approximation (hh-TDA)⁵⁵ implementation in TeraChem along with the semiempirical library, we can compute the corresponding excited states based on these methods in combination with the semiempirical tight-binding Hamiltonian. We compute the lowest four excited states of ethylene in combination with the LC-GFN1-xTB′(K_at) Hamiltonian and compare it to the first excited state from SI2-SA2-REKS(2,2). The results are listed in Table I. We also did the same comparison for the exchange-free GFN1-xTB′ Hamiltonian. We find that the visible ππ* state of ethylene, which is essentially a HOMO–LUMO transition, is the lowest excited state in SI2-SA2-REKS(2,2) but in no other method. The only exception is hh-TDA/GFN1-xTB′, where the excitations exactly correspond to the corresponding orbital energy differences since the Hamiltonian does not include exchange (see Ref. 55 for details). This indicates that the GFN1-xTB′ Hamiltonian erroneously places other electronic states in the minimal valence basis set (which cannot describe Rydberg states) to be lower than the ππ* state. The lowest excited state in all other schemes is in fact a σ(C–H)π* state. In the TDA scheme, the dark σ(C–C)π* state also lies lower (i.e., it is S₂) than the ππ* state, which is different in RPA, where the ππ* state lies lower. It is known that in the electrostatic monopole approximation, which is a consequence of the Mulliken-approximated two-electron integrals, the two-electron (ia|jb) terms vanish for i = σ and a = π* due to symmetry and, hence, such states are identical in RPA and TDA.⁶⁶ This is different for the ππ* states, which are strongly blue-shifted, predominantly due to the diagonal (ππ*|ππ*) term in the TDA/RPA A matrix. RPA somewhat reduces this shift, leading to lower excitation energies. This finding along with the strong blue-shift of the ππ* state in hh-TDA, if the coupling is activated due to the addition of exchange [LC-GFN1-xTB′(K_at)], leads to the conclusion that the Coulomb interactions in the GFN1-xTB′ Hamiltonian are somewhat overestimated or imbalanced with respect to the orbital energy differences, which are expected to have strong contributions from the extended Hückel term. Meanwhile, the parameters in the xTB-type Hamiltonian should offer enough flexibility to, in principle, provide a well-working Hamiltonian also for excited states, its original parametrization is specifically fitted (“overfitted”) to the original purpose of providing good ground state geometries and noncovalent interactions. The SI2-SA2-REKS(2,2) scheme still works well regardless, since the (2,2) active space prevents any other orbitals from contributing to the excited state. We conclude that a full reparameterization of the Hamiltonian is necessary to obtain a robust and accurate method for excited states. In this context, the presented framework allows the straightforward application of a broader range of excited state methods, just as we did now in this assessment for ethylene, and can now be used to check a re-parametrized Hamiltonian for consistency.

As a final point, we want to highlight the importance of the Fock exchange in Hamiltonians to be used for photochemical studies. It is known that vertical excitations are generally better described for approximate density functional and tight-binding models when nonlocal Fock exchange is included. This is also why we chose to augment the Hamiltonian with the Fock exchange since the GFN1-xTB Hamiltonian does not contain the Fock exchange or any spin-density dependent terms. To demonstrate the implications for photochemical applications beyond vertical excitations, which are long-term targets of the TeraChem/SQMBox framework, we plot the potential energy surfaces of the two lowest states of ethylene along the torsional and pyramidalization coordinates in Fig. 7. Here, the surfaces are computed within the SI2-SA2-REKS(2,2) framework in combination with the exchange-free GFN1-xTB′ and the exchange-inclusive LC-GFN1-xTB′(K_at) Hamiltonian (see Sec. II E). In Figs. 7(c) and 7(d), we see the direct consequence of employing a Hamiltonian, which is free from the exchange or spin-density components. Here, the S₀/S₁ intersection coincides with the S₁ minimum geometry and is found at a purely twisted and non-pyramidalized geometry. High-level reference calculations of complete active space second-order perturbation theory type indicate that both the S₁ minimum and the MECI geometry should be significantly pyramidalized (pyramidalization angles >70°).^55,67 This failure is immediately cured by the inclusion of Fock exchange as in the LC-GFN1-xTB′(K_at) Hamiltonian [Figs. 7(a) and 7(b)]. The MECI geometry is then found at 60 degrees of pyramidalization, i.e., much closer to the aforementioned reference, and the overall shape of the S₁ potential energy surface is significantly improved (cf. Fig. 1 in Ref. 68). We furthermore find that the potential energy surface at small pyramidalization angles and 90° torsion requires the addition of a short-range correction, which we introduce in the form of an atomic ab initio K contribution (see supplementary material for the plots with the LC-GFN1-xTB′ Hamiltonian). Hence, we conclude that a tight-binding approach with long-range exchange correction alone is not likely to be sufficient for photochemical studies.

We have presented a modular framework that facilitates the joint development of ab initio and semiempirical electronic structure workflows. In our approach, we group the semiempirical Hamiltonian terms by their orders in the one-particle density matrix and link them to the ab initio Hamiltonian quantities with an equivalent dependency on the density matrix (and nuclear coordinates). The semiempirical integral library then provides the Mulliken-approximated semiempirical Hamiltonian matrix equivalents that resemble the ab initio/Kohn–Sham Hamiltonian matrix quantities. In an existing ab initio electronic structure program, the semiempirical integral library is then called instead of the ab initio integral library, which enables the straightforward extension of electronic structure methods to semiempirical Hamiltonians. This removes any extra effort to re-derive and re-implement electronic structure methods for semiempirical Hamiltonians once an ab initio framework already exists. We showcase this by combining the recently presented extended tight-binding Hamiltonian GFN1-xTB (augmented with Fock exchange) with the state-interaction state-averaged spin-restricted ensemble-averaged Kohn–Sham (SI-SA-REKS) method. We apply this combined scheme to the optimization of conical intersections for a benchmark set from the literature and to the computation of the potential energy surfaces around the twisted-pyramidal S₁/S₀ minimum energy conical intersection of ethylene. The latter was also carried out with Fock exchange-augmented GFN1-xTB in combination with wave-function-based complete active space methods (state-averaged CASSCF and floating-occupation molecular orbital CASCI). This demonstrates that our framework gives full access to the existing functionality of the ab initio code, including gradients and nonadiabatic coupling vectors for the aforementioned method combinations. The derivation and implementation of these quantities, which can represent a considerable effort,^{49,50,63–65} is entirely avoided by leveraging the existing implementation in TeraChem with the modular integral library setup. Currently, the library is a proof-of-concept implemented within the TeraChem package and supports both GFN-xTB-type Hamiltonians. Current efforts aim at generalizing the workflow to set up and re-parametrize such semiempirical Hamiltonians, whereupon we hope to make the library more widely available. All Hamiltonian matrix and derivative expressions are thoroughly documented in the supplementary material.

We envisage the presented framework to enable faster method and workflow developments in which the ab initio and semiempirical electronic structure developments are no longer decoupled. Multi-level procedures, in which semiempirical electronic structure methods are used for exploratory calculations while the corresponding ab initio variant is then used for refinement, are naturally facilitated by this framework. Furthermore, testing in the course of implementing new electronic structure methods becomes simplified by resorting to the semiempirical integral library for this purpose: a standard consistency check in electronic structure developments is the comparison of analytic to numerical gradients, with the latter typically requiring a considerable amount of computation time. Given that the semiempirical integral library offers one-to-one correspondence to the terms from the ab initio library but is significantly faster (3–4 orders of magnitude), future testing and AO-direct electronic structure code developments can be accelerated as well by our framework.

See the supplementary material PDF file for detailed equations, derivations, and algorithms. Further benchmarking data for mixed precision Fock matrix construction schemes. ZIP file for molecular geometries used for benchmarking and testing.

This work was supported by the Office of Naval Research (Grant Nos. N00014-21-1-2151 and N00014-18-1-2659). C.B. thanks the German National Academy of Sciences Leopoldina for a Leopoldina Fellowship (Project No. LPDS 2018–09). We thank Jimmy K. Yu for providing a template script for plotting the ethylene potential energy surfaces.

T.J.M. is a co-founder of PetaChem, LLC.

Christoph Bannwarth: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Todd J. Martínez: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.