No need for external orthogonality in subsystem density-functional theory

Recent reports on the necessity of using externally orthogonal orbitals in subsystem density-functional theory (SDFT) [Annu. Rep. Comput. Chem., 8, 2012, 53; J. Phys. Chem. A, 118, 2014, 9182] are re-investigated. We show that in the basis-set limit, supermolecular Kohn–Sham-DFT (KS-DFT) densities can exactly be represented as a sum of subsystem densities, even if the subsystem orbitals are not externally orthogonal. This is illustrated using both an analytical example and in basis-set free numerical calculations for an atomic test case. We further show that even with finite basis sets, SDFT calculations using accurate reconstructed potentials can closely approach the supermolecular KS-DFT density, and that the deviations between SDFT and KS-DFT decrease as the basis-set limit is approached. Our results demonstrate that formally, there is no need to enforce external orthogonality in SDFT, even though this might be a useful strategy when developing projection-based DFT embedding schemes.


Introduction
Quantum-chemical subsystem and embedding methods (for recent reviews, see, e.g., ref. [1][2][3][4] allow for an efficient treatment of complex chemical systems. For approaches based on densityfunctional theory (DFT), one can distinguish two different strategies that allow for a formally exact subsystem treatment. Such formally exact approaches are particularly relevant because they can serve as a starting point for multiscale methods combining wavefunction-based quantum-chemical methods with DFT (wavefunction-in-DFT embedding methods). [5][6][7][8][9][10] Firstly, subsystem density-functional theory (SDFT) 11,12 and frozen-density embedding (FDE) 13 are based on a partitioning of the total electron density r(r) into (possibly overlapping) subsystem electron densities. The embedding potential obtained by reformulating Kohn-Sham-DFT (KS-DFT) in terms of these subsystem electron densities contains a contribution of the nonadditive kinetic energy, which can either be approximated 13,14 or treated exactly by reconstructing the embedding potential leading to a certain target density. [15][16][17][18] It is widely believed that with such an exact treatment of the nonadditive kinetic energy, SDFT and FDE will exactly reproduce the total electron density from supermolecular KS-DFT calculations. 2,4,19 Secondly, the partitioning can also be achieved by assigning the KS-like orbitals to different subsystems. To this end, nonlocal projection techniques are used to ensure orthogonality of the subsystem orbitals, so that the non-additive kinetic energy vanishes. 20 Hence, no approximate functionals for this energy component need to be introduced and the use of potential reconstruction techniques is not required. It can be demonstrated that such projection-based embedding calculations are exactly equivalent to supermolecular KS-DFT calculations. 20 The question of orthogonality between the different subsystem orbitals (''external orthogonality'') is recently causing some apparent confusion in the literature. Hoffmann and co-workers have questioned that SDFT is formally exact. 21,22 Instead, they argue that an exact treatment of the nonadditive kinetic energy alone does not suffice to make SDFT equivalent to KS-DFT, but that it is mandatory to enforce the external orthogonality of the subsystem orbitals. Also Chulhai and Jensen have picked up the argument, based on the work by Hoffmann and co-workers, that non-additive kinetic-energy potentials as used in SDFT and FDE may lead to incorrect results even for the exact kinetic-energy functional. 23 If these arguments were correct, the two apparent conclusions are: (i) if external orthogonality is not enforced, even the exact non-additive kinetic-energy potential will give incorrect results, and (ii) if external orthogonality is enforced, this potential is exactly zero and SDFT becomes equivalent to projection-based embedding schemes. Hence, there would be no reason whatsoever to attempt constructing approximate non-additive kinetic energy potentials or to develop numerical schemes for the reconstruction of accurate embedding potentials.
In this work, we re-investigate the question whether or not external orthogonality is mandatory for an exact SDFT treatment.
First, in Section 2 we reconsider the analysis given in ref. 21 and 22 to formally address this question. In a second step, we study this question on the basis an analytical example in Section 3 and on the basis of numerical results obtained with reconstructed potentials, which are derived from exact relations for the nonadditive kinetic-energy potential, in Section 4. We conclude from our results in Section 5.

Theory
In (supermolecular) KS-DFT, the ground-state electron density is obtained as, where n tot is the total number of electrons and where the occupied supermolecular KS-orbitals {f KS i } are obtained from the self-consistent solution of the KS equations, with the KS potential v KS eff [r](r) containing the usual contributions of the nuclear potential, the electronic Coulomb potential, and the exchange-correlation potential. Note that the solution of the KS equations yields a complete set of (occupied and virtual) ,N , which spans the one-electron Hilbert space.
In the following, we will consider the special case of SDFT for two subsystems (A and B), which already introduces all the basic questions that may arise in the more general case of arbitrarily many subsystems. The goal of SDFT is then to find electron densities for these two subsystems r A and r B which add up to the correct total electron density r, r tot (r) = r A (r) + r B (r). ( Each of the subsystem densities is expressed through its own set of KS-like orbitals, where n A and n B are the numbers of electrons in subsystems A and B, respectively. The orbitals within one set are internally orthogonal, but generally, orbitals of different subsystems in SDFT calculations are non-orthogonal, i.e., the subsystem orbitals are not externally orthogonal. In FDE, a frozen density r B (r) is chosen for subsystem B and the orbitals of subsystem A are then obtained from the solution of the KS-like equations, FDE is generally understood to be exact -in the sense that the sum of the subsystem densities, r A (r) + r B (r) is equal to the total electron density r tot (r) from a supermolecular KS-DFT calculation [cf. eqn (3)] -if the following conditions are fulfilled: (1) The contribution of the nonadditive kinetic energy T nadd s [r A ,r B ] to the embedding potential is treated exactly, for instance by using analytical potential reconstruction techniques. 24,25 Numerical schemes for the reconstruction of the embedding potential 15-18 might introduce errors due to the inaccuracies of optimized effective potential algorithms, in particular if finite basis sets are used. 26 (2) The exchange-correlation energy is treated consistently in all calculations by using the same approximate functional in the supermolecular calculation, the subsystem calculations, and when evaluating the embedding potential.
(3) The KS(-like) equations for the supermolecule [eqn (2)] and for subsystem A [eqn (5)] are solved exactly, i.e., without introducing a finite basis set for the orbitals. To which extent this condition can be relaxed will be discussed in the following.
(4) All appearing densities, most importantly r A (r) = r tot À r B (r), are noninteracting v s -representable. In particular, the frozen density r B (r) must be chosen such that r B (r) r r tot (r) holds at every point in space. This condition can be relaxed in SDFT, where the densities of both subsystems are updated iteratively in freeze-and-thaw cycles. 27 That is, they assume that it must be possible to expand the supermolecular KS orbitals in the basis of the subsystem orbitals. However, all that is required for SDFT to be exact is that the supermolecular density can be obtained, but not necessarily the supermolecular KS orbitals. Thus, the question is rather whether it is possible that the two different expansions of This journal is © the Owner Societies 2016 Phys. Chem. Chem. Phys. eqn (1) and of eqn (4) result in exactly the same total electron density, i.e., for the combined set of subsystem orbitals. Following the analysis in ref. 21 and 22, we proceed by creating an explicitly orthonormalized combined set of subsystem orbitals for subsystems A and B, {f orth where S ij = hf AB i |f AB j i is the overlap matrix of the combined set of subsystem orbitals. Note that so far there is no reason to assume that these orthonormalized subsystem orbitals {f orth i } span the same subspace as the supermolecular KS orbitals {f KS i }. The sum of the subsystem densities [right-hand side of eqn (7)] can now be expressed as, 21 We now express the orthonormalized subsystem orbitals {f AB i } i=1,n tot in terms of the occupied and virtual supermolecular KS orbitals {f KS p } p=1,N , which form a complete basis of the oneelectron Hilbert space, Note that at this point there is no reason to assume that only the occupied supermolecular KS orbitals are sufficient for expanding the subsystem orthonormalized subsystem orbitals. Expressing eqn (9) in the basis of the supermolecular KS orbitals, we arrive at U piSij U qj . For SDFT to be exact, this has to be equal to the left-hand side of eqn (7), that is, At this point, one might be tempted to conclude that this can only be the case if and only if the occupied-occupied block of S is the identity matrix (i.e., S ij = d ij for i, j = 1,. . .,n tot ) and all other blocks of S are zero. This is only possible if S = I, i.e., if and only if the subsystem orbitals are externally orthogonal. This step is at the heart of the arguments of Hoffmann and coworkers in ref. 21 and 22. However, this conclusion can only be drawn if the set of orbital products {f KS p f KS q } p,q=1,N is linearly independent. It is well known that for a complete basis set of the one-electron Hilbert space, the corresponding product basis set is linearly dependent. For a proof, we refer to the clear presentation by Görling et al. in Appendix A of ref. 28.
Specifically, since the supermolecular KS orbitals {f KS p } p=1,N form a complete basis, the set of orbital products {f KS p f KS q } p,q=1,N must be linearly dependent. Therefore, there is a non-trivial solution of the equation which in turn makes it possible that eqn (7) and (12) can be solved for S a I, i.e., for the case of subsystem orbitals that are not externally orthogonal. Thus, we conclude that the arguments for enforcing external orthogonality in SDFT given in ref. 21 and 22 do not hold in the basis-set limit. Even for incomplete basis sets, (near) linear dependencies in the orbital products may still occur. 28

An analytical example
It can be demonstrated that SDFT is able to exactly reproduce the supermolecular density with subsystem orbitals that are not externally orthogonal for a simple analytical model, which has been previously introduced by Savin and Wesolowski. 24 To this end, we consider a closed-shell four-electron system, in which the supermolecular electron density is given by, that is, there are two doubly occupied orbitals. As was already shown in ref. 24, the embedding potential can be calculated analytically in this case if the supermolecular KS potential is known. Moreover, it has been demonstrated that the supermolecular KS density can be exactly represented as a sum of subsystem densities obtained with this analytical embedding potential.
As an extreme case, we consider a decomposition of the total densities into equal subsystem densities, Each subsystem density then has to be composed from a single, doubly-occupied orbital, These two (identical) orbitals for subsystems A and B are not externally orthogonal, but have maximal overlap hf A 1 |f B 1 i = 1. At the same time, the sum of the subsystem densities r A (r) + r B (r) is by construction exactly equal to the supermolecular density r tot (r). Note that both the supermolecule and the two subsystems have been treated as closed-shell systems.
The potential in the KS-like equations for the subsystems A and B that leads to the required subsystem densities is analytically given by, 2 Dr A=B ðrÞ r A=B ðrÞ À 1 8 It is further obvious that the subsystem orbitals are not contained in the subspace spanned by only the occupied supermolecular orbitals. Instead, they can be expanded in the occupied and virtual supermolecular KS orbitals as with the expansion coefficients which are in general nonzero for all p. Note, however, that it is nevertheless possible to construct a finite basis set that will result in the exact orbitals and densities in both the supermolecular and in the subsystem calculations, for instance by using the basis set {f KS 1 ,f KS 2 ,f A/B 1 } or by orthogonalizing this set.

Basis-set limit
To investigate an example with more than two orbitals, we turn to the argon atom as a test case, because for atoms the KS equations can easily be solved numerically on a grid, which provides results close to the basis-set limit.
As total density r KS tot (r), we choose the one obtained in a KS-DFT calculation performed with the ADF program package 29,30 using the BP86 exchange-correlation functional 31,32 and a QZ4P Slater-type orbital basis set. 33 We partition this total density into two closed-shell subsystem densities by defining where {f KS i } are the occupied supersystem KS orbitals [see Fig. 1c and d for plots of the radial functions rf KS i (r)], and r B (r) = r KS tot (r) À r A (r).
By construction, these subsystem densities are both nonnegative at any point in space and each integrates to an integer number of electrons (n A = 10 and n B = 8). The radial densities r 2 r KS tot (r), r 2 r A (r), and r 2 r B (r) are shown in Fig. 1a.
To verify that these two subsystem densities can each be represented in terms of subsystem orbitals, we performed a reconstruction of the KS potentials v s [r A ] and v s [r B ] that yield the target densities r A (r) and r B (r), respectively, in a fully numerical solution of the KS equations on a logarithmic radial grid of 400 points. The modified van Leeuwen-Baerends algorithm 34 used for the potential reconstruction as well as the computational details of the numerical solution of the KS equations are described in ref. 26 and 35. It turns out to be possible to reconstruct the subsystem densities almost exactly. A comparison of the reference and the reconstructed subsystem densities as well as their difference on a logarithmic scale are shown in Fig. 1b is smaller than 10 À6 electrons in both cases. For both subsystems, the lowest-energy orbitals are occupied, i.e., the reconstructed densities correspond to the ground-state solution. The subsystem orbitals for subsystems A and B are shown in Fig. 1c and d, respectively, alongside the corresponding supersystem KS orbitals. For subsystem A, the density as defined by eqn (21) is mainly composed of the five lowest supersystem KS orbitals f KS 1s , and f KS 2s , f KS 2p (three degenerate orbitals), with a small admixture of the remaining orbitals. Consequently, the subsystem orbitals are very similar to the supersystem ones, with the deviations being caused by the admixture of small 3s-and 3p-orbital contributions in the subsystem density. On the other hand, the density of subsystem B is composed of the four orbitals f KS 3s and f KS 3p (three degenerate orbitals), with a small admixture of the remaining orbitals. In the subsystem calculation, this density needs to be reproduced by the ground-state orbitals f B 1s and f B 2p . Hence, in contrast to the corresponding supersystem KS orbitals, the radial parts of the subsystem orbitals are nodeless and large differences occur in the region from r = 0.0 bohr to r = 0.7 bohr. The fact that for both subsystem A and subsystem B the lowest s-and p-orbitals are each nodeless leads to a significant overlap between the subsystem orbitals, as can be seen from the overlap matrix S of the subsystem orbitals shown in Table 1. Thus, the subsystem orbitals are not externally orthogonal. Nevertheless, the sum of the subsystem densities reproduces the total density exactly. Of course, the s-and p-orbitals of the different subsystems are still mutually orthogonal because of the atomic symmetry of the specific test case considered here.
In summary, our results show that if the KS-like equations for the subsystems are solved numerically on a grid to avoid errors introduced by a finite basis set, a total density can be reproduced exactly as a sum of subsystem densities, even if the subsystem orbitals are not externally orthogonal. This confirms that SDFT is indeed exact if all the conditions listed in Section 2 are fulfilled, without the need to enforce external orthogonality.

Finite basis sets
Finally, we also want to investigate the case in which both the supermolecular KS-DFT calculation and the subsystem calculations are performed using a finite basis set for expanding the supermolecular and the subsystem orbitals. This is the scenario in which SDFT is usually applied in practical calculations. With finite basis sets, the linear dependence of the orbital products is not guaranteed. With small basis sets, the orbital product can in general be expected to be linearly independent. 28 As discussed in Section 2, in this case the arguments of ref. 21 and 22 hold and it will not be possible to express the supermolecular density as a sum of two subsystem densities unless the subsystem Table 1 Overlap matrix S of the subsystem orbitals obtained for the Ar atom test case View Article Online orbitals are externally orthogonal. However, as the basis set gets larger, (near) linear dependencies of the orbital products will appear. We consider one of the test cases investigated in ref. 22, an NH 3 Á Á ÁNH 3 complex (see Fig. 2), which can be partitioned into two NH 3 subsystems. All calculations have been performed using the ADF program package 29,30 and its FDE implementation, 36 in combination with the PyADF scripting framework. 37 In all of the subsystem calculations a supermolecular integration grid was used to ensure the transferability of the results. All subsystem calculations use the full supermolecular basis set [FDE(s)], i.e., ghost basis functions on the atoms of the frozen subsystem are included. The PW91 exchange-correlation functional 38,39 has been used in all calculations.
First, a supermolecular KS-DFT calculation using a finite Slater-type orbital basis set has been performed to obtain a supermolecular reference density r KS tot (r). Subsequently, SDFT using accurate reconstructed embedding potentials were performed to investigate to which extent it is possible to reproduce the supermolecular density as a sum of two subsystem densities when using the same finite basis set in the subsystem calculations.
These accurate SDFT calculations have been done as follows: in the initial iteration, the densities r (0) A (r) and r (0) B (r) calculated for the isolated subsystems A and B, respectively, are used as initial guess for the subsystem densities. In the n-th iteration, the kinetic-energy component of the embedding potentials for subsystem A and B are then reconstructed as 16,41 dT nadd respectively. Here, v s [r target ] is the local potential that yields the density r target (r) as its ground state density. This local potential is reconstructed using the algorithm of Wu and Yang 42 as described in ref. 16. Since in the context of this study, we are only interested in the resulting densities and not in obtaining high-quality potentials, no further regularization 26,43 is applied to the reconstructed potentials. With these accurate kineticenergy components of the embedding potential, the KS-like equations for subsystems A and B [eqn (5)] are solved to obtain updated subsystem densities r (n) A (r) and r (n) B (r), respectively. This procedure is repeated until the subsystem densities are converged. Typically, 10-20 of these iterations have been performed.
Such a freeze-and-thaw procedure ensures that both subsystem densities are noninteracting v s -representable and thus alleviates the shortcomings of the approach used in ref. 16 (see also ref. 44 and 45), where a frozen density constructed from localized orbitals was used. Since both subsystem densities are optimized, it should result in a pair of subsystem densities r A (r) and r B (r) that minimizes the total energy functional E[r A + r B ].
The accuracy of the sum of the resulting subsystem densities r A (r) + r B (r) compared to the supermolecular density r KS tot (r) can be measured as, The resulting density errors are summarized in Table 2. In addition, the difference densities Dr(r) = r KS tot (r) À r A (r) À r B (r) in a cut plane containing the N-HÁ Á ÁN atoms are shown in Fig. 3.
We find that in the SDFT calculations with reconstructed potentials (SDFT/RecPot, Fig. 3c), the sum of the subsystem densities is not identical to the supermolecular density, but approaches it closely. This is the case already with the small DZP Slater-type orbital basis set 33 and when using the ATZP Slater-type orbital basis set. 33,40 In fact, the integrated density error of 0.0046 electrons and 0.0047 electrons with DZP and ATZP, respectively, is one order of magnitude smaller that in an SDFT calculation in which the the nonadditive kinetic energy was approximated using the PW91k 46 kinetic-energy functional (SDFT/PW91k, Fig. 3b) and almost two orders of magnitude smaller than for the sum of the densities from KS-DFT calculations for the isolated subsystems (SumFrag, Fig. 3a). When repeating the calculations with a larger ET-pVQZ Slater-type orbital basis set, 47 the discrepancy is further reduced (Fig. 3d). This suggest that when approaching the basis-set limit, it becomes possible to reproduce the supermolecular density as a sum of subsystem densities.   (26) for subsystem calculations on a NH 3 Á Á ÁNH 3 complex compared to a supermolecular KS-DFT calculation. ''SumFrag'' refers to the sum of the isolated subsystem calculations, ''SDFT/PW91k'' to a freeze-and-thaw SDFT calculation using the approximate PW91k kinetic-energy functional, and ''SDFT/ PotRec'' to a freeze-and-thaw SDFT calculation using accurate reconstructed embedding potentials, as described in the text. Results are shown for both the small DZP basis set, the medium-sized ATZP basis set and the large ET-pVQZ basis sets

View Article Online
In order to evaluate the external orthogonality of the subsystem orbitals, the overlap matrix S of the subsystem orbitals from the SDFT/RecPot calculation with the ET-pVQZ basis set is shown in Table 3. Indeed, the overlap of the subsystem orbitals is non-zero, i.e., the subsystem orbitals that can closely reproduce the supermolecular density are not externally orthogonal.
In summary, with finite basis sets, it is in general not possible to exactly reproduce the supermolecular KS-DFT density as a sum of two subsystem densities. However, as the size of the basis set is increased and the basis-set limit is approached, the discrepancies decrease and with basis sets of reasonable size it is possible to reproduce the supermolecular density rather accurately. While most of the remaining difference in the densities is most likely due to the use of a finite basis set, a part could also be caused by inaccuracies of the numerical potential reconstruction algorithm, such as the use of a finite basis set for the potential in the Wu-Yang method.
It should also be noted that with finite basis sets, both the supermolecular KS-DFT and the SDFT calculation are affected by a basis set error. Even though both calculations minimize the total energy functional E[r tot ], the space of the possible densities is different since different expansions of the density are used. We would argue that it is not a priori possible to decide which expansion will result in the more accurate total density. If we assume that the supermolecular KS-DFT/ET-pVQZ calculation gives results close to the basis-set limit, both the supermolecular KS-DFT and the SDFT/RecPot calculation in the smaller ATZP basis set exhibit a similar basis set error: the integrated density error D abs err when comparing the supermolecular KS-DFT/ATZP calculation to the one in the larger ET-pVQZ basis set is 0.1315 electrons, whereas when comparing the sum of the SDFT/RecPot subsystem densities in the ATZP basis set to the supermolecular KS-DFT/ET-pVQZ density the integrated difference density D abs err is 0.1322 electrons. However, a more systematic investigation of the basis set errors will be necessary to decide whether in a finite basis set, supermolecular KS-DFT or SDFT calculations yield results closer to the basis-set limit.

Conclusion
Our results demonstrate that recent claims 21-23 that SDFT can only be exact if the external orthogonality of the subsystem orbitals is enforced are not justified. In fact, our analysis reveals that such claims are based on the assumption that the products of linearly independent orbitals are themselves linearly independent. This assumption is not true in the basis-set limit, and might also not hold for many finite basis sets used in practical calculations because (near) linear dependencies of orbital products appear.
We have illustrated our analysis by demonstrating for the analytically solvable case of a two-orbital system that the sum of subsystem densities arising from mutually non-orthogonal subsystem orbitals can exactly reproduce a supermolecular KS-DFT density. The same has been shown, within numerical accuracy, for a many-orbital atom when solving the KS-like equations for the subsystems in the basis-set limit. For calculations in a finite basis set, we could show that a supermolecular density can be approached closely by the sum of subsystem densities, with the differences decreasing as the basis-set limit is approached.
However, for finite basis sets in which the products of basis functions are indeed linearly independent, the arguments given in ref. 21 and 22 do hold and external orthogonality is needed to reproduce the supermolecular density exactly. Thus, the minimization of the total energy functional E[r tot ] is performed with respect to different spaces of possible densities in SDFT and in KS-DFT. However, it is not clear which one will result in the density that yields a lower total energy (and is thus closer to the basis-set limit) if the nonadditive kinetic energy is evaluated exactly. This question is complicated by the fact that in SDFT, the subsystem densities of externally orthogonal subsystem orbitals cannot be expected to be v s -representable. 16,24 Therefore, the supermolecular KS-DFT ground state density will not be accessible in SDFT, whereas the SDFT ground-state density might not be accessible in KS-DFT.
Alternatively, one can, of course, use subsystem approaches that are not based on a partitioning of the total density, but on a partitioning of the KS orbital space. 20 With such approaches, supermolecular KS-DFT and subsystem calculations will result in exactly the same density, even with finite basis sets. This can be achieved by enforcing the external orthogonality of the subsystem orbitals by using projection operators. It must be noted that in this case the nonadditive kinetic energy vanishes, and its contribution to the embedding potential is replaced by the projection operator enforcing the external orthogonality. Using both a projection operator and an (approximate) nonadditive kinetic-energy contribution in the embedding potential will introduce additional errors because of double counting. 22