# Examination of self interaction correction methods for Na clusters

## I Introduction

The local density approximation (LDA) provides a powerful practical technique to apply the Kohn-Sham framework to interacting many body problems. A problem of this method is unphysical self interaction, i.e. the interaction of a particle with itself. Perdew and Zunger proposed a prescription to remedy this shortcoming, which has been used for atoms, molecules, bulk systems and also metal clusters. Though there still remain some effects of self interaction, the major part of the problem is removed in this method. Compared with the Hartree-Fock theory, which is free from the self interaction problem, the local density approximation with the self interaction correction (SIC) has advantages such as, i) exchange and correlation energies can be relatively easily handled in the same manner, ii) the resultant single particle energies well approximate the physical removal energy from each orbit, iii) the numerical calculation is much lighter, especially for three dimensional calculations.

A characteristic feature of the SIC method of Perdew and Zunger is that the energy functional depends not only on the total density, but explicitly also on the density of each single particle orbital. In almost all calculations for closed shell atoms and metal clusters, the single particle densities in the energy functional are substituted by the spherically averaged densities. The central single particle potential is then deduced by taking functional derivative of the resultant energy functional with respect to the spherically averaged single particle density. Following more closely the original idea of Perdew and Zunger, on the other hand, Harrison  proposed a method of using the original single particle densities without introducing spherical averaging. Since the energy functional is not invariant under unitary transformation of single particle orbitals, Harrison represented the single particle orbitals by either spherical harmonics or Cartesian basis as two choices.

Though the method by Harrison works well for atoms, it has not been tested for metal clusters. We address this question in this paper by taking Na clusters as an example. We show that it does not work well, especially for large clusters which have single particle orbitals with large angular momentum. We confine our study to the exchange energy without referring to the correlation energy in order to make the argument clear and compare the results with those of Hartree-Fock calculations. We show that a better agreement with the Hartree-Fock calculations is obtained if one applies Harrison’s method only to the Hartree term.

In addition to the validity of Harrison’s method, we discuss in this paper the problem of non-diagonal Lagrange multipliers originating from the orthonormality of single particle orbitals in the self interaction correction method of Perdew and Zunger. We show that the non-diagonal property of the Lagrange multipliers introduces a negligible effect to single particle energies as well as the total energy.

The paper is organized as follows. In Sec.II the SIC method of Perdew and Zunger and Harrison’s approach are briefly explained. In Sec.III A the exchange energies calculated by several SIC methods are compared, and the effect of non-diagonal Lagrange multipliers is discussed in Sec.III B. Summary and conclusion are given in Sec.IV.

## II SIC formalism and Harrison’s method

The total energy in the self interaction corrected local density approximation (SIC-LDA) by Perdew and Zunger  is expressed as

$`E_{\mathrm{TOT}}=T+E_{\mathrm{ext}}+E_H+E_X^{\mathrm{SIC}},`$ (1)   
where

$`T={\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}{\displaystyle d^3r\psi _i^{}(𝐫)^2\psi _i(𝐫)},`$ (2)   
$`E_{\mathrm{ext}}[\rho ]={\displaystyle d^3rv_{\mathrm{ext}}(𝐫)\rho (𝐫)},`$ (3)   
$`E_H[\rho ]={\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{\rho (𝐫)\rho (𝐫^{})}{\left|𝐫𝐫^{}\right|}},`$ (4)   
$`E_X^{\mathrm{SIC}}=E_X^{\mathrm{LDA}}[\rho _{},\rho _{}]{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{E_H[\rho _i]+E_X^{\mathrm{LDA}}[\rho _i,0]\right\},`$ (5)   
$`E_X^{\mathrm{LDA}}[\rho _{},\rho _{}]={\displaystyle \frac{3}{2}}\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}{\displaystyle \underset{\sigma =,}{}}{\displaystyle d^3r\rho _\sigma (𝐫)^{4/3}}.`$ (6)   
Here all quantities are in Hartree atomic units, i.e. $`m=e^2=\mathrm{}=1`$. As mentioned in the introduction, the correlation energy has been neglected. We take the effects of ions into account in the spherical jellium model. The external potential is then given by

$`v_{\mathrm{ext}}(𝐫)`$ $`=`$ $`\{\begin{array}{cc}Z/(2R_{\mathrm{jell}})\left\{3\left(r/R_{\mathrm{jell}}\right)^2\right\}\hfill & rR_{\mathrm{jell}}\hfill \\ Z/r\hfill & r>R_{\mathrm{jell}},\hfill \end{array}`$ (9)   
where the jellium radius $`R_{\mathrm{jell}}`$ is related to the number of atoms in the cluster $`Z`$ by $`R_{\mathrm{jell}}=r_sZ^{1/3}`$, $`r_s`$ being the bulk Wigner-Seitz radius which is 4 a.u. for Na.

