A numerical Poisson solver with improved radial solutions for a self-consistent locally scaled self-interaction correction method

The universal applicability of density functional approximations is limited by self-interaction error made by these functionals. Recently, a novel one-electron self-interaction-correction (SIC) method that uses an iso-orbital indicator to apply the SIC at each point in space by scaling the exchange-correlation and Coulomb energy densities was proposed. The locally scaled SIC (LSIC) method is exact for the one-electron densities, and unlike the well-known Perdew–Zunger SIC (PZSIC) method recovers the uniform electron gas limit of the uncorrected density functional approximation, and reduces to PZSIC method as a special case when isoorbital indicator is set to the unity. Here, we present a numerical scheme that we have adopted to evaluate the Coulomb potential of the electron density scaled by the iso-orbital indicator required for the self-consistent LSIC calculations. After analyzing the behavior of the finite difference method (FDM) and the green function solution to the radial part of the Poisson equation, we adopt a hybrid approach that uses the FDM for the Coulomb potential due to the monopole and the GF for all higher-order terms. The performance of the resultant hybrid method is assessed using a variety of systems. The results show improved accuracy than earlier numerical schemes. We also find that, even with a generic set of radial grid parameters, accurate energy differences can be obtained using a numerical Coulomb solver in standard density functional studies.


