{'fraction_non_alphanumeric': 0.037941979774119734, 'fraction_numerical': 0.020890234000147632, 'mean_word_length': 4.448421475970239, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comprehensive Analysis of the Neglect of Diatomic Differential Overlap Approximation

Many modern semiempirical molecular orbital models are built on the neglect of diatomic differential overlap (NDDO) approximation. An in-depth understanding of this approximation is therefore indispensable to rationalize the success of these semiempirical molecular orbital models and to develop further improvements on them. The NDDO approximation provides a recipe to approximate electron-electron repulsion integrals (ERIs) in a symmetrically orthogonalized basis based on a far smaller number of ERIs in a locally orthogonalized basis. We first analyze the NDDO approximation by comparing ERIs in both bases for a selection of molecules and for a selection of basis sets. We find that the errors in Hartree-Fock and second-order Moller-Plesset perturbation theory energies grow roughly linearly with the number of basis functions. We then examine different approaches to correct for the errors caused by the NDDO approximation and propose a strategy to directly correct for them in the two-electron matrices that enter the Fock operator.


Introduction
The neglect of diatomic differential overlap (NDDO) approximation 1 is the foundation for many modern semiempirical molecular orbital (SEMO) models such as the PMx (x = 3, 6, 7) 2-4 and OMx (x = 1, 2, 3) 5-7 models. These SEMO models are chosen when, on the one hand, accurate ab initio electronic structure models are computationally unfeasible, but when, on the other hand, a calculation with an electronic structure model is favored over a classical force field to exploit the first-principles nature of the fundamental electrostatic interactions. Examples include the simulation of very large systems such as proteins, [8][9][10][11][12][13][14][15] virtual high-throughput screening schemes for materials and drug discovery, [16][17][18][19][20][21][22][23] or real-time interactive quantum chemical calculations. [24][25][26][27][28][29][30] Originally, the NDDO approximation was conceptualized as a means to reduce the number of electron-electron repulsion integrals (ERIs) which need to be explicitly calculated in the course of a Hartree-Fock (HF) calculation. 1 The NDDO approximation specifies how ERIs in a symmetrically orthogonalized basis may be approximated based on a far smaller number of ERIs in a locally orthogonalized basis. However, the NDDO approximation has not found acceptance in ab initio calculations due to the significant errors which it introduced in ERIs in the symmetrically orthogonalized basis. [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] Instead, the NDDO approximation has become popular in SEMO models where it is combined 50 with various other approximations made to the oneelectron matrix and to the nucleus-nucleus repulsion energy which benefit from mutual error compensation. Obviously, the dramatically reduced computational timings of NDDO-SEMO models come at a price: they are known to be notoriously unreliable. 51,52 For further improvements to NDDO-SEMO models, an in-depth analysis of the NDDO approximation is mandatory. De-spite decades of work on NDDO-SEMO models, a fully satisfactory analysis has not been provided yet. We intend to close this gap here and study the foundations of the NDDO approximation from a state-of-the-art perspective.
First, we determine how the NDDO approximation affects ERIs evaluated in a symmetrically orthogonalized basis. Previous analyses of the NDDO approximation were limited to few tens of molecules which consisted of few atoms (usually less than four heavy atoms). [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] In this work, we consider a diverse selection of molecules that are also much larger. As the errors in the ERIs will propagate to all quantities calculated in an NDDO framework, we study how the NDDO approximation affects the (absolute and relative) HF and second-order Møller-Plesset perturbation (MP2) theory energies. In this context, we examine different basis set choices.
In general, the NDDO approximation is only valid for a locally orthogonal basis set, [44][45][46]48 which appears to restrict contemporary NDDO-SEMO models to a minimal basis set. A minimal basis set, however, is generally unsuitable for the description of atoms in molecules because it does not yield reliable relative energies, force constants, electric dipole moments, static dipole polarizabilities, and other properties. [53][54][55][56][57][58] It was suggested to generalize NDDO-SEMO models to larger, e.g., double-zeta split-valence basis sets to obtain more accurate results. 59 Two studies examined 60, 61 the effects of the application of a double-zeta split-valence basis set in conjunction with the NDDO approximation, but came to the conclusion that, contrary to what one would expect, the results did not improve compared to the results obtained with a single-zeta basis set. In this work, we dissect in detail the origins of this counterintuitive observation and elaborate on a solution.
It does not come as a surprise that we find -in agreement with previous results 32-49 -the NDDO approximation to cause severe errors. We therefore examine how one can correct for these errors. We briefly review the error compensation strategy which contemporary NDDO-SEMO models apply and then propose a way to directly correct for errors in the two-electron matrices.
We show that our approach allows for rapid extended-basis set calculations invoking the NDDO approximation with error control.
2 Neglect of Diatomic Differential Overlap Approximation

Basic Notation
In a basis-set representation, each spatial molecular orbital ψ i = ψ i (r) is approximated as a linear combination of pre-defined basis functions. Following the well-established approach for finite systems in molecular physics, we choose the basis functions to be Gaussian-type atom-centered functions χ I µ = χ I µ (r). Our notation indicates that the µ-th basis function of type χ is centered on atom I. Additionally, we require the basis functions χ to be locally orthogonal which means that the overlap χ S µν between the basis functions χ I µ and χ J ν must fulfill, A molecular orbital ψ i is then given as the sum of the M basis functions χ I µ weighted with expansion coefficients χ C = { χ C µi }, Throughout this work, we denote the bases by a left superscript, i.e., ' χ C'.
In the χ-basis, the Roothaan-Hall equation in the spin-restricted formulation then reads 54 where χ F is the Fock matrix, which depends on χ C, χ S is the overlap matrix, and the diagonal matrix of orbital energies. As is invariant under unitary matrix transformations with which we may transform one basis into another one, it does not carry a left superscript. The evaluation of the Fock-matrix elements in the χ-basis, 54 requires the evaluation of one-electron integrals χ I µ |h|χ J ν , elements of the density matrix χ P ( χ C), and ERIs in the χ-basis ( χ ERIs). In a spin-restricted framework, the elements of χ P ( χ C) are given by where n is the number of electrons and real expansion coefficients, and the χ ERIs are calculated according to For the following discussion, it is convenient to divide the Fock matrix into a one-electron matrix χ H and a two-electron matrix χ G( χ C) which, in HF theory, consists of a Coulomb matrix χ J ( χ C) and an exchange matrix χ K( χ C), After reaching self-consistency, the total electronic HF energy of the system is calculated from χ P ( χ C), χ F ( χ C), and the nucleus-nucleus repulsion energy V , We need to introduce a second basis, the symmetrically orthogonalized 62 basis φ, to discuss the NDDO approximation. The symmetrically orthogonalized basis functions φ = {φ µ } and the locally orthogonal basis functions Consequently, we can calculate the ERIs in the φ-basis ( φ ERIs) by a transformation involving the χ ERIs, This is formally a four-index transformation which scales as O(M 5 ). 55
Throughout this work, we denote the error which arises from the application of Eq. (11) instead of Eq. (10) by E. The superscript to E indicates which quantity is affected; additional specifications are then given as subscripts. For example, the error introduced by the NDDO approximation for Obviously, M 4 different errors E φ ERI µνλσ need to be accounted for. First, we consider the effect of erroneous φ ERIs on the HF energy when the self-consistent solution φ P obtained from an exact HF calculation is applied. The errors in the φ ERIs affect the Coulomb matrix elements, and the exchange matrix elements, Interestingly, the matrix element φ J µν will always be exactly zero if χ I µ and χ J ν are centered on different atoms (I = J) irrespective of the number of atoms on which χ K λ and χ L σ are centered. By contrast, φ K µν will not be strictly zero in this case (see also Figure 1). By affecting the Coulomb and exchange matrices, the NDDO approximation will introduce an error compared to E HF el which we denote by E HF = E E HF el , in the closed-shell case.
If we apply the NDDO approximation, we, however, must iterate to self consistency. The self-consistently obtained χ P NDDO will likely not be the same as φ P . Hence, another error arises from the NDDO approximation by introducing errors in other quantities during the self-consistent-field (SCF) cycles, i.e., in the coefficient matrix and in the matrix of orbital energies. We denote this error by G and again indicate by a superscript which quantity is affected by the error (and give additional specifications as subscripts). By contrast, E denotes the error which is obtained when applying χ P NDDO . The difference of the two self-consistent solutions produces G HF = G E HF el , The electronic HF energy does, by definition, not contain effects from electron correlation. 55 Various electronic structure methods are available for calculating correlation energies. 55 The prevalent single-and multi-reference method require the calculation of ERIs in the molecular orbital basis ψ ( ψ ERIs). These ψ ERIs are obtained through a 4-index transformation of the φ ERIs (or the χ ERIs), This 4-index transformation is similar to the 4-index transformation with which the φ ERIs are determined from the χ ERIs (Eq. (10)). When applying Eq. (11), we may approximate the ψ ERIs as, The formal scaling of the ψ ERI evaluation step is therefore reduced from If we determine the coefficient matrix in a self-consistent field procedure, we will introduce an additional error G ψ ERI ijkl from applying a different coefficient matrix, In this work, we demonstrate how E ψ ERI ijkl and G ψ ERI ijkl affect the MP2 correlation energies by quantifying E ψ MP2 and G ψ MP2 , respectively.

Extension to Conventional Basis Sets
The χ-basis must fulfill the condition that it is locally orthogonal for the NDDO approximation to be justifiable. [44][45][46][47][48] This condition is, however, generally not fulfilled by ordinary basis sets. We denote ordinary Gaussian-type basis functions which are, in general, not locally orthogonal by τ I µ (µ-th basis function of type τ centered on atom I). The NDDO approximation is not straightforwardly applicable for an arbitrary τ -basis, because we observe E φ ERI µνλσ > 1.0 atomic units (a.u.) for a large number of φ ERIs at the example of the water molecule (see Figure 2). To cure this problem, we propose to transform {τ I µ } to a locally orthogonal basis {χ I µ } by the application of the transformation matrix T , so that Figure 2 illustrates that E φ ERI µνλσ 1.0 a.u. after local orthogonalization. We demonstrate in Section 4 that this statement holds true in general.
Red circles encode ERIs for which I = J and K = L in the χ-basis and in the τ -basis, and blue crosses encode ERIs for which I = J or K = L in the χ-basis and in the τ -basis.

Molecules in a χ-Basis
For the first part of our analysis of the NDDO approximations, we applied the ECP-3G basis set. 5, 63 We selected the ECP-3G basis set because its basis functions form a χ-basis and because it is applied (in slightly modified forms) in the OM1, OM2, and OM3 models. [5][6][7] The ECP-3G basis set specifies one

Effect on Electron-Electron Repulsion Integrals
For the simplest possible model system, H 2 , 2 4 = 16 ERIs arise for a given  is as large as φ 1 φ 2 |φ 1 φ 2 . Figure 3 shows that this is a poor approximation for R 12 for 1.6 Å < R 12 < 3.0 Å. For all R 12 , we find a significant error in at least one of the φ ERIs.
We also encounter nonnegligible E φ ERI µνλσ for all structures in the sets of molecules B, C, and D. The largest absolute error of an φ ERI in a given molecule is between 0.10 a.u. and 0.23 a.u. It is already obvious at this point that only an efficient error cancellation may yield useful observables based on these erroneous φ ERIs (or on density matrices obtained with them in SCF procedures). Figure 4 shows that the cumulative absolute error (CAE) in the φ ERIs, grows roughly linearly with the number of basis functions for reaction B ( Figure 4). We observe a linear growth not only in the overall cumulative absolute error, but also in individual contributions to it when we break down the corresponding χ ERIs in one-, two-, three-, and four-center ERIs. In agreement with previous studies, 40, 44-48 we find that assuming φ µ φ ν |φ λ φ σ to be as large as χ I µ χ I ν |χ J λ χ K σ for I = J = K and φ µ φ ν |φ λ φ σ to be as large as χ I µ χ I ν |χ I λ χ J σ for I = J introduces significant errors (purple and orange circles in Figure 4, respectively). These contributions account for 60-65% of the overall cumulative absolute error for B. Figure 4 shows that E φ ERI CAE also grows approximately linearly with the number of basis functions for C and that the spread of the individual E φ ERI CAE is large. The change of E φ ERI CAE with reaction progress for D in Figure 4 shows that the cumulative absolute error in the φ ERIs crucially depends on the arrangement of the atomic nuclei, i.e. on the underlying nuclear framework that generates the external potential. The cumulative absolute error E φ ERI CAE (black squares) is broken down into the corresponding χ ERIs for the molecule set B: χ I µ χ I ν |χ I λ χ I σ (red circles),

Error Propagation: the Hartree-Fock Energy
For H 2 , E HF depends linearly on the cumulative error of the φ ERIs because φ P can be determined analytically, 54  Figure 3). For larger or smaller R 12 , we encounter larger E HF because the errors in the φ ERIs do not cancel as efficiently.
For the calculation of E HF , the φ ERIs are contracted with elements of the density matrix. The elements of the density matrix are therefore a central ingredient for an efficient error cancellation. We see how differently the errors in the φ ERIs add up by comparing E HF and G HF for the molecule set B in Figure 5. Interestingly, |E HF | and G HF increase in an almost perfectly linear fashion with the number of basis functions. For a given number of basis functions, G HF is significantly larger than |E HF |. Not surprisingly, we observe large E HF and G HF for the vast majority of molecules in B and C. In line with previous studies, 33,40,[44][45][46][47][48][49] we see that the NDDO approximation is a rather crude approximation. We discuss strategies to overcome this situation in Section 5.

