Impact Factor 3.560 | CiteScore 3.1
More on impact ›

ORIGINAL RESEARCH article

Front. Phys., 11 March 2021 | https://doi.org/10.3389/fphy.2020.618142

Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities

www.frontiersin.orgMiriam Winkelmann1,2,3,4,5*, www.frontiersin.orgEdoardo Di Napoli1,2,3, www.frontiersin.orgDaniel Wortmann1,2,4 and www.frontiersin.orgStefan Blügel1,2,4
  • 1Institute for Advanced Simulation, Forschungszentrum Jülich, Jülich, Germany
  • 2JARA-CSD, Jülich, Germany
  • 3Jülich Supercomputing Centre, Forschungszentrum Jülich, Jülich, Germany
  • 4Peter Grünberg Institute, Forschungszentrum Jülich, Jülich, Germany
  • 5Physics Department, RWTH-Aachen University, Aachen, Germany

We present a general method for solving the modified Helmholtz equation without shape approximation for an arbitrary periodic charge distribution, whose solution is known as the Yukawa potential or the screened Coulomb potential. The method is an extension of Weinert’s pseudo-charge method [Weinert M, J Math Phys, 1981, 22:2433–2439] for solving the Poisson equation for the same class of charge density distributions. The inherent differences between the Poisson and the modified Helmholtz equation are in their respective radial solutions. These are polynomial functions, for the Poisson equation, and modified spherical Bessel functions, for the modified Helmholtz equation. This leads to a definition of a modified pseudo-charge density and modified multipole moments. We have shown that Weinert’s convergence analysis of an absolutely and uniformly convergent Fourier series of the pseudo-charge density is transferred to the modified pseudo-charge density. We conclude by illustrating the algorithmic changes necessary to turn an available implementation of the Poisson solver into a solver for the modified Helmholtz equation.

1 Introduction

A variety of problems in condensed matter physics require an efficient solution of the partial differential equation

(Δλ2)Vλ=4πρ,(1)

for a charge density ρ in a periodic domain. In real solids, the electronic charge density is a strongly oscillating function in the vicinity of the nuclei of atoms, making a solution in Fourier space, as anticipated by the periodicity, unfeasible due to slow convergence of the Fourier series of the charge density. This equation is frequently referred to as the modified Helmholtz equation or the Yukawa equation. The latter name derives from the Yukawa potential [1], Vλexp(λr)/r, in nuclear physics, which is the underlying free-space Green function of Eq. 1. In the field of condensed matter, e.g. in physics, chemistry, and biology, the Yukawa potential is also known as the screened Coulomb potential. It typically emerges in cases when a many-body system of charged particles is treated in terms of an effective single-particle theory applying a mean-field approximation. Then the many particles contribute to an effective screening of a Coulomb interaction generated by the single, representative charged particle when treated in linear response theory.

The relation between the bare Coulomb potential on the one hand and the screened Coulomb potential or the induced screening charge on the other hand is referred to as the dielectric constant or susceptibility, respectively. Depending on the context, such relations appear in the Debye–Hückel theory [2] in the form of a linearization of the Poisson–Boltzmann equation, where the Poisson equation describing the electrostatics of charged particles is a function of the charge density distribution obeying a Boltzmann statistics. Another example is the Thomas–Fermi model [3, 4] of the dielectric constant in metals, which describes the screening potential due to the linearized change of the electron distribution described by the Fermi–Dirac distribution with respect to the spatial variations of the electrostatic potential. In these theories, the constant λ represents the inverse of a typical length scale over which an individual charged particle exerts a notable effect.

The Thomas-Fermi theory can be regarded as a precursor of the density functional theory [5] (DFT). The latter is the most important theory and methodology for the modeling and simulation of material properties of a crystalline solid based on the quantum mechanical treatment of many electron systems. In addition, the Thomas-Fermi theory provides a rough but fast approximation of the common density functionals, which relate the electron density to the effective Kohn–Sham potentials [6]. Such a scheme makes the solution of Eq. 1 particularly valuable. For instance it could be used to obtain a good starting potential for the iterative solution of the Schrödinger-like Kohn–Sham equations, where the nuclear charge is included in ρ. Other examples are the attainment of an efficient approximate solution of the dielectric function, or the implementation of a hybrid functional [7] to DFT using the Yukawa screening of the Hartree–Fock exchange. In both cases the charge density ρ in Eq. 1 is replaced by an overlap charge density [8] obtained as a product of wave functions associated with different quantum numbers. Since the Thomas–Fermi model approximates the description of the response behavior of delocalized electrons in solids very well, Eq. 1 is frequently used to determine a preconditioner for the acceleration of the self-consistent solution of the DFT equation. So far the Kerker preconditioner [9] is applied to density functional methods with Fourier transformable charge density. By solving Eq. 1 for general periodic densities, the Kerker preconditioner [9] can be extended to density functional methods of general densities [10, 11].

Although most electronic structure methods implementing DFT applied to solid-state materials systems make explicit use of the underlying periodicity of the crystalline lattice, a straightforward solution of Eq. 1 using Fourier transformation techniques is in general not possible due to the strongly oscillating charge density close to the nuclei. This problem is well discussed for the solution of the Poisson equation, ΔV=4πρ, a limit of the modified Helmholtz equation for λ=0.

In a seminal work, Weinert [12] proposed an elegant and numerically efficient solution of the Poisson equation for periodic charges and corresponding electrostatic potentials without shape approximation. Weinert’s solution, to which we refer here as Weinert’s pseudo-charge method, is implemented (in several variants) in most full-potential all-electron DFT methods, such as the augmented spherical wave (ASW) method [13], the Korringa–Kohn–Rostoker Green function (KKR-GF) method [14], and the full-potential linearized augmented planewave (FLAPW) method [15], just to name a few.

