Development of NOTCH, an all-electron, beyond-NDDO semiempirical method: Application to diatomic molecules

Wang, Zikuan; Neese, Frank

doi:10.1063/5.0141686

In this work, we develop a new semiempirical method, dubbed NOTCH (Natural Orbital Tied Constructed Hamiltonian). Compared to existing semiempirical methods, NOTCH is less empirical in its functional form as well as parameterization. Specifically, in NOTCH, (1) the core electrons are treated explicitly; (2) the nuclear–nuclear repulsion term is calculated analytically, without any empirical parameterization; (3) the contraction coefficients of the atomic orbital (AO) basis depend on the coordinates of the neighboring atoms, which allows the size of AOs to depend on the molecular environment, despite the fact that a minimal basis set is used; (4) the one-center integrals of free atoms are derived from scalar relativistic multireference equation-of-motion coupled cluster calculations instead of empirical fitting, drastically reducing the number of necessary empirical parameters; (5) the (AA|AB) and (AB|AB)-type two-center integrals are explicitly included, going beyond the neglect of differential diatomic overlap approximation; and (6) the integrals depend on the atomic charges, effectively mimicking the “breathing” of AOs when the atomic charge varies. For this preliminary report, the model has been parameterized for the elements H–Ne, giving only 8 empirical global parameters. Preliminary results on the ionization potentials, electron affinities, and excitation energies of atoms and diatomic molecules, as well as the equilibrium geometries, vibrational frequencies dipole moments, and bond dissociation energies of diatomic molecules, show that the accuracy of NOTCH rivals or exceeds those of popular semiempirical methods (including PM3, PM7, OM2, OM3, GFN-xTB, and GFN2-xTB) as well as the cost-effective ab initio method Hartree–Fock-3c.

I. INTRODUCTION

Semiempirical methods are quantum chemistry methods with intermediate cost, accuracy, and degree of empiricism between density functional theory (DFT) and molecular mechanics (MM) methods and have found wide use in the study of large molecules (>1000 atoms) and the molecular dynamics/conformational search/pre-optimization of moderately sized molecules, to name a few.^1,2 As approximate, parameterized versions of minimal basis Hartree–Fock (HF) or DFT calculations, semiempirical methods are more robust and describe more of the underlying physics than classical MM methods while typically requiring only a few hundredths of the time of a DFT calculation with a double-zeta basis set, even when efficient approximations, such as resolution of the identity (RI)³ and chain-of-spheres exchange (COSX),^4–7 are applied to the latter. When properly parameterized, semiempirical methods can often achieve useful accuracy on a handful of properties of sufficiently “well-behaved” molecules, although like any other method with a fair amount of empiricism, it is notoriously difficult to transfer the accuracy to properties and types of molecules that are not used in the training procedure. Intriguingly, few methods are known whose timing and accuracy are between those of most semiempirical methods and RI-DFT. Thus, even if slightly more computational budget is available than what is required for a semiempirical calculation, one usually cannot take advantage of this to improve the accuracy of the computational results.

The most popular semiempirical methods mainly fall into two categories: the zero-differential overlap (ZDO) methods (e.g., the PMx^8–10 and OMx¹¹ families of methods) and tight-binding (TB) methods (e.g., the DFTB^12,13 and xTB^14–16 families of methods), respectively. In the ZDO family of methods,^17,18 the following approximations are made (among others): (1) Only the (AA|AA)- and (AA|BB)-type two-electron integrals are retained, i.e., $(μ_{A} ν_{B} | κ_{C} λ_{D}) = δ_{A B} δ_{C D} (μ_{A} ν_{A} | κ_{C} λ_{C})$ ⁠, which is the chief reason for the low cost of the ZDO methods; (2) consistent with the above approximation, the Roothaan–Hall equation FC = SCϵ—where F, C, S, and ϵ are the Fock, molecular orbital (MO) coefficient, overlap, and orbital energy matrices, respectively—is replaced by FC = Cϵ, i.e., the overlap matrix S is set to the unit matrix. On the other hand, the TB methods^12–14 usually do not approximate S as the unit matrix, while they describe the electrostatic (and exchange-correlation) interaction as a damped interaction between atomic charges (usually based on a Mulliken population analysis), thereby achieving a similar cost as the ZDO methods. However, most TB methods lack an exact exchange contribution, which introduces considerable self-interaction errors that can be difficult to remedy and makes them unable to describe finer electronic structure effects, such as electronic multiplets or spin state splittings.

The present work aims to develop a new semiempirical method, NOTCH (Natural Orbital Tied Constructed Hamiltonian), which is designed to fill up the timing and accuracy gap between the common semiempirical methods and RI-DFT. It makes fewer approximations and, at the same time, introduces less empiricism than both the ZDO and the TB methods. As its name suggests, the NOTCH atomic orbital (AO) basis functions are designed to mimic atomic natural orbitals (ANOs) as closely as possible. Even more importantly, the NOTCH one-center one-electron and two-electron integrals of free atoms, as well as those of monocationic and monoanionic atomic ions, are required to reproduce the ANO basis integrals of an effective ab initio Hamiltonian that implicitly includes dynamic correlation and scalar relativity, i.e., the second-order Douglas–Kroll–Hess (DKH2)^19,20 similarity-transformed multireference equation-of-motion coupled cluster (MR-EOM-CC)^21–24 Hamiltonian. In this way, the parameterization of atomic integrals is solved in a convenient, systematic, and completely nonempirical way. The nonempirical atomic integrals serve also as a starting point for building the one-electron and two-electron integrals in molecules.

A few other salient points of our construction include the following:

In NOTCH, the core electrons are treated explicitly, rather than approximated by absorbing their contribution into the nuclear attraction potential, as what is done in most semiempirical methods. This obviates the need for complicated treatments of the core–core repulsion energy, which have often been “abused” to absorb the errors of other parts of the semiempirical method and contributed significantly to the highly empirical nature of semiempirical methods and a partially severe lack of transferability.
In NOTCH, the radial distributions of the AOs depend on the molecular geometry, so that, e.g., an AO becomes more compact in a molecular environment than in the free atom. By this way, we mimic a double-zeta basis calculation within a minimal basis framework. Semiempirical methods employing geometry-dependent basis sets are rare, one prominent example being sTDA-xTB.²⁵
The radial distributions of the AOs also depend on the atomic charges in NOTCH, i.e., the more positive charge an atom carries, the more compact its AOs will be. A simple way to account for this is to add charge-dependent corrections to the one-center one-electron and two-electron integrals. Similar ideas have been explored in ZDO^26–28 and TB^15,16,29,30 methods.
The NOTCH AO two-electron integrals are first calculated analytically from the AO basis functions, which are subsequently corrected in order to incorporate correlation effects. This stands in contrast to the common practice of using crude approximations like the Dewar–Sabelli–Klopman formula^31–35 to approximate the integrals and then hope to absorb dynamic correlation effects and various errors by parameterization. Note that the left–right correlation (LRC) of bonds, which is inherently nonlocal, is also explicitly accounted for.
In NOTCH, we retain not only the AO basis (AA|AA) and (AA|BB)-type two-electron integrals but also the (AA|AB) and (AB|AB)-type ones.
In NOTCH, the dispersion interaction is not described as a post-self-consistent field (SCF) correction to the energy but as a correction to the (AA|BB)-type two-electron integrals. The main rationale for this choice is that dispersion is inherently a correlation effect and thus represents an interaction between electrons, not between atoms. Similar ideas have been pursued in some semiempirical methods³⁶ but have not gained wide attention.
In NOTCH, a classical charge polarization term is added as a correction to the one-electron integrals to account for the underestimation of polarizability by the minimal basis set formalism. The idea of adding a classical charge polarization term has been explored before, but as a post-SCF correction to the total energy instead of a correction to the integrals.^37,38 The same correction also mimics the charge polarization due to bonding and therefore improves the description of bonding as well.

It should be noted that the present work is preliminary, as only the elements from H to Ne are parameterized, and only results for diatomic molecules are reported. Nevertheless, our benchmark results on all the possible neutral diatomic molecules that can form between the elements H–Ne show that NOTCH provides equilibrium bond length, vibrational frequency, dipole moment, ionization potential (IP), electron affinity (EA), bond dissociation energy, and excitation energy results that are at least on par with, and frequently outperform, the best existing semiempirical methods. More importantly, NOTCH is characterized by a lack of severe outliers, which can at least partly be traced back to NOTCH’s greatly reduced empiricism compared to other semiempirical methods. The possible ways to extend NOTCH to multiatomic molecules will be briefly discussed at the end of the manuscript.

II. THEORY

A. Notations

Throughout the manuscript, the following notations are used: Atoms are denoted by A, B, C, …, AOs and orthogonalized AOs (OAOs) are denoted by μ, ν, κ, λ, …, occupied and virtual MOs are labeled by i, j, … and a, b, …, respectively; spin indices are denoted by σ, τ, …. The notation μ_A stands for a primitive AO, a contracted AO, or an OAO centered on atom A.

Integrals without any superscript, e.g., $(μ ν | κ λ)$ ⁠, h_μν, T_μν, or $(μ ν | A)$ ⁠, denote analytic, uncorrected AO integrals calculated over ano-pc-0 basis functions (vide infra). Integrals with overbars [e.g., $\bar{(μ ν | κ λ)}$ ] stand for corrected integrals that are explained in detail below.

As the present manuscript deals with monoatomic and diatomic systems only, it will sometimes be assumed that the number of atoms in the system is no larger than 2, in order to simplify the expressions. Likewise, for the sake of simplicity, we may sometimes assume that all the elements in the molecule belong to the first two periods, but it is trivial to generalize the conclusions to heavier elements.

B. Overview

The NOTCH SCF energy is given by the following formula, where all integrals are initially evaluated analytically at the HF/ano-pc-0 level of theory. However, the actually used integrals (indicated by the overbar) are corrected to effectively include correlation and relativity as well as further modifications to overcome the limitation of the minimal basis set. The energy expression reads as follows:

\begin{align} E^{NOTCH} & = E_{NN} + \sum_{μ ν σ} \bar{h_{μ ν}} [D] D_{μ ν σ} + \frac{1}{2} \sum_{μ ν σ} \sum_{κ λ τ} D_{μ ν σ} \\ \times (\bar{(μ ν | κ λ)} [D] - δ_{σ τ} \bar{(μ κ | ν λ)} [D]) D_{κ λ τ} + E^{LRC} [D], \end{align}

(1)

where D is the density matrix and E^LRC is the left–right correlation energy correction that will be detailed in Sec. II H. Herein ano-pc-0 is an all-electron minimal basis set with geometry-dependent contraction coefficients, whose construction will be detailed in Subsection II C. Thanks to the explicit description of the core electrons and thus the core-core repulsion, the E_NN in Eq. is given by the exact nuclear repulsion energy, without any empirical correction:

E_{NN} = \sum_{A > B} \frac{Z_{A} Z_{B}}{R_{A B}} .

(2)

Another special feature of NOTCH is that the integrals depend on the atomic charges (and therefore on the density matrix D). In this way we are able to mimic the expansion/shrinking of AOs as a result of the atom acquiring negative/positive charge, respectively (which otherwise requires at least a split valence basis to describe), despite using a minimal basis set.

Like in any other SCF method, the density matrix is required to minimize the electronic energy E^NOTCH, under the idempotency constraint. This leads to the familiar Roothaan–Hall equation,

F C = S C ϵ,

(3)

where

F_{μ ν σ}^{NOTCH} = \frac{\partial E^{NOTCH}}{\partial D_{μ ν σ}}

⁠. However, owing to the density matrix dependencies of the integrals, the Fock matrix consistent with Eq. is not

F_{μ ν σ}^{NOTCH} = \bar{h_{μ ν}} [D] + \sum_{κ λ τ} (\bar{(μ ν | κ λ)} [D] - δ_{σ τ} \bar{(μ κ | ν λ)} [D]) D_{κ λ τ}

(4)

but rather acquires additional terms arising from the density matrix element derivatives of the integrals, as a result of the variational principle, giving us

\begin{align} F_{μ ν σ}^{NOTCH} & = \bar{h_{μ ν}} [D] + \sum_{κ λ τ} \frac{\partial \bar{h_{κ λ}} [D]}{\partial D_{μ ν σ}} D_{κ λ τ} \\ + \sum_{κ λ τ} (\bar{(μ ν | κ λ)} [D] - δ_{σ τ} \bar{(μ κ | ν λ)} [D]) D_{κ λ τ} \\ + \frac{1}{2} \sum_{μ^{'} ν^{'} σ^{'} κ λ τ} (\frac{\partial \bar{(μ^{'} ν^{'} | κ λ)} [D]}{\partial D_{μ ν σ}} \\ - δ_{σ^{'} τ} \frac{\partial \bar{(μ^{'} κ | ν^{'} λ)} [D]}{\partial D_{μ ν σ}}) D_{μ^{'} ν^{'} σ^{'}} D_{κ λ τ} + \frac{\partial E^{LRC} [D]}{\partial D_{μ ν σ}} . \end{align}

(5)

For the sake of simplicity, in the remaining parts of the manuscript, we use $\bar{h_{μ ν}}$ and $\bar{(μ ν | κ λ)}$ to denote the “zero-charge” core Hamiltonian and two-electron integrals, respectively, i.e., the integrals one obtains in the case that all atomic charges are exactly zero.

C. Basis set

The basis set ano-pc-0 is obtained by a geometry-dependent re-contraction of the split valence, polarization function free basis set pc-0,^39,40

ψ_{μ_{A}}^{ano - p c - 0} (r; R_{A}, \{R_{B}\}) = \sum_{ν_{A}} c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\}) ϕ_{ν_{A}}^{pc - 0} (r; R_{A}),

(6)

where

c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\})

are the contraction coefficients which depend on the coordinates of the surrounding atoms

\{R_{B}\}

⁠, and

ϕ_{ν_{A}}^{pc - 0} (r; R_{A})

are the pc-0 primitive functions. For Li and Be, we also add the diffuse p function from the aug-pc-0 basis set⁴¹ to the ano-pc-0 primitive functions, to improve the description of the 2p orbitals. This yields a generally contracted minimal basis set with the following composition: 3s (uncontracted)/1s (contracted) for H and He, 5s2p (uncontracted)/2s1p (contracted) for Li and Be and 5s3p (uncontracted)/2s1p (contracted) for B–Ne. Note that the pc-0 primitive functions, which are optimized for partial general contraction, are more efficient than those of fully segmentally contracted basis sets like Slater-type orbitals (STO)-3G⁴² and MINIX,⁴³ that possess 6 primitive s functions instead of 5.

The contraction coefficients

c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\})

are determined via the following procedure:

A scalar relativistic second-order DKH2 state-averaged (SA)-complete active space self-consistent field (CASSCF) calculation with the decontracted pc-0 basis set is carried out for the neutral atoms of the elements Li–Ne, where the CAS active space includes all valence orbitals, and all states are averaged. (The DKH2 scalar relativistic treatment is of course unnecessary for the elements H–Ne, but will be necessary for the heavy elements that will be parameterized in the future; we have thus used DKH2 here as well, for consistency reasons.) For Li and Be, we consider the 2p orbitals as belonging to the valence orbitals, and furthermore we add the diffuse p function from the aug-pc-0 basis set to improve the description of the 2p orbitals. The converged inactive orbital is then taken as the 1s function of the ano-pc-0 basis set.
The valence functions of the ano-pc-0 basis set are also given as contractions of the pc-0 primitive functions (plus the aug-pc-0 diffuse p functions for Li and Be). We first variationally optimized the contraction coefficients of the valence functions to minimize the DKH2 state specific (SS)-CASSCF ground state energies of the homodimers of elements H–Ne at the bond lengths 3.0, 2.9, …, 0.1 Å, under the constraint that the 2s function (if it exists) is orthogonal to the 1s function. We then interpolate the optimized contraction coefficients as a function of the bond length, with the following formula:

\begin{align} c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\}) \\ = c_{μ_{A} ν_{A}}^{(0)} + \sum_{B \neq A} (c_{μ_{A} ν_{A}}^{(1)} \int ρ_{A}^{PModel} (|r - R_{A}|) V_{B}^{PModel} (|r - R_{B}|) d r \\ + c_{μ_{A} ν_{A}}^{(2)} \int {\frac{\partial ρ_{A}^{PModel} (r^{'})}{\partial r^{'}}|}_{r^{'} = |r - R_{A}|} V_{B}^{PModel} (|r - R_{B}|) d r), \end{align}

(7)

where

c_{μ_{A} ν_{A}}^{(0 - 2)}

are fitted coefficients, and

ρ_{A}^{PModel} (r)

is the spherically averaged atomic density of atom A, which is calculated at the DKH2 restricted open-shell Hartree-Fock (ROHF)/UGBS⁴⁴ level of theory, fitted by a sum of s-type Gaussian functions and tabulated in advance;

V_{B}^{PModel} (r)

is the Coulomb potential generated by the atomic density

ρ_{B}^{PModel} (r)

plus the nuclear charge of atom B. Since the atomic densities are expressed as a sum of Gaussians, the integrals in Eq. can be readily evaluated analytically. The functional form of Eq. is justified in Appendix A. Although

c_{μ_{A} ν_{A}}^{(0 - 2)}

are determined by fitting, they are fitted against variationally optimized contraction coefficients, instead of high-level reference data (which are usually used in existing semiempirical methods). They can be compared with the contraction coefficients of ab initio basis sets, which are often also determined by variationally minimizing the energies of atoms and simple molecules. Note also that, although

c_{μ_{A} ν_{A}}^{(0 - 2)}

are fitted against homodiatomic molecules only, they are found to be well transferable to heterodiatomic molecules.

The radial parts of some selected ano-pc-0 basis functions of atoms and diatomic molecules are plotted in Fig. 1 as a function of the bond length. As expected, the basis functions are generally more compact in covalently bonded diatomic molecules than in free atoms, and become progressively more compact when the bond shortens or (in the case of H 1s) when the other atom becomes more and more electronegative. Thus, ano-pc-0 mimics a double-zeta basis set despite being a minimal basis set, thanks to the geometry dependence of the contraction coefficients. It is also interesting to note that the radial distribution of the C 2s AO is nearly invariant with respect to the chemical environment, probably because it has to remain orthogonal to the 1s AO, which restricts its ability to shrink in a molecular environment. However, there are two important effects that ano-pc-0 cannot yet capture: (1) the p_x, p_y, and p_z functions of an atom can have different radial distribution profiles, a frequently cited example being linear molecules where the p orbital parallel to the bond axis (p_σ) is typically more compact than those perpendicular to the bond axis (p_π); and (2) the basis functions may deform by mixing with polarization functions. These effects are crucial for an accurate prediction of equilibrium bond lengths and vibrational frequencies, and will be discussed in Sec. II G.

FIG. 1.

View large Download slide

Radial parts of some ano-pc-0 basis functions in different diatomic molecules with different bond lengths. The C 2p (F-C, R = 1.3 Å) curve nearly coincides with the C 2p (C-C, R = 1.5 Å) one.

D. NOTCH integrals of free atoms and atomic ions

To keep the empiricism of NOTCH at the lowest possible level while still obtaining decent and uniform accuracy, the core Hamiltonian and two-electron NOTCH integrals of free atoms are determined from DKH2 multireference equation-of-motion coupled cluster singles and doubles (MR-EOM-CCSD)/ANO-RCC-Full^45–47 calculations of neutral atoms. We make full use of the fact that the ano-pc-0 AOs resemble the atomic natural orbitals (since they are optimized for atomic CASSCF calculations), so that for atomic calculations, the integrals of the former can be replaced by those of the latter to take advantage of the higher basis set quality and (more importantly) the implicit inclusion of dynamic correlation of the latter. This construction is at the heart of the NOTCH idea and is the reason why we insisted on using a generally contracted basis set. In addition, we require that the 1s and 2s functions of every element be orthogonal: a segmented contracted or non-orthogonal basis set may expand the correct s-space, but its basis functions generally do not resemble the exact atomic natural orbitals as orthogonal generally contracted basis sets can do.

Specifically, a DKH2 SA-CASSCF/ANO-RCC-Full calculation is first performed on the neutral atom, with the active space being all the valence orbitals, and all states being averaged. Then, the MR-EOM-T|T^†|SXD|U method²¹ is invoked, which consists of four sequential similarity transformations of the full Hamiltonian

\hat{H}

that approximately (but quite accurately) decouple the CAS configurations from their single- and double excitations in the first-order interacting space, until only the couplings with the one-particle (1p) and one-hole (1h) configurations remain,

\hat{\bar{H}} = e^{- \hat{T}} \hat{H} e^{\hat{T}},

(8)

\hat{\tilde{H}} = e^{- {\hat{T}}^{†}} {\hat{\bar{H}}}_{2} e^{{\hat{T}}^{†}},

(9)

\hat{\bar{F}} = {\{e^{\hat{S} + \hat{X} + \hat{D}}\}}^{- 1} {\hat{\tilde{H}}}_{2} \{e^{\hat{S} + \hat{X} + \hat{D}}\},

(10)

\hat{G} = e^{- \hat{U}} {\hat{\bar{F}}}_{2} e^{\hat{U}} .

(11)

Here, the curly braces denote the Kutzelnigg–Mukherjee normal ordering, and the subscript 2 denotes the truncation up to two-body operators. The definitions of the cluster operators

\hat{T}

⁠,

\hat{S}

⁠,

\hat{X}

⁠,

\hat{D}

and

\hat{U}

are given in the original references of the MR-EOM method.^21,22,24 The final transformed Hamiltonian

\hat{G}

for an atom A has the following form (here

{\hat{E}}_{ν_{A}}^{μ_{A}}

and

{\hat{E}}_{ν_{A} λ_{A}}^{μ_{A} κ_{A}}

are the usual spin-free excitation operators):

\hat{G} = g_{0} + \sum_{μ_{A} ν_{A}} g_{ν_{A}}^{μ_{A}} \{{\hat{E}}_{ν_{A}}^{μ_{A}}\} + \sum_{μ_{A} ν_{A} κ_{A} λ_{A}} g_{ν_{A} λ_{A}}^{μ_{A} κ_{A}} \{{\hat{E}}_{ν_{A} λ_{A}}^{μ_{A} κ_{A}}\} + \dots

(12)

The standard MR-EOM workflow proceeds to diagonalize

\hat{G}

in the space spanned by the CAS configurations plus the 1h and 1p configurations, giving the excitation energies of the system. In the parameterization of the NOTCH method, however, we stop here and set the atomic NOTCH one-electron and two-electron integrals as the one-body and two-body contributions of

\hat{G}

⁠, respectively [the superscript (A) stresses that the integrals are those of neutral atoms; the one-center integrals in charged atoms or molecules will generally have different values, see below]:

\bar{h_{μ_{A} ν_{A}}^{(A)}} = g_{ν_{A}}^{μ_{A}},

(13)

\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} = g_{ν_{A} λ_{A}}^{μ_{A} κ_{A}} .

(14)

The rationale of this construction is that, a multireference configurational interaction singles (MR-CIS) calculation with NOTCH integrals will then closely reproduce the MR-EOM excitation energies (the reproduction is not exact, because the NOTCH MR-CIS calculation does not include Rydberg orbitals). This is corroborated by our benchmark results of atomic excitation energies, as shown in Table I: the NOTCH MR-CIS excitation energies of atoms reproduce the MR-EOM ones (which in turn agrees with the experimental results⁴⁸ to about 0.1 eV) to within a few tenth of an eV, with only one exception (fluorine). Thus, the NOTCH two-electron integrals can be said to implicitly include dynamic correlation with close to MR-CCSD quality. While a slightly better accuracy is achieved by the INDO/S intermediate neglect of differential overlap/spectroscopy^49,50 method, this is mostly due to the fact that the INDO/S integrals were fitted against this exact set of experimental spectra, unlike the NOTCH integrals that were derived from higher-level calculations without any empiricism. Semiempirical methods that were not fitted to excitation energies at all, such as PM3⁸ and OM2,^51,52 perform very poorly and yield mean absolute deviations (MADs) around 2–3 eV. The largest outlier, the 1²P → 1²S excitation energy of the fluorine atom, is predicted to be 17.615 and 9.474 eV too low by PM3 and OM2, respectively. This is probably due to the fact that ground state, fluorine-containing molecules generally do not have a hole in the fluorine 2s orbital, making typical ground state molecule training sets incapable of properly constraining the integrals over 2s orbitals. While PM3 is well-known to be unsuitable for excitation energies in general,⁵³ the failure of OM2 is surprising since random phase approximation (RPA) and multireference configuration interaction (MRCI)⁵⁴ calculations with OM2 integrals are routinely used to calculate the excited states of molecules. This again highlights the importance of deriving the one-center integrals from the MR-EOM Hamiltonian, instead of fitting against reference data.

TABLE I.

Atomic excitation energies (eV) of elements Li–F computed by MR-CIS using NOTCH, INDO/S, PM3 and OM2 integrals, compared with DKH2 MR-EOM-CCSD/ANO-RCC-Full and experimental results.

Atom (ground state)	Excited state	NOTCH	INDO/S	PM3	OM2	MR-EOM	Expt.
Li [1s²2s¹ (²S)]	1s²2p¹ (²P)	1.846	1.850	2.407	N.A.a	1.846	1.848
Be [1s²2s² (¹S)]	1s²2s¹2p¹ (³P)	2.908	3.258	N.A.a	N.A.a	2.866	2.725
Be [1s²2s² (¹S)]	1s²2s¹2p¹ (¹P)	5.961	5.810	N.A.a	N.A.a	5.280	5.277
B [1s²2s²2p¹ (²P)]	1s²2s¹2p² (⁴P)	2.623	4.154	6.998	N.A.a	3.703	3.551
	1s²2s¹2p² (²D)	5.706	6.789	9.988	N.A.a	5.985	5.932
	1s²2s¹2p² (²S)	7.054	8.043	11.758	N.A.a	8.092	7.879
	1s²2s¹2p² (²P)	9.168	9.555	12.428	N.A.a	8.916	8.991
C [1s²2s²2p² (³P)]	1s²2s²2p² (¹D)	1.217	1.082	1.754	1.240	1.148	1.260
	1s²2s²2p² (¹S)	3.042	2.706	3.583	2.265	2.646	2.680
	1s²2s¹2p³ (⁵S)	3.114	5.168	3.847	4.516	4.310	4.179
	1s²2s¹2p³ (³D)	7.737	9.090	8.768	8.806	7.949	7.942
N [1s²2s²2p³ (⁴S)]	1s²2s²2p³ (²D)	2.278	2.325	1.421	2.085	2.251	2.384
	1s²2s²2p³ (²P)	2.934	3.876	2.268	2.891	3.519	3.576
	1s²2s¹2p⁴ (⁴P)	10.731	11.180	8.593	14.026	10.870	10.927
	1s²2s¹2p⁴ (²D)	15.624	15.717	10.677	18.556	14.952	15.027
O [1s²2s²2p⁴ (³P)]	1s²2s²2p⁴ (¹D)	1.808	1.657	1.248	1.540	1.867	1.958
	1s²2s²2p⁴ (¹S)	4.521	4.142	3.094	2.816	4.204	4.180
	1s²2s¹2p⁵ (³P)	16.385	15.896	18.336	20.669	15.546	15.650
F [1s²2s²2p⁵ (²P)]	1s²2s¹2p⁶ (²S)	23.104	20.257	3.506	11.647	21.121	21.039
MDb		0.036	0.289	−0.485 (0.585c)	0.056 (0.923c)
MADb		0.519	0.420	2.861 (1.938c)	2.117 (1.448c)

Atom (ground state)	Excited state	NOTCH	INDO/S	PM3	OM2	MR-EOM	Expt.
Li [1s²2s¹ (²S)]	1s²2p¹ (²P)	1.846	1.850	2.407	N.A.a	1.846	1.848
Be [1s²2s² (¹S)]	1s²2s¹2p¹ (³P)	2.908	3.258	N.A.a	N.A.a	2.866	2.725
Be [1s²2s² (¹S)]	1s²2s¹2p¹ (¹P)	5.961	5.810	N.A.a	N.A.a	5.280	5.277
B [1s²2s²2p¹ (²P)]	1s²2s¹2p² (⁴P)	2.623	4.154	6.998	N.A.a	3.703	3.551
	1s²2s¹2p² (²D)	5.706	6.789	9.988	N.A.a	5.985	5.932
	1s²2s¹2p² (²S)	7.054	8.043	11.758	N.A.a	8.092	7.879
	1s²2s¹2p² (²P)	9.168	9.555	12.428	N.A.a	8.916	8.991
C [1s²2s²2p² (³P)]	1s²2s²2p² (¹D)	1.217	1.082	1.754	1.240	1.148	1.260
	1s²2s²2p² (¹S)	3.042	2.706	3.583	2.265	2.646	2.680
	1s²2s¹2p³ (⁵S)	3.114	5.168	3.847	4.516	4.310	4.179
	1s²2s¹2p³ (³D)	7.737	9.090	8.768	8.806	7.949	7.942
N [1s²2s²2p³ (⁴S)]	1s²2s²2p³ (²D)	2.278	2.325	1.421	2.085	2.251	2.384
	1s²2s²2p³ (²P)	2.934	3.876	2.268	2.891	3.519	3.576
	1s²2s¹2p⁴ (⁴P)	10.731	11.180	8.593	14.026	10.870	10.927
	1s²2s¹2p⁴ (²D)	15.624	15.717	10.677	18.556	14.952	15.027
O [1s²2s²2p⁴ (³P)]	1s²2s²2p⁴ (¹D)	1.808	1.657	1.248	1.540	1.867	1.958
	1s²2s²2p⁴ (¹S)	4.521	4.142	3.094	2.816	4.204	4.180
	1s²2s¹2p⁵ (³P)	16.385	15.896	18.336	20.669	15.546	15.650
F [1s²2s²2p⁵ (²P)]	1s²2s¹2p⁶ (²S)	23.104	20.257	3.506	11.647	21.121	21.039
MDb		0.036	0.289	−0.485 (0.585c)	0.056 (0.923c)
MADb		0.519	0.420	2.861 (1.938c)	2.117 (1.448c)

^a

The semiempirical method is not parameterized for this element.

^b

With respect to the MR-EOM results.

^c

Statistics after the worst outlier, the 1²P → 1²S excitation energy of fluorine, is removed.

The NOTCH integrals of atomic monocations and monoanions,

\bar{h_{μ_{A} ν_{A}}^{(A^{\pm})}}

and

\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A^{\pm})}}

⁠, are similarly calculated and tabulated in advance. By this way, not only the total energies and excitation energies, but also the IPs and EAs of free atoms automatically acquire MR-EOM-like accuracy. The integrals of fractionally and multiply charged atoms can then be obtained by a quadratic interpolation/extrapolation from the neutral, cationic and anionic integrals:

\bar{h_{μ_{A} ν_{A}}^{(A)}} [q_{A}] = \bar{h_{μ_{A} ν_{A}}^{(A)}} + a_{μ_{A} ν_{A}}^{CDI, 1 e (1)} \tilde{q_{A}} + a_{μ_{A} ν_{A}}^{CDI, 1 e (2)} {\tilde{q_{A}}}^{2},

(15)

\tilde{q_{A}} = \frac{4}{π} \arctan q_{A},

(16)

a_{μ_{A} ν_{A}}^{CDI, 1 e (1)} = \frac{1}{2} (\bar{h_{μ_{A} ν_{A}}^{(A^{+})}} - \bar{h_{μ_{A} ν_{A}}^{(A^{-})}}),

(17)

a_{μ_{A} ν_{A}}^{CDI, 1 e (2)} = \frac{1}{2} (\bar{h_{μ_{A} ν_{A}}^{(A^{+})}} + \bar{h_{μ_{A} ν_{A}}^{(A^{-})}}) - \bar{h_{μ_{A} ν_{A}}^{(A)}} .

(18)

Here, the superscript “CDI” stands for “charge dependence of integrals,” and q_A is the atomic charge. The two-electron integrals of charged atoms are derived in a completely analogous manner:

\begin{align} \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} [q_{A}] & = \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} + a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (1)} \tilde{q_{A}} \\ + a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (2)} {\tilde{q_{A}}}^{2}, \end{align}

(19)

a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (1)} = \frac{1}{2} (\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A^{+})}} - \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A^{-})}}),

(20)

\begin{align} a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (2)} & = \frac{1}{2} (\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A^{+})}} + \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A^{-})}}) \\ - \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} . \end{align}

(21)

Note that q_A is “regularized” via Eq. before it is used in the quadratic interpolation/extrapolation [Eqs. and ], so that

\tilde{q_{A}}

is equal to q_A when

q_{A} \in \{- 1, 0, 1\}

⁠, and levels off to ±2 as

|q_{A}|

becomes very large, instead of increasing indefinitely in magnitude and potentially leading to unphysical integral values. In reality we have not observed such unphysical results even if we used q_A instead of

\tilde{q_{A}}

⁠, but we nevertheless chose to use

\tilde{q_{A}}

to ensure the robustness of NOTCH for exotic systems.

E. NOTCH integrals of diatomic molecules

For diatomic molecules, one needs to derive the NOTCH two-center integrals, in addition to the one-center integrals mentioned above. Less obvious is the fact that the NOTCH one-center integrals in a molecule must, in general, be different from the corresponding integral in the free atom. Consider a dihydrogen molecule whose H–H bond length is compressed to zero. The system then becomes a helium atom, whose 1 s orbitals are much tighter than those of the hydrogen atom. Therefore we have

\lim_{R_{A B} \to 0} {(s_{A} s_{A} | s_{A} s_{A})}^{(H_{2})} = {(s_{A} s_{A} | s_{A} s_{A})}^{(H e)} > {(s_{A} s_{A} | s_{A} s_{A})}^{(H)} .

(22)

The most straightforward approach of setting

\bar{{(s_{A} s_{A} | s_{A} s_{A})}^{(H_{2})}} \equiv \bar{{(s_{A} s_{A} | s_{A} s_{A})}^{(H)}}

can thus have an error up to

\bar{{(s_{A} s_{A} | s_{A} s_{A})}^{(H e)}} - \bar{{(s_{A} s_{A} | s_{A} s_{A})}^{(H)}} \approx 0.400 a . u

⁠. While R_AB = 0 is not a realistic situation, the huge error at the R_AB → 0 limit casts doubt on the validity of the approximation

\bar{{(s_{A} s_{A} | s_{A} s_{A})}^{(H_{2})}} \equiv {(s_{A} s_{A} | s_{A} s_{A})}^{(H)}

at non-zero bond lengths. The same situation occurs for the one-electron integrals. Interestingly, none of the most popular ZDO methods seem to have made the one-center integrals dependent on the chemical environment. This can be partly rationalized by the fact that the AOs of ZDO methods correspond more closely to the OAOs (rather than the AOs) of ab initio calculations,⁵⁵ and the molecular OAOs do not reduce to the atomic OAOs as R_AB → 0 (e.g., the OAOs of a dihydrogen molecule do not reduce to an OAO of the helium atom in the R_AB → 0 limit, but rather become linear combinations of the helium 1s and 2p orbitals). However, the environmental dependence of one-center integrals becomes a real issue in TB methods, where integrals are usually given in the AO basis, and indeed the one-center integrals were made to depend on the atomic coordination numbers in, e.g., the xTB family of methods.^14–16

A better alternative for the two-electron integrals of a diatomic molecule AB, would be to add a dynamic correlation correction to the HF/ano-pc-0 integrals of the molecule,

{(μ ν | κ λ)}^{(A B)}

[which differs from

{(μ ν | κ λ)}^{(A)}

due to the geometry dependence of the contraction coefficients of the ano-pc-0 basis set]. Since the ano-pc-0 basis set is generally more compact in the molecular environment than in the free atom, the molecular NOTCH one-center integrals calculated in this way will be generally larger in magnitude than the corresponding free atom NOTCH integrals, as expected. Without loss of generality, we may write

\begin{align} \bar{{(μ ν | κ λ)}^{(A B)}} & = {(μ ν | κ λ)}^{(A B)} + {(μ ν | κ λ)}^{corr, (A)} + {(μ ν | κ λ)}^{corr, (B)} \\ + {(μ ν | κ λ)}^{disp, (A B)} . \end{align}

(23)

Here,

{(μ ν | κ λ)}^{corr, (A)}

and

{(μ ν | κ λ)}^{corr, (B)}

are the local dynamic correlation contributions of atoms A and B to

\bar{{(μ ν | κ λ)}^{(A B)}}

⁠, respectively, and

{(μ ν | κ λ)}^{disp, (A B)}

is a nonlocal term that describes the dispersion effect, as will be detailed below.

The local dynamic correlation contributions

{(μ ν | κ λ)}^{corr, (A)}

⁠, where the AOs μ, ν, κ, λ can be on either atom A or atom B, can in principle be approximated by projecting the AOs onto the AO basis functions of atom A (taking advantage of the fact that the AO basis functions μ_A of a given atom A are orthonormal). Note that we have temporarily used the bra-ket notation

(μ | μ_{A})

for the overlap matrix elements

S_{μ μ_{A}}

to better illustrate how the projection is applied:

{(μ ν | κ λ)}^{corr, (A)} \approx \sum_{μ_{A}} (μ | μ_{A}) {(μ_{A} ν | κ λ)}^{corr, (A)}

\approx \dots \approx \sum_{μ_{A}} \sum_{ν_{A}} \sum_{κ_{A}} \sum_{λ_{A}} (μ | μ_{A}) (ν | ν_{A}) (κ | κ_{A}) (λ | λ_{A}) {(μ_{A} ν_{A} | κ_{A} λ_{A})}^{corr, (A)}

\equiv \sum_{μ_{A} ν_{A} κ_{A} λ_{A}} S_{μ μ_{A}} S_{ν ν_{A}} {(μ_{A} ν_{A} | κ_{A} λ_{A})}^{corr, (A)} S_{κ κ_{A}} S_{λ λ_{A}} .

(24)

Thus, we have reduced the task of finding the local dynamic correlation corrections

{(μ ν | κ λ)}^{corr, (A)}

to finding the one-center corrections

{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{corr, (A)}

⁠, which (without loss of generality) can be expressed as the correlation correction in the free atom times a factor,

\begin{align} {(μ_{A} ν_{A} | κ_{A} λ_{A})}^{corr, (A)} \\ = f_{μ_{A} ν_{A} κ_{A} λ_{A}}^{corr} (\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} - {(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}) . \end{align}

(25)

However, Eq. double-counts the corrections to two-center integrals in the large R_AB limit: since the

f_{μ_{A} ν_{A} κ_{A} λ_{A}}^{corr}

must approach 1 in this limit in order to reproduce the correct free atom results, both

{(μ ν | κ λ)}^{corr, (A)}

and

{(μ ν | κ λ)}^{corr, (B)}

now describe the full correlation, and the total local dynamic correlation correction to the two-center integrals is thus too large by a factor of 2. Thus, in the real working equations, the overlap matrix elements in Eq. are scaled by a factor that approaches 1 when R_AB is small, but smoothly tends to 1/2 when R_AB → +∞, in order to remove this double-counting while not spoiling the small R_AB behavior,

{(μ ν | κ λ)}^{corr, (A)} = \sum_{μ_{A} ν_{A} κ_{A} λ_{A}} {\tilde{S}}_{μ μ_{A}} {\tilde{S}}_{ν ν_{A}} {(μ_{A} ν_{A} | κ_{A} λ_{A})}^{corr, (A)} {\tilde{S}}_{κ κ_{A}} {\tilde{S}}_{λ λ_{A}},

(26)

{\tilde{S}}_{μ_{A} ν_{B}} = S_{μ_{A} ν_{B}} \frac{1 + e^{- α_{corr} {(R_{A B} / r_{A})}^{2}}}{2} .

(27)

Here, α_corr is an empirical constant. r_A is a nonempirical, element-specific parameter that correlates with the atomic radius of the atom A; its construction recipe will be defined later.

The correlation factors

f_{μ_{A} ν_{A} κ_{A} λ_{A}}^{corr}

account for the fact that an atom in a molecule generally feels less dynamic correlation than in the free atom, and are given by

f_{μ_{A} ν_{A} κ_{A} λ_{A}}^{corr} = {(f_{μ_{A}} f_{ν_{A}} f_{κ_{A}} f_{λ_{A}})}^{1 / 4} .

(28)

For the orbital specific correlation factors, we took inspiration from density functional theory and expressed them as,

f_{μ_{A}} = \frac{\int ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) \frac{ρ_{A}^{PModel} (|r - R_{A}|)}{ρ_{A B}^{PModel} (r)} f^{c, P W 92} (ρ_{A B}^{PModel} (r)) ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) d r}{\int ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) f^{c, P W 92} (ρ_{A}^{PModel} (|r - R_{A}|)) ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) d r},

