Efficient and accurate solver of the three-dimensional screened and unscreened Poisson’s equation with generic boundary conditions

Efficient and accurate solver of the three-dimensional screened and unscreened Poisson’s equation with generic boundary conditions

Alessandro Cerioni alessandro.cerioni@esrf.fr European Synchrotron Radiation Facility, 6 rue Horowitz, BP 220, 38043 Grenoble Cedex 9, France    Luigi Genovese luigi.genovese@cea.fr Laboratoire de simulation atomistique (L_Sim), SP2M, UMR-E CEA / UJF-Grenoble 1, INAC, Grenoble, F-38054, France    Alessandro Mirone    Vicente Armando Sole European Synchrotron Radiation Facility, 6 rue Horowitz, BP 220, 38043 Grenoble Cedex 9, France

We present an explicit solver of the three-dimensional screened and unscreened Poisson’s equation which combines accuracy, computational efficiency and versatility. The solver, based on a mixed plane-wave / interpolating scaling function representation, can deal with any kind of periodicity (along one, two, or three spatial axes) as well as with fully isolated boundary conditions. It can seamlessly accommodate a finite screening length, non-orthorhombic lattices and charged systems. This approach is particularly advantageous because convergence is attained by simply refining the real space grid, namely without any adjustable parameter. At the same time, the numerical method features scaling of the computational cost ( being the number of grid points) very much like plane-wave methods. The methodology, validated on model systems, is tailored for leading-edge computer simulations of materials (including ab initio electronic structure computations), but it might as well be beneficial for other research domains.


I Introduction

Poisson’s equation (screened or not) is involved in a large variety of problems in physics and chemistry as well as in engineering. Therefore, there is a quite strong motivation for developing efficient and accurate solving methods.

As far as electrostatics is concerned, the three-dimensional screened Poisson’s equation is written as follows (in Gaussian units):


where represents a continuous electric charge distribution (the input of the problem at hand), is the electrostatic potential (the output), and represents the reciprocal screening length as defined, for instance, in the Debye-Hückel or Thomas-Fermi approximations. In the special case , Eq. (1) reduces to the usual Poisson’s equation.

Any method aiming at providing a solution to Eq. (1) has to deal with boundary conditions (BC), which in general can be either periodic or free (otherwise referred to as “isolated” or “open”’) along each of the three directions . In the case of fully periodic BC, the most natural (and efficient) approach to the problem is that of the reciprocal space treatment. It amounts to expanding both the density and the potential as superpositions of plane waves (Fourier series), following which Eq. (1) becomes algebraic in the Fourier components of and . This equation is readily solved and the result is finally transformed back into real space. Forward and backward transformations are carried out via Fast Fourier Transform (FFT), hence the overall computational scaling of the method with respect to the number of grid points is a rather appealing .

Owing to the above mentioned advantages, many attempts have been made to tackle also free BC or mixed free/periodic BC within a nominally fully periodic framework. In the simplest implementations, the simulation box is artificially enlarged by vacuum padding, such as to suppress the spurious Coulombian interaction among periodic replicas (“super-cell” approximation). Yet, as a result of the the long-range nature of the Coulomb interaction, there are situations in which the simulation box ought to be made unfeasibly large, in particular for charge distributions exhibiting significant multipolar terms. Moreover, a non-zero net charge in the primitive cell would yield a divergent total electrostatic energy when infinitely replicated along every direction, unless a compensating artificial uniformly charged background (“jellium”) is introduced, see e.g. Refs. Leslie and Gillan (1985) and Makov and Payne (1995). Another option which has been put forward consists in cutting off the long tail of the Coulomb interaction beyond a spherical region in real space, the radius of which is adequately chosen. Correspondingly, the reciprocal space components of the bare Coulomb potential are multiplied by screening functions, known analytically for all types of periodicity Jarvis et al. (1997); Rozzi et al. (2006). In a series of papers Martyna and Tuckerman (1999); Mináry et al. (2002, 2004), the screening function formalism was combined with the explicit break-up of the short- and long-range components of the Coulomb interaction, as also done in the context of (smooth particle-mesh) Ewald summation techniques, whereas in Refs. Dabo et al. (2008, 2011) the errors induced by the periodic images are alleviated by introducing a corrective potential.