The Euler equation under the orthonormality condition

$`{\displaystyle \frac{\delta }{\delta \psi _i^{}(𝐫)}}\left\{E_{\mathrm{tot}}+{\displaystyle \underset{ij}{}}ϵ_{ij}\left(\delta _{ij}{\displaystyle d^3r\psi _j^{}(𝐫)\psi _i(𝐫)}\right)\right\}=0`$ (10)   
results in the following coupled equations for the single particle wave functions

$`\left\{{\displaystyle \frac{1}{2}}^2+v_{\mathrm{ext}}(𝐫)+{\displaystyle d^3r^{}\frac{\rho (𝐫^{})}{\left|𝐫𝐫^{}\right|}}+v_X^{\mathrm{S}IC(i)}(𝐫)\right\}\psi _i(𝐫)={\displaystyle \underset{j}{}}ϵ_{ij}\psi _j(𝐫),`$ (11)   
$`v_X^{\mathrm{S}IC(i)}(𝐫)=\left({\displaystyle \frac{3}{\pi }}\right)^{1/3}\rho ^{1/3}(𝐫)\left\{{\displaystyle d^3r^{}\frac{\rho _i(𝐫^{})}{\left|𝐫𝐫^{}\right|}}2\left({\displaystyle \frac{3}{4\pi }}\right)^{1/3}\rho _i^{1/3}(𝐫)\right\}.`$ (12)   
The Lagrange multiplier $`ϵ_{ij}`$ becomes non-diagonal because of the orbital dependence of the self interaction corrected exchange potential. In the following, we approximate it by the diagonal components and solve the following non-coupled equations

$`\left\{{\displaystyle \frac{1}{2}}^2+v_{\mathrm{ext}}(𝐫)+{\displaystyle d^3r^{}\frac{\rho (𝐫^{})}{\left|𝐫𝐫^{}\right|}}+v_X^{\mathrm{S}IC(i)}(𝐫)\right\}\psi _i(𝐫)=ϵ_i\psi _i(𝐫)`$ (13)   
We discuss the validity of this approximation in Sec. III B for Na clusters.

As we see in Eq.(5) the SIC method of Perdew and Zunger is characteristic in that the total energy functional depends explicitly on the individual orbital density. A consequence is that it loses invariance under the unitary transformation of single particle orbitals.

Another problem is that the numerical load is heavy, because one has to solve three-dimensional equations instead of the one-dimensional equations for the radial motion of electrons even for closed shell atoms and metal clusters. In applying this formalism to those systems, one usually replaces $`\rho _i(𝐫)`$ in the curly brackets in Eqs.(5) and (12) by the spherically averaged orbital density given by $`\stackrel{~}{\rho }_i(r)`$

$`\stackrel{~}{\rho }_i(r)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle 𝑑\widehat{𝐫}\rho _i(𝐫)}.`$ (14)   
This prescription certainly reduces the numerical load, because the resultant single particle potentials become central potentials. However, it does not optimize the SIC following the original scheme of Perdew and Zunger.

Harrison proposed an alternative procedure, which avoids replacing the single particle densities in the energy functional by spherically averaged ones. As mentioned in the introduction, he expressed them in the spherical harmonic basis

$`\rho _{nlm}^{\mathrm{SH}}(𝐫)=\left|{\displaystyle \frac{u_{nl}(r)}{r}}Y_l^m(\widehat{𝐫})\right|^2,`$ (15)   
or in Cartesian basis

$`\rho _{nlm}^\mathrm{C}(𝐫)=\{\begin{array}{cc}\left|{\displaystyle \frac{u_{nl}(r)}{r}}Y_l^0(\widehat{𝐫})\right|^2\hfill & m=0\hfill \\ \left|{\displaystyle \frac{u_{nl}(r)}{r}}{\displaystyle \frac{Y_l^{|m|}(\widehat{𝐫})\pm Y_l^{|m|}(\widehat{𝐫})}{\sqrt{2}}}\right|^2\hfill & m0.\hfill \end{array}`$ (18)   
He noticed that the resultant energy functional can be expressed in terms of the spherically averaged orbital density after the integration over angle,

