We present a robust strategy to numerically sample the Coulomb potential in reciprocal space for periodic Born–von Karman cells of general shape. Our approach tackles two common issues of plane-wave based implementations of Coulomb integrals under periodic boundary conditions: the treatment of the singularity at the Brillouin-zone center and discretization errors, which can cause severe convergence problems in anisotropic cells, necessary for the calculation of low-dimensional systems. We apply our strategy to the Hartree–Fock and coupled cluster (CC) theories and discuss the consequences of different sampling strategies on different theories. We show that sampling the Coulomb potential via the widely used probe-charge Ewald method is unsuitable for CC calculations in anisotropic cells. To demonstrate the applicability of our developed approach, we study two representative, low-dimensional use cases: the infinite carbon chain, for which we report the first periodic CCSD(T) potential energy surface, and a surface slab of lithium hydride, for which we demonstrate the impact of different sampling strategies for calculating surface energies. We find that our Coulomb sampling strategy serves as a vital solution, addressing the critical need for improved accuracy in plane-wave based CC calculations for low-dimensional systems.

## I. INTRODUCTION

Coulomb integrals manifest in numerous theories and models within the domains of many-body physics and quantum chemistry. However, under periodic boundary conditions, numerically modeling the periodic Coulomb potential is a non-trivial endeavor. While its well-known reciprocal form of 4*π*/*q*^{2} is beyond question, there is no ideal strategy to sample this function in convolutions with electron densities on a discretized grid. Such samplings are necessary for the implementation of computational methods such as the Hartree–Fock (HF) theory, Møller–Plesset perturbation theory,^{1} or coupled cluster (CC) theories.^{2–4}

In this study, our primary focus is on a specific contribution to the finite-size error stemming from the finite sampling of the Coulomb potential. While the integrable singularity at *q* = 0 has been extensively investigated in prior studies,^{5–12} we address a less-explored concern related to discretization errors in anisotropic Born–von Karman (BvK) cells,^{13,14} which can lead to non-physical or non-convergent behavior in HF and CC calculations. We show that this behavior comes from discretization errors that arise in anisotropic BvK cells, which are particularly useful in calculations of low-dimensional systems.

To this end, we present an effective sampling technique for the reciprocal Coulomb potential, which is designed for anisotropic unit cells and for mitigating the effects resulting from discretization errors. We show how this alleviates the above-identified problems that can arise when running many-electron calculations for these systems. The modified potential is universal in character and can just as well be applied to cells of general shape.

## II. THEORY

*V*, of a periodic system reads

^{11}

*φ*

_{i},

*φ*

_{j}, which are normalized in the unit cell, as well as two integrals over the first Brillouin zone with the volume Ω

_{BZ}.

**is a momentum transfer vector. The equivalent formulation in reciprocal space allows us to replace the two continuous real-space integrals with one discrete sum over all reciprocal lattice vectors,**

*q***,**

*G*In order to evaluate the exchange energy expression numerically, it is necessary to discretize the two Brillouin zone integrals. In this context, two prominent challenges manifest themselves. First, a discretization error is inevitably introduced, and second, the presence of the notorious, albeit integrable, singularity at ** G** +

**= 0 necessitates careful handling if the co-densities contain a monopole. In the exchange energy expression, this always happens for**

*q**i*=

*j*, as $\rho iik,q=0(G=0)=1$. We note that these issues originate solely from the used cell shape and discretization grids, which means that no specific class of material is affected, but rather all calculations in which Coulomb integrals are used are affected. While an extensive body of literature addresses the singularity issue,

^{5–12}less attention has been devoted to mitigating the discretization error,

^{13,14}especially in cases involving anisotropic cell shapes.

*i*,

*j*,

**, and**

*k***, one has to discretize the integral over**

*G***in the expression of the exchange energy in Eq. (2). This can be recast as**

*q**v*(

**) ≔ 4**

*q**π*/(

**+**

*G***)**

*q*^{2}and $f(q)\u2254\rho ijkq(G)2$. A mesh of