(29)

where

f^{c, P W 92} (ρ) = \frac{\partial^{2} E^{c, P W 92} (ρ)}{\partial ρ^{2}} = \frac{1}{2} (f_{α α}^{c, P W 92} (ρ) + f_{α β}^{c, P W 92} (ρ))

is the spin-unpolarized correlation kernel of the PW92 local density approximation (LDA) correlation functional,⁵⁶ and

ρ_{A B}^{PModel} (r) = ρ_{A}^{PModel} (|r - R_{A}|) + ρ_{B}^{PModel} (|r - R_{B}|)

is an approximation to the molecular density. The numerator of Eq. has a clear physical interpretation: it is the contribution of atom A (as defined via a Hirshfeld partitioning) to the LDA correlation kernel tensor element

f_{μ_{A} μ_{A} μ_{A} μ_{A}}^{c, P W 92} [ρ_{A B}^{PModel} (r)]

in the atom AB. The denominator similarly equals the LDA correlation kernel

f_{μ_{A} μ_{A} μ_{A} μ_{A}}^{c, P W 92} [ρ_{A}^{PModel} (|r - R_{A}|)]

in the free atom A.

Equation is further approximated by spherically averaging the basis functions, e.g.,

ψ_{p_{x}} (r)

is replaced by

\sqrt{ψ_{p_{x}} {(r)}^{2} + ψ_{p_{y}} {(r)}^{2} + ψ_{p_{z}} {(r)}^{2}}

⁠. Now the integrand in the denominator becomes spherically symmetric, and the integrand in the numerator becomes axially symmetric. They are then integrated on a product grid centered on atom A with 1 radial point, 4 theta points and 1 phi point, where the z axis of the product grid is aligned with the A–B bond axis; the latter ensures that the numerically integrated results are rotationally invariant. The radial point of the grid,

r_{μ_{A}}^{(A)}

⁠, is set to the unique value that exactly integrates the ratio

\frac{\int ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) f^{c, P W 92} (ρ_{A}^{PModel} (|r - R_{A}|)) ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) d r}{\int ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) ψ_{μ_{A}} (r) d r}

(30)

and is tabulated in advance. The value

r_{μ_{A}}^{(A)}

is unique because

f^{c, P W 92} (ρ_{A}^{PModel} (|r - R_{A}|))

is a monotonic function of

|r - R_{A}|

⁠; it roughly correlates with the radial extent of the AO

ψ_{μ_{A}} (r)

⁠. The theta points are obtained from the standard Gauss-Legendre quadrature. The use of only one radial point may seem like an extremely crude approximation, but in reality it suffices to converge the equilibrium geometries to about 0.01 Å compared to a more realistic grid comprising 100 radial points and 30 theta points, and it is actually preferable to the latter due to a lack of numerical noise. The use of only 4 theta points similarly suffices to converge the equilibrium geometries to about 0.001 Å.

Finally, the dispersion term

{(μ ν | κ λ)}^{disp}

has the following form:

\begin{align} {(μ ν | κ λ)}^{disp} & = - \sum_{μ_{A} ν_{A} κ_{B} λ_{B}} {\tilde{S}}_{μ μ_{A}} {\tilde{S}}_{ν ν_{A}} \sum_{n = 6, 8} s_{n} \\ \times \frac{η_{A B}^{(n)} {(r_{μ_{A}}^{(A)} r_{ν_{A}}^{(A)} r_{κ_{B}}^{(B)} r_{λ_{B}}^{(B)})}^{\frac{n}{4}}}{R_{A B}^{n - 1} + {(a_{1} {(r_{μ_{A}}^{(A)} r_{ν_{A}}^{(A)} r_{κ_{B}}^{(B)} r_{λ_{B}}^{(B)})}^{\frac{1}{4}} + a_{2})}^{n - 1}} \\ \times (μ_{A} ν_{A} | κ_{B} λ_{B}) {\tilde{S}}_{κ κ_{B}} {\tilde{S}}_{λ λ_{B}} . \end{align}

(31)

The n = 6 and n = 8 terms in Eq. decay asymptotically as

R_{A B}^{- 6}

and

R_{A B}^{- 8}

⁠, respectively [noting that

(μ_{A} ν_{A} | κ_{B} λ_{B})

decays as

R_{A B}^{- 1}

], and are damped to finite values when R_AB is small, in direct analogy with the well-known

R_{A B}^{- 6}

and

R_{A B}^{- 8}

-decaying terms of the Becke–Johnson (BJ) damped dispersion correction.^57–59 Unlike the conventional MM-type dispersion corrections (e.g., DFT-D3⁶⁰), however, Eq. automatically captures the dependence of the dispersion energy with respect to atomic charges, i.e., the magnitude of the dispersion energy between two atoms increases when one of the atoms acquire a negative charge, and decreases when it acquires a positive charge, since the contribution of Eq. to the total energy is proportional to the density matrix elements of the atoms [cf. Eq. ]:

E^{dispcorr} = \frac{1}{2} \sum_{μ ν σ} \sum_{κ λ τ} D_{μ ν σ} ({(μ ν | κ λ)}^{disp} - δ_{σ τ} {(μ κ | ν λ)}^{disp}) D_{κ λ τ} .

(32)

The charge-dependent dispersion correction Eq. can be compared with established charge-dependent dispersion correction schemes, such as DFT-D4.⁶¹ However, unlike the latter, Eq. has the property that when one of the atoms has no electrons, the atom’s contribution to the dispersion energy is exactly zero. The factor

{(r_{μ_{A}}^{(A)} r_{ν_{A}}^{(A)} r_{κ_{B}}^{(B)} r_{λ_{B}}^{(B)})}^{\frac{n}{4}}

in Eq. further makes the core electrons contribute much less to the dispersion energy than the same number of valence electrons do.

The parameters $η_{A B}^{(n)}$ are pre-computed for every pair of atoms AB, so that when the spherically averaged NOTCH ROHF atomic density matrices of atoms A and B are plugged into Eq. (31), the DFT-D3 $C_{A B}^{(6)}$ and $C_{A B}^{(8)}$ coefficients of the atom pair AB are reproduced. The remaining parameters are determined by fitting, except for s₆ which is fixed at the theoretical value of 1.0, as will be detailed in Sec. II I.

Finally, the one-electron integrals are corrected in a similar way,

\begin{align} \bar{h_{μ ν}^{(A B)}} & = h_{μ ν}^{(A B)} + \sum_{μ_{A} ν_{A}} S_{μ μ_{A}} (\bar{h_{μ_{A} ν_{A}}^{(A)}} - h_{μ_{A} ν_{A}}) S_{ν ν_{A}} \\ + \sum_{μ_{B} ν_{B}} S_{μ μ_{B}} (\bar{h_{μ_{B} ν_{B}}^{(B)}} - h_{μ_{B} ν_{B}}) S_{ν ν_{B}} + h_{μ ν}^{pol} . \end{align}

(33)

The terms

\bar{h_{μ_{A} ν_{A}}^{(A)}} - h_{μ_{A} ν_{A}}

and

\bar{h_{μ_{B} ν_{B}}^{(B)}} - h_{μ_{B} ν_{B}}

are included to recover the free atom NOTCH integrals in the R_AB → +∞ limit. The extra correction

h_{μ ν}^{pol}

is detailed in Sec. II G.

F. Charge dependence of the molecular NOTCH integrals

Now that we have the one-electron and two-electron integrals at zero atomic charge, the charge dependence of integrals, Eqs. –, needs to be generalized to the molecular case as well. The first thing to consider is the definition of the atomic charge q_A in molecules, which can, in principle, be obtained from any population analysis method. One choice that we find particularly robust and convenient is the Löwdin charges:

q_{A} = Z_{A} - \sum_{μ_{A} σ} D_{μ_{A} μ_{A} σ}^{OAO} .

(34)

Now we make the approximation that only one-center OAO basis integrals are affected by the atomic charges, while the two-center ones are not. Like the dispersion correction of the two-electron integrals, the free atom charge dependence corrections

\bar{h_{μ_{A} ν_{A}}^{(A)}} [q_{A}] - \bar{h_{μ_{A} ν_{A}}^{(A)}}

and

\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} [q_{A}] - \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}}

also need to be scaled down before they can be added onto the molecular integrals,

\begin{align} \bar{h_{μ_{A} ν_{A}}^{OAO}} [D] & = \bar{h_{μ_{A} ν_{A}}^{OAO}} + f_{A}^{CDI} (\bar{h_{μ_{A} ν_{A}}^{(A)}} [q_{A}] - \bar{h_{μ_{A} ν_{A}}^{(A)}}) \\ + (h_{μ_{A} ν_{A}}^{pol, O A O} [\{q_{B}\}] - h_{μ_{A} ν_{A}}^{pol, O A O}), \end{align}

(35)

\begin{align} \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{OAO}} [D] & = \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{OAO}} \\ + f_{A}^{CDI} (\bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}} [q_{A}] \\ - \bar{{(μ_{A} ν_{A} | κ_{A} λ_{A})}^{(A)}}) . \end{align}

(36)

The last term of Eq. make the one-center integrals depend on the charges of other atoms (⁠

\{q_{B} | B \neq A\}

⁠) too, and will be explained in Sec. II G. The damping factor has the simple form

f_{A}^{CDI} = \prod_{B \neq A} (1 - e^{- β^{CDI} R_{A B} / r_{B}})

(37)

where β^CDI is an empirical parameter, and

r_{A} = \max_{μ_{A} \in A} r_{μ_{A}}^{(A)}

is the radial extent of the most diffuse AO of atom A, which should correlate with the size of the atom itself; note that r_A is also used in Eq. . As discussed before, the charge dependence corrections are designed to mimic the expansion and shrinking of AOs due to non-zero atomic charges, and Eq. accounts for the fact that, when other atoms are bonded to atom A, the AOs of those atoms can provide some of the variational flexibility required to expand or shrink the AOs, so that there is less need to correct the integrals.

G. Basis set incompleteness correction

Even with the above constructions, the NOTCH method still suffers from a systematic overestimation of equilibrium bond length (by 0.03–0.1 Å), as well as IPs and EAs (by 1.0–4.0 eV). This was found to be due to fundamental restrictions of the minimal basis set that are not remedied by the geometry-dependent re-contraction, i.e., (1) the constraint that the radial distributions of the p_x, p_y and p_z orbitals of a given atom must be the same, and (2) the lack of polarization functions. These effects prove to be very difficult to capture without adding more basis functions to the minimal basis set.

To overcome these remaining basis set deficiencies without deviating from a minimal basis set formalism, we first consider the case where the atoms are sufficiently far away from each other, in which case these basis set limitations manifest themselves primarily as a severe underestimation of atomic polarizabilities, a pitfall that is well-known in other semiempirical methods as well. We can thus define a “basis set correction” to the total energy as the difference between the atomic polarization energies with a nearly complete basis set and with the ano-pc-0 basis set,

\begin{gathered} Δ E^{pol} [\{q_{B}\}] = \sum_{A} Δ E_{A}^{pol} [\{q_{B}\}], \\ Δ E_{A}^{pol} [\{q_{B}\}] = - \frac{1}{2} (α_{A}^{def 2 - Q Z V P P D} [\tilde{q_{A}}] - α_{A}^{ano - p c - 0} [\tilde{q_{A}}]) {|F_{A} [\{q_{B}\}]|}^{2} . \end{gathered}

(38)

The charge-dependent polarizabilities

α_{A}^{def 2 - Q Z V P P D} [\tilde{q_{A}}]

and

α_{A}^{ano - p c - 0} [\tilde{q_{A}}]

are obtained from ground state polarizabilities of the neutral and monocationic states of the atom A calculated at the DKH2 SA-CASSCF/def2-QZVPPD^62,63 and DKH2 SA-CASSCF/ano-pc-0 levels of theory, respectively, and interpolated as a linear function of the regularized charge

\tilde{q_{A}}

⁠.

F_{A} [\{q_{B}\}]

is the electric field that all atoms other than A exert on the outermost orbital shell [whose density contribution is approximated by the Fukui function

ρ_{A}^{Fukui} (|r - R_{A}|) = ρ_{A}^{PModel} (|r - R_{A}|) - ρ_{A^{+}}^{PModel} (|r - R_{A}|)

] of atom A,

\begin{align} F_{A} [\{q_{B}\}] & = \sum_{B \neq A} f_{A B}^{pol} \int ρ_{A}^{Fukui} (|r - R_{A}|) \frac{r - r^{'}}{{|r - r^{'}|}^{3}} \\ \times (ρ_{B}^{PModel} [\tilde{q_{B}}] (|r^{'} - R_{B}|) - Z_{B} δ (r^{'} - R_{B})) d r d r^{'} . \end{align}

(39)

Here,

ρ_{B}^{PModel} [\tilde{q_{B}}]

is obtained from a linear interpolation between the fitted densities of the spherically averaged neutral atom (⁠

ρ_{B}^{PModel}

⁠) and monocation (⁠

ρ_{B^{+}}^{PModel}

⁠) of B. Thanks to the ρ^PModel densities being a sum of s-type Gaussian functions, the above integral can be easily evaluated analytically.

To make Eq. applicable to short interatomic distances as well, we have introduced the empirical factor

f_{A B}^{pol}

⁠, which damps

F_{A} [\{q_{B}\}]

when R_AB is small, since the in this case the atom B starts to probe the inner region of the electron cloud of atom A, which is less polarizable than what would be expected from the total atomic polarizability,

f_{A B}^{pol} = 1 - e^{- β^{pol} R_{A B} / r_{A}},

(40)

where β^pol is an empirical parameter. The energy correction is then equally distributed among the valence electrons of A, as a correction to the diagonal one-electron integrals (where

N_{A}^{val}

is the number of valence electrons in the neutral atom of A),

\begin{align} h_{μ_{A} ν_{A}}^{pol, O A O} [\{q_{B}\}] \\ = \{\begin{matrix} δ_{μ_{A} ν_{A}} \frac{Δ E_{A}^{pol} [\{q_{B}\}]}{N_{A}^{val}}, when both μ_{A} and ν_{A} are valence AOs, \\ 0, otherwise . \end{matrix} \end{align}

(41)

H. Left–right correlation correction

The only remaining important contribution to the total energy is the left–right correlation of the bonded electrons, also known as the non-dynamic correlation. It is well known that semilocal DFT correlation functionals mainly describe dynamic correlation and severely underestimate (if not completely miss) the left–right correlation contribution,⁶⁴ especially in elongated bonds. This is one of the major reasons why DFT functionals with 100% HF exchange typically perform poorly. As the NOTCH correlation correction described in Sec. II E uses an LDA-inspired form, it suffers from the same lack of left–right correlation, which is primarily manifested as a systematic overestimation of vibrational frequencies and IPs for certain molecules, like N₂, O₂, and F₂, as well as an erroneous prediction of dipole moments. Therefore, we now add a final correction, $E^{LRC} [D]$ ⁠, to approximately account for left–right correlation.

Recall that the quantity

4 {|D_{μ_{A} ν_{B} σ}^{OAO}|}^{2}

measures the expected number of shared spin-σ electrons between the OAOs μ_A and ν_B, which is the basis of the Wiberg bond order analysis.⁶⁵ Similarly, we argue that the analogous quantity constructed from the uncorrelated OAO two-body density matrix (2-RDM) element,

16 {|Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO}|}^{2}

⁠, where

Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO} = D_{μ_{A} κ_{B} σ}^{OAO} D_{ν_{A} λ_{B} τ}^{OAO} - δ_{σ τ} D_{μ_{A} λ_{B} σ}^{OAO} D_{ν_{A} κ_{B} τ}^{OAO},

(42)

for single reference methods like NOTCH, should represent a generalization of the Wiberg bond order analysis to 2-RDMs and measure the expected number of pairs of shared electrons [note that the pair of electrons, as defined in this section, need not occupy the same spatial orbital; thus, for example, there are 6 × (6 − 1)/2 = 15 pairs of shared electrons in N₂, instead of 3]. Specifically, it is easy to see from Eq. that

16 {|Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO}|}^{2}

is close to 1 when either (1) there is a spin-σ electron shared between OAOs μ_A and κ_B, and another spin-τ electron shared between OAOs ν_A and λ_B (in which case

|D_{μ_{A} κ_{B} σ}^{OAO}| \approx |D_{ν_{A} λ_{B} τ}^{OAO}| \approx \frac{1}{2}

and

D_{μ_{A} λ_{B} σ}^{OAO} = D_{ν_{A} κ_{B} τ}^{OAO} = 0

⁠); (2) there is a spin-σ electron shared between OAOs μ_A and λ_B, and another spin-σ electron shared between OAOs ν_A and κ_B, and σ = τ (in which case,

|D_{μ_{A} λ_{B} σ}^{OAO}| \approx |D_{ν_{A} κ_{B} τ}^{OAO}| \approx \frac{1}{2}

and

D_{μ_{A} κ_{B} σ}^{OAO} = D_{ν_{A} λ_{B} τ}^{OAO} = 0

⁠).

Since the left–right correlation is primarily due to the double excitations of shared electrons, it is natural to approximate the left–right correlation (LRC) energy as proportional to

{|Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO}|}^{2}

⁠, i.e., there is a contribution from every pair of shared electrons (the factor 4 instead of 16 stems from the double-counting of the OAOs in the sum),

E^{LRC} [D] = \sum_{A < B} \sum_{σ τ} \sum_{\begin{array}{c} μ_{A} ν_{B} κ_{A} λ_{B} \\ \in \{v a l e n c e O A O s\} \end{array}} 4 {|Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO}|}^{2} E_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC, (A B)},

(43)

where

E_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC, (A B)}

is given by the configuration interaction energy of a 2 × 2 Hamiltonian including a doubly excited determinant and the ground state determinant,

E_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC, (A B)} = Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC} - \sqrt{{(Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC})}^{2} + {(ξ_{A B}^{LRC})}^{2}} .

(44)

The quantity

Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC}

describes the energy gap between the doubly excited determinant and the ground state determinant. It should reduce to the AO orbital energy difference in the R_AB → +∞ limit,

\begin{align} Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC} (R_{A B} \to + \infty) & = \max (|ϵ_{μ_{A}} + ϵ_{ν_{B}} - ϵ_{κ_{A}} - ϵ_{λ_{B}}|, \\ |ϵ_{μ_{A}} + ϵ_{κ_{A}} - ϵ_{ν_{B}} - ϵ_{λ_{B}}|, |ϵ_{μ_{A}} + ϵ_{λ_{B}} - ϵ_{ν_{B}} - ϵ_{κ_{A}}|), \end{align}

(45)

i.e., the sum of the two higher AO orbital energies minus the sum of the two lower AO orbital energies. The orbital energies

ϵ_{μ_{A}}

⁠,

ϵ_{ν_{B}}

⁠,

ϵ_{κ_{A}}

⁠, and

ϵ_{λ_{B}}

are pre-computed at the ROHF-NOTCH level for spherically averaged atoms. At finite R_AB, the energy gap

Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}