Although we acknowledge that much remarkable work has been done on the subject, providing a thorough review would go beyond the scope of the present study. We therefore refer the reader to the original literature and proceed by presenting our approach, which differs from those mentioned above in that convergence is attained with no adjustable parameter and thus it aims at being fully generic.

This paper builds upon Refs. Genovese et al. (2007, 2006), where a novel method for solving the unscreened Poisson’s equation with free and surface-like BC was first presented. Such a method is direct (rather than iterative) in that the solution along the isolated directions is found in its integral form by using the Green’s function method. For instance, in the case of a fully isolated system (or “cluster-like”),


where . Homogeneous Dirichlet BC ( at ) along the isolated directions are explicitly enforced by the selection of the Green’s function.

The method has been in use for a few years in a number of ab-initio codes, namely ABINIT Gonze et al. (2005, 2009), BigDFT Genovese et al. (2008, 2011), CP2K VandeVondele et al. (2005) (see also Ref. Vácha et al. (2012) for a recent application thereof), Octopus Marques et al. (2003); Castro et al. (2006) and has proven to be highly efficient and accurate in every application attempted to date. It is based on a mixed plane-wave / interpolating scaling function (ISF) representation of and which allows to model any sort of periodicity in the most natural, clean and mathematically rigorous way. Clearly, periodic (isolated) directions are represented in terms of plane waves (interpolating scaling functions). ISFs - arising in the wavelet theory Daubechies (1992); Goedecker (1998) - enjoy several properties which make them superior to other basis sets. For instance, the representation in terms of -th order ISFs make the first moments of the continuous and discrete charge distributions coincide Genovese et al. (2006). As a consequence the representation is definitely faithful (other than handy), since the different moments of the charge distribution capture the major features of the potential. Moreover, ISFs are genuinely localized due to their compact support (the length of which is equal to ) and endowed with the so-called “refinement relations” which easily allow to switch from a representation on a grid with spacing to a doubly refined grid with spacing .

We have extended the previous implementation to account also for screening, for the case of periodicity along only one direction (“wire-like BC”) and for non-orthorhombic cells, this investigation providing a detailed account of such improvements. The inclusion of such new functionalities is motivated by the strong theoretical, experimental and technological interest in the characterization of nanostructured materials (among which polar nanorods, see e.g. Ref. Hine et al. (2012) and references therein), since solving Poisson’s equation is only one of the many steps involved in state-of-the-art computer simulations and is repeated several times. Moreover, in the context of Kohn-Sham (KS) density functional theory (DFT) and extensions thereof, there are quantities which are computed via convolution integrals very similar to that in Eq. (2): for instance, the exact exchange term arising within those generalizations of KS-DFT employing orbital-dependent or hybrid functionals (see Kronik et al. (2012) and references therein), or the coupling-matrix in time-dependent DFT Natarajan et al. (2012). In this respect, the electrostatic problem of concern here provides the paradigm for many other computations, even well beyond the scope of electrostatics.

Regarding the possibility of accounting also for screening, we note that this novel feature might for example be used to solve the Schrödinger equation iteratively (see e.g. Ref. Harrison et al. (2004)). In fact, one can exploit the formal analogy between Eq. (1) and the Schrödinger equation which becomes apparent if the latter is written in the following fashion:


where represents a bound eigenstate with negative energy ().

The next sections are structured as follows: we first present our solution method for free, wire-like and surface-like BC. We then discuss the accuracy of the proposed solver by reporting a collection of numerical benchmarks. We conclude by highlighting the benefits of using the present methodology.

Ii Free Boundary Conditions

In the case of free BC, the Green’s function which has to be plugged into Eq. (2) is




where . Along the same lines as Ref. Genovese et al. (2006), both the charge distribution and the electrostatic potential are expanded in terms of ISFs, here denoted by :