*N*q-points, {

*q*_{n}}, is commonly introduced in order to discretize Eq. (4) as

^{15}or uniform

**Γ**-centered meshes are considered, where

**Γ**is the center of the Brillouin zone (BZ). However, it should be noted that

*v*(

**) can vary significantly stronger between the q-points than the co-densities in**

*q**f*(

**), in particular if |**

*q***| is small or if the mesh is not equally dense in all spatial dimensions (e.g., in anisotropic cells). This is only the case for the integral over the so-called momentum transfer**

*G***, not for the**

*q***integral. The following**

*k**exact*reformulation of Eq. (4), where we simply split the integration range into

*N*distinct sub-spaces $\Omega BZ=\u222anN\Omega BZn$, suggests an improved strategy to discretize the

**integral:**

*q*

*q*_{n}than they are to any other points

*q*_{m≠n}. If we assume

*f*(

**) to be roughly constant within the subspace $\Omega BZn$, and hence**

*q**f*(

**) ≈**

*q**f*(

*q*_{n}), we can approximate

*n*and

**, as a**

*G**mean*Coulomb potential of the subspace $\Omega BZn$,

*mean”*in this manuscript, decouples the integration meshes for

*f*(

**) and**

*q**v*(

**). While conventional q-meshes can be used to evaluate Eq. (7), more sophisticated numerical integration techniques have to be applied to evaluate the deceptively simple integral in Eq. (8). In Sec. III B, we describe our refined numerical integration technique. Since $v\u0304n(G)$ depends only on**

*q**n*and

**, it can be pre-calculated and stored in memory once the lattice and the q-mesh are defined. Using**

*G**mean*strategy from Eq. (7), the final expression for the exchange energy for a mesh of

*N*k-points and q-points finally reads

The *mean* strategy can readily be applied to Coulomb integrals of any post-HF method. In addition to the Coulomb integrals of the HF method itself, we apply it to the Coulomb integrals for the equations of coupled cluster with single, double, and perturbative triple particle–hole excitation operators [CCSD(T)]. For more details on the implemented equations in CCSD(T) theory, we refer to Ref. 16.

*mean*strategy with two other common strategies, which are summarized in Table I. With “

*probe*,” we abbreviate the so-called probe-charge Ewald method.

^{7}In the

*probe*strategy, a usual discretization is performed, as formulated in Eq. (5), but for the singularity at

**+**

*G***= 0, the difference between an analytically calculated self-energy of a**

*q**probe*-charge and the corresponding numerically computed self-energy using the given q-grid is employed,

*α*> 0 but very small, such that

*v*

_{0}is almost independent of

*α*. The

*probe*strategy is the default setting in the vasp code and provides numerically comparable results as the auxiliary function strategy by Gygi and Baldereschi.

^{5,10}

Strategy . | Coulomb potential . |
---|---|

Mean | $1\Omega BZn\u222b\Omega BZnd3q4\pi (G+q)2$ |

Probe | $4\pi (G+qn)2,G+qn\u22600,v0,G+qn=0$ |

Disr | $4\pi (G+qn)2,G+qn\u22600,0,G+qn=0$ |

Strategy . | Coulomb potential . |
---|---|

Mean | $1\Omega BZn\u222b\Omega BZnd3q4\pi (G+q)2$ |

Probe | $4\pi (G+qn)2,G+qn\u22600,v0,G+qn=0$ |

Disr | $4\pi (G+qn)2,G+qn\u22600,0,G+qn=0$ |

Second, with “*disr*,” we refer to the strategy to simply neglect (“disregard”) the grid point at the singularity at ** G** +

**= 0 in the discretization (5). The**

*q**disr*strategy is the easiest way to achieve a functional periodic implementation of the Coulomb integral and is used in some studies.

^{12,17}Note that

*disr*only addresses the integrable divergence but not the discretization error.