Error Propagation: the MP2 Energy
The error in the φ ERIs also propagates to the ψ ERIs (Eq. (18)). We show the effect of the NDDO approximation on selected ψ ERIs, ψ 1 ψ 1 |ψ 1 ψ 1 , ψ 1 ψ 1 |ψ 2 ψ 2 , ψ 2 ψ 2 |ψ 2 ψ 2 , and ψ 1 ψ 2 |ψ 1 ψ 2 , in H 2 in Figure 6. We choose to study these ψ ERIs because they are applied in the calculation of the MP2 correlation energy ( ψ 1 ψ 1 |ψ 2 ψ 2 ) and of the full-configuration interaction (FCI) correlation energy 54 (all of these ψ ERIs). All of these ψ ERIs are significantly affected by the NDDO approximation. Hence, it comes as no surprise that the MP2 and FCI correlation energies are deteriorated by the NDDO approximation. The general shape of E MP2 follows that of E ψ ERI 1122 and has a minimum at R 12 = 0.86 Å where E MP2 = −0.009 a.u. As a consequence, instead of E MP2 = −0.016 a.u. a E MP2 = −0.007 a.u. is obtained when the NDDO approximation is invoked. The MP2 correlation energy is therefore significantly underestimated. Also the FCI correlation energy is also significantly underestimated for R 12 < 0.72 Å. An unfortunate addition of the errors in ψ 2 ψ 2 |ψ 2 ψ 2 and ψ 1 ψ 1 |ψ 1 ψ 1 , however, leads to an overestimation of E FCI for R 12 ≥ 0.72 Å. The errors in the ψ ERIs and, hence, in E FCI , only vanish in the limit of very large R 12 (R 12 > 25.0 Å).