Typical to these all-electron DFT-methods is the domain decomposition into atomic spheres around the atoms and an interstitial region in-between. Weinert’s pseudo-charge method is based on the observation that the relation between the charge density inside a sphere and its multipole expansion outside the sphere is not unique. A smooth Fourier transformable pseudo-charge density with the same multipole moments as the true density is constructed. The latter provides the true potential through Fourier transformation of the Poisson equation and a subsequent solution of a Dirichlet boundary value problem on the sphere boundary.

In this article, we extend Weinert’s pseudo-charge method to the modified Helmholtz Eq. 1 for values of λ>0 in order to determine the Yukawa potential Vλ for general periodic charge densities without shape approximation. We formulate the new method for general charge densities, including continuous charge densities as for electron densities, discrete charge densities as for nuclear charges or more abstract densities that arise of products of wave functions. Such an approach is consistent with the real-space representation of the charge density and potential in all-electron methods.

As a matter of choice, and motivated by the original work of Weinert [12], we demonstrate this extension explicitly for the FLAPW method [15] as implemented in the Fleur code [16]. We provide i) complete derivations of the modified multipole expansion of the charge density inside the atomic spheres using a Green function method, as well as the interstitial charge density’s modified multipole moments in the atomic spheres using some Bessel function integration properties, which yields the coefficients of the pseudo-charge density, and ii) a mathematical analysis of the convergence of the pseudo-charge density’s Fourier series. Our results confirm all expressions for all the quantities provided by Tran et al. [7] relevant in an actual computation. Beyond Tran et al. [7] we provide detailed derivations and justify the use of identities such as Proposition 1 for the modified spherical Bessel function i and the transferability of Weinert’s convergence analysis to the case λ>0. We also discuss and provide further insight into the convergence properties of the Fourier series, and in addition restructure and simplify Weinert’s approach, hopefully making it more comprehensible. We point out the algorithmic changes required to extend the solution of the Poisson equation to a solution of Eq. 1, which can then be straightforwardly transferred to other all-electron full-potential band structure methods.

This paper is organized as follows: In Section 2.1, we introduce the muffin-tin and interstitial region typical of the FLAPW method and the corresponding domain decomposition for the charge density and the potential. We summarize the main statements of Weinert’s pseudo-charge method and give a definition of the pseudo-charge density. Since we know the true charge density inside the muffin-tin spheres and with the assumption that we would know the interstitial potential, we construct in Section 2.2 the Yukawa potential inside the sphere by solving the Dirichlet boundary value problem. We develop two radial Green functions that are products of two linearly independent fundamental set solutions of the homogeneous radial modified Helmholtz equation. These Green functions are set apart by the boundary conditions they fulfill either at the muffin-tin sphere or in free-space. In Section 2.2.1, the radial free-space Green function is used to define the modified multipole expansion of the Yukawa potential. In Section 2.3, we construct a pseudo-charge density in reciprocal space consistent with the modified Helmholtz equation by making use of the definition of the modified multiple moments put forward in Section 2.2.1. We obtain the Yukawa potential for the interstitial region by solving the modified Helmholtz equation in Fourier space for the pseudo-charge density—the solution is a simple algebraic expression. This is followed by an analysis of the convergence properties of the Fourier series of the pseudo-charge density. The entire algorithm that solves the modified Helmholtz equation is summarized in Section 2.4 together with the minimal modifications necessary to change Weinert’s original algorithm. The conclusions and the outlook are presented in Section 3.

2 Yukawa Potential for a Muffin-Tin Decomposition of a 3D-Periodic Domain

2.1 Weinert’s Pseudo-Charge Method

In order to deal with the 1/r singularities of the Coulomb potential due to the point-like charge of the nucleus and the associated rapid oscillations of the charge density in the vicinity of the singularity, in all-electron electronic structure methods the space is typically partitioned into muffin-tin spheres BRα(τα) of radius Rα centered around the atoms α—the union of those is called the muffin-tin (MT) region—and the interstitial region (I) between the atoms. In FLAPW both charge densities