We do not consider the spherical truncation scheme^{9} since it is inadequate for anisotropic BvK cells. The Wigner–Seitz truncation strategy^{11} is also not considered, as no implementation was available for this work.

While preparing this manuscript, it was brought to our attention that an analogous approach was used for GW calculations of excited states in Ref. 18. While it would be of interest to make a comparison between the method presented here and Ref. 18, currently, no technical capability exists to do so because the method developed here focuses on CC ground-state calculations.

## III. METHODS AND IMPLEMENTATION

### A. Software and methods

We implemented the *mean* Coulomb potential strategy in the Vienna *ab initio* simulation package (vasp)^{19} and present the implementation in Sec. III B. We refer to the associated INCAR tag in vasp as LHFMEANPOT = T.

Through the interface between vasp and the cc4s code,^{20} we can directly compare *mean*, *probe*, and *disr* in periodic HF and CC calculations. All calculations were performed using a plane-wave cutoff of 500 eV. The following projector augmented wave (PAW) pseudopotentials with frozen-core potentials ($POTCAR$ files) were employed for the test systems: $PAW_PBEC_GW_new$, $Li_GW$, and $H_GW$. The unoccupied space for the CC calculations was compressed to 15/12 natural orbitals per occupied orbital for the carbon chain/LiH system. Natural orbitals are obtained on one-particle reduced density matrices obtained from the RPA^{21} or approximate MP2 level^{22} of theory. For the carbon chain, a highly effective basis-set correction scheme^{23} was employed to approach the complete-basis set limit of the CC correlation energies. Using *mean* potential, we employed BvK cells with a vacuum of $8A\u030a$ and 20/10 replicas of the two-atomic unit cell to calculate the potential energy surface of the carbon chain at the HF/CCSD correlation level. The used BvK cells for the LiH slab are discussed in Sec. IV. Furthermore, in this work, we exploit the equivalence between an *N* × *N* × *N* supercell employing Gamma-only sampling of the BZ and a single unit cell employing *N* × *N* × *N* sampling of the BZ. Here, all BvK cells are modeled as supercells using Gamma-only sampling of the BZ.

### B. Implementation of the mean potential strategy

*mean*potential of Eq. (8) is described in the following, addressing the perils of its seemingly simple form. We partition the integral into a term designed for a radial grid and a term designed for a uniform grid by means of the identity $1=e\u2212\gamma q2+(1\u2212e\u2212\gamma q2)$ so that we can write

_{n}(

*q*,

*θ*,

*φ*) is 1 if the point (

*q*,

*θ*,

*φ*) is in the integration volume $\Omega BZn$, and 0, otherwise.

The first integral on the right-hand side of Eq. (11) is numerically implemented on a radial grid, whereas the second integral is implemented on a uniform grid along the axes of the reciprocal vectors. Note that no singularity occurs in either integrands. While the Jacobian determinant of the spherical coordinates cancels the singularity in the first integrand, the second integrand approaches the finite value $(1\u2212e\u2212\gamma q2)/q2\u2192\gamma $ for *q* → 0. In our implementation, we set $\gamma =8ln(10)\u22c5Lmax2/(2\pi )2$, where *L*_{max} is the largest edge length of the BvK cell. This choice ensures that the factor $e\u2212\gamma q2$ reaches virtually zero inside the neighboring integration regions, $\Omega BZn$, of the region containing the origin ** G** +

**= 0. For all other integration regions, the first integral is virtually zero, and its numerical evaluation can be dropped.**

*q*The second integral, using the uniform grid, varies little even at the origin because of the numerator $1\u2212e\u2212\gamma q2$ and by the choice of *γ*. Thus, a relatively small number of grid points can be used. For the integration region containing the origin ** G** +

**= 0, the number of grid points in the direction**

*q**i*is chosen according to $Ni\u22611000A\u030a/Li$, where

*L*

_{i}is the corresponding edge length of the BvK cell in units of Ångström. Based on this number, we halve the number of grid points for each neighboring integration domain, $\Omega BZn$, that moves away from the origin.

