ORIGINAL RESEARCH article
Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities
- 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.
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. 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 ,
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  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  (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 . 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  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  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  is applied to density functional methods with Fourier transformable charge density. By solving Eq. 1 for general periodic densities, the Kerker preconditioner  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,
In a seminal work, Weinert  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 , the Korringa–Kohn–Rostoker Green function (KKR-GF) method , and the full-potential linearized augmented planewave (FLAPW) method , 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
As a matter of choice, and motivated by the original work of Weinert , we demonstrate this extension explicitly for the FLAPW method  as implemented in the Fleur code . 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.  relevant in an actual computation. Beyond Tran et al.  we provide detailed derivations and justify the use of identities such as Proposition 1 for the modified spherical Bessel function
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
are represented in plane waves
Weinert’s pseudo-charge method for the Poisson equation is based on the crucial observation that several charge densities ρ inside a sphere
and thus, the same potential
outside the sphere. Here,
• It has the same multipole moments
• It is equal to the true charge density
• It has a fast convergent Fourier expansion.
The Fourier Components of
where G is a Green function associated with the solution of the Poisson equation,
can also be evaluated for
are charge densities localized in the atomic spheres
The approach above is generally the same for the Yukawa potential, i.e., for
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
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
Step 1. We solve the homogeneous modified Helmholtz equation,
and its fundamental set of solutions are for each
Step 2. A function
subject to the Dirichlet boundary condition
The radial Green function takes the form of the product of the two linearly independent solutions with the proper boundary conditions,
Since the Wronskian W is linear, the addition of
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
is expanded in the form
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
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
using the plane-wave expansion
Furthermore, the normal derivative of G on the sphere boundary is
We recall that
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
and the radial part of the spherical harmonics expansion of the total potential,
Due to the kink of
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
The modified spherical Bessel function
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
The summation of
2.3 Interstitial Yukawa Potential
Suppose we had found a Fourier transformable pseudo-charge density
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
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. , but since the modified multipole moments are different from the known multipole moments for the Coulomb potential, we go through this step of deriving
The modified multipole moments of
Both equations above can be derived by partial integration in two different ways and applying the identities in Ref. ,
2.3.2 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
it follows from the ansatz Eq. 51
The Fourier transform of this expression is then given by
which is only different from 0 for
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
2.3.3 Smoothness of the Pseudo-Charge Density and Convergence of Its Fourier Series
In Section 2.3.2, we have set
itself and all its first
Note that this includes localization of
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
For an explicit rule on how to choose n, Weinert  discussed the factors
His reasoning is based on the zeros of the K-dependent function
Note that in this method
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
1: Modified multipole moments
2: Modified multipole moments
3: Modified Multipole Moments
5: Pseudo-charge density
6: Interstitial potential
7: Boundary terms
8: Radial parts
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,
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 . 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
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.
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.
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.
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. .
Proposition 1. Let
holds for all
Proof by mathematical induction on n. The base case
The induction step uses the property for
be the left-hand side of Eq. 1. Then
where we used the induction hypothesis in the last equation, and thus the proposition follows. ∎Proposition 2. For nonnegative λ and R, and
holds, with equality if and only if
and it remains to show that
The first term of the sum for
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
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
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
16.FLEUR. The Jülich FLAPW code family. Available from: http://www.flapw.de/ (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, firstname.lastname@example.org