where and represent the (uniform) grid spacing and the number of grid points along each direction, respectively. Since the function is, by construction, such that , , the expansion coefficients are readily found to be:


In other words, the expansion coefficients coincide exactly with the values of on a uniform grid. In this respect, the ISF representation appears genuinely tailored for numerical studies, where the whole available information reduces to knowing the values at the grid points. Clearly, the grid spacing has to be chosen adequately, depending on the typical spatial scales over which the density distribution exhibits significant variations. The underlying mathematics then assures that the moments built upon the discrete charge distribution coincide with those of the continuous charge distribution up to order , where is the order of the ISF (see Ref. Genovese et al. (2006) for the proof):


if . A representation analogous to that in Eq. (6) can also be given for the potential , where


replaces . As a consequence of the chosen representation, the convolution integral in Eq. (2) is more conveniently expressed in Cartesian coordinates, and upon plugging Eq. (6) into Eq. (2) we obtain the following equation in discrete form:




is the convolution kernel.

The numerical evaluation of Eq. (11) would be too onerous if performed directly. The computational effort can be drastically reduced by expressing the Green’s function as a linear combination of Gaussian functions, as the integral would become separable along and the resulting 1D integrals can be evaluated very efficiently.

We proceed by approximating the Green’s function as


where and are determined so as to minimize the error on a given range of . More specifically, we found out the following approximation,


with satisfactory accuracy for any (see Fig. 1). The best fit was performed using the Levenberg-Marquardt algorithm Levenberg (1944); Marquardt (1963) with 136 Gaussian functions, the -values ranging between and . The actual Green’s function can therefore be written as in Eq. (13) with


In the unscreened case () we use the Gaussian fit of the function which was already proposed in Ref. Genovese et al. (2006) and that we recall as being affected by an error for any .

Figure 1: Accuracy of the approximation of the function with 136 Gaussians used in the solution of the screened Poisson’s equation for the case of free BC. The range of the independent variable is . We plot both the absolute and the relative error because the latter is a better indicator close to the origin (where the fitted function takes on very large values), while the former is a reliable signature of the goodness of the fit towards the opposite end. Note that at the function is already smaller than the machine precision.

The convolution kernel can be written as follows:




and can be computed by evaluating 1D integrals - being the number of Gaussian functions - hence at a much lower cost than Eq. (11), which would require the computation of 3D integrals, instead.

We point out that the numerical evaluation of Eq. (16) is performed using the same method described in Ref. Genovese et al. (2006), which exploits the refinement relations fulfilled by the ISFs and yields an accuracy as high as the machine precision, even for the narrowest Gaussians. To this effect, ISFs show their superiority over other basis sets (cf. e.g. the explicit method laid out in Ref. Lee and Tuckerman (2008), where a Gaussian approximation similar to ours is carried out within a discrete variable representation approach).

Iii Wire-like Boundary Conditions

We now consider a system which is periodic along the direction (with period equal to ) and isolated over the -plane. We can hence expand the continuous charge density distribution as a sum over its Fourier components along :


noting that in this case the Fourier coefficients are solely functions of and . After expanding the electrostatic potential in a similar manner, the screened Poisson’s equation yields the following relation between the potential’s reciprocal space components and those of the density:


where . The symmetry of the problem suggests writing the Green’s function of Eq. (18) in cylindrical coordinates:


where and . The solution of Eq. (19) is given by


where is the zero-th order modified Bessel function of the second kind. We then express the 2D Fourier components of both density and potential in terms of ISFs, thereby completing the required steps towards the mixed plane-wave / ISF representation for the case investigated here:


where and are the grid spacings along the non-periodic directions. Combining Eq. (18) with Eq. (17) and Eq. (21) one obtains

where the kernel is very similar to that in Eq. (11), except that in (III) the integral is restricted to the non-periodic directions. Within the approximation (12) the kernel elements can be evaluated as in Eq. (15) with .