Finally, we note that smearing the step function Θ_{n}(*q*, *θ*, *φ*) improves the convergence of the numerical integration on the radial grid significantly. In our implementation, we multiply Fermi–Dirac-like functions, $exp(8d/\Delta q)+1\u22121$, in each direction, to smear the sharp step of Θ_{n} at the boundaries of the integration regions $\Omega BZn$. Here, *d* represents the signed distance between an integration boundary in one direction to the grid point (*q*, *θ*, *φ*) and Δ*q* is the grid spacing of the *q* coordinate. We use the following radial grid in our implementation: *N*_{q} = *q*_{max}/Δ*q* ≡ 20 with *q*_{max} = 2*π*/*L*_{max}, and *N*_{φ} = ⌊2*π* × *N*_{q}⌉ = 126, as well as *N*_{θ} = ⌊*π* × *N*_{q}⌉ = 63. The symbol ⌊⋯⌉ denotes rounding to the next integer.

The described implementation of the *mean* Coulomb potential strategy allows us to achieve sub-meV accuracy of the Hartree–Fock energy per unit cell with respect to the number of grid points.

## IV. RESULTS AND DISCUSSION

### A. Infinite carbon chain

We now turn to our first benchmark system, the infinite chain of carbon atoms in vacuum, where we aim for the calculation of the lattice constant and the bond length alternation, which goes back to a Peierls distortion.^{24} To our knowledge, no canonical periodic CCSD(T) result has been reported in the literature so far. Based on an extrapolation of finite hydrogen-terminated chains, a CCSD(T) result of $0.125A\u030a$ was estimated for the infinite chain by Wanko *et al.*^{25} However, this is in disagreement with the bond length alternation (BLA) of $0.136A\u030a$ and the lattice constant of $2.582A\u030a$ from a periodic Diffusion Monte Carlo (DMC) calculation reported in Ref. 26. A further many-electron correlation calculation under periodic boundary conditions was reported by Ramberger and Kresse^{24} using the random phase approximation (RPA), resulting in a BLA of $0.129A\u030a$, but using a fixed lattice constant of $2.575A\u030a$. Here, we report the first canonical and periodic CCSD(T) result for the BLA and the lattice constant. We calculate five CCSD(T) energies each for BLAs in the range of $(0.075,0.175A\u030a)$ and lattice constants in the range of $(2.425,2.725A\u030a)$. These 25 data points are then fitted using the function $E(b,a)=\u2211n,m=03cnm(b\u2212b0)n(a\u2212a0)m$, where *b* and *a* are the BLA and the lattice constant, respectively, and we explicitly neglect the coefficients *c*_{11} = *c*_{31} = *c*_{32} = *c*_{13} = *c*_{23} = *c*_{33} = 0. Optimizing the 12 fitting parameters, we find a BLA of $b0=0.128A\u030a$ and a lattice constant of $a0=2.578A\u030a$ with a standard deviation of the residuals of 0.5 meV per unit cell. The CCSD(T) data points and a heatmap of the potential energy surface can be found in the supplementary material.

These calculations require strongly anisotropic BvK cells, which prohibited previous HF and CCSD(T) calculations, the reasons for which will be evident below. For growing BvK cells in the chain direction, *z*, the differences of the three strategies *mean*, *probe*, and *disr* are evident for absolute HF energies, as can be observed in Fig. 2(a). While *disr* diverges, both *mean* and *probe* seem to converge to finite, although perhaps slightly different, limits. This is not surprising, as with fixed dimensions in *x* and *y*, we are not approaching the TDL, which should indeed be unique.

The divergence of *disr* is due to the fact that we approach the non-integrable singularity of 1/*z*^{2} at *z* = 0. More precisely, this non-integrable singularity is approached if only the *z* component of the BvK cell is increased, turning only ∑_{z} → ∫d*z*, while ∑_{x,y} remains untouched. Note that in three dimensions, d*x*d*y*d*z*/(*x*^{2} + *y*^{2} + *z*^{2}) = d*r* sin(*θ*)d*θ*d*φ* is free from a singularity, while d*z*/(*x*^{2} + *y*^{2} + *z*^{2}) possesses a non-integrable singularity for the discrete sampling points *x* = *y* = 0 at *z* = 0.