Introduction
The Kohn-Sham formulation of density functional theory (DFT)is the workhorse method for material simulations when a quantum mechanical description is required [1].Among the several factors, such as efficient numerical implementations in easy-to-use codes and the rapid advancement of computing power, one important factor that contributed to this explosive growth in the use of DFT is the development of sufficiently accurate density functional approximations (DFAs) to the exact exchange-correlation functional [2].The majority of these approximations, however, suffer from the self-interaction error (SIE).Fortunately, SIE gets cancelled in the computation of many equilibrium properties but can be problematic in dissociation limits, systems with fractional charges, or when the properties are sensitive to the asymptotic behavior of the effective potential (e.g.atomic anions or polarizabilities).The most commonly used approach that removes the SIE on an orbital-by-orbital basis is the well-known Perdew-Zunger self-interaction correction (PZSIC) method [3].The performance of the PZSIC method is however paradoxical [4,5].The PZSIC violates the uniform gas limit of the uncorrected DFA.Recently introduced locally scaled self-interaction correction (LSIC) method recovers the uniform gas limit of the uncorrected functional and like PZSIC method is exact for one electron densities [6].Conceptually, in LSIC method, an iso-orbital indicator is used to identify the nature of the charge density at each point in space, distinguishing one-electron-like regions from regions where the density is slowly varying and many-electron-like.The self-Coulomb and self-Coulomb energy densities of each orbital are then scaled locally such that full SIC is maintained in one-electron regions but in slowly varying regions, the SIC is scaled down.This results in significant improvements over PZSIC.LSIC method has provided accurate results for properties such as atomization energies [6,7], electron affinities, magnetic properties, dipole moments, and polarizabilities [8][9][10], vertical electron detachment energies [11,12], exchange-coupling constants [13] and spin-state gaps of octahedral complexes [14].It also significantly reduces fractional electron or delocalization errors [15].
In LSIC method, the total energy is given by expressed as Here, and The z σ (⃗ r) is the iso-orbital indicator which is used to identify the one-electron regions and to scale the energy densities in the SIC terms.We refer the interested to [6,7] for more details.The functional derivative of E LSIC with respect to variation in density is written as follows, δE LSIC δρ (r) This equation can be used to construct the Fock matrix needed for the self-consistent LSIC calculations.The second term in this equation requires the computation of the Coulomb potential of the electron density scaled at each point.The present manuscript focuses on the efficient evaluation of the second term.However, before we detail our numerical approach, we provide a brief overview of the typical methods employed in mainstream quantum chemistry to compute the Coulomb potential or energy.Quantum chemistry is typically carried out using Gaussian-type orbitals (GTOs) due to the advantage of the analyticity of Gaussian product rules which allows efficient analytic evaluation of integrals like the basis overlap integrals, Coulomb integral, exchange integrals, dipole integrals etc [16].Over the years, many efficient schemes have been developed around GTOs to tackle computationally costly two-electron integral in Coulomb problems.References [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] Starting from the Boys function for Coulomb potential due to a spherically symmetric Gaussian charge distribution, the potential due to charge involving higher angular momentum components can be efficiently constructed recursively.The most common recursive schemes are the ones proposed by McMurchie and Davidson [17] and Obara and Saika [18].These common GTO integral evaluations have been implemented in integrals libraries such as libint and libcint [33,34] which are at the heart of many software packages.Despite the fact that recursive algorithms considerably enhance speed, the poor scalability of 4-center integrals remains an issue.In order to obtain better scaling, Termath and Handy [24] proposed to analytically evaluate the Coulomb potential on a numerical grid, building on the idea of combining a numerical grid and a finite basis set [19].This strategy not only significantly reduces complexity, but it may also be paired with either normal far-field expansion [35] or an efficient and accurate algorithm such as the faster multipole method (FMM) and tree codes [25][26][27][28][29][30].The advantage becomes more apparent when greater angular momentum basis functions are involved.
A pure numerical basis-set-free approach to the Coulomb problem in polyatomic systems was presented by Becke and Dickson [36].Their method utilizes a multi-center numerical grid [31], that partitions the space into multiple regions.Such partitioning of space facilitates the integration within each region to be performed independently by solving the Poisson equations for the charge density confined in each region.By combining with multipolar expansions, one can then obtain a set of radial parts of the Poisson equations.Once the radial equations are solved, the solutions can be used to construct Coulomb potential in the form of spherical harmonic expansion (SHE).This approach is appealing but unfortunately, errors accumulate as the system size increases [24] An extensive literature search showed that only a handful of implementations [24,[37][38][39] and the majority of these studies concentrated on small molecules [38][39][40][41][42].
Shortly after, based on the same theoretical framework, Delley used the integral form of the Green's function (GF) solution to the Laplacian [32] as an alternative to solve the radial part of the Poisson equations.This approach removes the problem of error accumulation with system size, thereby opening a gateway to broader numerical Poisson applications.The method is now more widely used and has been implemented in a few DFT codes [43,44] for large scale calculations.
The numerical Poisson solver based on multi-center grid and multipolar expansion has been discussed in several earlier works and shown promises.But it is mostly adopted in grid [39], Slater-type orbital (STO) [45,46] or numerical atomic orbital (NAO) based codes [44,47], and rarely applied to Gaussian-based DFT codes as the Coulomb potential can already be calculated analytically [17,18].There have been only a handful of earlier works that attempted at applying the numerical Poisson solver to GTO based first-principles calculations on rather small systems (i.e.mostly less than 10 atoms).As discussed earlier in this section, our interest in the numerical Poisson solver is due to the need to evaluate the Coulomb potential of electron density scaled by an arbitrary function (iso-orbital indicator) for self-consistent LSIC method.
The numerical evaluation of Coulomb potential can also be advantageous for large systems as each center is treated independently, and therefore can be parallelized easily for large-scale calculations [43,46].
In the majority of instances, Dunlap correction demonstrates quadratic convergence in energy, and combining electrostatic energies is a common practice to mitigate error accumulation.However, in specialized cases like ours in LSIC method [6,7], implementing additional correction techniques can be impractical or cost-prohibitive.In response to this challenge, we propose a straightforward method to directly enhance the accuracy of radial solutions, thereby reducing errors.As the method adopts two existing methods for different components of spherical harmonic expansion (SHE), it does not imposes additional limitation beyond the original methods.This approach is particularly helpful because it addresses the challenges associated with finding suitable screening charges for computing the Coulomb potential and energy of scaled orbital charge densities of arbitrary shapes, which are necessary for calculating self-interaction corrections (SIC) in the one-electron SIC approaches [3,[5][6][7].Moreover, for the standard Coulomb problem, our implementation remains independent of Dunlap correction and charge screening.As a result, it presents an additional option for enhancing accuracy at essentially no extra cost.
We first provide a brief overview of two prominent numerical methods [32,36].We then perform numerical tests of these methods with various parameters and techniques, with the goal of determining the source of the strengths and weaknesses of each method.By analyzing the results, we suggest a hybrid method to further enhance the efficiency of the numerical Poisson solver.

Multi-center grid and single-center decomposition
The numerical Poisson solver utilizes so-called multi-center grid [31], or often collectively referred to as Becke mesh, that was originally designed for 3D molecular integrals for functionals of the form F(ρ(⃗ r), ∇ρ(⃗ r), ∇ 2 ρ(⃗ r)).To properly described the cusps at the nuclear cores, the grid itself is constructed as a superposition of multiple spherical integration grids where each spherical grid is constructed by multiplying a radial quadrature onto a spherical mesh to form concentric spherical shells of mesh centered on each atomic site.
While Lebedev quadrature is generally believed to be the most efficient for spherical mesh [48], there have been a wide variety of options for choosing radial grids [49][50][51][52][53][54][55][56][57], and each has its own advantages.In this study, we adopt the radial quadrature proposed by Mura and Knowles [55], as it has been proven to be numerically efficient [52] and has a simple form where and N is the number of radial points.Both α and m are empirical parameters controlling how the points are distributed.The recommended values are m = 3 and α = 5.0 (7.0 for alkali and rare-earth metals) [54].This simple expression makes it easier for analytically evaluating the coefficients of the finite-difference method (FDM) operator for solving 1D radial Poisson equations.
The overlap between the spherical meshes centered on different atomic sites can be avoided by scaling down the integrated function ρ(r) with the partition weight function w n as and This trick, similar to Voronoi charge construction, divides the entire space into multiple independent cells but separated by smoother transition boundaries.In this case, each cell contains only one nucleus and can therefore be efficiently expressed in spherical coordinates.The partition weight functions w n are constrained to satisfy the condition and defined as where P n (⃗ r) is the cell function.
Over the years, many cell function generating schemes have been proposed.A prominent one is due to Becke [31] which requires predetermined empirical parameters to accommodate different atomic species of different sizes.Delley [32] provides some easier atomic density dependent implementations that are convenient for the molecules that contain different atomic species.Some weight schemes are also proposed to achieve better scaling for finite systems [58] and for periodic systems [46].There are also variants based on Becke's original scheme, which are designed for certain properties of interest [59,60].
Although the mesh generated with these cell function schemes can be applied to both general three dimensional molecular integration and numerical Poisson solver, they are not necessarily the same [32].The weight functions created with different cell functions are generally in good agreement; nevertheless, our experience indicates that the approach described by Stratmann et al [58] produces the most accurate results in particular geometries, confirming the earlier findings [44].As a result, Stratmann's method is used for all of the calculations reported in this paper.

Multipolar expansion
The starting point is to partition the source charge density in the Poisson equation using equations (7) and (8).This results in significant simplification as one only needs to deal with the equation containing one atomic center (n) at a time.
As the mesh for each center is constructed based on spherical mesh, we can efficiently expand both ρ n and V n with spherical harmonic functions [35] respectively as and where and V (n) lm (r) are to be solved on the radial grid r i .

Radial Poisson
For each lm component in SHE, the angular and radial degrees of freedom are separable, the 3D Poisson equation can then be converted into a set of 1D problems.There are two major approaches to solving for V lm which is also where Delley and Becke's methods differ.
The following only provides a summary that is relevant to the discussion.For more details, readers are directed to the original papers for both methods.
In Becke's approach, V lm (r) are obtained by solving a set of 1D differential equations.With the substitution one arrives at the following general expression with the coefficients ∂x ∂r 2 , ( 16) and Equation ( 15) can be evaluated on the radial grid r = r i using FDM.
The boundary conditions (BCs) are imposed as U lm (r) = 0 for r = 0 (r 0 ) and r → ∞ (r N+1 ) for all lm's except for the monopole l = 0 where U 00 (r → ∞) = √ 4π q n .On the other hand, Delley applied the GF solution to the Laplacian [61] by integrating the radial charge directly for I lm (r i ) [32,43,44] in the following form where s i lm is the fitted radial charge density evaluated at mesh point r i .In this study, we also try to further smooth out the fitted curve s i lm (r ′ ) by rearranging the integral form of both terms in equation ( 19) into the following expression where p = 0, 1 or 2 and we perform cubic spline interpolation on r p s i lm (r ′ ) instead.The solution of the radial Poisson defined in equation ( 13) can be expressed as This approach is in principle straightforward; however, the integration for equations (19) and ( 20) is tricky, as the integration needs to be done for every r i and can no longer benefit from the quadrature weights.The numerical instability could potentially occur in both the interpolation of the charge and the evaluation of equation (19) as the term involves higher order polynomial terms.In some of the earlier works, such as [44], the integration is first evaluated by spline-interpolating the multipolar charge density onto a denser mesh, and then the integration is done numerically.Here, we adopt the most recent Franchini's integration scheme [43] instead, where the multipolar charge density, s i lm (r) at interval [i, i + 1], is expressed as piece-wise polynomials then the integrals in equation ( 19) can be evaluated analytically for each segment.This approach should, in principle, yield the most accurate results for a given spline compared to the former and is also computationally more efficient.

Reconstructing the Coulomb potential
Once the 1D radial solutions V lm are calculated, the total Coulomb potential due to the charge in a particular cell can be reconstructed on any integration mesh using equation (13), where V lm = r −1 U lm (r i ) for Becke's FDM and V lm = 4π 2l+1 I lm (r i ) for Delley's GF integral.For V lm of any arbitrary given r ′ that does not coincide any radial grid point r i , the evaluation is done using cubic spline interpolation

Screening charge
As the charge density ρ(r) is normally rapid varying in the proximity of nuclei, it is often useful to introduce some type of screening charge [32,44].The ρ screen is usually chosen to smooth out the curve around the cusps in ρ(r) in such a way that its corresponding Coulomb potential can be evaluated analytically and efficiently.One then only needs to consider the Poisson problem due to the variation from the screening charge.A convenient choice for screening density is the superposition of spherically symmetric neutral atomic charge densities [32,44].
and the final total Coulomb potential can be restored by adding the analytically evaluated potential V atoms due to the screening charge,

Computational details
We use Pederson-Porezag basis set [62] in all our calculations.A Lebedev unit sphere quadrature of order l = 47 with 590 points is used for the integration mesh and the numerical Poisson solver.To preserve the banded shape of the operator matrix, we employ an 11-point central difference formula for the middle radial points and a non-central difference formula for the points involving the boundary instead of the 7-point formula suggested by Becke and Dickson [36].However, we found no discernible difference beyond the 7-point formula.The Matlab function [63] is used to generate the coefficients for FDM operators.We employed Franchini's density fitting [43] with a little modification as introduced in equation ( 20) for the GF approach to solve radial Coulomb potentials.We use a piece-wise cubic spline interpolation subroutine modified from the subroutines in numerical recipes [64] to interpolate V lm onto any given radial point r.
As the treatment for the angular part of the solutions is well-defined within SHE, the main purpose of this study will be to focus on improving the radial solutions, we choose a generic set of parameters α = 6 and m = 3, similar to the recommended values in [55], for the radial quadrature in all our calculations unless explicitly stated.

Atomic system (single-center)
Since this study aims to directly improve the quality of the radial solutions to the Poisson, the common techniques like Dunlap correction, although very efficient in improving the accuracy of Hartree energy by reducing the error due to charge deviation, will not be considered.We begin with single atoms because space partitioning is not required for atoms, allowing us to better understand how efficiently the two approaches work for a given radial quadrature.Figure 1 shows the error in Hartree energy for Mn atom (a) without and (b) with charge screening compared to the reference value.In the reference calculation, the Coulomb potential is analytically evaluated [24].It is apparent from figure 1(a) that in their primitive form, FDM performs better than GF (p = 0), as it requires less than 120 radial points to achieve 10 −5 Hartree accuracy, while GF provides barely a mHatree accuracy.However, once the additional polynomial r p (p = 1 and 2) as defined in equation ( 19) is introduced to smooth out the density peak close to the origin, the performance of the GF improves drastically.Particularly with p = 2, the GF result is almost comparable with FDM.This  suggests that the charge density variation in the nucleus region is too large for the cubic spline polynomials to efficiently describe.Smoothing out the curve becomes essential to reduce the interpolation error.Figure 1(b) shows the results when a screening charge is applied.Unsurprisingly, with the inclusion of screening charge that naturally reduces the curvature of charge density near the nucleus, the error for both methods is greatly reduced to about 10 −7 Hartree, with only about 100 radial points.It is obvious that FDM still outperforms all the GF cases with different p's, even though the accuracy gap has become much smaller.It is also interesting to note that in this particular case, p = 1 or 2 does not seem to improve the result at all.
Figure 2 presents the same analysis as figure 1(b) but on a single Copper atom.Please note that the black dashed line for GF (p = 0) follows the axis to the right, while the other three follow the left axis.By applying screening charge, all four cases offer adequate accuracy with 100 radial points (i.e.≲10 −5 Hartree).Noticeably, the error of FDM is still smaller than all three GF cases, while the result of GF improves as p increases, especially when compared to its original form (p = 0) where the error is at least an order of magnitude larger than the rest.
In this case, unlike in Mn atom, the polynomial factor improves the accuracy even in the presence of charge screening.Although both screening charge and r p can improve the interpolation for the GF method, it is still possible that the piece-wise cubic spline is not sufficient, as it incorporates only two adjacent data points (i.e.radial quadrature) for each segment, while FDM naturally incorporates multiple data points during the 2 N+1-point (N = 5) difference operator construction.As a result, FDM is always more accurate with a more stable convergence even without charge screening.
In principle, the result for GF can be further improved by optimizing radial parameters α and m.However, the purpose of this analysis is to demonstrate the effectiveness of the FDM, where the method generally describes the monopole better and is less sensitive to the radial quadrature parameters.This is particularly useful for solving for the potentials due to orbital densities needed in the self-consistent one  (a) FDM  electron self-interaction methods [3,6,65,66] where suitable screening charges are either difficult or expensive to find.

Molecular system (multi-center)
Although FDM outperforms GF for single-atom systems with only one center, the outcome is reversed in multi-center systems (i.e.molecules).The table 1 displays the error per atom for various radial Poisson solvers on a variety of molecules ranging from 10 to 33 atoms.The FDM becomes problematic in this case.Not only is the error large in general, but it also scales with system size.As previous studies have already pointed out [24], the error appears to be accumulative and becomes worse with increasing system size.The increasing error with system size has limited the method from being used widely.Furthermore, we also noticed that the result is very sensitive to the radial mesh scaling factor α.An unreasonably large α often yields better accuracy than a small one with larger N rpt .This is rather counter-intuitive.Since FDM performs so well in single-atom systems which means the solver is capable of describing the monopole that contains a sharp peak.Naively, one would expect the charge partitioning should mostly affect the multipolar expansion (i.e.requires larger l).Since the spherical symmetry is broken by the partitioning.
Table 2 shows the convergence of the l-resolved Hartree energy contributions for Glycine defined as , where  is evaluated from the potential defined as The data reveal something rather intriguing.Particularly in FDM, a perfect convergence can be seen in all terms except for l = 1 which exhibits errors that are 3 ∼ 4 orders of magnitude larger than other l components.The result in table 2(a) for FDM is also plotted in figure 3 for better visualization.From figure 3(a), it becomes obvious that the error is entirely coming from l = 1 while all other terms converge rapidly to 0, even the dominant l = 0 term is well converged with around N rpt = 100.Figure 3(b) compares the error of only l = 1 and the error of the entire Hartree energy.The two curves coincide nearly perfectly which further confirms not only that l = 1 is the only major source of error but also that the contributions from all the higher orders terms converge perfectly.
To better understand this rather peculiar cause of error, we analyze the radial part of the potential, U lm , for the Oxygen site in the Glycine molecule of the lowest four orders as shown in figure 4 using four distinctively varied sets of radial mesh parameters.Since not all terms are nonzero due to symmetry consideration, only the largest contribution is shown for every l.The parameters (α, N rpt ) for the four quadratures are chosen as (6.0, 90), (6.0, 300), (30.0, 90) and (30.0, 300).A large number of radial points N rpt = 300, is chosen here intended as an accurate reference.A large scaling factor α simply stretches out the radial mesh to cover wider space range and as a result will make the mesh grid sparse.At the first glance, all four radial meshes agree nicely at l = 0 and l = 3 (and the same for all the higher order terms not shown here).On the other hand, a large deviation in l = 1 and a very small but noticeable difference in l = 2 can be seen, which indicates that these two l−orders are more sensitive to the radial grid.Upon closer inspection, only the ones with the larger α have converged for l = 1 regardless what N rpt is.This means that N rpt = 90 is sufficient to describe the region around the nucleus and the maximum range of the quadrature r max = r Nrpt which scaled directly with α ( see equation ( 5)) has a significant influence on the accuracy.
As expected, (6, 300) gives a better result than (6,90), but this is partly because the increase of N rpt also increases the r max i = 20.5395for (6, 90) and 27.6709 for (6,300), rather than the fineness of the mesh.It is now clear that the source of the discrepancy is about how the asymptotic behavior is described by FDM rather than the description of the nucleus core region.By looking at the asymptotic tails at a larger distance r = 20, U max lm , both l = 0 and l = 4 (and above) have reached their asymptotic limit √ 4πq n and 0 respectively.On the other hand for l = 1, U max 1m (r = 20) is still significant and decreasing to 0 in an extremely slow pace (roughly 1/r), hence introduces a large deviation.For the l = 2 term, there is a visible but much smaller deviation at r = 20, therefore leaving smaller room for discrepancy compared to l = 1.
Intuitively speaking, if U lm (r) at the largest radial point r = r Nrpt has a significant deviation from the asymptotic value (i.e.√ 4πq n for l = 0 or 0 for l ̸ = 0), a significant error is expected, as FDM only knows the boundary condition at the infinity, the behavior beyond the last radial point r Nrpt can only be extrapolated from the last few radial points.This explains why a significant error only emerges in the lower order terms but not l = 0, as it reaches √ 4πq n rapidly within few Bohr's, while higher order terms are less likely to suffer from this since they are shorter in range and decay to 0 much faster.
This also explains why the method still works well for smaller molecules, as the deviation in the l = 1 term in figure 4 is not obvious until r is large enough.It is also worth mentioning that in the original attempt, Becke uses Gauss-Chebyshev formula of the second kind for the radial quadrature, which is known to emphasize too much in the extended region while not putting enough points into the chemical bonding region [52].However, despite not being integrationally optimal, using this quadrature could potentially mitigate the boundary issue that appears in U 1m terms described above, particularly for smaller size systems.

Hybrid radial solver
As discussed in the previous section, the FDM method struggles for l = 1 and possibly l = 2. On the other hand, GF method only requires the evaluation of the charge integrals of the form r k(l) ρ lm (r) (see equations (19) and ( 20)), where k is just some integer depends on l chosen manually and does not suffer from the same 'boundary' issue.GF generally converges rather nicely with the number of radial points for all terms except for l = 0 as shown in the table 2. This is because the charge density in each center is more localized in space after the truncation imposed by the weight function.
By considering the strengths of both methods, we propose a hybrid approach where FDM is used for the 'near-field' l = 0 and the rest are calculated using GF integration.The detail of the procedure is outlined in figure 5.For a large system, the first term of equation ( 20) which corresponds to the far-field multipole expansion, is employed to ensure the asymptotic behavior is accurately preserved.In the rest of the discussion, the screening charge is applied to all the calculations unless explicitly stated.
The test used in the analysis of FDM and GF method, is also used to analyze the performance of the proposed hybrid approach.The results are shown in table 1.It is evident from the table that using the same generic set of parameters (α, N rpt ) = (6.0,90), the accuracy of the proposed hybrid approach is improved by several times to an order of magnitude compared to the GF method.More importantly, the error per atom in the proposed hybrid approach does not scale with the system size.According to our analysis, both FDM and GF describe higher order terms (i.e.l > 3) equally efficient and accurate.Therefore, one would expect to see the same improvement for the FDM-based method as long as l = 1 terms are calculated using GF.

Energy difference between different spin states
As an assessment of accuracy of the proposed hybrid in practical DFT calculations, we compute the energy difference between different spin states, namely high spin (HS) and broken symmetry (BS), for two systems p-C 8 H 8 and Cu 2 Cl 2− 6 .The results are summarized in table 3. The error in the energy difference between different spin states in both cases is much smaller than the error in total energy alone by 1 ∼ 2 orders of magnitude, with the error percentage only about 0.01%, as it is often easier to achieve higher accuracy in energy difference than in total energy alone.This makes the method particularly useful for studying properties involving energy differences, such as exchange coupling or   magnetic anisotropy energy (MAE).The similar pattern for the energy difference is also observed in [44].
An additional convergence test with respect to the highest l included is shown in figure 6 for the total energy error compared to the reference energies for (a) Glycine and (b) Cu 2 Cl 2− 6 of both HS and BS states and the energy difference between the two states.While all total energies require around l = 16 to reach convergence, it requires only up to l = 6 (see figure 6(b)) to converge the energy difference.Similar patterns can also be found in many different contexts, such as Brillouin zone integration in MAE.

Conclusions
In summary, we have examined the behavior of the two main numerical techniques (GF and FDM) for solving Poisson equations in order to enhance their efficiency and/or modify them to meet our needs for calculating SIC energy and potentials in the one-electron SIC methods.We observed that, given the same set of radial parameters, FDM appears to be more accurate in describing the monopole (l = 0) but it inherits a serious issue coming from the multipolar contribution of mainly l = 2 and sometimes l = 2.We propose a hybrid scheme by combining the two approaches, where the radial potential of monopole (i.e.U 00 ) in spherical harmonic expansion is calculated using FDM while the rest of the higher order terms are calculated using GF's function integration.We then performed a series of tests on different systems and the convergence with difference numerical parameters to demonstrate the effect of the implementation.The overall results suggest that the accuracy is improved in all the cases.More importantly, our analysis shows that like GF method, the proposed hybrid approach is free from error accumulation with system size and is more accurate for the monopole than the GF method.This observation and its favorable scaling especially when large basis sets are employed makes it an attractive method for studies on large complexes.The hybrid approach presented in this manuscript has allowed us to perform self-consistent LSIC calculations for one electron SIE-free density functional calculations [7,15].We have found that the hybrid approach presented here, while generally useful, is especially advantageous for the computation of self-interaction corrections in one electron SIC method (scaled), wherein the Coulomb potential and energy of arbitrary-shaped (noded and rapidly varying) orbital charge densities need to be computed.

Figure 1 .
Figure 1.Hartree energy EHart for a single Mn atom is calculated using both FDM and GF (a) without and (b) with screening charge.Different modifications of charge interpolation are also performed as a comparison as shown for p = 0, 1 and 2.

Figure 2 .
Figure 2. The error in Hartree energy for a single Cu atom versus the number of radial points Nrpt of modified charge interpolations (i.e.r p s i lm with different shown.Generally, the best performance can be achieved with n = 2.

Figure 3 .
Figure 3.The deviation of Hartree energy of Glycine molecule of each l-component from the same quantity but evaluated with Nrpt = 300 as a convergence test.(a) The Hartree energy contributions from the lowest six orders (l = 0 ∼ 5).(b) the largest error contribution (l = 1) compared to the error of the total Hartree energy.

Figure 4 .
Figure 4. sets of radial parameters are used to test the convergence of the radial solutions.The maximum component of the radial solutions (U max lm ) of the four lowest expansion orders (up to l = 3) for the oxygen site in Glycine molecule are shown as (a)-(d).For (b)-(d) the insets present a closer look of the convergence of the tails.

P-H Chang et alFigure 5 .
Figure 5.The flow chart of the numerical Poisson solver.

Figure 6 .
Figure6.The energy convergence versus multipolar expansion order l for (a). the total energy error of Glycine compared to the reference value and (b).total energy of Cu2Cl 2− 6 of both HS and BS spin states and the energy difference between the two.

Table 1 .
Error per atom (∆E) in the Hartree energy (in Hartree atomic unit) for a set of testing molecules of different sizes range from 10 to 60 atoms.

Table 2 .
The l-resolved energy convergence for the Glycine molecule of (a) FDM and (b) GF method versus number of radial points.The reference energies are taken as the energy evaluated with Nrpt = 300 for both methods.

Table 3 .
Total energy and the total energy difference between high spin and broken symmetry states for p-C8H8, and Cu2Cl 2− 6 .All energies are in unit of Hartree.