Clearly, the accuracy of the whole method depends on the accuracy of the Gaussian approximation of the Green’s function. Note that close to the origin the function behaves like , whereas it decreases exponentially for large arguments:


We were able to approximate as a sum of Gaussians (we denote the fitting parameters by , cf. Eq. (13)) with an accuracy better than in the range , where the upper bound for was chosen by realizing that the value attained by at is , i.e. already comparable with the machine precision. The function fitting was carried out using the same algorithm as before, but with 144 Gaussians, with -values ranging between and . We used the very same 144 Gaussians () but different relative weights () also to approximate , achieving an accuracy of in the range (see Fig. 2).

Figure 2: Accuracy of the approximation of the Green’s function with 144 Gaussians as used in the solution of the screened Poisson’s equation for the case of wire-like BC. The range of the independent variable is for the function and for .

Therefore, for any the convolution kernel is the same as in Eq. (15) with , and , whereas for we exploit the scaling properties of the logarithm to adapt the best fit obtained for to any . We write it explicitly for sake of clarity:

where is deliberately chosen so that . In fact, is the box size along , i.e. the maximum range covered by the charge distribution along each axis; is the extent of the ISF definition domain. The convolution of any charge distribution function with an ISF can be non-zero at most over a rectangular domain of diagonal length .

In order to find the solution for the electrostatic potential in real space, we first compute the Fourier coefficients of the density () through a 1D FFT along the periodic direction . The corresponding quantities for the potential are then obtained by calculating the convolutions in Eq. (III) via a zero-padded FFT procedure Goedecker (1993). Finally, the potential is transformed back into real space along the direction. Let us notice that real-to-complex FFTs can be used instead of complex-to-complex FFTs, since all the quantities are real and the kernel is symmetric.

Iv Surface-like Boundary Conditions

While referring the reader to Ref. Genovese et al. (2007) for a more extensive description of our treatment of surface-like BC, in the following we just discuss our methodological improvements . In particular, we are now able to allow also for the screening and for monoclinic lattices, simply by redefining the in the equation relating the 2D Fourier components of density and potential,


as follows:


where (i.e. the periodic directions with period equal to , ) and is the 2D contra-variant metric tensor,


being the angle between the and unit vectors.

We point out that in the case of surface-like BC there is no Gaussian approximation of the Green’s function involved in the procedure. Consequently, in this case the accuracy of the method is limited only by the machine precision.

V Numerical results

In order to measure the accuracy of our method, we used several test charge distributions for which the Poisson’s equation is exactly solvable, and compared the approximate numerical solution against the exact one. In the following the accuracy is given in terms of the infinity norm:


Without loss of generality, all tests were run on a cubic simulation box (), with a uniform grid and an equal number of points along each direction.

In particular, in the case of free BC we used a Gaussian density distribution, . As our procedure relies on a convolution - see Eq. (2) - from which no divergence can arise (as long as is regular), we expect that the numerical solution is regular everywhere. The latter statement offers an unambiguous prescription for fixing the integration constants of the analytic exact solution. The latter is eventually found to be as follows:


The results of our tests - with such that - are reported in Fig. 3 for all the ISF supported in our code but no screening and in Fig. 4 for 16th order ISF and selected values of covering four orders of magnitude. We point out that, owing to the chosen value of , the value of the density on the simulation box faces is , hence smaller than our solver’s accuracy. In order to study the influence of the box size on accuracy, we performed several runs with different box sizes, computed the Hartree energy, and compared the result to the exact Hartree energy corresponding to an infinite box size and unitary monopole ():


Throughout this second set of tests the value of was chosen so as to keep . The results are reported in Fig. 5, together with the value of the charge density distribution at the box faces. We can observe that there is absolutely no need to enlarge the box size beyond the reference size, as the latter is already large enough to capture all the features within the reach of our solver. On the other hand, accuracy decreases on reducing the box size, as the ideal free BC becomes increasingly violated. Nevertheless, the attained accuracy remains interesting over a broad range of box sizes even smaller than the reference one, at variance with plane-wave methods which would fail our test, firstly because of the presence of a non-zero monopole and secondly because of the aliasing due to insufficiently large box.