Molecules in the τ -Basis
When applying ordinary τ -basis sets, large errors in the φ ERIs arise. For H 2 and R 12 = 0.74 Å described in a def2-QZVP basis set, 64 for example, the largest absolute error in an φ ERI amounts to 1.12 a.u. The application of these erroneous φ ERIs even leads to a positive E HF = 3.48 a.u. By contrast, the largest absolute error in an φ ERI is reduced to 0.22 a.u. when we locally orthogonalize the basis (Eq. (23)). After the local orthogonalization of the basis set, we obtain G HF = 0.01 a.u. A prior local orthogonalization led to significantly smaller largest absolute errors in the φ ERIs and cumulative errors in the φ ERIs for all molecules (B, C, and D) by up to two orders of magnitude.
We found a fundamental limitation of the NDDO approximation in the course of our analysis of different τ -basis sets (see Table 1): Electronic energies calculated with the NDDO approximation do not converge with respect to the basis set size. The HF limit for H 2 for R 12 = 0.74 Å was determined to be E HF el = −1.133629 a.u. 70 No convergence for E HF el of H 2 with increasing basis-set size is observed when invoking the NDDO approximation (see Table 1). We also found this behavior for 50 randomly selected molecules out of molecule set C.

Improving on the NDDO Approximation
The NDDO approximation introduces significant errors in the φ ERIs, observed here for the simplest possible neutral molecule, H 2 , and for a diverse selection of medium-sized organic compounds. Obviously, these errors are too large and too unsystematic for the NDDO approximation to be useful in purely ab initio electronic structure models.