$`E_X^{\mathrm{SIC}}`$ $`=`$ $`E_X^{\mathrm{LDA}}[\rho _{},\rho _{}]{\displaystyle \underset{n,l}{}}2(2l+1)\left\{{\displaystyle \underset{k=0}{\overset{2l}{}}}E_H^{k,l}[\stackrel{~}{\rho }_{nl}]+c_X^lE_X^{\mathrm{LDA}}[\stackrel{~}{\rho }_{nl},0]\right\},`$ (19)   
$`E_H^{k,l}[\stackrel{~}{\rho }_{nl}]`$ $`=`$ $`{\displaystyle \frac{1}{2}}c_H^{k,l}{\displaystyle 𝑑r𝑑r^{}\left|u_{nl}(r)\right|^2\left|u_{nl}(r^{})\right|^2\frac{r_<^k}{r_>^{k+1}}}.`$ (20)   
The coefficients $`c_H^{k,l}`$ and $`c_X^l`$ for each of the spherical harmonic and Cartesian representations are listed in Tables I and II. Harrison then calculated the corresponding central potential for each set of quantum nembers $`n,l`$ by taking the functional derivative with respect to the spherically averaged orbital density $`\stackrel{~}{\rho }_{nl}(r)`$

$`v_X^{\mathrm{S}IC(nl)}(r)`$ $`=`$ $`{\displaystyle \frac{\delta E_X^{\mathrm{SIC}}[\stackrel{~}{\rho }_{nl}]}{\delta \stackrel{~}{\rho }_{nl}(r)}}.`$ (21)   
Though this method restricts the variational space smaller than that in the original scheme of Perdew and Zunger since it presumes a spherically symmetric potential from the beginning, Harrison showed that his method still improves both the exchange energy and the total energy for atoms compared with the simple procedure where $`\rho _i(𝐫)`$ is replaced by the spherically averaged orbital density.

It is an interesting question to see whether Harrison’s method can be applied to metal clusters. A simple minded consideration would suggest that Harrison’s method becomes more powerful in metal clusters. This is because Harrison’s treatment should have a large effect on high angular momentum orbitals which play more important roles in metal clusters than in atoms where the main contribution to the SIC originates from the $`1s`$ state . In the next section, we apply Harrison’s method to Na clusters, and show that it does not work well contrary to the simple expectation.

## III Results and discussion

### A Exchange and total energy

We compare in Table III the exchange energy per electron for Na clusters calculated by several methods. In this table and in what follows, the abbreviations SA-SICX, SH-SICX and C-SICX stand for the SICX calculation using the spherically averaged, spherical harmonic and Cartesian orbital densities, respectively. The difference between the HF and the other calculations represents the error of each method since the HF calculation provides the exact exchange energy. Strictly speaking, the exact exchange energy in the Kohn-Sham formalism is the one given by the optimized effective potential method . It is known, however, that it is nearly the same as that given by the HF calculation for atoms.

We first compare the results of HF, LDAX and SA-SICX for four different Na clusters. The order of the estimated exchange energies, which are negative, is LDAX $`>`$ SA-SICX $`>`$ HF irrespectively of the size of the Na cluster. This is different from the order for atoms, where LDAX $`>`$ HF $`>`$ SA-SICX . The deviation of the result of LDAX from that of HF gets smaller with the size of the Na cluster, while that of SA-SICX is almost constant.

We next compare the results of HF and three SIC methods. The absolute value of the exchange energy calculated by using the spherical harmonic and Cartesian orbital densities becomes smaller than that estimated by the SA-SICX. Including the negative sign, the order is C-SICX $`>`$ SH-SICX $`>`$ SA-SICX. Consequently, the deviation of the results of SH-SICX and C-SICX from that of the HF gets larger than that for SA-SICX. These two calculations are even worse than LDAX for large systems. Their deviation from the HF calculation gets larger with increasing size of the cluster.

For atoms, the order of the exchange energies C-SICX $`>`$ SH-SICX $`>`$ SA-SICX is the same as that for Na clusters. However, the result of C-SICX is still below that of HF. This means that both of SH-SICX and C-SICX have smaller deviation from the HF result than SA-SICX and provide better prescriptions for atoms.

We now investigate the reason why the SH-SICX and C-SICX methods are inferior to SA-SICX for Na clusters. Since the self-consistent density in the three SICX methods is almost the same, the first term of Eq.(19) is nearly the same for the three methods. The major difference among them is therefore associated with the second term of Eq.(19), which consists of contributions from orbitals with various $`n,l`$ quantum numbers. In order to see the physics clearly, we consider the following self interaction averaged over the azimuthal quantum number for each set of quantum numbers $`n,l`$,

$`{\displaystyle \frac{1}{2l+1}}{\displaystyle \underset{m}{}}\left\{E_H[\rho _{lnm}]+E_X^{\mathrm{LDA}}[\rho _{nlm},0]\right\}.`$ (22)   
Its value evaluated by inserting the orbital densities given by Eqs.(14), (15) and (18) is compared in Fig.1 for the occupied orbitals in Na<sub>40</sub>.