We remark that, as opposed to what claimed in Ref. Hine et al. (2011) in relation to Refs. Genovese et al. (2007, 2006) (where the method deployed here was first proposed), our solver is actually reliable without any resort to the so-called “minimum image convention”, namely without rendering the input density charge distribution periodic within the simulation box.

Figure 3: Accuracy test for the case of free/isolated boundary conditions in the absence of screening ( is the order of the ISF, the grid spacing).
Figure 4: Accuracy test for the case of free/isolated boundary conditions in the presence of screening ( is the order of the ISF, the grid spacing).
Figure 5: Influence of the (cubic) simulation box size on the accuracy of the Hartree energy. stands for the box size, whereas stands for the box size for which .

In the case of surface-like and wire-like BC, we consider a charge density distribution obtained by applying explicitly the screened Poisson’s operator on an exact potential written as


where each of the ad hoc functions ’s entering the product is either periodic,


or localized,


depending on the intended BC. Results are shown in Figs. 6, 7, 8, 9 and indicate a very good overall convergence rate.

Figure 6: Accuracy test for the case of surface-like boundary conditions in the absence of screening ( is the order of the ISF, the grid spacing).
Figure 7: Accuracy test for the case of surface-like boundary conditions in the presence of screening ( is the order of the ISF, the grid spacing).
Figure 8: Accuracy test for the case of wire-like boundary conditions in the absence of screening ( is the order of the ISF, the grid spacing).
Figure 9: Accuracy test for the case of wire-like boundary conditions in the presence of screening ( is the order of the ISF, the grid spacing).

In the case of wire-like BC, we furthered our tests by choosing a 2D Gaussian charge distribution,


where the charge density along is implicitly set to unity. The corresponding exact potential to be used as reference is


where is the exponential integral function. On deriving Eq. (35), integration constants were fixed unambiguously by following the same criteria which led to Eq. (29). In particular, one integration constant guarantees the regularity of the potential at the origin,


where is the Euler-Mascheroni constant, while the second integration constant (an additive constant) is set to zero so that as , in accordance with the behavior of the Green’s function. This test case is relevant to our studies for a two-fold reason. Firstly, it is typically out-of-reach for plane-wave methods, since the density distribution in Eq. (34) exhibits a non-zero monopole. Secondly, it allows to probe the goodness of the Gaussian approximation of the function involved in Eq. (20) within the approximation (12). The charge distribution being constant along , the only non-trivial term in the Fourier expansion along is the zero-mode () and, in case the screening is absent (), only the -branch of the Green’s function (20) plays a role in the computation. We have also analyzed the rate of convergence towards the exact Hartree linear energy density,


so as to have a further confirmation that the great accuracy in the numerically evaluated electrostatic potential ensures the reliable computation of derived physical quantities. The results are shown in Figs. 10-11 and are definitely good. We also include in Fig. 12 a plot of the charge density distribution (magnified by a factor 10) and the corresponding potentials obtained at different ’s. We observe that for the potential does not fall to zero for increasing , whereas it tends to be more and more localized around the origin as increases, as expected.

Figure 10: Accuracy test for the case of wire-like boundary conditions (WBC) with monopolar charge density distribution ( is the order of the ISF, the grid spacing).
Figure 11: Accuracy in the computation of the Hartree linear energy density - see Eq. (37) - in the case of wire-like boundary conditions (WBC) with monopolar charge density distribution ( is the order of the ISF, the grid spacing). In our setup .
Figure 12: Gaussian density charge distribution (red, dashed mesh) as a function of the isolated directions (, in our notation) and the corresponding electrostatic potential (black, solid mesh) evaluated for different values of the screening, namely upon increasing boundary thickness. The charge density is implicitly periodic along (wire-like BC). The amplitude of the plotted density charge distribution is multiplied by a factor with respect to the actual value in order to improve the readability of the picture. Only half of the solution is drawn to highlight its profile.