The *probe* strategy cancels this divergence of the non-integrable singularity and converges to a finite HF result, as *v*_{0} from Eq. (10) approaches −*∞* if only the *z* dimension of the BvK cell is increased. Both the value of *v*_{0} and the HF solution of the *disr* strategy approach −*∞* for the same reason, as explained above. A negative value of *v*_{0} can be observed in Fig. 1. There, the ratio of longest to shortest cell length is about 4.3. As discussed, *v*_{0} approaches −*∞* with increasing anisotropy. This *probe* correction provides reasonable results at the HF level of theory but fails miserably at the CCSD level, as shown in Fig. 2(d). Instead of converging to a limit, the CCSD correlation energy falls into the negative, and finally, at a length ratio of about 4.7 (i.e., at a chain length of $28.3A\u030a$), the solver fails to converge the CCSD energy after 200 iterations if the *probe* strategy is used. This is in contrast to *mean* and *disr*, for which the CCSD solver successfully converges for each length of the BvK cell.

We believe that negative values of *v*_{0} are responsible for the rapid drop of the CCSD correlation energy and the failure of the CCSD solver, as observable in Fig. 2(d). Note that the Coulomb potentials of the *probe* and *disr* strategy are equivalent, except for the value of *v*_{0}. With increasing anisotropy, the *probe* strategy turns the Coulomb potential into an unphysically pseudized potential. While this leads only to a simple energy correction at the HF level, the impact of negative *v*_{0} on the solution of the non-linear amplitude equations of the CCSD theory is unknown. What can be stated, however, is that *v*_{0} enters the CCSD amplitude equations only in the so-called ladder contributions and via an energy shift of the HF orbital energies. A deeper analysis of the failure of the CCSD iterations is beyond the scope of this communication.

In contrast to the absolute HF energies, the HF energy differences converge only for the *mean* strategy, as is apparent in Fig. 2(b). Here, an energy difference between the BLAs $0.075$ and $0.125A\u030a$ is considered. Furthermore, it becomes obvious that both strategies, *probe* and *disr*, are equivalent for energy differences within the same BvK cell at the HF level. This underlines that *v*_{0} simply introduces a constant energy shift in the HF theory while the discretization error remains.

This is no longer true for CCSD theory, as evident in Fig. 2(e), which shows the CCSD correlation contribution to this energy difference between the two BLAs. While *mean* provides relatively constant results with mild but tolerable fluctuations, the failure of the *probe* strategy is markedly visible. Except for one strong outlier, the *disr* strategy here appears to be only shifted to *mean*.

We close the discussion of the convergence behavior with the dependence on the size of the vacuum, i.e., the *x* and *y* dimensions of the BvK cell. In Fig. 2(c), one can see that all strategies seemingly approach the same TDL at the HF level, although at significantly different rates. Note that the *z* dimension of the BvK cell is already quite large with $41.2A\u030a$. This is also true at the CCSD level of theory, as depicted in Fig. 2(f), since we leave the strong anisotropic regime. Here, a reduced length of $18.0A\u030a$ was used. In the small vacuum of $4A\u030a$, we observed another failure of the CCSD solver using *probe* due to the strong anisotropy.

### B. Lithium hydride surface

In our second benchmark system, we address the impact of different sampling strategies on the electronic ground state of a surface slab model, as necessary for the calculation of surface formation energies or applications in heterogeneous catalysis. The electronic structure of a material’s surface determines its physical and chemical properties, i.e., its functionality in scientific and industrial applications. Under periodic boundary conditions, the BvK cells used to model the surface slab can exhibit anisotropic shapes, primarily due to the necessary vacuum above the slab to prevent interactions between periodic images along the normal direction to the surface. Here, we study the impact of the different sampling techniques *mean*, *probe*, and *disr* on the convergence of the per-atom energy of a (001) lithium hydride (LiH) surface slab with respect to the number of surface layers. This energy is a key ingredient for the calculation of surface free energies or the simulation of surface reconstructions. Note that we do not report a surface energy value for LiH in this study, as this requires calculations outside the primary concern of this work, including the determination of the bulk thermodynamic limit of LiH^{27} and achieving the basis-set limit.