should in general increase as a result of the resonance integral, β_AB, which couples the AOs of the two atoms,

Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC} = \sqrt{{(Δ_{μ_{A} ν_{B} κ_{A} λ_{B}}^{LRC} (R_{A B} \to + \infty))}^{2} + 4 β_{A B}^{2}} .

(46)

The resonance integrals, while usually approximated as proportional to the overlap matrix elements in semiempirical methods, are shown in the OMx methods¹¹ to be more accurately approximated by the following functional form:

β_{A B} = β^{(0)} \sqrt{R_{A B}} e^{- α^{(0)} R_{A B}^{2}} .

(47)

We have thus adopted the same functional form in NOTCH. Here,

β^{(0)}

and

α^{(0)}

are parameters that traditionally depend on the element types of A and B, as well as the types of AOs involved in the bonding, but are treated as global parameters in NOTCH.

Meanwhile, the term

ξ_{A B}^{LRC}

in Eq. describes the off-diagonal matrix element between the ground state and the doubly excited determinant, whose R_AB → +∞ asymptotic behavior is determined by the asymptotic behavior of the MO two-electron integrals (where

ξ_{A B}^{0}

and ζ_AB are constants),

ξ_{A B}^{LRC} = ξ_{A B}^{0} - \frac{ζ_{A B}}{R_{A B}} + O (\frac{1}{R_{A B}^{2}}) .

(48)

Here, the constant term

ξ_{A B}^{0}

and the first-order term

\frac{ζ_{A B}}{R_{A B}}

come from the (AA|AA) and monopolar (AA|BB) integrals, respectively. On the other hand, we would like

ξ_{A B}^{LRC}

to vanish in the R_AB → 0 limit, since in this limit all correlation effects become completely local and are thus already well described by the local correlation corrections described in Sec. II E. These suggest the following form:

ξ_{A B}^{LRC} = ξ_{A B}^{0} - \frac{ζ_{A B}}{R_{A B}} (1 - e^{- ξ_{A B}^{0} R_{A B} / ζ_{A B}}) .

(49)

For homodiatomic molecules, the parameters

ξ_{A B}^{0}

and ζ_AB can be determined by matching the asymptotic behaviors of their restricted HF (RHF) NOTCH and unrestricted HF (UHF) NOTCH dissociation curves up to

O (\frac{1}{R_{A B}})

(for He₂, Be₂, and Ne₂, where there is no shared electron in the R_AB → 0 limit, we set

ξ_{A B}^{0} = ζ_{A B} = 0

⁠). The parameters for a heterodiatomic molecule AB are then taken as the maximum of those of the homodiatomic molecules AA and BB.

We conclude by a few remarks of the functional form Eq. (43):

In the DFT community, the left–right correlation is usually modeled by using a less than 100% amount of HF exchange; the difference of DFT exchange and HF exchange approximates the left–right correlation energy because it makes the exchange–correlation hole more local (as it usually should) than if the functional comprises 100% HF exchange plus a semilocal correlation functional. However, this strategy introduces the notorious self-interaction error problem, since the so-modeled left–right correlation is nonzero even for single-electron densities. The same problem obviously also exists in all self-consistent TB methods. In NOTCH, however, $E^{LRC} [D]$ by construction vanishes for single-electron densities, since the 2-RDM matrix elements $Γ_{μ_{A} κ_{B} σ ν_{A} λ_{B} τ}^{OAO}$ are zero when there is only one electron. The theory thus remains self-interaction free.
Equation (43) is a quartic functional of the density matrix, so it cannot be accurately absorbed into any ansatz where the energy only depends on the density matrix up to third order. The energy expressions of all mainstream ZDO methods are quadratic functionals of the density matrix, while some TB methods include terms cubic in the density matrix (i.e., cubic charge dependence terms of the on-site energy).^15,16,29,30 This in part explains why heavily parameterized core–core repulsion potentials are ubiquitous in existing semiempirical methods because the employed functional forms are not flexible enough to describe the left–right correlation, so that many more empirical parameters have to be introduced to fit the reference data.
The parameterization procedure of $ξ_{A B}^{0}$ and ζ_AB implies that NOTCH can dissociate homodiatomic molecules correctly even under a restricted HF formalism, which even many ab initio methods struggle to achieve. The dissociation is in general not exactly correct for heterodiatomic molecules or diatomic molecular ions but in general still better than if $E^{LRC} [D]$ was not present. While this may not be directly relevant in real-world calculations (since it is usually acceptable to use broken-symmetry wavefunctions, which naturally lead to the correct dissociation limits), it is responsible for curing the systematic overestimation of energy at long but finite bond lengths, leading to better vibrational frequencies at equilibrium geometries (and potentially also the systematic overestimation of barrier height, which is common for methods with 100% HF exchange).
The use of the uncorrelated 2-RDM [Eq. (42)] is crucial here. Take H₂ as an example: The exact correlated 2-RDM element decays to zero in the R_AB → +∞ limit, giving zero correlation correction, while the uncorrelated RHF 2-RDM shows that there is still one pair of shared electrons, and thus correctly yields a correlation correction that restores the UHF dissociation limit (see Fig. S1 for the uncorrelated and correlated 2-RDM elements of H₂ as a function of the bond length). Thus, Eq. (42) should not be understood as an approximation of the correlated 2-RDM. This is reminiscent of the fact that, the Wiberg bond order of, e.g., H₂ calculated from the correlated density matrix goes to zero in the R_AB → +∞ limit and therefore does not reflect the number of shared electrons in the RHF wavefunction.

To illustrate the effect of the LRC correction (especially near the dissociation region), we have plotted the potential energy curves of six representative diatomic molecules (H₂, N₂, F₂, CH, HF, LiF) with NOTCH, both with and without the LRC correction (Fig. 2). As LRC strongly stabilizes the RHF wavefunction in the dissociation region but less so near the equilibrium geometry, it corrects the energy overestimation of the other terms of the NOTCH method in the dissociation region, including the highly ionic molecule LiF, although the correction is sometimes too large (especially for F₂ and CH). Perhaps more importantly, the preferential stabilization of the dissociation region also increases the equilibrium bond length and reduces the vibrational frequency at the equilibrium geometry, which is especially obvious for F₂ and CH. Therefore, despite the fact that the LRC correction does not consistently improve the bond energy and UHF potential energy curves in the dissociation region, we expect that it should in general be beneficial for geometries and vibrational frequencies, as will be confirmed in Secs. III B–III C.

FIG. 2.

View large Download slide

Potential energy curves of selected diatomic molecules, calculated by NOTCH with and without the LRC correction, compared against DKH2 UCCSD(T)/ANO-RCC-TZP results.

Another property of the LRC correction is that the RHF instability point of the NOTCH method occurs later than if the LRC correction is not present, especially in F₂ and CH; this can be compared to the well-known fact that the ground state wavefunction instabilities of DFT methods generally occur at longer bond lengths than the HF method. For H₂, the restricted NOTCH wavefunction even remains stable all the way until dissociation, so that the RHF NOTCH and UHF NOTCH curves coincide. For N₂ and F₂, the LRC correction also makes the RHF NOTCH dissociation limit coincide with the UHF one (by construction), but an RHF instability is still observed when the molecules are not fully dissociated; thus, while the RHF NOTCH curves exhibit unphysical bumps, the UHF NOTCH ones remain qualitatively correct. In the case of the HF molecule, the RHF NOTCH dissociation limit is wrong, but again the UHF NOTCH curve is at least as well behaved as when the LRC correction is not present.

Overall, our preliminary version of NOTCH is a semiempirical method featuring minimal empiricism and explicitly describes much more physics than existing semiempirical methods: (1) the local dynamic correlation and scalar relativistic effects are included through nonempirical MR-EOM calculations (Sec. II D) and the former is properly scaled down in a molecular environment (Sec. II E); (2) the nonlocal component of dynamic correlation, i.e., dispersion, is accounted for as a correction to the (AA|BB)-type two-electron integrals (Sec. II E); (3) the non-dynamic or left–right correlation is described by a separate term inspired by Wiberg bond orders (Sec. II H); (4) the radial distortions of AOs as a function of geometry (Sec. II C) and atomic charge (Sec. II F) are accounted for by corrections to the AO contraction coefficients and one-center integrals, respectively; and (5) the angular distortions of AOs are approximated by adding a classical charge polarization term (Sec. II G). It can be seen that although our method is built upon a HF/minimal basis framework, all relativity, correlation, and basis set incompleteness effects are included in a physically transparent way, and no attempt is made to absorb any effect into terms whose functional forms are not designed to capture that effect.

I. Parameterization

The present form of NOTCH contains eight global empirical parameters: (1) the parameter α_corr responsible for removing the double-counting of dynamic correlation [Eq. (27)]; (2) the dispersion correction parameters s₈, a₁, and a₂ [Eq. (31)]; (3) the parameter β^CDI used for damping the charge dependence correction of integrals [Eq. (37)]; (4) the parameter β^pol used for damping the basis set incompleteness correction [Eq. (40)]; and (5) the parameters $β^{(0)}$ and $α^{(0)}$ involved in the resonance integral in the LRC correction [Eq. (47)]. They were fitted via the following procedure:

The parameters α_corr, s₈, a₁, and a₂ were fitted against the domain-based local pair natural orbital (DLPNO)-coupled-cluster singles and doubles with perturbative triples [CCSD(T)]^66,67 local energy decomposition (LED)^68,69 dispersion energies of all 55 diatomic molecules constructible from the elements H–Ne, at bond distances that are 0.9, 1.0, 1.2, 1.5, and 2.0 times the sum of the respective van der Waals radii (for details, see Appendix B).
The parameters β^pol, $β^{(0)}$ ⁠, and $α^{(0)}$ are obtained by fitting the RHF NOTCH potential energy curve of H₂ against DKH2 CCSD(T)/ANO-RCC-TZP.^45–47
The parameter β^CDI is fitted against the DKH2 HF/ANO-RCC-TZP dissociation energy of H₂⁺.

Since there are some couplings between the parameters, the above procedure is iterated a few times until all the parameters converge. Although the parameters β^pol, $β^{(0)}$ ⁠, $α^{(0)}$ ⁠, and β^CDI are fitted against H₂ and H₂⁺ only, the benchmark results (Sec. III) do not improve noticeably when we add more molecules (such as N₂ and F₂) to the training set, suggesting that the parameters are well transferable among different diatomic molecules.

In addition, NOTCH has many parameters that are not derived from fitting against high-level reference data, including (1) the parameters $c_{μ_{A} ν_{A}}^{(0 - 2)}$ used for calculating the geometry-dependent contraction coefficients of the ano-pc-0 basis set (Sec. II C), fitted from CASSCF calculations of ground state homodiatomic molecules; (2) the Gaussian exponents and coefficients used for expanding the spherically averaged density of atoms, $ρ_{A}^{PModel}$ (Sec. II C), and atomic monocations, $ρ_{A^{+}}^{PModel}$ (Sec. II G), fitted against spherically averaged ROHF atomic densities; (3) the MR-EOM one-electron and two-electron integrals (Sec. II D), generated from atomic MR-EOM calculations; (4) the location of the radial grid point in the numerical integration of Eq. (29) (Sec. II E), fitted from numerically integrating Eq. (30) for spherically averaged atoms; (5) the atomic polarizability correction $α_{A}^{def 2 - Q Z V P P D} [\tilde{q_{A}}] - α_{A}^{ano - p c - 0} [\tilde{q_{A}}]$ (Sec. II G), obtained from polarizability calculations of atoms; and (6) the element-specific parameters $ξ_{A A}^{0}$ and ζ_AA (Sec. II H) involved in the LRC correction, determined by matching the RHF and UHF NOTCH energies for homodiatomic molecules at the dissociation limit. These parameters, as well as detailed descriptions for the procedures to generate them, are tabulated in the supplementary material for convenience. Note that all the nonempirical parameters are derived from the calculations of either atoms or homodiatomic molecules; thus, it is easy to extend the parameterization to heavier elements. While the global, empirical parameters (Table II) may need to be re-fitted when we parameterize more elements, we expect that the reference data of a very few number of heavy element-containing molecules should suffice for the fitting.

TABLE II.

The global NOTCH parameters. All parameters are given in atomic units.

Parameter	α_corr	s₈	a₁	a₂	β^CDI	β^pol	$β^{(0)}$	$α^{(0)}$
Value	0.080	0.711	0.752	5.773	0.793	0.240	0.130	0.094

Parameter	α_corr	s₈	a₁	a₂	β^CDI	β^pol	$β^{(0)}$	$α^{(0)}$
Value	0.080	0.711	0.752	5.773	0.793	0.240	0.130	0.094

Finally, based on the comments of a referee to this paper, we would like to briefly discuss the meaning of the term “empirical” in semiempirical theories. Empiricism, as a deviation from pure ab initio theories, may enter a semiempirical method in various ways:

by fitting physically more or less well-defined parameters that are atom or even atom pair specific to a set of molecular reference data;
by fitting a set of element-specific parameters to experimental data obtained on free atoms and ions;
by replacing well-defined quantities in ab initio theory, for example, analytical AO integrals, with simplified and possibly parameterized expressions—the latter may or may not respect the known asymptotic behavior of the corresponding ab initio quantities; or
by inventing correction terms that have no counterpart in ab initio theory and are based on plausibility arguments.

Based on these levels of empiricism, we wish to comment on the level of empiricism in NOTCH:

In NOTCH, the one-center integrals are derived from high-level calculations on atoms and ions. There is no fitting involved in this process. The parameters are taken directly from the twice similarity-transformed Hamiltonian of MR-EOM calculations in an ANO basis. This is manifestly less empirical than fitting to experimental data of free atoms and ions, let alone the experimental data of molecules. We therefore regard our current approach as nonempirical.
In NOTCH, the “breathing” of basis functions due to the chemical environment is modeled by fitting them to basis functions optimized in an ab initio fashion for homodiatomic molecules. Importantly, the latter basis functions are not optimized to reproduce high-level data obtained either experimentally or computationally (like the case of most semiempirical methods, whose basis set exponents are treated as free parameters and fitted against experimental data) but are rather optimized to minimize the diatomic CASSCF energy computed using that basis set (like the case of any ab initio basis set whose basis functions are optimized variationally on atoms and possibly molecules). Thus, we regard this form of fitting also as nonempirical, even though molecular calculations are involved.
In NOTCH, there are no element-specific empirical parameters that are derived from fitting to atomic or molecular reference data. There are a handful of semiempirical parameters that are fitted to molecular reference data, which is indeed an empirical element of the method. However, the total number of parameters determined in this way (8) is much smaller as for example, in PM7⁹ (225) or GFN2-xTB¹⁶ (194) for the elements H–Ne.
In NOTCH, physically motivated, semiempirical expressions are used to approximately incorporate well-defined physical phenomena that all have a counterpart in ab initio theory. This obviously is a form of empiricism, but here we stress the close ties to ab initio theory. Most of the NOTCH correction terms can be turned on or off with well-defined and physically meaningful changes in the results of the methods (Sec. III); this is in stark contrast with many other semiempirical methods or even some density functionals, where individual energy contributions do not have rigorous ties to any physical quantity, since only the sum of the contributions was constrained by the parameterization process. Our expressions reproduce most of the known asymptotic behaviors of their ab initio counterparts to the extent that semiempirical methods can possibly achieve.

Based on these criteria, we believe that NOTCH is significantly less empirical than practically all semiempirical methods in use, while at least up to this point in the development, NOTCH is characterized by a far more uniform accuracy across a large array of molecular properties, as will be shown in Sec. III. We are hopeful that this behavior will translate to larger molecules.

J. Implementation

The NOTCH method has been implemented in a development version of ORCA.⁷⁰ The program workflow of a NOTCH single point calculation closely resembles a minimal basis HF calculation and is given as follows:

The ano-pc-0 basis set contraction coefficients of the molecule in question, $c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\})$ ⁠, are generated on the fly, via the procedure detailed in Sec. II C.
The zero-charge one-electron and two-electron integrals, $\bar{h_{μ ν}^{(A B)}}$ and $\bar{{(μ ν | κ λ)}^{(A B)}}$ ⁠, are calculated under the AO basis (Sec. II E), transformed to the OAO basis, and stored on disk. The transformation to the OAO basis is not necessary for diatomic calculations (since the Fock builds can be done under the AO basis instead) but may become essential for multiatomic molecules (see Sec. IV).
A standard integral conventional HF calculation is then carried out. During each SCF iteration, the usual one-electron and two-electron components of the Fock matrix, F^OAO, are assembled under the OAO basis using the pre-calculated, zero-charge integrals and the OAO basis density matrix D^OAO. The contributions from the $\sum_{κ λ τ} \frac{\partial \bar{h_{κ λ}} [D]}{\partial D_{μ ν σ}} D_{κ λ τ}$ ⁠, $\sum_{μ^{'} ν^{'} σ^{'} κ λ τ} (\frac{\partial \bar{(μ^{'} ν^{'} | κ λ)} [D]}{\partial D_{μ ν σ}} - δ_{σ^{'} τ} \frac{\partial \bar{(μ^{'} κ | ν^{'} λ)} [D]}{\partial D_{μ ν σ}}) D_{μ^{'} ν^{'} σ^{'}} D_{κ λ τ}$ ⁠, and $\frac{\partial E^{LRC} [D]}{\partial D_{μ ν σ}}$ terms in Eq. (5) are then calculated under the OAO basis and added onto F^OAO. The derivatives of the integrals with respect to the density matrix elements are given by the chain rule, e.g.,
$\frac{\partial \bar{h_{κ λ}^{OAO}} [D^{OAO}]}{\partial D_{μ ν σ}^{OAO}} = \sum_{A} \frac{\partial \bar{h_{κ λ}^{OAO}} [D^{OAO}]}{\partial \tilde{q_{A}}} \frac{\partial \tilde{q_{A}}}{\partial q_{A}} \frac{\partial q_{A}}{\partial D_{μ ν σ}^{OAO}} .$
(50)
It is then elementary to derive the working formulas of the derivatives of the integrals, which are given in the supplementary material.
The contribution of the charge dependence of integrals (Sec. II G) and the LRC correction (Sec. II H) to the energy are then calculated and added to the total energy [Eq. (1)].
The final Fock matrix can be transformed back to the AO basis to yield the Roothaan–Hall equation (3), which is then solved in the standard way to give the molecular orbitals of the next SCF iteration. For large systems, however, it may be advantageous to directly diagonalize the Fock matrix under the OAO basis to save a few matrix operations.

III. RESULTS AND DISCUSSIONS

In this section, we calculate a few properties of atoms and diatomic molecules using NOTCH and some other semiempirical methods (PM3,⁸ PM7,⁹ OM2,^51,52 OM3,⁷¹ GFN-xTB,¹⁵ GFN2-xTB¹⁶), as well as two cost-effective ab initio methods (HF-3c⁴³ and B3LYP-D3^60,72,73/def2-SVP⁶²). All calculations at the NOTCH, PM3, GFN-xTB, GFN2-xTB, HF-3c, and B3LYP-D3/def2-SVP levels were conducted by a development version of ORCA,⁷⁰ while the PM7 and OM2/OM3 calculations were performed with the MOPAC⁷⁴ and MNDO⁷⁵ programs, respectively. The reference values for all ground state properties were generated at the DKH2 CCSD(T)/ANO-RCC-Full (for atomic calculations) and DKH2 CCSD(T)/ANO-RCC-TZP (for molecular calculations) levels of theories. The reference excitation energies were generated at the DKH2 multireference average coupled-pair functional (MRACPF)⁷⁶/ANO-RCC-TZP level of theory; herein MRACPF was used instead of the more rigorous MR-EOM-CCSD method due to the better numerical stability of the former.

A. Atomic IPs and EAs

As a first check, we report the atomic IPs and EAs calculated at the NOTCH level and compare them against the reference results. All IPs and EAs are calculated as differences of single-reference unrestricted SCF energies. As can be seen from Tables III and IV, the MADs of NOTCH atomic IPs and EAs are on the order of 0.4 eV and do not exhibit systematic bias. However, if the neutral atom integrals are used in the calculations of the atomic cations and anions, the MADs increase to 1.63 and 2.42 eV for IPs and EAs, respectively, with the IPs being systematically overestimated and the EAs being systematically underestimated; in particular, all the atomic monoanions are predicted to be unbound. This clearly illustrates the importance of introducing the charge dependence to the integrals. It is interesting to note that none of the other semiempirical methods outperform NOTCH for either the IPs or the EAs, and while the B3LYP-D3/def2-SVP IPs are more accurate than the NOTCH ones, the corresponding EAs are inferior to the NOTCH results. The latter shows that the charge dependence of NOTCH integrals introduces a radial flexibility in the basis set that goes beyond that of even def2-SVP. While the GFN-xTB and GFN2-xTB methods do effectively account for the charge dependence of integrals owing to an energy term that depends cubically on the atomic charges, in practice GFN-xTB and GFN2-xTB perform very poorly for the IPs and EAs of atoms. Indeed, it has been recommended to use a specially parameterized version of GFN-xTB to calculate the IPs and EAs of molecules; even then, an empirical shift on the order of 5 eV must be applied to the IPs and EAs.¹⁵ Also of note is the presence of several severe outliers in various semiempirical methods [e.g., the IPs of Ne (PM7), He, O (GFN-xTB), and the EAs of Li (PM7) and Be (GFN2-xTB)], which clearly represent corner cases that are not found in typical training sets. The nonempirical atomic integrals of NOTCH, however, ensure that NOTCH remains reliable even in these cases. It must be noted, however, that atoms (and many of the diatomic molecules that will be mentioned in the following sections) are not the target of most semiempirical methods.

TABLE III.

Atomic IPs (eV) of elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-Full results.

Element	NOTCH	NOTCHa	PM3	PM7	OM2	OM3	GFN-xTB	GFN2-xTB	HF-3c	B3LYP-D3/def2-SVP	CCSD(T)
H	13.61	13.61	13.07	11.07	12.65	12.46	17.32	17.10	12.68	13.55	13.61
He	24.56	26.09	N.A.b	22.33	N.A.b	N.A.b	43.09	34.27	24.59	24.78	24.56
Li	5.42	5.40	5.13	4.80	N.A.b	N.A.b	11.00	9.42	5.26	5.46	5.39
Be	8.52	8.73	N.A.b	7.84	N.A.b	N.A.b	14.09	N.A.c	7.74	9.00	9.31
B	8.41	8.62	5.32	4.04	N.A.b	N.A.b	13.54	N.A.c	6.88	8.73	8.26
C	11.51	12.24	9.86	8.22	10.01	9.87	17.23	N.A.c	9.87	11.51	11.23
N	14.90	16.41	13.28	14.09	13.63	13.42	19.56	18.91	13.79	14.58	14.52
O	12.17	16.67	12.76	10.54	13.43	13.41	26.17	22.39	12.87	13.91	13.51
F	17.99	20.91	16.39	17.25	17.06	17.29	28.58	24.84	16.90	17.44	17.34
Ne	22.35	25.70	N.A.b	5.30	N.A.b	N.A.b	24.76	31.86	21.47	21.35	21.53
MD	0.02	1.51	−2.17	−3.38	−0.68	−0.75	7.61	6.90	−0.72	0.11
MAD	0.45	1.63	2.17	3.38	0.68	0.75	7.61	6.90	0.72	0.22

Element	NOTCH	NOTCHa	PM3	PM7	OM2	OM3	GFN-xTB	GFN2-xTB	HF-3c	B3LYP-D3/def2-SVP	CCSD(T)
H	13.61	13.61	13.07	11.07	12.65	12.46	17.32	17.10	12.68	13.55	13.61
He	24.56	26.09	N.A.b	22.33	N.A.b	N.A.b	43.09	34.27	24.59	24.78	24.56
Li	5.42	5.40	5.13	4.80	N.A.b	N.A.b	11.00	9.42	5.26	5.46	5.39
Be	8.52	8.73	N.A.b	7.84	N.A.b	N.A.b	14.09	N.A.c	7.74	9.00	9.31
B	8.41	8.62	5.32	4.04	N.A.b	N.A.b	13.54	N.A.c	6.88	8.73	8.26
C	11.51	12.24	9.86	8.22	10.01	9.87	17.23	N.A.c	9.87	11.51	11.23
N	14.90	16.41	13.28	14.09	13.63	13.42	19.56	18.91	13.79	14.58	14.52
O	12.17	16.67	12.76	10.54	13.43	13.41	26.17	22.39	12.87	13.91	13.51
F	17.99	20.91	16.39	17.25	17.06	17.29	28.58	24.84	16.90	17.44	17.34
Ne	22.35	25.70	N.A.b	5.30	N.A.b	N.A.b	24.76	31.86	21.47	21.35	21.53
MD	0.02	1.51	−2.17	−3.38	−0.68	−0.75	7.61	6.90	−0.72	0.11
MAD	0.45	1.63	2.17	3.38	0.68	0.75	7.61	6.90	0.72	0.22

^a

The charge dependence corrections $a_{μ_{A} ν_{A}}^{CDI, 1 e (1)}$ ⁠, $a_{μ_{A} ν_{A}}^{CDI, 1 e (2)}$ ⁠, $a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (1)}$ ⁠, and $a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (2)}$ are set to zero.

^b

The element is not parameterized.

^c

The SCF calculation fails to converge.

TABLE IV.

Atomic EAs (eV) of elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-Full results. The EAs of He and Ne are omitted since their anions are either not calculable (due to the use of a minimal basis set) or unbound at all levels of theories.

Element	NOTCH	NOTCHa	PM3	PM7	OM2	OM3	GFN-xTB	GFN2-xTB	HF-3c	B3LYP-D3/def2-SVP	CCSD(T)
H	0.72	−2.32	−1.72	−3.08	−0.20	−0.39	4.53	5.91	−8.60	−0.42	0.72
Li	0.27	−0.72	−2.17	−14.10	N.A.b	N.A.b	5.41	2.75	−1.64	0.41	0.60
Be	−0.61	−1.79	N.A.b	−3.95	N.A.b	N.A.b	2.54	−14.87	−7.05	−0.82	−0.73
B	0.32	−1.54	−1.89	2.39	N.A.b	N.A.b	4.07	−1.93	−4.70	−0.30	0.14
C	1.43	−0.74	1.70	−0.63	0.79	0.65	4.74	3.37	−4.53	0.38	1.24
N	−1.27	−3.07	0.58	−0.03	−0.74	−0.95	6.30	5.86	−8.21	−1.71	−0.24
O	1.87	−1.64	0.98	0.99	1.22	1.20	9.65	8.20	−6.23	−0.39	1.40
F	4.03	−0.14	2.17	3.36	3.05	3.28	8.85	7.90	−4.78	1.21	3.37
MD	−0.06	−2.42	−0.86	−2.69	−0.47	−0.54	4.95	1.33	−6.53	−1.02
MAD	0.38	2.42	1.36	3.31	0.47	0.54	4.95	6.31	6.53	1.02

Element	NOTCH	NOTCHa	PM3	PM7	OM2	OM3	GFN-xTB	GFN2-xTB	HF-3c	B3LYP-D3/def2-SVP	CCSD(T)
H	0.72	−2.32	−1.72	−3.08	−0.20	−0.39	4.53	5.91	−8.60	−0.42	0.72
Li	0.27	−0.72	−2.17	−14.10	N.A.b	N.A.b	5.41	2.75	−1.64	0.41	0.60
Be	−0.61	−1.79	N.A.b	−3.95	N.A.b	N.A.b	2.54	−14.87	−7.05	−0.82	−0.73
B	0.32	−1.54	−1.89	2.39	N.A.b	N.A.b	4.07	−1.93	−4.70	−0.30	0.14
C	1.43	−0.74	1.70	−0.63	0.79	0.65	4.74	3.37	−4.53	0.38	1.24
N	−1.27	−3.07	0.58	−0.03	−0.74	−0.95	6.30	5.86	−8.21	−1.71	−0.24
O	1.87	−1.64	0.98	0.99	1.22	1.20	9.65	8.20	−6.23	−0.39	1.40
F	4.03	−0.14	2.17	3.36	3.05	3.28	8.85	7.90	−4.78	1.21	3.37
MD	−0.06	−2.42	−0.86	−2.69	−0.47	−0.54	4.95	1.33	−6.53	−1.02
MAD	0.38	2.42	1.36	3.31	0.47	0.54	4.95	6.31	6.53	1.02

^a

The charge dependence corrections $a_{μ_{A} ν_{A}}^{CDI, 1 e (1)}$ ⁠, $a_{μ_{A} ν_{A}}^{CDI, 1 e (2)}$ ⁠, $a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (1)}$ ⁠, and $a_{μ_{A} ν_{A} κ_{A} λ_{A}}^{CDI, 2 e (2)}$ are set to zero.

^b

The element is not parameterized.

B. Diatomic equilibrium bond lengths

Encouraged by the success with atoms, we then move on to diatomic molecules. The first thing to look at is the ground state equilibrium bond lengths and vibrational frequencies of the molecules, which are of paramount importance since geometry optimization and molecular dynamics are perhaps the most common applications of semiempirical methods. As only 55 different diatomic molecules can possibly form between the elements H–Ne, we can do an exhaustive analysis of all these molecules to rule out most qualitative failures of NOTCH on the diatomic level. To the best of our knowledge, such a systematic analysis of diatomic molecules has never been pursued for semiempirical methods.

As can be seen from Fig. 3, the NOTCH equilibrium bond lengths agree very well with the reference values. The largest absolute deviations are on the order of 1 Å and occur in the molecules HeLi and LiNe, but since they are van der Waals complexes with very long interatomic distances [6.043 and 5.112 Å at the CCSD(T) level, respectively; see Table S2 for the complete set of equilibrium geometry benchmark results], such deviations should be tolerable for most purposes. The largest negative deviations are similarly found in van der Waals complexes (BeNe, BNe, CNe), whose equilibrium bond lengths are very long. The largest deviations (>0.05 Å) seen in covalent and ionic bonds are the multiply bonded molecules BeN, NO, and O₂ (which have appreciable multireference character), and the lithium compounds LiH, LiBe, LiB, LiN, LiO, and LiF (which are highly ionic).

FIG. 3.

View large Download slide

Errors of equilibrium bond lengths (Å) of all neutral, ground state diatomic molecules formed between the elements H–Ne (left), as well as excluding Be₂ and all molecules containing He and/or Ne (right), computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. Important outliers are labeled. Method abbreviations: NOTCH a—NOTCH where the geometry dependence of the basis contraction coefficients are neglected, i.e., $c_{μ_{A} ν_{A}}^{(1)} = c_{μ_{A} ν_{A}}^{(2)} = 0$ ⁠; NOTCH b—NOTCH where α_corr is set to 0; NOTCH c—NOTCH where the dispersion correction ${(μ ν | κ λ)}^{disp}$ is set to 0; NOTCH d—NOTCH where the charge dependence of integrals is neglected, i.e., $f_{A}^{CDI}$ is set to 0; NOTCH e—NOTCH where the basis set incompleteness correction is neglected, i.e., $f_{A B}^{pol}$ is set to 0; NOTCH f—NOTCH where the LRC correction $E^{LRC} [D]$ is set to 0; B3LYP—B3LYP-D3/def2-SVP.

To illustrate the effects of some key ingredients of NOTCH, we next perform benchmark studies where certain terms in NOTCH are turned off or certain parameters are set to zero. When the geometry dependence of the ano-pc-0 contraction coefficients is neglected (“NOTCH a”), we observe a systematic overestimation of bond lengths, especially for hydrogen-containing molecules; this is consistent with the well-known fact that the optimal hydrogen 1s exponents of minimal basis sets for molecular calculations are significantly larger than the free atom value 1.0 (typically around 1.2⁴²), while the optimal exponents of other elements differ significantly less from the corresponding free atom values. Meanwhile, the use of a scaled-down version of the overlap matrix elements in Eq. (27) appears essential for the correct description of a few non-covalent complexes containing He and Ne, as otherwise their bond lengths are severely underestimated and sometimes even end up at covalent bond-like values (e.g., the Li–Ne equilibrium bond length is 1.826 Å when α_corr is set to zero, as in “NOTCH b” in Fig. 3, compared to the reference value 5.112 Å). The dispersion correction Eq. (31) is of course essential for the qualitative description of most He- and Ne-containing molecules as well as Be₂, which are under-bound or even unbound when the dispersion correction is turned off (“NOTCH c”); the effect is notable even for some weak covalent bonds, such as Li₂ (Fig. 3). The charge dependence of integrals does not in general affect the bond lengths appreciably but is responsible for preventing the overbinding of some He- and Ne-containing molecules when the dependence is neglected (“NOTCH d”), which can undergo charge transfers on the 0.1e magnitude at compressed geometries and are thus susceptible to the treatment of the charge dependence of integrals. Neglecting the basis set incompleteness correction (“NOTCH e”) results in a systematic overestimation of bond lengths, especially for hydrogen-containing molecules and molecules containing the elements N, O and F. Finally, neglecting the LRC correction (“NOTCH f”) gives systematically underestimated bond lengths for hydrogen-containing molecules and F₂, consistent with the well-known fact that HF tends to underestimate bond lengths due to the lack of left–right correlation. In general, neglecting the various terms in NOTCH results in either worse outliers, wider error spread or increased systematic error, which illustrates the importance of the terms (Fig. 3).

Subsequently, we test the accuracy of some other representative semiempirical methods, as well as HF-3c and B3LYP-D3/def2-SVP (Table S3). All methods, except for the OM2 and OM3 methods parameterized for the elements H, C, N, O, and F only, have severe outliers with qualitatively wrong bond lengths, as shown in Fig. 3 and Table S2. For example, the PM3 method systematically overestimates the bond lengths of Li- and B-containing molecules by as much as 0.452 and 0.210 Å, respectively. The PM7, GFN-xTB, and GFN2-xTB methods, as well as HF-3c and B3LYP-D3/def2-SVP, severely underestimate the bond lengths of several He- and Ne-containing molecules. In particular, the PM7 and GFN2-xTB methods predict that the bond lengths of some He- and Ne-containing molecules are even shorter than typical covalent bonds, for example, HeO [R(GFN2-xTB) = 0.808 Å], BNe [R(GFN2-xTB) = 1.213 Å], HNe [R(PM7) = 0.970 Å], and BeNe [R(PM7) = 0.743 Å]. Even the B3LYP-D3/def2-SVP method, which is widely used in production level geometry optimizations, fails dramatically for LiNe, predicting a bond length of 2.677 Å, which is even shorter than the Li–Li bond length obtained at the same level, 2.766 Å. The widespread failure of semiempirical and cost-effective ab initio methods (especially B3LYP-D3/def2-SVP) for these molecules can be partly attributed to basis set superposition error (BSSE), although the aforementioned failures of GFN2-xTB and PM7 for HeO, BNe, HNe, and BeNe must be primarily due to other deficiencies of the methods, given the magnitudes of the deviations. While the latter deficiencies presumably do not show up in non-covalent complexes in the training sets of the methods, they lead to spurious underestimations of energy near the typical interatomic distance of covalent bonds, a region that is scarcely probed by most He- and/or Ne-containing systems in typical training sets. By contrast, the use of less empirical and more physically grounded constructions in the NOTCH method helped us to avoid unprecedented failures of the method at extremely short bond distances, and the explicit consideration of the charge dependence of integrals (which reflects the radial flexibility of the basis set) makes the NOTCH AOs behave like a much more complete basis set than it would otherwise be (as evidenced by the failure of NOTCH for He- and Ne-containing molecules when $f_{A}^{CDI}$ is set to zero), therefore reducing the BSSE. Finally, even when non-covalent molecules are excluded from the statistics, NOTCH continues to have a narrower error spread than PM3, PM7, GFN-xTB, GFN2-xTB, and HF-3c (Fig. 3, right), which shows that NOTCH is systematically superior for covalently bonded molecules as well. The excellent performance of OM2 and OM3 can be explained by the fact that they are not parameterized for the elements He, Li, Be, B, and Ne, which are generally “difficult” elements for bond lengths, as evidenced by the outliers shown in Fig. 3.

C. Diatomic vibrational frequencies

The errors of harmonic vibrational frequencies (Fig. 4; for the complete set of data, see Table S4) generally mirror the trend of the equilibrium bond lengths, with the frequencies tending to be overestimated whenever the equilibrium bond length is underestimated, and vice versa. However, as NOTCH has a systematic tendency to overestimate vibrational frequencies, it may appear that some of the terms in NOTCH, such as the basis set incompleteness correction $h_{μ ν}^{pol}$ ⁠, actually worsen the agreement of the vibrational frequencies with the reference values (cf. the “NOTCH e” result in Fig. 4, where $h_{μ ν}^{pol}$ is set to 0); nevertheless, the effect is not very large (on the order of 100 cm⁻¹) and may thus be viewed as a necessary compromise for obtaining satisfactory equilibrium bond lengths. Other terms generally improve the accuracy of the vibrational frequencies in most cases. Neglecting the geometry dependence of the AO contraction coefficients (“NOTCH a”) leads to too small vibrational frequencies for H₂ and increases the overestimation of vibrational frequencies for period 2 elements. Setting α_corr to zero (“NOTCH b”) worsens the results of Li- and Be-containing molecules, again reflecting the role of Eq. (27) in preventing the overbinding of weakly bound molecules. Neglecting the charge dependence of integrals (“NOTCH d”) results in too large vibrational frequencies for highly polar molecules like HF and LiF. Finally, the LRC correction (“NOTCH f”) is found to play a major role in the correct prediction of vibrational frequencies of hydrogen-containing molecules as well as O₂ and F₂.

FIG. 4.

Errors of harmonic vibrational frequencies at equilibrium geometries (cm−1) of all neutral, ground state diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. Important outliers are labeled. For the meanings of the method abbreviations, please refer to Fig. 3.

View large Download slide

Errors of harmonic vibrational frequencies at equilibrium geometries (cm⁻¹) of all neutral, ground state diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. Important outliers are labeled. For the meanings of the method abbreviations, please refer to Fig. 3.

Similar to the case of equilibrium geometries, the vibrational frequencies computed by other semiempirical methods and HF-3c (Table S5) invariably show severe outliers. In addition to the outliers that are already expectable from the bad equilibrium geometries (such as the HeO, BNe, HNe, and BeNe molecules mentioned above), there are a few more outliers, like the BN and O₂ molecules at the PM3 level and the H₂, HO, and F₂ molecules at the OM2 level. One may notice that many of these molecules have significant multireference character, such that single reference methods are prone to yield erroneous results. However, the much smaller errors of NOTCH compared to other semiempirical methods as well as the fact that the B3LYP-D3/def2-SVP frequencies of these multireference molecules are not much worse than the single-reference ones suggest that the failure of the existing semiempirical methods is not an intrinsic problem of the single-reference theoretical framework but rather a problem of the functional forms and parameterization procedures of the methods.

D. Diatomic dipole moments

The dipole moments are then benchmarked as a measure of the quality of the semiempirical ground state wavefunctions. We first note that there is some ambiguity in defining the NOTCH dipole moment of a molecule, since the classical polarization term Eq. introduces a contribution to the dipole moment, in addition to the usual nuclear and electronic contributions,

Δ μ^{pol} = \sum_{A} (α_{A}^{def 2 - Q Z V P P D} [\tilde{q_{A}}] - α_{A}^{ano - p c - 0} [\tilde{q_{A}}]) F_{A} [\{q_{B}\}] .

(51)

We have thus computed the NOTCH dipole moment both including and excluding Δμ^pol. As shown in Fig. 5 and Table S6, when Δμ^pol is included (“NOTCH + p”), the dipole moments of most hydrogen-containing molecules are systematically overestimated, while those of Li-containing molecules are underestimated. When Δμ^pol is removed, the errors of H- and Li-containing molecules are noticeably reduced, but those of multiply bonded molecules between B, C, N, and O are somewhat increased; in particular, the dipole moment of CO (−0.355 D) is now predicted to have the wrong sign (C^δ+ − O^δ−), unlike when Δμ^pol is included (0.206 D, compared with the reference value 0.288 D). Whether one should include Δμ^pol in computing the NOTCH dipole moments is as yet unclear and is the subject of future research.