Having checked that our solver performs quite well in all the above mentioned cases, it seemed tempting to further probe its capabilities with some other charge density distributions. In particular, we modelled a planar capacitor (hence involving surface-like BC), and a cylindrical capacitor (wire-like BC). In the latter case, the input charge distribution was mimicked by a positive 2D Gaussian distribution sharply peaked at the origin, together with a ring of negative Gaussian distributions centered at . The relative amplitudes of the central and of the peripheral Gaussians were chosen to yield zero total charge. The solutions obtained are shown in Figs. 13-14, are definitely consistent with the intuitively expected behavior.

Figure 13: Electrostatic potential generated by a pair of planar charge distributions of opposite sign (positive in red/solid; negative in black/dashed) modelling a planar capacitor unlimited in the periodic directions. The piecewise planar behavior corresponds to the case with no screening (), whereas the other solutions are obtained with upon increasing the boundary thickness. The potential is more and more localized around the capacitor’s plates and falls rapidly to zero as is increased. Each curve is normalised to one to improve readability.
Figure 14: Electrostatic potential generated by a cylindrical capacitor, periodic in the vertical direction. The different solutions correspond to upon increasing the boundary thickness. Each curve is normalised to one for sake of readability. Only half of the solution is drawn so as to highlight the potential profile.

Vi Conclusion

We have presented a numerical method for the solution of Poisson’s equation which can tackle any type of periodicity, the presence of screening, non-orthorhombic geometries and charged systems, showing that convergence to highly accurate results is attained with no adjustable parameter other than the (unavoidable) grid spacing.

The charge density distribution and the electrostatic potential are both expressed in terms of plane waves along the periodic directions, and in terms of interpolating scaling functions along the non-periodic directions. The latter representation proves to be very handy, in that the expansion coefficients of any continuous quantity are simply given by its values on a uniform grid. Moreover, -th order interpolating scaling functions preserve the matching between the moments built upon the discretized quantity with those of the continuous one. This occurrence is particularly important in electrostatics, because of the interest in resolving the multipolar features of the electrostatic potential and other derived quantities.

In our approach, the solution is obtained via the Green’s function method. The (in principle) most expensive operation, namely the convolution of the Green’s function with the input charge density distribution, becomes affordable by making use of highly-optimized FFT routines (zero-padded along the isolated directions), while the evaluation of the convolution kernel is carried out separately along the three spatial directions by approximating the Green’s function as a sum of Gaussian functions.

Owing to the above mentioned advantages, our solver is suitable for intensive computer simulations of electronic structure and molecular dynamics, where the Poisson’s equation has to be solved several times and it is important to limit the growth and propagation of numerical errors as much as possible, especially because other physical quantities are computed starting from the solution of the Poisson’s equation. There are other contexts in computational physics and chemistry in which relevant quantities are obtainable as convolutions involving the same Green’s functions found in electrostatics. This is the case, for instance, of the exact exchange term within the generalizations of Kohn-Sham DFT employing orbital-dependent or hybrid functionals. Actually, our methodology features a level of generality which allows to address also problems well beyond electrostatics.

As a possible outlook, we are working towards enabling the computation of range-separated hybrid functionals, where the Coulomb potential is split into a long- and a short-range component. In the basic range-separated approach the Coulomb interaction is written as , being an adjustable length scale, although other options have been put forward (in which, for instance, the long- and short-range parts can be weighted differently). We are thus planning to model the range-separated Coulomb interaction within our framework.

Our solver is already designed for taking full advantage of multi-core CPUs and GPUs, and is currently integrated in BigDFT, the sources of which are freely downloadable from http://inac.cea.fr/L_Sim/BigDFT/. The release of a stand-alone package is also envisaged for the near future. The details on the GPU acceleration and the performance of the solver in the context of massively parallel electronic structure computations will be described in a forthcoming paper.

The authors thank Thierry Deutsch for valuable suggestions on the manuscript, and Claudio Ferrero for the critical proofreading. A.C. acknowledges the financial support of the French National Research Agency in the frame of the “NEWCASTLE” project.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description