We model the LiH surface using a 2 × 2 × *n* conventional slab, where we increase the number of layers *n* from 2 to 6. Each layer consists of 16 atoms, and a vacuum of about $30A\u030a$ is added in the normal direction to the slab. The BvK cell of the *n* = 6 layer slab is provided in the supplementary material.

Figure 3 shows the relative CCSD energy per atom of the LiH surface slabs in dependence of the number of surface layers for different sampling techniques, as well as density function theory results using the Perdew–Burke–Ernzerhof exchange–correlation functional (DFT-PBE).^{28} For a better comparison of the convergence behavior, all energies are shifted to match the two-layer slab. For both *probe* and *disr*, no convergence of the per-atom CCSD energy can be observed. A differentiated view reveals that in the case of *probe*, the HF energy converges but not the CCSD correlation contribution, and using *disr*, it is exactly the opposite. This is documented in the supplementary material. In contrast, employing the *mean* strategy, the CCSD results converge at nearly the same rate as DFT-PBE. Both DFT-PBE and CCSD-mean exhibit a smooth 1/*n* behavior. Note that periodic DFT-PBE implementations do not rely on two-electron, four-orbital Coulomb integrals and, therefore, do not suffer from the discussed sampling issues. Thus, the agreement is a strong indication that the *mean* strategy effectively eliminates the discretization errors that lead to problems in HF and CCSD calculations of finite anisotropic cells.

## V. SUMMARY AND CONCLUSION

In this manuscript, we have introduced the *mean* Coulomb potential as a sampling technique to enable many-electron ground-state calculations at the level of Hartree–Fock and coupled cluster theories in anisotropic Born–von Karman cells. By examining two illustrative cases, the infinite carbon chain and the (001) LiH surface, we identified and effectively resolved deficiencies inherent in existing sampling techniques for the reciprocal Coulomb potential. The presented implementation for the evaluation of the *mean* Coulomb potential is specifically designed for sub-meV accuracy of ground states and can be directly adopted in other codes employing the reciprocal form of the Coulomb potential. This advancement not only enables predictive *ab initio* calculations for low-dimensional systems using highly accurate methods, such as CCSD(T), but can also be readily applied to Born–von Karman cells of general shape. Furthermore, the *mean* potential technique can just as well be applied without further ado to modified reciprocal Coulomb potentials, such as range-separated potentials^{29,30} or truncated potentials,^{9,11,31} thereby both improving their sampling in reciprocal space and circumventing any special treatment of singularities.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the CCSD(T) data to calculate the potential energy surface of the carbon chain, a decomposition of the CCSD convergence of the (001) LiH surface slab energies into the HF and correlation contribution, and the atomic structure of the six-layer surface slab (POSCAR files).

## ACKNOWLEDGMENTS

T.S. acknowledges support from the Austrian Science Fund (FWF) [10.55776/ESP335]. The computational results presented have been largely achieved using the Vienna Scientific Cluster (VSC). The research presented here was funded in part by the National Science Foundation under Grant No. NSF CHE-2045046 (J.J.S. and W.Z.V.B.).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Tobias Schäfer**: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (lead). **William Z. Van Benschoten**: Conceptualization (equal); Investigation (equal). **James J. Shepherd**: Conceptualization (equal). **Andreas Grüneis**: Conceptualization (equal); Methodology (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

## REFERENCES

*Many-Body Methods in Chemistry and Physics*

*ab initio*and force-field-based calculations in clusters

*ab initio*tool for excited state calculations