FIG. 5.

View large Download slide

Errors of signed dipole moments (Debye) of all neutral, ground state heterodiatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results (evaluated on the linearized coupled cluster densities). The sign convention is chosen such that the signs of all DKH2 CCSD(T)/ANO-RCC-TZP dipole moments are positive. Important outliers are labeled. “NOTCH + p” stands for NOTCH dipole moments that include the contribution from Eq. (51); other results do not include the contribution from Eq. (51). For the meanings of the other method abbreviations, please refer to Fig. 3.

Despite this, many things can still be said about the influence of the various terms on the dipole moment. The error spread of the dipole moments is worsened when the geometry dependence of the ano-pc-0 contraction coefficients is neglected (“NOTCH a”), suggesting that a correct description of the radial extent of the AOs is crucial for the accuracy of dipole moments (Fig. 5). Neglecting the charge dependence of integrals (“NOTCH d”) improves the dipole moments of some boron compounds but worsens the dipole moments of lithium compounds; as the boron compounds generally have strong multireference character, we believe that the former improvement is fortuitous. The basis set incompleteness correction $h_{μ ν}^{pol}$ is likewise found to slightly deteriorate the dipole moment results for the molecules formed between the elements B, C, and N but improve the results of the hydrides of B, C, N, O, and F (Table S6). Finally, removing the LRC correction (“NOTCH f”) deteriorates the results of molecules like NO, NF, and OF but also improves the dipole moments of some other molecules, notably LiF (Table S6). Overall, many of the terms in NOTCH are not generally beneficial for the accuracy of dipole moments, unlike the case of geometries and vibrational frequencies.

Nevertheless, the NOTCH dipole moments still prove much more robust than the dipole moments calculated by many semiempirical methods and HF-3c (Table S7). The dipole moments of GFN-xTB, GFN2-xTB, and HF-3c have rather large outliers for uncommon bonds, especially Li- and Be-containing molecules; PM7 additionally performs exceptionally poorly for Ne-containing molecules, which can be traced to its strong underestimation of the IP of Ne (Table III). The PM3, OM2, and OM3 methods give accuracies that are similar to or better than NOTCH, but this may partly be attributed to the fact that beryllium is not parameterized in these methods, which seems to be a difficult case for methods that do parameterize it. Finally, B3LYP-D3/def2-SVP gives excellent predictions for most molecules, with the exception of LiBe, LiC, LiN, LiO, and LiF whose dipole moments are significantly underestimated as well as LiB and BeC whose dipole moments are overestimated (obviously owing to their multireference character). The former suggests that the errors of the NOTCH dipole moments of lithium-containing molecules (Table S6), while large in magnitude, are already remarkably small when compared with other methods. Overall, NOTCH has a smaller error range than all the studied semiempirical methods and HF-3c, except for OM2 and OM3, which are only parameterized for H, C, N, O, and F (Fig. 5).

E. Diatomic IPs

Our next target is the IPs of diatomic molecules. As can be seen from Fig. 6 and Table S8, the errors of the NOTCH IPs generally lie within 1.5 eV, but there are a few positive outliers, namely a few molecules formed between B, C, N, and O, as well as a single negative outlier, F₂. All these molecules have significant multireference character, and while one may expect that the LRC correction may play a significant role here, it is evident from Fig. 6 that the LRC correction overcorrects F₂ and under-corrects most of the other multireference molecules. In general, however, the LRC correction proves to be beneficial for reducing the largest IP errors. The geometry dependence of the AO contraction coefficient also helps to reduce the systematic overestimation of the IPs, especially for molecules containing carbon and/or nitrogen. The parameter α_corr and the charge dependence of integrals do not contribute as much to the IPs although the former reduces the IPs of lithium-containing molecules by about 0.1–0.5 eV and thus helps to reduce the error.

FIG. 6.

View large Download slide

Errors of vertical IPs (eV) of all neutral, ground state diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. The DKH2 CCSD(T)/ANO-RCC-TZP equilibrium geometries are used. Important outliers are labeled. For the meanings of the method abbreviations, please refer to Fig. 3.

The performance of other methods is then briefly discussed (Table S9). The GFN-xTB and GFN2-xTB methods drastically overestimate all the IPs by about 2–10 eV with a large scatter, consistent with what is observed in atoms (Table III). The PM7 and HF-3c methods, on the other hand, systematically underestimate the IPs; while the former is well expected from the large negative mean deviation (MD) of the atomic IPs of PM7 (Table III), the latter is better explained by the tendency of HF-3c to significantly underestimate atomic EAs (Table IV), which make the molecular IPs too low when one of the atoms in the molecules carries significant negative charge (e.g., the molecules LiN, LiO and LiF in Table S9). The PM3, OM2, and OM3 methods perform reasonably well and, in particular, are free from the systematic overestimation of the IPs in B-, C-, N-, and O-containing molecules as found in NOTCH. Whether the latter is physically meaningful or due to error cancellation is unclear at the moment. Finally, we note that the B3LYP-D3/def2-SVP method, while performing excellently for all covalently bonded molecules, underestimates the IPs of He₂ and Ne₂ by large margins, despite NOTCH performing very well for these systems; this is clearly due to the self-interaction error of B3LYP, which spuriously stabilizes stretched radical ions like He₂⁺ and Ne₂⁺ and illustrates the importance of using 100% HF exchange in the NOTCH method.

F. Diatomic EAs

The EAs of diatomic molecules computed under different levels of theories are summarized in Fig. 7 and Table S10; note that molecules whose anions are unbound at the reference level of theory are excluded from the analysis. The accuracy of the NOTCH EAs is similar to that of the IPs; there are a few outliers (Li₂, F₂, and boron-containing molecules) that can probably be attributed to some problems of the LRC correction. A few other contributions, such as the geometry dependence of the AO contraction coefficients, the parameter α_corr, and the charge dependence of integrals, prove notably beneficial for the accuracy of the EAs, indicating the crucial role of the radial flexibility of the AOs (and to a lesser extent, the correct treatment of correlation) in the correct prediction of the EAs.

FIG. 7.

View large Download slide

Errors of vertical EAs (eV) of all neutral, ground state diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. The DKH2 CCSD(T)/ANO-RCC-TZP equilibrium geometries are used. Molecules whose DKH2 CCSD(T)/ANO-RCC-TZP EAs are negative are omitted. Important outliers are labeled. For the meanings of the method abbreviations, please refer to Fig. 3.

Most of the other methods studied herein (Table S11) have difficulties predicting the EAs of some of the diatomic molecules. For example, PM3 predicts much too low EAs for LiB, B₂, C₂, and F₂, while PM7 gives systematically too low EAs for most hydrogen- and lithium-containing molecules, except for Li₂ where the EA is significantly overestimated. The OM2 and OM3 methods perform well, with just a single outlier, F₂; however this may be partially due to the fact that the more difficult elements Li, Be, and B are not parameterized in these methods. The EAs are systematically overestimated at the GFN-xTB and GFN2-xTB levels, except for some beryllium-containing molecules, whose EAs are predicted to be too negative at the GFN2-xTB level, in accord with the drastically underestimated EA of the beryllium atom (Table IV). The HF-3c and, to a much less extent, the B3LYP-D3/def2-SVP methods predict systematically too low EAs, which can be obviously attributed to the use of a too small basis set and the lack of corrections for the lack of diffuse functions.

G. Diatomic bond dissociation energies

The final ground state property to look at is the bond dissociation energy. As shown in Fig. 8 and Table S12, NOTCH gives bond dissociation energies that are usually accurate to about 1 eV; the worse outliers (BeF, F₂, BF, N₂) generally possess some multireference character. The terms that affect the NOTCH bond dissociation energies the most are the scaled overlap matrix elements [Eq. (27)] and the LRC correction [Eq. (43)]. When α_corr is set to zero (“NOTCH b”), the bond dissociation energies are overestimated by about 1 eV on average and feature a broader error distribution due to overestimation of dynamic correlation. Interestingly, removing the LRC correction (“NOTCH f”) actually improves the bond dissociation energies (the MAD is reduced from 0.813 to 0.485 eV; see Table S12). The latter suggests that the LRC correction may still have room for improvement.

FIG. 8.

View large Download slide

Errors of bond dissociation energies (eV) of all neutral, ground state diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical and cost-effective ab initio methods, compared with DKH2 CCSD(T)/ANO-RCC-TZP results. The equilibrium geometries optimized at the respective levels of theory are used. Important outliers are labeled. For the meanings of the method abbreviations, please refer to Fig. 3.

The bond dissociation energies of other semiempirical methods and HF-3c (Table S13) have broader error distributions than NOTCH does, again with the exception of OM2 and OM3, which are parameterized for H, C, N, O, and F only. Naturally, the biggest outliers are typically molecules possessing multireference character, like BO, C₂, and N₂. The PM3 and HF-3c methods, bearing 100% HF exchange, tend to underbind these multiply bonded, multireferential molecules, while the GFN-xTB and GFN2-xTB methods, which do not have HF exchange, typically overbind them. The PM7 method manages to achieve less systematic bias throughout the set of molecules studied, probably due to it being more extensively parameterized; however, like what we have seen in Sec. III B, this comes at a price that some non-covalent complexes (HeB, BeNe) are predicted to be extremely stable, covalently bonded molecules, with dissociation energies (9.246 and 8.594 eV, respectively) even exceeding that of N₂ at the same level of theory (8.369 eV). B3LYP-D3/def2-SVP performs excellently for the bond dissociation energies, as expected.

H. Diatomic excitation energies

Finally, we examine the excitation energies of the diatomic molecules. Since the charge dependence of integrals (Sec. II F), the basis set incompleteness correction (Sec. II G) and the LRC correction (Sec. II H) introduce additional response contributions, a straightforward configuration interaction singles (CISs) or RPA calculation using NOTCH integrals derived from ground state atomic charges is not appropriate due to neglecting, e.g., the change of the atomic charges during excitation. Therefore, we calculate the NOTCH excitation energies via the second-order response of the NOTCH energy, giving the Casida equation,⁷⁷

(\begin{matrix} A & B \\ B & A \end{matrix}) (\begin{matrix} X \\ Y \end{matrix}) = ω (\begin{matrix} I & 0 \\ 0 & - I \end{matrix}) (\begin{matrix} X \\ Y \end{matrix}),

(52)

where the orbital Hessians A and B are defined as

A_{i a σ j b τ} = δ_{i j} δ_{a b} δ_{σ τ} (ϵ_{a σ} - ϵ_{i σ}) + K_{i a σ j b τ},

(53)

B_{i a σ j b τ} = K_{i a σ b j τ},

(54)

K_{p q σ r s τ} = \frac{\partial^{2} E^{NOTCH}}{\partial D_{p q σ} \partial D_{r s τ}} .

(55)

The working expressions of the second-order derivatives

\frac{\partial^{2} E^{NOTCH}}{\partial D_{p q σ} \partial D_{r s τ}}

are given in the supplementary material. Note that all the molecular excitation energies in this section are calculated in a single-reference manner since it is nontrivial to extend Eq. to multireference calculations.

The excitation energy benchmark results of NOTCH, INDO/S, PM3, and OM2 are listed in Table V. While NOTCH performs very well for π → π* excitations (e.g., N₂ and CO) and some σ → π* excitations [e.g., Li₂ (³Π_u), HN (³Π), BC (⁴Π), BF, LiC (⁴Π), and BeF (²Π)], giving an error of within 1 eV in most cases, it systematically overestimates the energies of excitations into σ* orbitals [e.g., H₂ (σ → σ*), F₂ (π* → σ*), HN (⁵Σ⁻, σ → σ*), HF (π → σ*), LiF (π* → σ*), BeF (²Σ⁺, σ → σ*), BO (²Π, π → σ*, and ²Σ⁺, σ → σ*), and CN (²Π, π → σ*, and ²Σ⁺, σ → σ*)]. Neglecting the geometry dependence of the ano-pc-0 basis set improves the excitation energies of hydrogen-containing molecules by as much as 5 eV, suggesting that optimizing the AO basis set for the ground state leads to H 1s orbitals that are too compact for expanding the σ* orbital. However, the fact that a sizable error still remains even if the (very diffuse) free atom AOs are used and that the excitation energies of F₂, LiF, BeF, and BO are not much affected by this change suggests that there is probably no minimal basis set that can provide a reasonable compromise between describing the ground state properties and the σ* orbital. Turning off the charge dependence of the integrals does not in general lead to appreciable changes of the excitation energies but changes the excitation energies of highly polar molecules (e.g., HF, LiF) by up to 1.6 eV; while the changes do not always improve the results, it can be said that the effect of the charge dependence of integrals on the excitation energies is large. The LRC correction is occasionally useful in curing the underestimation of triplet excitation energies (N₂), which in extreme cases can lead to triplet instabilities (O₂, F₂), as already discussed in Sec. II H.

TABLE V.

The lowest two vertical excitation energies (eV) of neutral, diatomic molecules formed between the elements H–Ne, computed by NOTCH as well as other semiempirical methods via RPA, compared with DKH2 MRACPF/ANO-RCC-TZP results. Only those molecules that do not contain rare gas atoms, whose ground states are spatially nondegenerate, and whose ground state NOTCH wavefunctions do not have instabilities are listed. The DKH2 CCSD(T)/ANO-RCC-TZP equilibrium geometries are used.

Molecule	State	NOTCH	NOTCH a	NOTCH b	NOTCH c	INDO/S	PM3	OM2	MRACPF
H₂	³Σ_u⁺	18.305	15.458	18.305	15.782	13.700	6.926	5.066	10.585
H₂	¹Σ_u⁺	22.371	17.463	24.670	21.518	18.317	10.006	7.236	12.732
Li₂	³Σ_u⁺	1.432	1.600	1.432	1.212	d	d	e	1.304
Li₂	³Π_u	1.833	1.153	1.833	1.041	d	d	e	1.384
N₂	³Σ_u⁺	7.613	7.820	7.613	6.219	3.266	5.261	5.951	7.569
N₂	³Π_g	7.940	7.220	7.940	7.425	4.979	5.615	6.500	7.935
O₂	³Δ_u	6.503	6.673	6.503	d	2.311	4.641	6.009f	6.145
O₂	³Σ_u⁺	4.045	3.171	5.757	d	2.311	4.694	6.053f	6.292
F₂	³Π_u	6.178	6.134	6.178	d	7.978	3.916	2.298	3.152
F₂	¹Π_u	8.405	8.353	8.405	5.641	9.171	4.206	3.300	4.347
HLi	³Σ⁺	5.394	5.501	5.403	3.936	4.937	3.297	e	3.263
HLi	¹Σ⁺	5.667	5.585	5.740	4.814	5.426	4.863	e	3.629
HN	³Π	4.207	3.632	4.208	3.968	3.976	2.876	3.959f	3.701
HN	⁵Σ^-	17.950	13.783	17.912	16.182	12.266	6.207	7.137f	8.416
HF	³Π	16.357	15.197	17.324	19.024	14.519	8.266	9.667	10.200
HF	¹Π	17.133	16.320	18.790	20.258	15.737	8.632	10.355	10.569
LiC	⁴Π	4.445	4.763	4.311	4.462	3.442	4.102	e	3.254
LiC	⁴Σ^-	5.156	5.490	5.041	5.172	3.718	4.293	e	3.778
LiF	³Π	8.570	8.956	7.557	9.381	9.502	7.854	e	6.734
LiF	¹Π	8.790	9.100	7.809	9.542	9.568	7.987	e	6.778
BeF	²Π	5.508	5.335	5.540	5.179	3.888	e	e	4.119
BeF	²Σ⁺	11.929	12.632	12.761	12.878	10.416	e	e	8.218
BC	⁴Π	0.360	0.316	d	0.577	2.437	3.997	e	1.509
BC	⁴Σ^-	4.430	4.541	d	4.110	2.566	4.419	e	2.260
BO	²Π	8.172	8.003	7.684	8.131	3.339	2.981	e	3.496
BO	²Σ⁺	11.040	9.395	8.181	8.699	4.777	6.366	e	5.561
BF	³Π	2.638	2.521	2.464	2.134	3.216	5.073	e	3.652
BF	¹Π	6.892	6.663	6.825	6.312	6.102	7.383	e	6.521
CN	²Π	4.596	4.649	4.295	4.696	0.886	1.392	2.209f	1.329
CN	²Σ⁺	5.800	5.665	5.454	5.597	1.962	2.170	4.541f	3.195
CO	³Π	5.789	5.269	5.965	5.389	4.975	4.771	5.040	6.191
CO	³Σ⁺	8.808	8.388	8.874	8.400	6.053	6.626	7.165	8.326
NF	³Π	8.664	8.631	8.668	6.996	8.247	5.004	3.803f	7.107
NF	³Π	10.564	10.396	10.586	9.358	10.793	5.085	6.031f	7.264
MDg		2.440	1.919	2.603	2.229	0.842	−0.800	−1.388
MADg		2.723	2.354	2.724	2.608	2.238	1.672	1.649
MDg,h		1.450	1.299	1.348	1.005	0.200	−0.101	−0.711
MADg,h		1.807	1.845	1.504	1.494	1.804	1.228	1.125

Molecule	State	NOTCH	NOTCH a	NOTCH b	NOTCH c	INDO/S	PM3	OM2	MRACPF
H₂	³Σ_u⁺	18.305	15.458	18.305	15.782	13.700	6.926	5.066	10.585
H₂	¹Σ_u⁺	22.371	17.463	24.670	21.518	18.317	10.006	7.236	12.732
Li₂	³Σ_u⁺	1.432	1.600	1.432	1.212	d	d	e	1.304
Li₂	³Π_u	1.833	1.153	1.833	1.041	d	d	e	1.384
N₂	³Σ_u⁺	7.613	7.820	7.613	6.219	3.266	5.261	5.951	7.569
N₂	³Π_g	7.940	7.220	7.940	7.425	4.979	5.615	6.500	7.935
O₂	³Δ_u	6.503	6.673	6.503	d	2.311	4.641	6.009f	6.145
O₂	³Σ_u⁺	4.045	3.171	5.757	d	2.311	4.694	6.053f	6.292
F₂	³Π_u	6.178	6.134	6.178	d	7.978	3.916	2.298	3.152
F₂	¹Π_u	8.405	8.353	8.405	5.641	9.171	4.206	3.300	4.347
HLi	³Σ⁺	5.394	5.501	5.403	3.936	4.937	3.297	e	3.263
HLi	¹Σ⁺	5.667	5.585	5.740	4.814	5.426	4.863	e	3.629
HN	³Π	4.207	3.632	4.208	3.968	3.976	2.876	3.959f	3.701
HN	⁵Σ^-	17.950	13.783	17.912	16.182	12.266	6.207	7.137f	8.416
HF	³Π	16.357	15.197	17.324	19.024	14.519	8.266	9.667	10.200
HF	¹Π	17.133	16.320	18.790	20.258	15.737	8.632	10.355	10.569
LiC	⁴Π	4.445	4.763	4.311	4.462	3.442	4.102	e	3.254
LiC	⁴Σ^-	5.156	5.490	5.041	5.172	3.718	4.293	e	3.778
LiF	³Π	8.570	8.956	7.557	9.381	9.502	7.854	e	6.734
LiF	¹Π	8.790	9.100	7.809	9.542	9.568	7.987	e	6.778
BeF	²Π	5.508	5.335	5.540	5.179	3.888	e	e	4.119
BeF	²Σ⁺	11.929	12.632	12.761	12.878	10.416	e	e	8.218
BC	⁴Π	0.360	0.316	d	0.577	2.437	3.997	e	1.509
BC	⁴Σ^-	4.430	4.541	d	4.110	2.566	4.419	e	2.260
BO	²Π	8.172	8.003	7.684	8.131	3.339	2.981	e	3.496
BO	²Σ⁺	11.040	9.395	8.181	8.699	4.777	6.366	e	5.561
BF	³Π	2.638	2.521	2.464	2.134	3.216	5.073	e	3.652
BF	¹Π	6.892	6.663	6.825	6.312	6.102	7.383	e	6.521
CN	²Π	4.596	4.649	4.295	4.696	0.886	1.392	2.209f	1.329
CN	²Σ⁺	5.800	5.665	5.454	5.597	1.962	2.170	4.541f	3.195
CO	³Π	5.789	5.269	5.965	5.389	4.975	4.771	5.040	6.191
CO	³Σ⁺	8.808	8.388	8.874	8.400	6.053	6.626	7.165	8.326
NF	³Π	8.664	8.631	8.668	6.996	8.247	5.004	3.803f	7.107
NF	³Π	10.564	10.396	10.586	9.358	10.793	5.085	6.031f	7.264
MDg		2.440	1.919	2.603	2.229	0.842	−0.800	−1.388
MADg		2.723	2.354	2.724	2.608	2.238	1.672	1.649
MDg,h		1.450	1.299	1.348	1.005	0.200	−0.101	−0.711
MADg,h		1.807	1.845	1.504	1.494	1.804	1.228	1.125