The three methods agree well for the $`s`$ orbitals. The small differences originate from a subtle difference of the total densities in the three methods. On the other hand, a large difference appears among them for the finite angular momentum states, i.e. for the $`p`$, $`d`$ and $`f`$ orbitals. Two specific features can be remarked, i) always SA-SICX $`>`$ SH-SICX $`>`$ C-SICX, and ii) the higher the orbital angular momentum is, the larger the differences among them are.

The key point to understand these features is the degree of localization of single particle density adopted in the three SICX methods. Clearly they differ in their angular treatment of the orbital density. C-SICX has the most distinct localization in angle, while SA-SICX has, of course, no angular localization. Symborically, we express this situation of the degree of localization as C-SICX $`>`$ SH-SICX $`>`$ SA-SICX. This difference gets more prominent for higher angular momentum.

The problem is that the angular localization property affects the self Hartree term, i.e. the first term in the curly brackets in Eq.(5), and the self exchange terms, i.e. the second term, in different way. Since the self Hartree term uses the exact long range Coulomb interaction between electrons, it does not depend so much on the density profile of electrons including the angular localization. On the other hand, the self exchange term is evaluated based on the LDA. This is equivalent to assuming a short range $`\delta `$ interaction, which leads to a strong sensitivity of the exchange energy on the details of the density profile of electrons. The more the orbital localizes, the larger the absolute value of the self exchange energy is. In C-SICX, which has the most prominent angular localization, the negative self exchange energy increases with angular momentum and eventually overwhelms the self Hartree term resulting in the sign change of the total self interaction energy. This can be clearly seen in Fig.1. A less prominent, but a similar, trend can be seen also for SH-SICX. However, this strong sensitivity of the self correction energy to the angular localization of electronic density may be unphysical because the LDA for the exchange term cannot be justified for the localized density in SH-SICX and C-SICX. The use of spherically averaged density in the SA-SICX moderates to some extent the overamplification of the angular localization dependendence leading to a better approximation.

The localization of radial wave functions in atoms is very different from that in Na clusters when one compares the states with the same nodal quantum number. For example, $`1s`$ and $`1p`$ orbitals have very similar radial distibution in Na clusters, while the latter ($`2p`$ in the atomic notation) is much more extended than the former in atoms. Consequently, the self interaction correction mainly originates from the most localized $`1s`$ state in atoms. Moreover, only small angular momentum orbitals appear in atoms. Thus, the problem stated above, i.e. the unphysical strong angular localization dependence does not cause so serious trouble in atoms.

To remedy the problem for Na clusters, we propose to utilize spherical harmonic and Cartesian orbital densities only for the self Hartree term. This is a natural consequence of the considerations mentioned above. The exchange energy calculated in this way is given in the lower side in the rows for SH-SICX and C-SICX in Table III. They are designated by the label “(H only)”. All of the H-only exchange energies better agree with that in HF than the corresponding results of the SH-SICX and C-SICX which have been calculated using the spherical harmonic or Cartesian orbital densities both for the self Hartree and exchange terms (the upper side in each row). Especially, H-only SH-SICX is superior to the SA-SICX for most systems. We remark that a similar improvement has been obtained by the H-only method concerning the total energy.

### B Effects of non-diagonal Lagrange multipliers

Before we conclude the paper we comment on the effects of off-diagonal Lagrange multipliers in Eq.(11). We performed calculations by keeping off-diagonal Lagrange multipliers for Na<sub>20</sub> and Na<sub>40</sub> with SA-SICX, SH-SICX and C-SICX. We found that the differences of the total energy and single particle energy obtained by self-consistently solving Eq.(11) and Eq.(13) are less than 10<sup>-1</sup> eV and 10<sup>-2</sup> eV, respectively. We therefore conjecture that one can safely ignore the off-diagonal components of the Lagrange multipliers for Na clusters. This result is consistent with that in Ref.  for atoms.

## IV Summary and conclusion

We have calculated the self interaction corrected exchange energy for Na clusters by using spherically averaged, spherical harmonic and Cartesian orbital densities. We found that both calculations using the spherical harmonic and Cartesian orbital densities deviate from the HF results more than the calculations with spherically averaged orbital densities. The deviation is especially large for systems with large angular momentum orbitals. We attribute this problem to the LDA to evaluate the self exchange term. From this consideration, we propose to use the spherical harmonic and Cartesian orbital deinsities only for the self Hartree term and to use spherically averaged orbital densities for the self exchange term. We have shown that this treatment improves indeed the exchange energy to well reproduce the results of HF calculations. We expect that the remaining errors can be diminished by the SIC using GGA(generalized gradient approximation) .

## Acknowledgement

We are thankful to David M. Brink for enlightening discussions.