Development of a Kohn-Sham like potential in the Self-Consistent Atomic Deformation Model M. J. Mehl, L. L. Boyer and H. T. Stokes 6 9 Complex Systems Theory Branch 9 Naval Research Laboratory 1 Washington, D.C. 20375-5345 n a February 8, 2008 J 9 2 Abstract 1 v This is a brief description of how to derive the local “atomic” potentials from the Self- 3 Consistent Atomic Deformation (SCAD) model density function. Particular attention is paid 0 to the spherically averagedcase. 0 1 0 6 1 The SCAD Model 9 / h The Self-Consistent Atomic Deformation (SCAD) model [1] has been described as an extension of -t the work of Gordon and Kim [2]. In can also be viewed as an approximation to density functional l r theory, which is the point of view taken in this paper. t m The SCAD model approximates the exact electronic density functional [3] by the form: : v Xi E[n(r)] = T0[vi(r)]+ TTF[n(r)]− TTF[ni(r)] +F[n(r)] . (1) ( ) i i r X X a Here: 1. T0[ni(r)] is the kinetic energy of non-interacting electrons centered about the site Ri, T0[ni(r)] = − ψi∗α(r)∇2ψiα(r)d3r , (2) α Z X where the ψ are an orthonormal set of wave functions such that iα n (r) = |ψ (r)|2 . (3) i iα α X 2 (We use atomic units, with h¯ = 2m = e /2 = 1, where m is the mass of the electron and −e is its charge.) 1 2. We assume that the charge associated with each site, N = n (r)d3r , (4) i i Z is fixed. That is, there is no charge transfer between sites during the iteration of the problem to self-consistency. 3. The total density of the system is given by n(r) = n (r+R ) . (5) i i i X 4. The Thomas-Fermi kinetic energy is given (in Rydberg units) by TTF[n(r)] = 3(3π2)23 n35(r)d3r . (6) 5 Z The sum over atomic Thomas-Fermi kinetic energies in (1) is to cancel this contribution to E[n] from the individual atoms, where the kinetic energy is given by (2). The Thomas-Fermi terms in parenthesis in (1) thus represent the overlap kinetic energy, and vanish when the atoms are separated by large distances. 5. The functional F[n(r)] contains the remaining density functional terms[3, 4], including the Coulomb interaction terms (electron-electron, electron-nucleus, and nucleus-nucleus), which is a simple functional of the total density (5), and the exchange-correlation term, which is localfunctionalofthedensitywithintheLocalDensityApproximation(LDA).Thefunctional derivative of F with respect to the total density is just the Kohn-Sham potential [4], δF[n(r)] v [n(r);r] = . (7) KS δn(r) This potential includes both the Coulomb and the exchange-correlation potentials. 2 Finding the potential In this section we wish to find a way to determine the individual charge densities n (r). We begin i by noting that if n is, as we assume, “v-representable”, then there is a one-to-one correspondence i between the density n and a potential v , such that the wave functions described in (2) and (3) i i satisfy the Schro¨dinger equation −∇2ψ (r)+v (r)ψ (r) = ǫ ψ (r) , (8) iα i iα iα iα This is equivalent to rewriting the SCAD density Functional (1) in the form E[{vi(r)}] = T0[vi(r)]+ TTF[n(r)]− TTF[ni(r)] +F[n(r)] . (9) ( ) i i X X 2 Now consider changing the potential in (9) v (r) → v (r)+δv (r) . (10) i i i This change in the potential at the site i changes the associated wave functions ψ (r) → ψ (r)+δψ (r) , (11) iα iα iα where, to first order in δv , i [ǫ +∇2−v (r)]δψ (r) = (δv (r)−δǫ )ψ (r) . (12) iα i iα i iα iα Without loss of generality we may take each δψ to be orthogonal to the corresponding ψ: δψ∗ (r) ψ (r) =0. (13) iα iα Z Multiplying both sides of (12) by ψ∗ (r) and integrating, we obtain the first order change in the iα eigenvalue: δǫ = |ψ (r)|2δv (r)d3r . (14) iα iα i Z The electron density on site i then changes by an amount δn (r) = 2ℜ δψ∗ (r) ψ (r) , (15) i iα iα α X which is a change of O[δv ]. (ℜz is the real part of z.) Note that the densities n on the other sites i j will not change, since by supposition we are only changing the spherical potential on site i, and there is no charge transfer. Thus δn(r) = δn (r+R ) . (16) i i The kinetic energy T0 (2) changes by an amount δT0 = −2ℜ δψi∗α(r)∇2ψiα(r)d3r . (17) α Z X By (8) this becomes δT0 = 2ℜ δψiα(r)[ǫiα−vi(r)]ψiα(r)d3r . (18) α Z X Applying (13) and (15), we find δT0 = − vi(r)δni(r) . (19) Z The other parts of the density functional change in a straight-forward way: δF[n(r)] = v [n(r);r]δn (r)d3r , (20) KS i Z 3 δT [n(r)] = v [n(r);r]δn (r)d3r , (21) TF TF i Z and δT [n (r)] = v [n(r);r]δn (r)d3r , (22) TF i TF i Z where vTF[n(r),r] = δTTδFn([nr)(r)] = (3π2)32n32(r) (23) is the Thomas-Fermi potential. Substitutingequations(19-23)intotheenergyformula(9), wefindthatitchangesbyanamount δE[{v (r)}] = {v [n(r);r]+v [n(r);r]−v [n (r);r]−v (r)} δn (r)d3r . (24) i KS TF TF i i i Z Equation (24) must hold for an arbitrary infinitesimal change in the potential δv . By (15) and i (12), this produces a corresponding change in δn . Hence (24) must hold for arbitrary infinitesimal i changes in δn , which can only occur if the term in curly brackets vanishes, i.e., if i v (r) = v [n(r);r]+v [n(r);r]−v [n (r);r] . (25) i KS TF TF i 3 Spherical Averages When using the SCAD model for crystals with atoms at high symmetry sites, it is often useful to assume that the densities n are spherically symmetric. This can be achieved by assuming that the i n representaclosed shell system, and thateach potential v is spherical. Undertheseassumptions, i i the formalism developed in Section 2 remains the same until equation (24), except that we replace the arbitrary potential v (r) by the spherical potential v (r), and restrict the potential changes to i i the spherical δv (r). At that point we note that because of the assumption that n represents a i i closed shell, δn (r) must also be spherical. Thus (24) can be written i δE[{v (r)}] = {v [n(r);r]+v [n(r);r]−v [n (r);r]−v (r)} δn (r)d3r . (26) i KS TF TF i i i Z This holds for arbitrary spherical potential changes δv (r), whence v (r) must obey i i 1 v (r)= {v [n(r);r]+v [n(r);r]−v [n (r);r]} dΩ , (27) i KS TF TF i 4π Z where the integral is over the surface of a sphere of radius r. Note that this is just the spherical average of the full potential (25): 1 v (r)= v (r)dΩ . (28) i i 4π Z 4 4 Discussion The essential difference between the SCAD model and the more familiar Kohn-Sham LDA formu- lation is that the SCAD model approximates kinetic energy in low-density regions (outside atomic cores) with theaid of a local functional. While the total density in the SCADmodelis expressed as a sum of localized densities, this should not be a major source of error, provided that the localized densities are given adequate flexibility. Solving the atomic calculations self consistently, with the potential for each site determined from the densities and potentials of the rest of the system, puts the atomic energy levels on a common scale. Just as in the Kohn-Sham treatment, the lowest energy levels are occupied to achieve the minimum total energy. The SCAD model offers some advantages over the Kohn-Sham LDA approach. It is generally more efficient, particularly for large systems, because the computational labor increases linearly withthenumberofatomsinthesystem. Moreover, itisideallysuitedforapplication onamassively parallel computer. Finally, we note that the SCAD model has the correct behavior in the limit of infinitely separated atoms, whereas, each electron in the Kohn-Sham model is delocalized and remains so in this limit. 5 Acknowledgment This work was supported by the Office of Naval Research. References [1] L. L. Boyer and M. J. Mehl, Ferroelectrics, 150, 13-24 (1993). [2] R. G. Gordon and Y. S. Kim, J. Chem. Phys. 56, 3122 (1972). [3] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). [4] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 5