^a

The geometry dependence of the AO contraction coefficients is neglected, i.e., $c_{μ_{A} ν_{A}}^{(1)} = c_{μ_{A} ν_{A}}^{(2)} = 0$ ⁠.

^b

The charge dependence of integrals is neglected, i.e., $f_{A}^{CDI} = 0$ ⁠.

^c

The LRC correction $E^{LRC} [D]$ is neglected.

^d

The excitation energy cannot be calculated because the molecule has a triplet instability.

^e

The molecule contains an unparameterized element.

^f

A valence space CASSCF calculation is performed instead of RPA since the MNDO program does not support RPA calculations from an unrestricted reference.

^g

Molecules that possess triplet instabilities or unparameterized elements are not included in the statistics.

^h

States whose excitation energies are above 8 eV are excluded in the statistics.

By comparison, the other semiempirical methods INDO/S, PM3 and (to a lesser extent) OM2 generally underestimate the π → π* excitation energies of N₂, O₂, and CO, which is a well-known problem of the ZDO formalism itself. In ZDO methods, the bonding and antibonding MOs formed between two AOs are split symmetrically, while in reality the virtual orbitals are higher than what would be expected from a symmetric splitting;⁵¹ while this is partially remedied by the orthogonalization corrections in the OMx methods,⁵⁴ the corrections are only carried out to low order. Thus, while these methods generally perform better than NOTCH for excitations involving σ* orbitals, it is reasonable to believe that the better performance is due to a fortunate cancellation of the negative errors due to the ZDO formalism and the positive errors due to using a minimal basis set. In some cases, the INDO/S, PM3, and OM2 methods significantly underestimate the excitation energies (e.g., H₂ at the OM2 level, and NF at the PM3 and OM2 levels), sometimes even resulting in triplet instabilities (e.g., Li₂ at the INDO/S and PM3 levels), when one would predict an overestimation of excitation energies due to the use of a minimal basis set. Consequently, despite that the largest outliers occur in NOTCH instead of the other semiempirical methods, we believe that NOTCH should probably be more robust than the other methods for realistically sized organic and organometallic molecules, where the virtual orbitals are less diffuse and much easier to describe by a minimal basis than in diatomic molecules. Moreover, if one desires, it should be relatively easy to improve NOTCH for excited states involving σ* orbitals, since the cause of the overestimation of the excitation energies is well understood.

IV. SUMMARY AND OUTLOOK

The present work describes the preliminary formulation and benchmark results of a new semiempirical method, NOTCH. The method can be considered as a parameterized all-electron HF/minimal basis set calculation, where (contrary to most ZDO methods) all two-center two-electron integrals are retained and all scalar relativity, correlation (local dynamic correlation, dispersion, left–right correlation) and basis set incompleteness (the lack of radial and angular flexibility of the basis set) effects are accounted for by individual, physically motivated corrections.

In our opinion, the theory is based on well-designed and physically appealing approximations and only introduce eight semiempirical parameters which are all global rather than element or even element pair specific. This level of empiricism is even lower than that of many density functionals. Preliminary equilibrium geometry, vibrational frequency, dipole moment, IP, EA, bond dissociation energy, and excitation energy calculations on the atoms and diatomic molecules of the elements H–Ne show that NOTCH is exceptionally robust, with a lack of serious and practically relevant outliers that plague almost all the other semiempirical methods as well as a few cost-effective ab initio methods (e.g., HF-3c), such as the formation of unreasonably short covalent bonds in non-covalent complexes or IPs and EAs that are in error in excess of 5 eV. Meanwhile, the accuracy of NOTCH for the more “conventional” systems, where all methods perform reasonably well, is at least comparable with (and sometimes better than) other semiempirical methods. Moreover, in contrast to other semiempirical methods, NOTCH is balanced in the sense that the quality of the results is consistent across many different properties including the ground state as well as excited states. This feature makes us optimistic that in NOTCH one obtains the right answer for the right reasons. Taken together, the numerical results provide support for the hope that, once fully developed for polyatomic molecules, NOTCH will be useful as a robust, reliable, general-purpose semiempirical method.

While multiple ideas have contributed to the development of NOTCH, we believe that one of the most important ingredients is the parameterization of the one-center one-electron and two-electron integrals on the basis of MR-EOM reference calculations on free atoms and ions, without involving any fitting. The leading idea to approach the design of a semiempirical method from this angle is rooted in the hypothesis that, starting from a good description of the atoms and their excitation and ionization events would be a solid basis for embedding these atoms in a reasonable molecular potential, thus leading to systematically correct electronic structure descriptions. The MR-EOM construction appeared to be particularly promising in this respect given its excellent accuracy relative to experimental data and the fact that the sequence of similarity transformations leads to a form of the theory that can be expressed in terms of “dressed” integrals that can be directly taken as a basis for the parameterization. Thus, MR-EOM does not only provide reliable integrals for atomic calculations, but it also helped us tremendously to avoid large-scale fitting procedures that might have led to overfitting and fortuitously good-looking molecular results.

Insisting on not further manipulating the atomic integrals, we were then forced to identify the real physical reasons for the discrepancies of the calculated results with the reference data, and design physically motivated corrections whose functional forms are genuinely original and have never been proposed or used in existing semiempirical methods. In particular, the geometry dependence of the basis set (Sec. II C), the basis set incompleteness correction (Sec. II G), and the left–right correlation correction (Sec. II H) describe effects that most semiempirical methods do not explicitly consider. In our opinion, this is far more physically satisfying than trying to absorb all of these effects into a large number of parameters that were originally designed to capture other effects. For example, modifying the parameters that define the core–core repulsion potential to absorb, e.g., correlation errors and basis set limitations is particularly dissatisfying to us. While many semiempirical methods incorporate more and more flexible forms of the core–core repulsion potential as they are developed, we have taken the opposite route, i.e., to describe the core–core repulsion exactly using an all-electron approach. Consequently, there is no excuse for adding corrections to the core–core repulsion energy (because it is already exact), and we were forced to develop our other approximations in a way that retains the balance with this “electrostatic cornerstone.” This is clearly a more difficult task, but to us it appears to be more physically appealing than the traditional alternatives.

The obvious next step in the development is the extension of NOTCH to multiatomic molecules. In multiatomic molecules, there are a few issues that are not present in diatomic molecules and that must be properly taken care of:

We would like to omit, to the largest extent possible, three- and four-center two-electron integrals in order to keep the computational performance of NOTCH close to existing semiempirical methods. While simply neglecting these integrals is an unacceptably crude approximation under the AO basis, it is far better justified in the OAO basis.⁷⁸ If all surviving two-electron integrals can be conveniently derived from one-center and two-center two-electron integrals, then for sufficiently large systems, NOTCH will be at most, say, 2–5 times more expensive than existing semiempirical methods. The additional cost stems partly from the inclusion of the (AA|AB)- and (AB|AB)-type integrals and partly because of the explicit treatment of core electrons.
The $C_{A B}^{(6)}$ and $C_{A B}^{(8)}$ coefficients of a pair of atoms depend on the chemical environments of the atoms. It remains to be seen if ideas from established dispersion correction methods, such as expressing $C_{A B}^{(6)}$ and $C_{A B}^{(8)}$ as functions of the atomic coordination numbers (as used in the DFT-D3⁶⁰ and DFT-D4⁶¹ methods), can give useful accuracy in conjunction with our newly proposed dispersion correction formalism, Eq. (31).

Of course, analytic gradients and perhaps also analytic Hessians must be implemented such that geometry optimization, frequency analysis, and molecular dynamics simulations of realistically sized molecules become feasible.

Our hope is to be able to extend NOTCH to heavier elements as well. In particular, open shell transition metals appear to be a major challenge for most existing semiempirical methods. Since NOTCH is based on a formulation in which all terms are in balance by construction, and the exchange is treated in a detailed and physically correct fashion without introducing any self-interaction error, we believe that it holds significant promise in the area. A particularly exciting possibility is to combine NOTCH with a minimal active space CASSCF or MR-CIS approach, which is essential for correctly describing the spin multiplet structure of transition metal complexes. When active spaces spanning more than one atom have to be used, one can turn off the contribution of these atoms to the LRC correction in the multireference NOTCH calculation, to avoid double-counting the non-dynamic correlation between these atoms while still retaining the dynamic correlation contribution. This is a key advantage compared to most existing semiempirical methods, which do not separate the non-dynamic and dynamic contributions to correlation.

Eventually, we believe that NOTCH will become a robust, accurate, general-purpose semiempirical method, which offers intermediate cost and accuracy between existing semiempirical methods and DFT methods that allows for a very wide range of applications in (bio)chemistry, material science, and spectroscopy/photochemistry. Our efforts along this way will be reported in due course.

DEDICATION

This paper is dedicated to the memory of Professor Walter Thiel. Walter was a wonderful and kind friend and mentor. One of us (FN) had the privilege to discuss the very early stages of this project with him shortly before his untimely passing. We deeply regret that he could not see this work to fruition, but we hope that this work may help to reinvigorate interest in semiempirical quantum chemistry, a subject that always remained close to Walter’s heart.

SUPPLEMENTARY MATERIAL

See the supplementary material for further computational details of the benchmark studies, derivations of the NOTCH Fock matrix elements and orbital Hessian matrix elements, raw benchmark data, and lists of the nonempirical NOTCH parameters.

ACKNOWLEDGMENTS

The authors gratefully acknowledge generous financial support of this work by the Max Planck society.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Zikuan Wang: Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Frank Neese: Conceptualization (lead); Funding acquisition (lead); Methodology (equal); Software (supporting); Supervision (lead); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

APPENDIX A: DERIVATION OF THE ano-pc-0 CONTRACTION COEFFICIENTS AS A FUNCTION OF THE MOLECULAR GEOMETRY

As mentioned in the Introduction, the AOs in molecules tend to be more compact than those in the corresponding free atoms. Two possibly predominant reasons are (a) when the length of a bond A–B is sufficiently short, the AOs of atom A can enjoy more nuclear attraction from atom B by shrinking and (b) the AOs of atom A may overlap more favorably with the AOs of atom B after shrinking, resulting in stronger bonding. While (a) does not necessarily dominate over (b), it is much easier to derive a formula for the contribution (a) than for (b). We have therefore decided to use a functional form that is consistent with (a) and contains a few adjustable parameters and fit the parameters against optimized molecular basis set contraction coefficients so that the effects of contribution (b) are effectively accounted for.

Assuming that the density of atom A shrinks uniformly, the new atomic density can be written in terms of the free atom density

ρ_{A}^{PModel} (|r - R_{A}|)

as

ρ_{A} (|r - R_{A}|) = λ^{3} ρ_{A}^{PModel} (λ |r - R_{A}|)

(A1)

for some scale factor λ. We can then find the optimal scale factor λ_opt by variationally optimizing λ to minimize the total energy,

λ_{opt} = \arg \min_{λ} E (λ) .

(A2)

By definition, λ_opt is 1 if A is a free atom. Thus, the deviation of λ_opt from 1 in a molecular environment, which measures the extent of shrinking of the AOs, should be proportional to the derivative of

E (λ)

with respect to λ,

λ_{opt} - 1 \propto {\frac{\partial E (λ)}{\partial λ}|}_{λ = 1} .

(A3)

If we only consider the electrostatic interactions between atoms, then

E (λ)

can be given as

E (λ) = E (λ, f r e e a t o m) + \sum_{B \neq A} \int ρ_{A} (|r - R_{A}|) V_{B}^{PModel} (|r - R_{B}|) d r .

(A4)

Taking the derivative with respect to λ yields

