Solution to the modified Helmholtz equation for arbitrary periodic charge densities

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 [M. Weinert, J. Math. Phys. 22, 2433 (1981)] 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.


I. INTRODUCTION
A variety of problems in condensed matter physics require an efficient solution of the partial differential equation for a charge density ρ in a periodic domain. 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 (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 apa) Corresponding author: m.hinzen@fz-juelich.de pear 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 equation (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 (1) is replaced by an overlap charge density 8 obtained as a product of wave functions associated with different quantum numbers.
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 (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 9 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 fullpotential all-electron DFT methods, such as the augmented spherical wave (ASW) method 10 , the Korringa-Kohn-Rostoker Green function (KKR-GF) method 11 , and the full-potential linearized augmented planewave (FLAPW) method 12 , 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 equation (1) for values of λ > 0, and for general periodic charge densities without shape approximation. We formulate our 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 9 , we demonstrate this extension explicitly for the FLAPW method 12 as implemented in the FLEUR code 13 . We provide a complete derivation of the modified multipole expansion using a Green function method, and the derivation of 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. We also discuss the convergence of the Fourier series of the pseudo-charge density. We point out which are the algorithmic changes required to extend the solution of the Poisson equation to a solution of (1), which can then be straightforwardly transferred to other all-electron full-potential band structure methods. This paper is organized as follows: In Sect. II A, 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 pseudocharge 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 Sect II B 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 Sect. II B 1, the radial free-space Green function is used to define the modified multipole expansion of the Yukawa potential. In Sect. II C, 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 Sect. II B 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 II D together with the minimal modifications necessary to change Weinert's original algorithm. The conclusions and the outlook are presented in Sect. III.

A. 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 B Rα (τ α ) 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 and potentials are represented in plane waves e iK·r , where K defines the reciprocal lattice vector dual to the lattice vectors defining the periodic domain, and in spherical harmonics, Y L , 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, withr α = r−τα |r−τα| its unit vector. The precision of the representation is determined by the cut-off parameters K max for the wave vectors, K, with length K ≤ K max , 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 B Rα (τ α ) can generate the same multipole moments and thus, the same potential outside the sphere. Here, Y * L denotes the complex conjugate of Y L . The pseudo-charge density,ρ, defined by Weinert in Ref. 9 is such a charge density. It fulfills the following three conditions: • It has the same multipole moments q α • It is equal to the true charge density ρ I in the interstitial region.
• It has a fast convergent Fourier expansion.
The Fourier components of V I are then simply while V I (0) will be set to a constant. Once the interstitial potential V I has been calculated, the muffin-tin potential can be obtained by solving the Dirichlet boundary value problem on the sphere 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 α + τ α ∈ B Rα (τ α ). 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 V α S due to the charge density distribution inside the sphere and a boundary term V α B due to the interstitial potential at the boundary of the sphere. The Fourier coefficients of the pseudo-charge density basically have the form 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 can also be evaluated for r ∈ ∪ α B Rα (τ α ), then the true charge ρ can also be written as whereρ are charge densities localized in the atomic spheres B Rα (τ α ). If these localized densities are now substituted by other localized densitiesρ α , then the charge density is still correct in I. We now continue with the derivation of the muffin-tin and interstitial potentials for λ > 0.

B. 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 B Rα (τ α ) with centre τ α through the solution of the boundary value problem 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 V α L (r α ), the expansion coefficients to the spherical harmonics expansion inside the sphere (see (3)).
Step 1. We solve the homogeneous modified Helmholtz equation, (∆ − λ 2 )U = 0, in spherical coordinates. Following the solution 14 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 Y L (r), U (r,r) = L u (r)Y L (r). The term −λ 2 U in the homogeneous modified Helmholtz equation only has an effect on the radial solution. The spherical harmonics, Y L , 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 14,15 , (14) and its fundamental set of solutions are for each the two modified spherical Bessel functions 14,15 i (λr) and k (λr), the first of which is the regular solution welldefined at the origin, but grows fast with growing radius r and the second is the irregular solution that goes to infinity for r → 0. To realize the proper boundary condition for the radial Green functions two conditions have to be fulfilled: The first solution, u 1 , must be finite at r = 0. We conclude that since k (λr) → ∞ for r → 0. The second solution, u 2 , must be 0 at r = R α . This is achieved by a linear combination of the two modified spherical Bessel functions, Step 2.
subject to the Dirichlet boundary condition where (∆ r − λ 2 ) is the linear radial differential operator in (14) and δ denotes the radial Dirac delta function, for which The radial Green function takes the form of the product of the two linearly independent solutions with the proper boundary conditions, where r < = min(r, r ), r > = max(r, r ) and Since the Wronskian W is linear, the addition of cu 1 (c =const) onto k (λr) to suffice the boundary condition u 2 (R α ) = 0 has no influence on the Wronskian: W (u 1 , cu 1 ) = 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 r → 0 or the asymptotic values for r → ∞ to find Therefore, C = 4πλ, and the radial Green function finally reads Step 3. Considering the standard expression of Dirac's delta function separated according to the radial and angular coordinates the three-dimensional (3D) Green function G(r, r ) solving is expanded in the form and the solution to the inhomogeneous equation (12) is given by (7). For the derivation of the 3D Green function and the 3D inhomogeneous solution V α = V α S + V α B (7) we refer the reader to Ref. 16. Both the integral over the sphere B Rα (0) and the boundary integral simplify by exploiting the orthonormality relation of the spherical harmonics, ∂B1(0) Y * L (r)Y L (r) dω = δ LL . The integral over the sphere B Rα (0) provides the source contribution to the muffin-tin potential In order to obtain the boundary contribution to the muffin-tin potential, we evaluate the boundary integral in (7) by expanding the interstitial potential V I (K) (see Sect. II C) on the sphere boundaries ∂B Rα (τ α ) r in spherical coordinates using the plane-wave expansion where (31) Furthermore, the normal derivative of G on the sphere boundary is Since r α < R α = r α and since G takes the form (20), we obtain We recall that u 2 (R α ) = 0 and reuse (22) to obtain yielding ∂G (r α , r α ) ∂r α With this and the knowledge of the interstitial potential at the sphere boundary from (31), the boundary contribution to the muffin-tin potential becomes and the radial part of the spherical harmonics expansion of the total potential, V α S + V α B , in the sphere B Rα (0) becomes 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 α :

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 , The modified spherical Bessel function k already contains the proper boundary condition for r → ∞. A charge densityρ α localized in a sphere B Rα (τ α ) embedded in free space, produces a Yukawa potential outside the sphere, i.e. r = r α + τ α / ∈ B Rα (τ α ), which can be expressed analogously to the Coulomb potential in terms of the modified multipole expansion with the modified multipole moments An analogous definition holds true for the modified multiple moments q α L [ρ α ] of the true charge ρ α in the sphere. With the standard expansion of the charge density inside the sphere into spherical harmonics (3), ρ α (r α + τ α ) = L ρ α L (r α ) Y L (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 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 (40), since for the free-space Green function (39) the boundary value term disappears for r → ∞.

C. Interstitial Yukawa Potential
Suppose we had found a Fourier transformable pseudocharge 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 (8) and the solution of the modified Helmholtz equation (1) through Fourier transformation yields an algebraic equation from which we calculate the interstitial Yukawa potential, In the previous subsection II B 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, , 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, . From (41) we conclude that this is fulfilled if the modified multipole moments (42) of both charge densities, q α L [ρ α ] = q α L [ρ α ] − q α L [ρ I ] and q α L [ρ α ], are equal. The modified multiple moments q α L [ρ α ] are already known through (43). Next we determine the modified multiple moments q α L [ρ I ] of the interstitial charge extended into the muffin-tin spheres and then construct the pseudo-charge density.

Modified Multiple 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. 9, but since the modified multipole moments are different from the known multipole moments for the Coulomb potential, we go through this step of deriving q α L [ρ I ] in detail. We write ρ I relative to the sphere centre τ α and employ the Rayleigh expansion (30) to e iK·rα , which yields The modified multipole moments of ρ I in the sphere B Rα (τ α ), defined analogously to (42), are For K = 0, the latter integral becomes If K = 0, observe that j (0) = δ 0 (Kronecker δ) and so the integral becomes Both equations above can be derived by partial integration in two different ways and applying the identities in Ref. 15, for f = i with the plus sign and for f = j with the minus sign, and for both f = i and f = j . In conclusion, this yields

Construction of Pseudo-Charge Density
We construct the pseudo-charge density by following the Ansatz of Weinert, 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 Sect. II C 3, it is beneficial to choose a η = (−1) n−η R 2(n−η) α n η a n for η = 0, . . . , n (53) and ν η = + 2η, where n is yet to be determined. As will become apparent later when we derive the coefficients Q α Ln in (62), the coefficient a n cancels out in any relevant equation, like (55) or (56). With these choices of parameters and with the binomials theorem applied to it follows from the ansatz (52) The Fourier transform of this expression is then given bȳ |Ω| is the volume of the periodic domain and we used the Rayleigh expansion (30) and the orthonormality relation of the spherical harmonics. With Prop. A.1 in Appendix A A α n (K) finally reduces to A α n (K) = a n (−2) n n!R +n+2 In the same way we derive the coefficients Q α Ln : We insert (55) in the definition of the modified multipole mo-ments (42) and use (A.1) for f = i and κ = λ to obtain and thus, λ +n+1 a n (−2) n n!R +n+2 .
(62) Since A α n (K) and Q α Ln share the term a n (−2) n n!R +n+2 α , the term cancels out in the product Q α Ln A α n (K) entering (56), and setting ν = + n + 1 (in accordance with Sect. II C 3) this leads tō Forρ α (0) we take the limit which is only different from 0 for = 0 and thus yields 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 local- calculated from the modified multipole moments q α L [ρ α ] (43) and q α L [ρ I ] (51) of ρ α and ρ I respectively, in the sphere B Rα (τ α ).

Smoothness of the Pseudo-Charge Density and Convergence of Its Fourier Series
In Sect. II C 2, we have set a η by (53), ν η = + 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 (52) and derived the much simpler form (55). The function itself and all its first n − 1 derivatives with respect to r α are equal to zero at r α = R α . Consequently, we ensure smoothness on the boundary of the sphere, Note that this includes localization ofρ α in B Rα (τ α ) and thus the pseudo-charge density equals the charge density in I (condition 2 in Sect. II A). 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 n−1-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 1 K n−1 , 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 9 discussed the factors Q α Ln A α n (K)/q α L [ρ α ], where in the Coulomb case (69) His reasoning is based on the zeros of the K-dependent function Q α Ln A α n (K)/q α L [ρ α ]. Viewed as a function of K, however, our factor Q α Ln A α n (K)/q α L [ρ α ] differs from Weinert's one only by a multiplicative constant. Thus Weinert's arguments apply here as well. With Prop. A.2 in Appendix A it follows that this multiplicative constant is smaller than 1 for λ > 0. We see this confirmed in Fig. 1, which shows the KR α -dependence of the factor Q α Ln A α n (K)/q α L [ρ α ] in the Yukawa (63) and Coulomb (69) cases for several combinations of n and angular degree -it reveals a smaller amplitude in the Yukawa case. In agreement with Ref. 9, in Fig. 1 we make two observations: (i) Q α Ln A α n /q α L [ρ α ] 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 K max , the factor Q α Ln A α n /q α L [ρ α ] must be small for K > K max , 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 Fig. 1. From this arises Weinert's criterion for choosing n( ), which we adopt here: • Choose ν ∈ N such that the first zero of j ν (z) is approximately equal to (KR α ) max .
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 (64) depending on ν need to be computed just once.

D. Algorithm: Construction of Yukawa Potential
Algorithm 1 summarizes the construction of the Yukawa potential derived in this paper.
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, q α L [ρ I ] (51) and q α L [ρ α ] (43) respectively, and in the Fourier components of the pseudo-charge densityρ α (K) (64) and (66). 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 (44) undergoes changes indirectly through the pseudo-charge density and directly by the prefactor 4π K 2 +λ 2 that substitutes 4π K 2 . Since the K = 0-term is well-defined, it is not

Algorithm 1 Bulk-Case Yukawa Potential
Input: charge density ρ, integer ν chosen as described on page 8 and preconditioning parameter λ. Output: Yukawa potential V λ solving the modified Helmholtz equation (1)  set to a constant as in the original method. The muffintin potential is affected only in its radial dependence in both the boundary and the source contribution of (38). 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, V I L (R α ; τ α ) (31), only changes indirectly through the changed values of V I λ (K).

III. 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 9 . 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 allelectron 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.

IV. ACKNOWLEDGMENTS
We would like to thank Gregor Michalicek, Jan Winkelmann and Rudolf Zeller for fruitful discussions and help with conceptual and computational matters. This work has been supported by a JARA-HPC seedfund project and by the MaX Center of Excellence funded by the EU through the H2020-EINFRA-2015-1 project: GA 676598. 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 17 at Forschungszentrum Jülich.  (r 2 −R 2 ) n r +2 f (κr) dr = (−2) n n!R +n+2 f +n+1 (κR) κ n+1 holds for all n ∈ N 0 and all ∈ N 0 .
Proof by mathematical induction on n. The base case n = 0 follows immediately from the differentiation property (50) of the functions j and i , by The induction step uses the property for n − 1 and + 1, and partial integration with u(r) = (r 2 − R 2 ) n and v (r) = r +2 f (κr), to derive the statement for n and . Let Proof. Since i ν (x) = i −ν j ν (ix) and j ν has the expansion 15 where (·) s is the Pochhammer symbol defined for general a ∈ R and s ∈ N 0 by (a) 0 = 1 , (a) s = a(a + 1) · · · (a + s − 1) , i ν can be expanded by i ν (x) = x ν (2ν + 1)!!  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.

DATA AVAILIBILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.