Modifications
In ab initio electronic structure calculations, it is customary to apply screening techniques to determine which χ ERIs are negligibly small. 55 Errors are controlled by thresholds with respect to which the χ ERIs are neglected. 55 By contrast, the NDDO approximation cannot be applied as a screening tool for χ ERIs because it can only predict whether the corresponding ERIs in the φ-basis are negligibly small or not. Each φ ERI, however, encodes information on all χ ERIs (Eq. (10). It is, however, computationally difficult to improve on the NDDO approximation by an explicit transformation of the ERIs from the χ-to the φ-basis because this is a 4-index transformation that scales as O(M 5 ). Some attempts 49,71 were made to improve on the approximation of the φ ERIs for which the corresponding χ ERIs are one-center ERIs ( χ I µ χ I ν |χ I λ χ I σ ) and two-center ERIs ( χ I µ χ I ν |χ J λ χ J σ for I = J). However, these suggestions did not find widespread use. Moreover, our results suggest that a correction of these φ ERIs does not suffice to obtain a reliable HF energy E HF el . A large portion of the error originates from the φ ERIs which are assumed to be zero, but are not exactly zero. If we had to estimate these φ ERIs, we would again need close to O(M 4 ) operations. These estimates would have to be accurate due to the plethora of small φ ERIs which would again compromise the computational efficiency of the NDDO approximation.

Capitalizing on Error Cancellation
In SEMO models, further approximations (beyond NDDO) allow for error compensation to a certain degree. All successful NDDO-SEMO models introduce various parametric expressions to assemble the one-electron matrix and to calculate the nucleus-nucleus repulsion energy. [2][3][4][5][6][7]50 For contemporary NDDO-SEMO models, the parametric expressions were designed such that the result of the SCF optimization yields a result close to a reference energy. Overall, the results obtained with NDDO-SEMO models achieve a remarkably high accuracy with respect to these reference data. 51,52 At the same time, NDDO-SEMO models are notoriously unreliable for molecules not considered in the parametrization procedure. 51,52 It is virtually impossible to rationalize why errors occur due to the number and the diversity of the approximations invoked in an NDDO-SEMO model. At least some of these errors are likely to be due to the NDDO approximation. In general, we may anticipate that the improvement of the parametric expressions in NDDO-SEMO models is as complicated as the direct correction for the error introduced by the NDDO approximation as discussed in Section 5.1.
Therefore, the only viable use of NDDO-SEMO models appears to be their combination with system-focused rigorous error estimation schemes as proposed in Refs. 72, 73.

Correcting the Two-Electron Matrices
In 1969, Roby and Sinanoǧlu made an attempt to accelerate single-point HF calculations for a diverse set of small molecules. 32 They suggested to scale χ G NDDO with a matrix Γ to obtain a better estimate for φ G, Their attempt to define general rules to assemble Γ turned out to be impossible. 32 In this Section, we reconsider and build upon the Roby-Sinanoǧlu approach. We can exactly determine Γ({R n I }) for a given structure {R n I } from a reference HF, KS-DFT, or multi-configurational SCF calculation (yielding Eq. (29) defines a system-focused NDDO model that can be applied in connection with any Fock operator. It allows us to maintain complete error control on the resulting model because we can determine Γ {R The original Roby-Sinanoǧlu approach is not the only one conceivable for the construction of correction matrices. In fact, a multiplicative correction matrix changes matrix elements through a combination of elements of the original matrix. An additive correction appears easier and more straightforward to achieve the goal of re-adjusting individual matrix elements. We may therefore define separate additive correction with the matrices Γ J and Γ K to χ J NDDO and to χ K NDDO , or a unification of Γ J and Γ K as a total additive correction. For reaction A, this additive approach leads to smaller errors in general (see Figure 8) for the first part of the reaction.
Both correction approaches suffer from limitations. The nuclear coordinates and also the density matrices obviously differ in a sequence of structures. Eq. (29) will only yield sensible results when the change of both quantities remains small. We are currently exploring the possibility to apply this strategy in practice for sequences of structures as they occur during structure optimization, in Born-Oppenheimer molecular-dynamics trajectories, or in interactive reactivity studies.

Conclusions
The NDDO approximation is a central ingredient for many modern semiempirical molecular orbital models. We studied the effect of the NDDO approximation on the ERIs in the symmetrically orthogonalized basis for the simplest possible model system, H 2 , and for a diverse set of molecules. As expected and in agreement with previous results, [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] we find that the NDDO approximation leads to significant errors for molecules in their equilibrium structure. The errors do only vanish in the atomization limit where the overlap between different basis functions vanishes. The errors in the ERIs in the symmetrically orthogonalized basis increase roughly linearly with the number of basis functions. Additionally, we found that the errors in the ERIs in the symmetrically orthogonalized basis depend strongly on the arrangement of the atomic nuclei. These nonnegligible errors introduced by NDDO may only be alleviated by an efficient error cancellation which is in operation in SEMO models, but not in HF theory. For HF calculations, error cancellation is unlikely to occur because of the fact that the ERIs in the symmetrically orthogonalized basis are contracted with elements of the density matrix.
We were then able to dissect how the NDDO approximation affects ERIs in the molecular orbital basis and, hence, correlation energies. We found that MP2 correlation energies are underestimated and the underestimation increases with the number of basis functions. This finding explains previous reports 61, 68 that correlation energies obtained with respect to an NDDO-SEMO reference wave function are far too small.
We proposed a local orthogonalization that allowed us to transgress the domain of minimal basis sets and to apply ordinary basis sets in conjunction with the NDDO approximation. While we observed a drastic reduction in the largest absolute errors in the ERIs in the symmetrically orthogonalized basis, we discovered another limitation of the NDDO approximation. Electronic energies calculated with the NDDO approximation do not converge with respect to the basis set size so that this solution to the small basis-set restriction of NDDO-SEMO models does not pay off.
We then continued to propose how one could still capitalize on the efficiency enabled by the NDDO approximation without significant loss of accuracy in a system-focused manner for similar structures. Specifically, we proposed a strategy to correct for the errors introduced by the NDDO approximation in the two-electron matrices which was inspired by work of Roby and Sinanoǧlu. 32 The two-electron matrix obtained within the NDDO approximation is modified with a correction matrix obtained from a reference HF, KS-DFT, or multi-configurational SCF calculation. These correction matrices are transferable to a certain degree within sequences of related structures.

Appendix: Computational Methodology
All calculations in this work were carried out with our modified version of PySCF (version 1.4). 76,77 The ERIs in the τ -basis were transformed to the corresponding ERIs in the χ-basis or in the φ-basis with the ao2mo integral transformation module of PySCF. We applied the ECP-3G, 5, 63 STO-3G, 78 3-21G, 79 6-311G, 80 and def2 64 basis sets in calculations.
Raymond and co-workers assembled a database considered to be representative of chemical space. 65 We randomly chose 5000 entries of the QM9 data set 65,66 (set C) to study the error of the NDDO approximation across a large set of molecules. Additionally, we worked with linearly growing alkane chains (set B) with the stoichiometry C x H 2x+2 , (x = 1, 2, ..., 15). We include the optimized structures as Supplementary material. Finally, we selected three reactions (set D) which we had considered 29 for interactive reactivity explorations in the framework of real-time quantum chemistry; 25-30, 81, 82 they can be found in the supplementary material in Ref. 29.