ρ(r)={KρI(K)eiKrrILρLα(rα)YL(r^α)r=τα+rαBRα(τα)(2)

and potentials

V(r)={KVI(K)eiKrrILVLα(rα)YL(r^α)r=τα+rαBRα(τα)(3)

are represented in plane waves eiKr, where K defines the reciprocal lattice vector dual to the lattice vectors defining the periodic domain, and in spherical harmonics, YL, of degree and order m, where L is defined as L:=(,m). rαRα is the length of the vector rα=rτα, measured from the center of the atom α placed at position τα in the periodic domain, with r^α=rτα/|rτα| its unit vector. The precision of the representation is determined by the cut-off parameters Kmax for the wave vectors, K, with length KKmax, and max for the degree in the angular-momentum expansion. max sets also a natural cut-off of the angular-momentum expansions of all other charge densities or multipole moments throughout the paper.

Weinert’s pseudo-charge method for the Poisson equation is based on the crucial observation that several charge densities ρ inside a sphere BRα(τα) can generate the same multipole moments

qLα[ρ]=BRα(0)ρ(rα+τα)rαYL*(r^α)drα(4)

and thus, the same potential

VI[ρ](r)=L4π2+1qLα[ρ]1rα+1YL(r^α)(5)

outside the sphere. Here, YL* denotes the complex conjugate of YL. The pseudo-charge density, ρ˜, defined by Weinert in Ref. [12] is such a charge density. It fulfills the following three conditions:

• It has the same multipole moments qLα[ρ˜]=qLα[ρ] in every sphere BRα(τα).

• It is equal to the true charge density ρI in the interstitial region.

• It has a fast convergent Fourier expansion.

The Fourier Components of VI are then simply

VI(K)=4πK2ρ˜(K)forK0,(6)

while VI(0) will be set to a constant. Once the interstitial potential VI has been calculated, the muffin-tin potential can be obtained by solving the Dirichlet boundary value problem on the sphere

Vα(rα+τα)=VSα(rα+τα)+VBα(rα+τα)=BRα(0)G(rα,rα)ρ(rα+τα)drαRα24πBRα(0)VI(rα+τα)Gn(rα,rα)dω,(7)

where G is a Green function associated with the solution of the Poisson equation, dω=sinθdθdϕ denotes the solid angle element and r=rα+ταBRα(τα). Although the Green function depends on the muffin-tin radius Rα, for simplicity we drop the index α in the Green function and in related quantities. The muffin-tin potential Vα is fed by two terms, a source term VSα due to the charge density distribution inside the sphere and a boundary term VBα due to the interstitial potential at the boundary of the sphere. The Fourier coefficients of the pseudo-charge density basically have the form

ρ˜(K)=ρI(K)+αρ¯α(K),(8)

where ρ¯α is a Fourier transformable pseudo-charge density inside the muffin-tin sphere. The idea behind this is the following: if the domain of definition of ρI is formally expanded to the full space, i.e. including the muffin-tin spheres, such that

ρI(r)=KρI(K)eiKr(9)

can also be evaluated for rαBRα(τα), then the true charge ρ can also be written as

ρ=ρI+αρ˘α,(10)

where

ρ˘α:={0inIραρIinBRα(τα)(11)

are charge densities localized in the atomic spheres BRα(τα). If these localized densities are now substituted by other localized densities ρ¯α, then the charge density is still correct in I.

The approach above is generally the same for the Yukawa potential, i.e., for λ>0, only Eqs. 46 and the pseudo-charge density’s Fourier coefficients are slightly different. We now continue with the derivation of the muffin-tin and interstitial potentials in the case λ>0.

2.2 Muffin-Tin Yukawa Potential

Assume the interstitial potential is obtained, then the Green function method is used to determine the screened Coulomb potential Vλα inside the muffin-tin sphere BRα(τα) with centre τα through the solution of the boundary value problem

(Δλ2)Vλα=4πρin BRα(τα)(12)
Vλα=VλIonBRα(τα).(13)

The solution is divided into three steps, which we describe in the following. The derivation focuses on the construction of the Green function and its application. For simplicity, we leave out the index λ in Green functions and potentials in this subsection. The solution is given in terms of radial functions VLα(rα), the expansion coefficients to the spherical harmonics expansion inside the sphere (see Eq. 3).

Step 1. We solve the homogeneous modified Helmholtz equation, (Δλ2)U=0, in spherical coordinates. Following the solution [17] of the Laplace equation, Δψ=0, in spherical coordinates (r,r^), the homogeneous potential can be factorized into products of radial functions u(r) and angular functions YL(r^), U(r,r^)=Lu(r)YL(r^). The term λ2U in the homogeneous modified Helmholtz equation only has an effect on the radial solution. The spherical harmonics, YL, are the eigensolutions of the angular part of the Laplace equation with eigenvalues (+1). The radial part of the homogeneous modified Helmholtz equation is known as the modified spherical Bessel differential equation [17, 18],

d2u(r)dr2+2rdu(r)dr((+1)r2+λ2)u(r)=0,(14)

and its fundamental set of solutions are for each the two modified spherical Bessel functions [17, 18] i(λr) and k(λr), the first of which is the regular solution well-defined at the origin, but grows fast with growing radius r and the second is the irregular solution that goes to infinity for r0. To realize the proper boundary condition for the radial Green functions two conditions have to be fulfilled: The first solution, u1, must be finite at r=0. We conclude that

u1(r)=i(λr),(15)

since k(λr) for r0. The second solution, u2, must be 0 at r=Rα. This is achieved by a linear combination of the two modified spherical Bessel functions,

u2(r)=k(λr)i(λr)k(λRα)i(λRα).(16)

Step 2. A function GC0([0,Rα]×[0,Rα])C2([0,Rα]×[0,Rα]\{(r,r)|r[0,Rα]}) is called a radial Green function, if it is the solution to

(Δrλ2)G(r,r)=4πr2δ(rr)(17)

subject to the Dirichlet boundary condition

G(r,r)=0,ifr=Rαorr=Rα,(18)

where (Δrλ2) is the linear radial differential operator in Eq. 14 and δ denotes the radial Dirac delta function, for which

abδ(rr)f(r)dr={f(r),if r[a,b]0,otherwise. (19)

The radial Green function takes the form of the product of the two linearly independent solutions with the proper boundary conditions,

G(r,r)=Cu1(r<)u2(r>),(20)

where r<=min(r,r), r>=max(r,r) and

C=4πr2W1[u1(r),u2(r)].(21)

Since the Wronskian W is linear, the addition of cu1 (c = const) onto k(λr) to suffice the boundary condition u2(Rα)=0 has no influence on the Wronskian: W(u1,cu1)=0, with W(i,k) remaining. To calculate the Wronskian of i(λr) and k(λr) one can either Taylor expand the two functions or simply take the limiting values for r0 or the asymptotic values for r to find

W[u1(r),u2(r)]=u1(r)du2(r)drdu1(r)dru2(r)=1λr2.(22)

Therefore, C=4πλ, and the radial Green function finally reads

G(r,r)=4πλi(λr<)k(λr>)4πλi(λr)i(λr)k(λRα)i(λRα).(23)

Step 3. Considering the standard expression of Dirac’s delta function separated according to the radial and angular coordinates

δ(rr)=1r2δ(rr)δ(r^r^)=1r2δ(rr)LYL*(r^)YL(r^),(24)

the three-dimensional (3D) Green function G(r,r) solving

(Δλ2)G(r,r)=4πδ(rr)in BRα(0),(25)
G(r,r)=0on BRα(0),(26)

is expanded in the form

G(r,r)=LG(r,r)YL*(r^)YL(r^)(27)

and the solution to the inhomogeneous Eq. 12 is given by Eq. 7. For the derivation of the 3D Green function and the 3D inhomogeneous solution Vα=VSα+VBαEq. 7 we refer the reader to Ref. [19]. Both the integral over the sphere BRα(0) and the boundary integral simplify by exploiting the orthonormality relation of the spherical harmonics, B1(0)YL*(r^)YL(r^)dω=δLL. The integral over the sphere BRα(0) provides the source contribution to the muffin-tin potential

VSα(rα+τα)=BRα(0)G(rα,rα)ρ(rα+τα)drα=L[0RαG(rα,rα)ρLα(rα)rα2drα]YL(r^α).(28)

In order to obtain the boundary contribution to the muffin-tin potential, we evaluate the boundary integral in Eq. 7 by expanding the interstitial potential VI(K) (see Section 2.3) on the sphere boundaries BRα(τα)r in spherical coordinates

VI(rα+τα)=KVI(K)eiKταeiKrα=LVLI(Rα;τα)YL(r^α)(29)

using the plane-wave expansion

eiKr=L4πij(Kr)YL*(K^)YL(r^),(30)

where

VLI(Rα;τα)=4πiKVI(K)eiKταj(KRα)YL*(K^).(31)

Furthermore, the normal derivative of G on the sphere boundary is

Gn(rα,rα)=G(rα,rα)rα|rα=Rα=LG(rα,rα)rα|rα=RαYL*(r^α)YL(r^α).

Since rα<Rα=rα and since G takes the form Eq. 20, we obtain

G(rα,rα)rα|rα=Rα=4πλu1(rα)u2(Rα).(33)

We recall that u2(Rα)=0 and reuse Eq. 22 to obtain

u2(Rα)=1λRα2i(λRα),(34)

yielding

G(rα,rα)rα|rα=Rα=4πRα2i(λrα)i(λRα).(35)

With this and the knowledge of the interstitial potential at the sphere boundary from Eq. 31, the boundary contribution to the muffin-tin potential becomes

VBα(rα+τα)=Rα24πBRα(0)VI(rα+τα)Gn(rα,rα)dω=LVLI(Rα;τα)i(λrα)i(λRα)YL(r^α)(36)

and the radial part of the spherical harmonics expansion of the total potential, VSα+VBα, in the sphere BRα(0) becomes

VLα(rα)=0RαGα(rα,rα)ρLα(rα)rα2drα+VLI(Rα;τα)i(λrα)i(λRα).(37)

Due to the kink of G at rα=rα, for practical calculations the integral is split in a part where rα<rα, a part where rα>rα and a third part where the integrand is symmetric in rα and rα:

VLα(rα)=4πλ([0rαρLα(rα)i(λrα)rα2drα]k(λrα)+[rαRαρLα(rα)k(λrα)rα2drα]i(λrα)[0RαρLα(rα)i(λrα)rα2drα]i(λrα)k(λRα)i(λRα))+VLI(Rα;τα)i(λrα)i(λRα).(38)

2.2.1 Modified Multipole Expansion

In the same way we obtain the radial representation of the free-space Green function, well-known as the Yukawa potential for a Dirac test charge at r,

eλ|rr||rr|=4πλLi(λr<)k(λr>)YL*(r^)YL(r^).(39)

The modified spherical Bessel function k already contains the proper boundary condition for r. A charge density ρ˘α localized in a sphere BRα(τα) embedded in free space, produces a Yukawa potential outside the sphere, i.e., r=rα+ταBRα(τα),

VI[ρ˘α](r)=3G(rα,rα)ρ˘α(rα+τα)drα,(40)

which can be expressed analogously to the Coulomb potential in terms of the modified multipole expansion

VI[ρ˘α](r)=L4πλ+1(2+1)!!qLα[ρ˘α]k(λrα)YL(r^α),(41)

with the modified multipole moments

qLα[ρ˘α]=(2+1)!!λBRα(0)ρ˘α(rα+τα)i(λrα)YL*(r^α)drα.(42)

An analogous definition holds true for the modified multiple moments qLα[ρα] of the true charge ρα in the sphere. With the standard expansion of the charge density inside the sphere into spherical harmonics Eq. 3, ρα(rα+τα)=LρLα(rα)YL(r^α) and the application of their orthonormality relation, the calculation of the modified multipole moments inside the muffin-tin spheres is straightforward and results in

qLα[ρα]=(2+1)!!λ0RαρLα(rα)i(λrα)rα2drα.(43)

The summation of VλI[ρ˘α] over all spheres α finally provides the contributions of the charges of all spheres to the interstitial potential. Equation 7 reduces to Eq. 40, since for the free-space Green function Eq. 39 the boundary value term disappears for r.

2.3 Interstitial Yukawa Potential

Suppose we had found a Fourier transformable pseudo-charge density ρ¯α inside the sphere consistent with the Yukawa potential produced outside the sphere, with coefficients ρ¯α(K), and the Fourier series would converge rapidly throughout the periodic domain. Then we can find the Fourier coefficients of the pseudo-charge density, ρ˜(K), by Eq. 8 and the solution of the modified Helmholtz Eq. 1 through Fourier transformation yields an algebraic equation from which we calculate the interstitial Yukawa potential,

VλI(K)=4πK2+λ2ρ˜(K).(44)

In the previous Section 2.2 the interstitial Yukawa potential is used as boundary values for the Yukawa potential in the atomic spheres.

This subsection is concerned with the construction of the Fourier transformable pseudo-charge density ρ¯α that replaces the true local charge density ρ˘α [see (11)] inside the muffin-tin sphere such that the Yukawa potential in the interstitial region, VλI[ρ˘α]=VλI[ρα]VλI[ρI], due to the true charge density inside the sphere, is equal to the Yukawa potential in the interstitial region produced by the a priori unknown pseudo-charge density, VλI[ρ¯α]=VλI[ρ˘α]. From Eq. 41 we conclude that this is fulfilled if the modified multipole moments Eq. 42 of both charge densities, qLα[ρ˘α]=qLα[ρα]qLα[ρI] and qLα[ρ¯α], are equal. The modified multiple moments qLα[ρα] are already known through Eq. 43. Next we determine the modified multiple moments qLα[ρI] of the interstitial charge extended into the muffin-tin spheres and then construct the pseudo-charge density.

2.3.1 Modified Multipole Moments of Interstitial Charge Density Extended Into Sphere

The determination of the modified multiple moments of the interstitial charge density is in principle the same as in Ref. [12], but since the modified multipole moments are different from the known multipole moments for the Coulomb potential, we go through this step of deriving qLα[ρI] in detail. We write ρI relative to the sphere centre τα and employ the Rayleigh expansion Eq. 30 to eiKrα, which yields

ρI(r)=KρI(K)eiKrαeiKτα=KρI(K)eiKταL4πij(Krα)YL*(K^)YL(r^α).(45)

The modified multipole moments of ρI in the sphere BRα(τα), defined analogously to Eq. 42, are

qLα[ρI]=(2+1)!!λBRα(0)YL*(r^α)i(λrα)ρI(rα+τα)drα=(2+1)!!λ4πiKρI(K)eiKταYL*(K^)0Rαi(λrα)j(Krα)rα2drα.

For K0, the latter integral becomes

0Rαi(λrα)j(Krα)rα2drα=Rα2K2+λ2(Ki(λRα)j+1(KRα)+λi+1(λRα)j(KRα))=Rα2K2+λ2(λi1(λRα)j(KRα)Ki(λRα)j1(KRα)).(46)

If K=0, observe that j(0)=δ0 (Kronecker δ) and so the integral becomes

δ00Rαi0(λrα)rα2drα=Rα3i1(λRα)λRαδ0.(47)

Both equations above can be derived by partial integration in two different ways and applying the identities in Ref. [18],

ddr(rf(r))=±rf+1(r)(48)

for f=i with the plus sign and for f=j with the minus sign, and

ddr(r+2f+1(r))=r+2f(r)(49)

for both f=i and f=j. In conclusion, this yields

qLα[ρI]=δ04πRα2i1(λRα)λρI(0)+K0(2+1)!!λ4πiρI(K)eiKταYL*(K^)Rα2λ2+K2(Ki(λRα)j+1(KRα)+λi+1(λRα)j(KRα)).(50)

2.3.2 Construction of Pseudo-Charge Density

We construct the pseudo-charge density by following the ansatz of Weinert,

ρ¯α(rα+τα)=Lρ¯Lα(rα)YL(r^α)=LQLnα(η=0naηrανη)YL(r^α),(51)

in which the radial dependence of the charge density is expressed in terms of a polynomial expansion up to degree νn, which depends on atom α and angular degree , and otherwise use spherical harmonics for the angular part—this being the usual representation for charge densities in the muffin-tin region. As we will discuss in Section 2.3.3, it is beneficial to choose

aη=(1)nηRα2(nη)(nη)anfor η=0,,n(52)

and νη=+2η, where n is yet to be determined. As will become apparent later when we derive the coefficients QLnα in Eq. 61, the coefficient an cancels out in any relevant equation, like Eq. 54 or Eq. 55. With these choices of parameters and with the binomials theorem applied to

η=0n(1)nη(nη)(rαRα)2η=((rαRα)21)n(53)

it follows from the ansatz Eq. 51

ρ¯α(rα+τα)=an(rα2Rα2)nLQLnαrαYL(r^α).(54)

The Fourier transform of this expression is then given by

ρ¯α(K)=1|Ω|eiKταBRα(0)ρ¯α(rα+τα)eiKrαdrα=4π|Ω|eiKταL(i)QLnαAnα(K)YL(K^),(55)

where

Anα(K)=an0Rα(rα2Rα2)nrα+2j(Krα)drα,(56)

|Ω| is the volume of the periodic domain and we used the Rayleigh expansion Eq. 30 and the orthonormality relation of the spherical harmonics. With Proposition 1 in the appendix Anα(K) finally reduces to

Anα(K)=an(2)nn!Rα+n+2j+n+1(KRα)Kn+1.(57)

In the same way we derive the coefficients QLnα: We insert Eq. 54 in the definition of the modified multipole moments Eq. 42 and use Eq. 1 for f=i and κ=λ to obtain

qLα[ρ¯α]=(2+1)!!λBRα(0)ρ¯α(rα+τα)i(λrα)YL*(r^α)drα(58)
=(2+1)!!λQLnαan0Rα(rα2Rα2)nrα+2i(λrα)drα(59)
=(2+1)!!λQLnαan(2)nn!Rα+n+2i+n+1(λRα)λn+1,(60)

and thus,

QLnα=qLα[ρ¯α]λ+n+1an(2)nn!Rα+n+2(2+1)!!i+n+1(λRα).(61)

Since Anα(K) and QLnα share the term an(2)nn!Rα+n+2, the term cancels out in the product QLnαAnα(K) entering Eq. 55,

QLnαAnα(K)=j+n+1(KRα)Kn+1(2+1)!!λ+n+1i+n+1(λRα)qLα[ρ¯α],(62)

and setting ν=+n+1 (in accordance with Section 2.3.3) this leads to

ρ¯α(K)=4π|Ω|eiKταL(i)jν(KRα)Kν(2+1)!!λνiν(λRα)qLα[ρ¯α]YL(K^).(63)

For ρ¯α(0) we take the limit

limK0jν(KRα)Kν=limK0KRαν(2ν+1)!!,(64)

which is only different from 0 for =0 and thus yields

ρ¯α(0)=4π|Ω|(λRα)ν(2ν+1)!!iν(λRα)q00α[ρ¯α].(65)

Due to the condition that the pseudo-charge density has the correct modified multipole moments, the modified multipole moments used for the actual computation of the Fourier coefficients are the ones of the true localized charge density ρ˘α, qLα[ρ˘α]=qLα[ρα]qLα[ρI] calculated from the modified multipole moments qLα[ρα]Eq. 43 and qLα[ρI]Eq. 50 of ρα and ρI respectively, in the sphere BRα(τα).

2.3.3 Smoothness of the Pseudo-Charge Density and Convergence of Its Fourier Series

In Section 2.3.2, we have set aη by Eq. 52, νη=+2η and ν=+n+1 without having determined n yet. Here we motivate our choices and finally determine a proper n.

With our choices for aη and νη we have eradicated the sum in ansatz Eq. 51 and derived the much simpler form Eq. 54. The function

rα(rα2Rα2)n(66)

itself and all its first n1 derivatives with respect to rα are equal to zero at rα=Rα. Consequently, we ensure smoothness on the boundary of the sphere,

dkdrαkρ¯Lα(rα)|rα=Rα=0k=0,,n1.(67)

Note that this includes localization of ρ¯α in BRα(τα) and thus the pseudo-charge density equals the charge density in I (condition 2 in Section 2.1).

The smoothness of the pseudo-charge density is connected to the convergence properties of its Fourier series. Applying the Riemann–Lebesgue lemma for Fourier series to a n1-fold differentiable function in combination with the differentiation rule for a Fourier transform, one can show that the Fourier coefficients for K go faster to zero than 1Kn1, i.e., we obtain the fastest convergence of the Fourier series for large K and the convergence becomes the better the larger n is. If we choose n too large, however, we are left with the small-K Fourier coefficients only, and thus the Fourier series is unbalanced, in the sense that smaller coefficients have a larger weight. So, ideally, our choice of n is guided by the cut-off of the Fourier series.

For an explicit rule on how to choose n, Weinert [12] discussed the factors QLnαAnα(K)/qLα[ρ¯α], where in the Coulomb case

QLnαAnα(K)=j+n+1(KRα)Kn+1(2+1)!!(2+2n+3)!!Rα+n+1qLα[ρ¯α].(68)

His reasoning is based on the zeros of the K-dependent function QLnαAnα(K)/qLα(ρ¯α). Viewed as a function of K, however, our factor QLnαAnα(K)/qLα[ρ¯α] differs from Weinert’s one only by a multiplicative constant. Thus Weinert’s arguments apply here as well. With Proposition 2 in the Appendix it follows that this multiplicative constant is smaller than 1 for λ>0. We see this confirmed in Figure 1, which shows the KRα-dependence of the factor QLnαAnα(K)/qLα[ρ¯α] in the Yukawa Eq. 62 and Coulomb Eq. 68 cases for several combinations of n and angular degree it reveals a smaller amplitude in the Yukawa case. In agreement with Ref. [12], in Figure 1 we make two observations: (i) QLnαAnα/qLα[ρ¯α] as a function of KRα has larger oscillations for smaller n, and (ii) for fixed n, the largest contribution to the Fourier series comes from KRα less than the first zero of the Bessel function j+n+1. Since we deal with a finite number of K vectors, we would like to reduce the oscillations mentioned in (i) by choosing a large n. On the other hand, due to the cut-off of the Fourier series at some Kmax, the factor QLnαAnα/qLα[ρ¯α] must be small for K>Kmax, which limits +n+1 to a certain value, since the first zero is pushed towards infinity for growing +n+1, as can be seen from a comparison between the pink and yellow, or the blue and purple lines in Figure 1. From this arises Weinert’s criterion for choosing n(), which we adopt here:

• Choose ν such that the first zero of jν(z) is approximately equal to (KRα)max.

• Then n() is fixed by the relation ν=+n+1.

FIGURE 1
www.frontiersin.org

FIGURE 1. The factor QLnαAnα/qLα[ρ¯α] for several (,n) in the Yukawa and Coulomb cases with Rα=1 and λ=2.

Note that in this method is compensated by n in such a way that ν is de facto not depending on . Since the discretization of K vectors and the muffin-tin radii usually do not change over the course of the self-consistent-field iteration, the terms in Eq. 63 depending on ν need to be computed just once.

2.4 Algorithm: Construction of Yukawa Potential

Algorithm 1 summarizes the construction of the Yukawa potential derived in this paper.

Algorithm 1Bulk-Case Yukawa Potential

Input: charge density ρ, integer ν chosen as described on page 8 and preconditioning parameter λ.

Output: Yukawa potential Vλ solving the modified Helmholtz Eq. 1 with periodic boundary conditions.

Pseudo-Charge Densityρ˜ρ

1: Modified multipole moments qLα[ρI] of the interstitial charge density in BRα(τα). Eq. 50

2: Modified multipole moments qLα[ρα] of the muffin-tin charge density in BRα(τα). Eq. 43

3: Modified Multipole Moments qLα[ρ¯α]=qLα[ρα]qLα[ρI] of ρ¯α

4: Sphere-localized part ρ¯α(K) of the pseudo-charge density. Eqs. 63 and 65

5: Pseudo-charge density ρ˜(K). Eq. 8Interstitial potentialVλIρ˜

6: Interstitial potential VλI(K). Eq. 44Muffin-tin potentialVλMTρMT,VλI

7: Boundary terms VLI(Rα;τα) of muffin-tin potential. Eq. 31

8: Radial parts VLα(rα) of muffin-tin potential. Eq. 38

In the case that Weinert’s method is available as an implemented algorithm, then only relatively few changes are necessary to make it available for the solution of the modified Helmholtz equation. The changes to be made in practice are limited to the following: The slightly different radial behavior of the Green function leads to small changes in the multipole moments of the interstitial and muffin-tin charge densities in each sphere, qLα[ρI]Eq. 50 and qLα[ρα]Eq. 43 respectively, and in the Fourier components of the pseudo-charge density ρ¯α(K)Eqs. 63 and 65. The integer ν in the formula for the pseudo-charge density’s Fourier components, which determines the convergence of the Fourier series, is chosen exactly the same as in Weinert’s original method. The interstitial potential Eq. 44 undergoes changes indirectly through the pseudo-charge density and directly by the prefactor 4π/K2+λ2) that substitutes 4π/K2. Since the K = 0-term is well-defined, it is not set to a constant as in the original method. The muffin-tin potential is affected only in its radial dependence in both the boundary and the source contribution of Eq. 2.2. Basically, the polynomials r and 1/r+1 in these quantities are substituted by the modified spherical Bessel functions i(λr) and k(λr), respectively, and the prefactor 4π/(2+1 is substituted by 4πλ. The interstitial potential on the boundary of the spheres, VLI(Rα;τα)Eq. 31, only changes indirectly through the changed values of VλI(K).

3 Conclusion and Outlook

We have presented a general method for solving the modified Helmholtz equation for a 3D-periodic system of charge densities not restricted by any shape approximation of three-dimensional volume. The three-dimensional domain is typically decomposed into non-overlapping atom-centered spheres (the muffin-tin region) and the space between these spheres (the interstitial region). The solution, the Yukawa potential, suffices 3D-periodic boundary conditions. Since the Yukawa differential equation is similar to the Poisson equation, we leveraged our derivations on the work of Weinert [12]. Our work can be considered an extension of Weinert’s work determining instead of the bare Coulomb potential with zero screening, the screened Coulomb potential with finite screening length 1/λ. Like Weinert’s pseudo-charge method, our extension is based on the concept of the non-uniqueness of the multipole expansion as well as on the Dirichlet boundary value problem applied to a sphere. The difference between the modified Helmholtz equation and the Poisson equation lies solely in the radial behavior and thus the homogeneous solutions to the radial part of the differential equation are modified spherical Bessel functions instead of polynomial functions. The consequence is a different radial behavior of the Green function, resulting in the screening of the potential, which is now expanded in modified multipole moments, and this in turn affects the pseudo-charge density. Furthermore, the modification of the multipole moments implies that the modified monopole is not connected anymore to the total charge. We have shown that Weinert’s convergence analysis of an absolutely and uniformly convergent Fourier series of the pseudo-charge density is transferred to the modified pseudo-charge density and thus we can therefore best choose the same integer parameters for convergence. Finally we layed out the minor changes necessary to change an implemented method for solving the Poisson equation available to an implementation for solving the modified Helmholtz equation.

Considering that Weinert’s pseudo-charge method has become the standard method for calculating the electrostatic potential without shape approximation in all-electron band structure methods for applications of periodic solids, and since we have extended it to a modified pseudo-charge method with only minor modifications involving some radial integrals, this now allows to treat the screened Coulomb potentials without shape approximation described by the modified Helmholtz equation with all-electron methods. The screened Coulomb or Yukawa potential typically occurs in single particle or mean-field theories to the problem of many charged particles, where all charged particles contribute effectively to the screening of the bare Coulomb potential.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

MW did all the derivations, provided the proof and the figure. EN and DW supervised the work. SB conceptualized the work. All authors contributed through discussions and to the writing of the paper.

Funding

This work has been supported by a JARA-HPC seed-fund project and by the MaX Center of Excellence funded by the EU through the H2020-INFRAEDI-2018-1 (Grant No. 824143). The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC partition’s part of the supercomputer JURECA 20 at Forschungszentrum Jülich, as well as the funds received for open access publication fees from the Forschungszentrum Jülich GmbH.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Gregor Michalicek, Jan Winkelmann, and Rudolf Zeller for fruitful discussions and help with conceptual and computational matters. This manuscript has been released as a pre-print at Hinzen et al. [21].

Appendix

Proposition 1. LetR>0,κ>0andfbe either the spherical Bessel function of the first kind,f=j, or the modified spherical Bessel function of the first kind,f=i. Then

0R(r2R2)nr+2f(κr)dr=(2)nn!R+n+2f+n+1(κR)κn+1

holds for alln0and all0.

Proof by mathematical induction on n. The base case n=0 follows immediately from the differentiation property Eq. 49 of the functions j and i, by

0Rr+2f(κr)dr=κ30κRr+2f(r)dr=κ1R+2f+1(κR).

The induction step uses the property for n1 and +1 and partial integration with u(r)=(r2R2)n and v(r)=r+2f(κr), to derive the statement for n and . Let

F(n,)=0R(r2R2)nr+2f(κr)dr

be the left-hand side of Eq. 1. Then

F(n,)=[(r2R2)nκ1r+2f+1(κr)]0R0Rn(r2R2)n12rκ1r+2f+1(κr)dr=2nκ1F(n1,+1)=(2)nn!R+n+2f+n+1(κR)κn+1,

where we used the induction hypothesis in the last equation, and thus the proposition follows. ∎Proposition 2. For nonnegative λ and R, and ν

iν(λR)λνRν(2ν+1)!!

holds, with equality if and only ifλ=0orR=0. Proof. Since iν(x)=ivjν(ix) and jν has the expansion [18]

jν(x)=xν(2ν+1)!!s=0(1)ss!(ν+32)s(x2)2s,

where ()s is the Pochhammer symbol defined for general a and s0 by

(a)0=1,(a)s=a(a+1)(a+s1),

iν can be expanded by

iν(x)=xν(2ν+1)!!s=01s!(ν+32)s(x2)2s.

iν(λR)/λν thus becomes

iν(λR)λν=Rν(2ν+1)!!s=01s!(ν+32)s(λR2)2s

and it remains to show that

s=01s!(ν+32)s(λR2)2s1.

The first term of the sum for s=0 is always 1 irrespective of the values of λ and R. The other summed terms for s>0 are positive, if both λ and R are positive, and zero, if one of the two variables is zero. ∎

References

1. Yukawa H. On the interaction of elementary particles. I. Proc Physico-Math Soc Jpn (1935) 33:48–57. doi:10.11429/ppmsj1919.17.0_48

CrossRef Full Text | Google Scholar

2. Hückel E, Debye P. The theory of electrolytes. I. lowering of freezing point and related phenomena. Phys Z (1923) 24:185–206.

Google Scholar

3. Thomas L. The calculation of atomic fields. Math Proc Cambridge Philos Soc (1927) 23:542–8. doi:10.1017/S0305004100011683

CrossRef Full Text | Google Scholar

4. Fermi E. Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente. Zeitschrift für Physik (1927) 48:73–9. doi:10.1007/BF01351576

CrossRef Full Text | Google Scholar

5. Hohenberg P, Kohn W. Inhomogeneous electron gas. Phys Rev (1964) 136:B864–71. doi:10.1103/PhysRev.136.B864

CrossRef Full Text | Google Scholar

6. Kohn W, Sham L. Density function theory. J Phys Rev (1965) 140:A1133–8.

Google Scholar

7. Tran F, Blaha P. Implementation of screened hybrid functionals based on the Yukawa potential within the LAPW basis set. Phys Rev B (2011) 83:235118. doi:10.1103/PhysRevB.83.235118

CrossRef Full Text | Google Scholar

8. Massidda S, Posternak M, Baldereschi A. Hartree-Fock lapw approach to the electronic properties of periodic systems. Phys Rev B (1993) 48:5058–68. doi:10.1103/PhysRevB.48.5058

CrossRef Full Text | Google Scholar

9. Kerker G. Efficient iteration scheme for self-consistent pseudopotential calculations. Phys Rev B (1981) 23:3082.

CrossRef Full Text | Google Scholar

10. Winkelmann M, Di Napoli E, Wortmann D, Blügel S. Kerker mixing scheme for self-consistent muffin-tin based all-electron electronic structure calculations. Phys Rev B (2020) 102:195138. doi:10.1103/PhysRevB.102.195138

CrossRef Full Text | Google Scholar

11. Kim J, Gulans A, Draxl C. Robust mixing in self-consistent linearized augmented planewave calculations. Electron Struct (2020) 2:037001. doi:10.1088/2516-1075/ababde

CrossRef Full Text | Google Scholar

12. Weinert M. Solution of Poisson’s equation: beyond Ewald-type methods. J Math Phys (1981) 22:2433–9. doi:10.1063/1.524800

CrossRef Full Text | Google Scholar

13. Eyert V. The plane-wave based full-potential ASW method. Berlin, Heidelberg, Germany: Springer (2013). p. 113–72. doi:10.1007/978-3-642-25864-0˙4

CrossRef Full Text

14. Drittler B. KKR-Greensche Funktionsmethode für das volle Zellpotential. [PhD thesis]. Aachen (Germany): RWTH Aachen (1991).

15. Wimmer E, Krakauer H, Weinert M, Freeman AJ. Full-potential self-consistent linearized-augmented-plane-wave method for calculating the electronic structure of molecules and surfaces: O2 molecule. Phys Rev B (1981) 24:864–75. doi:10.1103/PhysRevB.24.864

CrossRef Full Text | Google Scholar

16.FLEUR. The Jülich FLAPW code family. Available from: http://www.flapw.de/ (2020).

17. Abramowitz M, Stegun IA. “Handbook of mathematical functions with formulas, graphs, and mathematical tables,” in Ninth dover printing, tenth gpo printing. New York, NY: Dover (1964).

18. Arfken G, Weber H, Harris F. Mathematical methods for physicists: a comprehensive guide. Amsterdam, Netherlands: Elsevier Science (2013).

19. Jackson JD. Classical electrodynamics. 3rd ed. New York, NY: Wiley (1999).

20.Jülich Supercomputing Centre. JURECA: modular supercomputer at Jülich supercomputing centre. J Large-Scale Res Facil (2018) 4:33. doi:10.17815/jlsrf-4-121-1

CrossRef Full Text | Google Scholar

21. Hinzen M, Napoli ED, Wortmann D, Blügel S. Solution to the modified helmholtz equation for arbitrary periodic charge densities (2020).

Keywords: partial differential equations, density functional theory, electronic structure methods, Green functions technique, materials science, electrostatics, Fourier analysis, muffin-tin approximation

Citation: Winkelmann M, Di Napoli E, Wortmann D and Blügel S (2021) Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities. Front. Phys. 8:618142. doi: 10.3389/fphy.2020.618142

Received: 16 October 2020; Accepted: 16 December 2020;
Published: 11 March 2021.

Edited by:

Manuel Asorey, University of Zaragoza, Spain

Reviewed by:

Kazuharu Bamba, Fukushima University, Japan
Babak Shiri, Neijiang Normal University, China

Copyright © 2021 Winkelmann, Di Napoli, Wortmann and Blügel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Miriam Winkelmann, m.winkelmann@fz-juelich.de