\begin{align} {\frac{\partial E (λ)}{\partial λ}|}_{λ = 1} & = \sum_{B \neq A} (3 \int ρ_{A}^{PModel} (|r - R_{A}|) V_{B}^{PModel} (|r - R_{B}|) d r \\ + \int {\frac{\partial ρ_{A}^{PModel} (r^{'})}{\partial r^{'}}|}_{r^{'} = |r - R_{A}|} V_{B}^{PModel} (|r - R_{B}|) d r) . \end{align}

(A5)

Here, the condition

{\frac{\partial E (λ, f r e e a t o m)}{\partial λ}|}_{λ = 1} = 0

is used, since

E (λ)

is minimized at λ = 1 for a free atom. Finally, we assume that the desired contraction coefficients

c_{μ_{A} ν_{A}}^{ano - p c - 0} (R_{A}, \{R_{B}\})

are linear functions of the scale factor λ_opt, which immediately leads to Eq. . Note that while Eq. implies that

c_{μ_{A} ν_{A}}^{(1)} = 3 c_{μ_{A} ν_{A}}^{(2)}

⁠, in reality the fitted values of

c_{μ_{A} ν_{A}}^{(1)}

and

c_{μ_{A} ν_{A}}^{(2)}

usually do not satisfy this condition, probably due to the influence of the contribution (b) mentioned above. In addition, as we are scaling the contraction coefficients rather than the exponents, the AOs actually do not shrink uniformly as would be implied by Eq. .

APPENDIX B: FITTING THE NOTCH PARAMETERS α_corr, s₈, a₁, AND a₂ AGAINST DLPNO-CCSD(T) LED DISPERSION ENERGIES

The dispersion correction parameters of existing semiempirical methods are usually fitted against non-covalent complexation energies or other properties related to the total energy, such as the equilibrium geometries. In these approaches, the dispersion correction parameters will absorb some of the errors of the intermolecular repulsion term and other inaccuracies of the non-dispersion contribution to the interaction energy. While this obviously helps to improve the accuracy of the method on the training set, it has a few drawbacks: (1) It will hide the fundamental deficiencies of the non-dispersion terms, which may be masked by error compensation in the training set but show up in systems not seen in the training set; and (2) one cannot easily extract the dispersion contribution from the total energy, since only the total energy but not its components are meaningful. Therefore, in NOTCH we would like to fit the dispersion correction parameters to reproduce the dispersion energy itself rather than the total energy.

The LED method^68,69 is a robust and accurate approach for decomposing the DLPNO coupled cluster correlation energy into contributions from dispersion and various non-dispersion contributions. The method takes advantage of the fact that the DLPNO coupled cluster correlation energy can be expressed as a sum of contributions from occupied localized MO (LMO) pairs,

E^{corr} = \sum_{i j} ε_{i j} .

(B1)

Each of the pair energies ɛ_ij can be further decomposed according to which virtual orbitals the electrons in LMOs i and j excite into

ε_{i j} = \sum_{{\tilde{a}}_{i j} {\tilde{b}}_{i j}} (i {\tilde{a}}_{i j} | j {\tilde{b}}_{i j}) {\tilde{τ}}_{{\tilde{a}}_{i j} {\tilde{b}}_{i j}}^{i j},

(B2)

where

{\tilde{τ}}_{{\tilde{a}}_{i j} {\tilde{b}}_{i j}}^{i j}

are contravariant coupled cluster amplitudes. In DLPNO methods, pair natural orbitals (PNOs) are used as the virtual orbitals and the subscripts of

{\tilde{a}}_{i j}

and

{\tilde{b}}_{i j}

stress that different occupied orbital pairs ij have different sets of PNOs.

It is well-known that the correlation energy can be decomposed into a sum of contributions from all double excitations, and since any given LMO or PNO of the molecule can often be unambiguously assigned to a single molecular fragment as long as the fragments are not very strongly interacting, the sum can be regrouped according to in which fragments the occupied and virtual orbitals involved in the double excitations are localized. For example, the dispersion energy between two fragments X and Y can be defined as

\begin{align} E_{X Y}^{disp} & = \sum_{i \in X, j \in Y} (\sum_{{\tilde{a}}_{i j} \in X, {\tilde{b}}_{i j} \in Y} (i {\tilde{a}}_{i j} | j {\tilde{b}}_{i j}) {\tilde{τ}}_{{\tilde{a}}_{i j} {\tilde{b}}_{i j}}^{i j} \\ + \sum_{{\tilde{a}}_{i j} \in Y, {\tilde{b}}_{i j} \in X} (i {\tilde{a}}_{i j} | j {\tilde{b}}_{i j}) {\tilde{τ}}_{{\tilde{a}}_{i j} {\tilde{b}}_{i j}}^{i j}), \end{align}

(B3)

i.e., the contribution of double excitations where one occupied orbital and one virtual orbital are localized on X, and the other two orbitals are localized on Y. Note that the energy partition scheme is more complicated in real DLPNO-CCSD(T) calculations, since (1) the counterpoise correction is applied, (2) some occupied LMO pairs (the so-called “weak pairs”) are treated perturbatively instead of at the coupled cluster level, and (3) the total correlation energy contains contributions from the perturbative triples correction (T). The interested reader is referred to the original references for the details of the DLPNO-CCSD(T) LED approach.^68,69

The DLPNO-CCSD(T) LED dispersion energy $E_{X Y}^{disp}$ represents an accurate metric of the dispersion contribution to the total interaction energy that is not contaminated by non-dispersive interactions. However, its physical meaning is different from the contribution of the NOTCH dispersion correction to the total energy, E^dispcorr [Eq. (32)], because the local dynamic correlation correction (Sec. II E) also describes some of the short-range dispersion. Thus, a naïve approach of fitting E^dispcorr against the reference $E_{X Y}^{disp}$ data would result in a drastic overestimate of the magnitude of dispersion. Instead, we have to define a NOTCH dispersion energy that is consistent with the physical picture of Eq. (B3). Herein, we have chosen the following approach:

Perform a usual NOTCH calculation of the dimer XY.
Perform another NOTCH calculation of XY, but all OAO-basis NOTCH two-electron integrals $\bar{{(μ ν | κ λ)}^{OAO}}$ where μ, ν ∈ X and κ, λ ∈ Y are replaced by their HF values ${(μ ν | κ λ)}^{OAO}$ ⁠. This roughly amounts to selectively turning off the correlation corrections that are responsible for dispersion.
The NOTCH dispersion energy E^disp between the fragments X and Y is defined as the difference between the aforementioned two energies.

There are a few possible variations of this approach, for example, (a) instead of doing a full SCF calculation in step (2), one may also do a single SCF energy evaluation on the converged wavefunction of step (1); and/or (b) one can turn off the correlation under the LMO basis, instead of the OAO basis. The differences between these variations are expected to be small when the interfragment distance is not much shorter than the van der Waals distance.

Thus, we fitted the NOTCH dispersion parameters s₈, a₁, and a₂, as well as the parameter α_corr [which controls the long-range behavior of the local dynamic correlation correction, Eq. (27)], against the diatomic LED dispersion energies at the DKH2 DLPNO-CCSD(T)/ANO-RCC-Full level of theory. When both atoms are open shell, we align their spins in parallel to each other, so as to minimize bonding interactions between the atoms. The resolution of the identity (RI) auxiliary basis sets were generated by the AutoAux method.⁷⁹ All the core electrons were included in the correlation treatment and the TightPNO setting was enforced to obtain more tightly converged dispersion energies. The sum of strong pair and weak pair contributions to the dispersion energy was used as the reference data $E_{X Y}^{disp}$ in the fitting.

REFERENCES

1.

W.

Thiel

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

4

(

2

),

145

–

157

(

2014

).

https://doi.org/10.1002/wcms.1161

Google Scholar

Crossref

2.

N. D.

Yilmazer

and

M.

Korth

,

Comput. Struct. Biotechnol. J.

13

,

169

–

175

(

2015

).

https://doi.org/10.1016/j.csbj.2015.02.004

Google Scholar

Crossref

PubMed

3.

K.

Eichkorn

,

O.

Treutler

,

H.

Öhm

,

M.

Häser

, and

R.

Ahlrichs

,

Chem. Phys. Lett.

240

(

4

),

283

–

290

(

1995

).

https://doi.org/10.1016/0009-2614(95)00621-a

Google Scholar

Crossref

4.

F.

Neese

,

F.

Wennmohs

,

A.

Hansen

, and

U.

Becker

,

Chem. Phys.

356

(

1–3

),

98

–

109

(

2009

).

https://doi.org/10.1016/j.chemphys.2008.10.036

Google Scholar

Crossref

5.

R.

Izsák

and

F.

Neese

,

J. Chem. Phys.

135

(

14

),

144105

(

2011

).

https://doi.org/10.1063/1.3646921

Google Scholar

Crossref

PubMed

6.

R.

Izsák

,

F.

Neese

, and

W.

Klopper

,

J. Chem. Phys.

139

(

9

),

094111

(

2013

).

https://doi.org/10.1063/1.4819264

Google Scholar

Crossref

PubMed

7.

B.

Helmich-Paris

,

B.

de Souza

,

F.

Neese

, and

R.

Izsák

,

J. Chem. Phys.

155

(

10

),

104109

(

2021

).

https://doi.org/10.1063/5.0058766

Google Scholar

Crossref

PubMed

8.

J. J. P.

Stewart

,

J. Comput. Chem.

10

(

2

),

209

–

220

(

1989

).

https://doi.org/10.1002/jcc.540100208

Google Scholar

Crossref

9.

J. J. P.

Stewart

,

J. Mol. Model.

13

(

12

),

1173

–

1213

(

2007

).

https://doi.org/10.1007/s00894-007-0233-4

Google Scholar

Crossref

PubMed

10.

J. J. P.

Stewart

,

J. Mol. Model.

19

(

1

),

1

–

32

(

2013

).

https://doi.org/10.1007/s00894-012-1667-x

Google Scholar

Crossref

PubMed

11.

P. O.

Dral

,

X.

Wu

,

L.

Spörkel

,

A.

Koslowski

,

W.

Weber

,

R.

Steiger

,

M.

Scholten

, and

W.

Thiel

,

J. Chem. Theory Comput.

12

(

3

),

1082

–

1096

(

2016

).

https://doi.org/10.1021/acs.jctc.5b01046

Google Scholar

Crossref

PubMed

12.

Q.

Cui

and

M.

Elstner

,

Phys. Chem. Chem. Phys.

16

(

28

),

14368

–

14377

(

2014

).

https://doi.org/10.1039/c4cp00908h

Google Scholar

Crossref

PubMed

13.

G.

Seifert

and

J.-O.

Joswig

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

2

(

3

),

456

–

465

(

2012

).

https://doi.org/10.1002/wcms.1094

Google Scholar

Crossref

14.

C.

Bannwarth

,

E.

Caldeweyher

,

S.

Ehlert

,

A.

Hansen

,

P.

Pracht

,

J.

Seibert

,

S.

Spicher

, and

S.

Grimme

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

11

(

2

),

e1493

(

2021

).

https://doi.org/10.1002/wcms.1493

Google Scholar

Crossref

15.

S.

Grimme

,

C.

Bannwarth

, and

P.

Shushkov

,

J. Chem. Theory Comput.

13

(

5

),

1989

–

2009

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00118

Google Scholar

Crossref

PubMed

16.

C.

Bannwarth

,

S.

Ehlert

, and

S.

Grimme

,

J. Chem. Theory Comput.

15

(

3

),

1652

–

1671

(

2019

).

https://doi.org/10.1021/acs.jctc.8b01176

Google Scholar

Crossref

PubMed

17.

T.

Husch

,

A. C.

Vaucher

, and

M.

Reiher

,

Int. J. Quantum Chem.

118

(

24

),

e25799

(

2018

).

https://doi.org/10.1002/qua.25799

Google Scholar

Crossref

18.

A. A.

Voityuk

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

3

(

5

),

515

–

527

(

2013

).

https://doi.org/10.1002/wcms.1141

Google Scholar

Crossref

19.

M.

Douglas

and

N. M.

Kroll

,

Ann. Phys.

82

(

1

),

89

–

155

(

1974

).

https://doi.org/10.1016/0003-4916(74)90333-9

Google Scholar

Crossref

20.

G.

Jansen

and

B. A.

Hess

,

Phys. Rev. A

39

(

11

),

6016

–

6017

(

1989

).

https://doi.org/10.1103/physreva.39.6016

Google Scholar

Crossref

21.

L. M. J.

Huntington

,

O.

Demel

, and

M.

Nooijen

,

J. Chem. Theory Comput.

12

(

1

),

114

–

132

(

2016

).

https://doi.org/10.1021/acs.jctc.5b00799

Google Scholar

Crossref

22.

D.

Datta

and

M.

Nooijen

,

J. Chem. Phys.

137

(

20

),

204107

(

2012

).

https://doi.org/10.1063/1.4766361

Google Scholar

Crossref

23.

M.

Nooijen

,

O.

Demel

,

D.

Datta

,

L.

Kong

,

K. R.

Shamasundar

,

V.

Lotrich

,

L. M.

Huntington

, and

F.

Neese

,

J. Chem. Phys.

140

(

8

),

081102

(

2014

).

https://doi.org/10.1063/1.4866795

Google Scholar

Crossref

24.

O.

Demel

,

D.

Datta

, and

M.

Nooijen

,

J. Chem. Phys.

138

(

13

),

134108

(

2013

).

https://doi.org/10.1063/1.4796523

Google Scholar

Crossref

25.

S.

Grimme

and

C.

Bannwarth

,

J. Chem. Phys.

145

(

5

),

054103

(

2016

).

https://doi.org/10.1063/1.4959605

Google Scholar

Crossref

26.

G. M.

Loubriel

and

R. G.

Selsby

,

Int. J. Quantum Chem.

8

(

4

),

547

–

557

(

1974

).

https://doi.org/10.1002/qua.560080406

Google Scholar

Crossref

27.

R.

Iffert

and

K.

Jug

,

Theor. Chim. Acta

72

(

5

),

373

–

378

(

1987

).

https://doi.org/10.1007/bf01192229

Google Scholar

Crossref

28.

J. D.

Baker

and

M. C.

Zerner

,

J. Phys. Chem.

94

(

7

),

2866

–

2872

(

1990

).

https://doi.org/10.1021/j100370a026

Google Scholar

Crossref

29.

M.

Gaus

,

Q.

Cui

, and

M.

Elstner

,

J. Chem. Theory Comput.

7

(

4

),

931

–

948

(

2011

).

https://doi.org/10.1021/ct100684s

Google Scholar

Crossref

30.

Y.

Yang

,

H.

Yu

,

D.

York

,

Q.

Cui

, and

M.

Elstner

,

J. Phys. Chem. A

111

(

42

),

10861

–

10873

(

2007

).

https://doi.org/10.1021/jp074167r

Google Scholar

Crossref

31.

M. J. S.

Dewar

and

N. L.

Hojvat

,

J. Chem. Phys.

34

(

4

),

1232

–

1236

(

1961

).

https://doi.org/10.1063/1.1731723

Google Scholar

Crossref

32.

M. J. S.

Dewar

and

N. L.

Hojvat

,

Proc. R. Soc. Lond. A

264

(

1319

),

431

–

444

(

1961

).

https://doi.org/10.1098/rspa.1961.0209

Google Scholar

Crossref

33.

M. J. S.

Dewar

and

N. L.

Sabelli

,

J. Phys. Chem. A

66

(

12

),

2310

–

2316

(

1962

).

https://doi.org/10.1021/j100818a007

Google Scholar

Crossref

34.

G.

Klopman

,

J. Am. Chem. Soc.

86

(

21

),

4550

–

4557

(

1964

).

https://doi.org/10.1021/ja01075a008

Google Scholar

Crossref

35.

M. J. S.

Dewar

and

W.

Thiel

,

J. Am. Chem. Soc.

99

(

15

),

4899

–

4907

(

1977

).

https://doi.org/10.1021/ja00457a004

Google Scholar

Crossref

36.

D. N.

Laikov

,

J. Chem. Phys.

135

(

13

),

134120

(

2011

).

https://doi.org/10.1063/1.3646498

Google Scholar

Crossref

37.

T. J.

Giese

and

D. M.

York

,

J. Chem. Phys.

123

(

16

),

164108

(

2005

).

https://doi.org/10.1063/1.2080007

Google Scholar

Crossref

38.

T. J.

Giese

and

D. M.

York

,

J. Chem. Phys.

127

(

19

),

194101

(

2007

).

https://doi.org/10.1063/1.2778428

Google Scholar

Crossref

39.

F.

Jensen

,

J. Chem. Phys.

115

(

20

),

9113

–

9125

(

2001

).

https://doi.org/10.1063/1.1413524

Google Scholar

Crossref

40.

F.

Jensen

,

J. Phys. Chem. A

111

(

44

),

11198

–

11204

(

2007

).

https://doi.org/10.1021/jp068677h

Google Scholar

Crossref

41.

F.

Jensen

,

J. Chem. Phys.

117

(

20

),

9234

–

9240

(

2002

).

https://doi.org/10.1063/1.1515484

Google Scholar

Crossref

42.

W. J.

Hehre

,

R. F.

Stewart

, and

J. A.

Pople

,

J. Chem. Phys.

51

(

6

),

2657

–

2664

(

1969

).

https://doi.org/10.1063/1.1672392

Google Scholar

Crossref

43.

R.

Sure

and

S.

Grimme

,

J. Comput. Chem.

34

(

19

),

1672

–

1685

(

2013

).

https://doi.org/10.1002/jcc.23317

Google Scholar

Crossref

44.

E. V. R.

de Castro

and

F. E.

Jorge

,

J. Chem. Phys.

108

(

13

),

5225

–

5229

(

1998

).

https://doi.org/10.1063/1.475959

Google Scholar

Crossref

45.

P.-O.

Widmark

,

P.-Å.

Malmqvist

, and

B. O.

Roos

,

Theor. Chim. Acta

77

(

5

),

291

–

306

(

1990

).

https://doi.org/10.1007/bf01120130

Google Scholar

Crossref

46.

B. O.

Roos

,

V.

Veryazov

, and

P.-O.

Widmark

,

Theor. Chem. Acc.

111

(

2

),

345

–

351

(

2004

).

https://doi.org/10.1007/s00214-003-0537-0

Google Scholar

Crossref

47.

B. O.

Roos

,

R.

Lindh

,

P.-Å.

Malmqvist

,

V.

Veryazov

, and

P.-O.

Widmark

,

J. Phys. Chem. A

108

(

15

),

2851

–

2858

(

2004

).

https://doi.org/10.1021/jp031064+

Google Scholar

Crossref

48.

A.

Kramida

,

Y.

Ralchenko

,

J.

Reader

, and

NIST ASD Team

, NIST Atomic Spectra Database (version 5.10),

National Institute of Standards and Technology

,

Gaithersburg, MD

,

2022

, available at: https://physics.nist.gov/asd.

49.

J. A.

Pople

,

D. L.

Beveridge

, and

P. A.

Dobosh

,

J. Chem. Phys.

47

(

6

),

2026

–

2033

(

1967

).

https://doi.org/10.1063/1.1712233

Google Scholar

Crossref

50.

J.

Ridley

and

M.

Zerner

,

Theor. Chim. Acta

32

(

2

),

111

–

134

(

1973

).

https://doi.org/10.1007/bf00528484

Google Scholar

Crossref

51.

W.

Weber

and

W.

Thiel

,

Theor. Chem. Acc.

103

(

6

),

495

–

506

(

2000

).

https://doi.org/10.1007/s002149900083

Google Scholar

Crossref

52.

W.

Weber

, “

Ein neues semiempirisches NDDO-Verfahren mit Orthogonalisierungskorrekturen: Entwicklung des Modells, Parametrisierung und Anwendungen

,” Ph.D. thesis,

Universität Zürich

,

1996

.

Google Scholar

53.

M. R.

Silva-Junior

and

W.

Thiel

,

J. Chem. Theory Comput.

6

(

5

),

1546

–

1564

(

2010

).

https://doi.org/10.1021/ct100030j

Google Scholar

Crossref

54.

D.

Tuna

,

Y.

Lu

,

A.

Koslowski

, and

W.

Thiel

,

J. Chem. Theory Comput.

12

(

9

),

4400

–

4422

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00403

Google Scholar

Crossref

55.

X.

Wu

,

P. O.

Dral

,

A.

Koslowski

, and

W.

Thiel

,

J. Comput. Chem.

40

(

4

),

638

–

649

(

2019

).

https://doi.org/10.1002/jcc.25748

Google Scholar

Crossref

56.

J. P.

Perdew

and

Y.

Wang

,

Phys. Rev. B

45

(

23

),

13244

(

1992

).

https://doi.org/10.1103/physrevb.45.13244

Google Scholar

Crossref

57.

A. D.

Becke

and

E. R.

Johnson

,

J. Chem. Phys.

123

(

15

),

154101

(

2005

).

https://doi.org/10.1063/1.2065267

Google Scholar

Crossref

58.

E. R.

Johnson

and

A. D.

Becke

,

J. Chem. Phys.

123

(

2

),

024101

(

2005

).

https://doi.org/10.1063/1.1949201

Google Scholar

Crossref

59.

E. R.

Johnson

and

A. D.

Becke

,

J. Chem. Phys.

124

(

17

),

174104

(

2006

).

https://doi.org/10.1063/1.2190220

Google Scholar

Crossref

60.

S.

Grimme

,

J.

Antony

,

S.

Ehrlich

, and

H.

Krieg

,

J. Chem. Phys.

132

(

15

),

154104

(

2010

).

https://doi.org/10.1063/1.3382344

Google Scholar

Crossref

61.

E.

Caldeweyher

,

C.

Bannwarth

, and

S.

Grimme

,

J. Chem. Phys.

147

(

3

),

034112

(

2017

).

https://doi.org/10.1063/1.4993215

Google Scholar

Crossref

62.

F.

Weigend

and

R.

Ahlrichs

,

Phys. Chem. Chem. Phys.

7

(

18

),

3297

–

3305

(

2005

).

https://doi.org/10.1039/b508541a

Google Scholar

Crossref

63.

D.

Rappoport

and

F.

Furche

,

J. Chem. Phys.

133

(

13

),

134105

(

2010

).

https://doi.org/10.1063/1.3484283

Google Scholar

Crossref

64.

N. C.

Handy

and

A. J.

Cohen

,

Mol. Phys.

99

(

5

),

403

–

412

(

2001

).

https://doi.org/10.1080/00268970010018431

Google Scholar

Crossref

65.

K. B.

Wiberg

,

Tetrahedron

24

(

3

),

1083

–

1096

(

1968

).

https://doi.org/10.1016/0040-4020(68)88057-3

Google Scholar

Crossref

66.

C.

Riplinger

,

B.

Sandhoefer

,

A.

Hansen

, and

F.

Neese

,

J. Chem. Phys.

139

(

13

),

134101

(

2013

).

https://doi.org/10.1063/1.4821834

Google Scholar

Crossref

67.

C.

Riplinger

,

P.

Pinski

,

U.

Becker

,

E. F.

Valeev

, and

F.

Neese

,

J. Chem. Phys.

144

(

2

),

024109

(

2016

).

https://doi.org/10.1063/1.4939030

Google Scholar

Crossref

68.

A.

Altun

,

F.

Neese

, and

G.

Bistoni

,

J. Chem. Theory Comput.

15

(

1

),

215

–

228

(

2019

).

https://doi.org/10.1021/acs.jctc.8b00915

Google Scholar

Crossref

69.

G.

Bistoni

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

10

(

3

),

e1442

(

2020

).

https://doi.org/10.1002/wcms.1442

Google Scholar

Crossref

70.

F.

Neese

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

12

(

5

),

e1606

(

2022

).

https://doi.org/10.1002/wcms.1606

Google Scholar

Crossref

71.

M.

Scholten

, “

Semiempirische Verfahren mit Orthogonalisierungskorrekturen: Die OM3 Methode

,” Ph.D. thesis,

Universität Düsseldorf

,

2003

.

Google Scholar

72.

A. D.

Becke

,

Phys. Rev. A

38

(

6

),

3098

(

1988

).

https://doi.org/10.1103/physreva.38.3098

Google Scholar

Crossref

73.

C.

Lee

,

W.

Yang

, and

R. G.

Parr

,

Phys. Rev. B

37

(

2

),

785

(

1988

).

https://doi.org/10.1103/physrevb.37.785

Google Scholar

Crossref

74.

J. J. P.

Stewart

,

A.

Klamt

,

W.

Thiel

,

D.

Danovich

,

G. B.

Rocha

,

R. L.

Gieseking

,

J. E.

Moussa

,

H. A.

Kurtz

,

P.

Korambath

,

K. M.

Merz

, and

B.

Wang

(

2022

). “

MOPAC (22.0.0)

,” Zenodo. https://doi.org/10.5281/zenodo.6511959

Google Scholar

75.

W.

Thiel

, Program MNDO2020,

Muelheim

,

2020

, https://mndo.kofo.mpg.de/.

76.

R. J.

Gdanitz

and

R.

Ahlrichs

,

Chem. Phys. Lett.

143

(

5

),

413

–

420

(

1988

).

https://doi.org/10.1016/0009-2614(88)87388-3

Google Scholar

Crossref

77.

M. E.

Casida

,

Recent Advances in Density Functional Methods: (Part I)

(

World Scientific

,

1995

), pp.

155

–

192

.

Google Scholar

Crossref

78.

T.

Husch

and

M.

Reiher

,

J. Chem. Theory Comput.

14

(

10

),

5169

–

5179

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00601

Google Scholar

Crossref

79.

G. L.

Stoychev

,

A. A.

Auer

, and

F.

Neese

,

J. Chem. Theory Comput.

13

(

2

),

554

–

562

(

2017

).

https://doi.org/10.1021/acs.jctc.6b01041

Google Scholar

Crossref

Development of NOTCH, an all-electron, beyond-NDDO semiempirical method: Application to diatomic molecules

I. INTRODUCTION

II. THEORY

A. Notations

B. Overview

C. Basis set

D. NOTCH integrals of free atoms and atomic ions

E. NOTCH integrals of diatomic molecules

F. Charge dependence of the molecular NOTCH integrals

G. Basis set incompleteness correction

H. Left–right correlation correction

I. Parameterization

J. Implementation

III. RESULTS AND DISCUSSIONS

A. Atomic IPs and EAs

B. Diatomic equilibrium bond lengths

C. Diatomic vibrational frequencies

D. Diatomic dipole moments

E. Diatomic IPs

F. Diatomic EAs

G. Diatomic bond dissociation energies

H. Diatomic excitation energies

IV. SUMMARY AND OUTLOOK

DEDICATION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DERIVATION OF THE ano-pc-0 CONTRACTION COEFFICIENTS AS A FUNCTION OF THE MOLECULAR GEOMETRY

APPENDIX B: FITTING THE NOTCH PARAMETERS αcorr, s8, a1, AND a2 AGAINST DLPNO-CCSD(T) LED DISPERSION ENERGIES

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX B: FITTING THE NOTCH PARAMETERS α_corr, s₈, a₁, AND a₂ AGAINST DLPNO-CCSD(T) LED DISPERSION ENERGIES