A boundary integral method for computing electron repulsion integrals in heterogeneous dielectric environments

Introduction

The four-index electron repulsion integral (ERI) tensor is a fundamental computational primitive of electronic structure theory [1, 2], describing electron-electron interactions between single-particle orbitals, i.e., the wavefunctions of individual electrons, which are spatially localized in our setting. Given a set of $N_o$ real-valued orbitals $\varphi:= \{\varphi_i\}_{i=1}^{N_o}$, the rank-four ERI tensor $V := \{V_{ijkl}\}_{i,j,k,l=1}^{N_o}$ has entries

where $W(x,x')$ is a two-body interaction kernel, $\rho:= \{\rho_{ij}\}_{i,j=1}^{N_o}$ is a set of $N_o^2$ pair densities, each defined as

and $\phi:= \{\phi_{ij}\}_{i,j=1}^{N_o}$ are the potentials generated by each pair density through the kernel $W$, thus defined as

The ERI tensor appears in many-body electronic structure methods from Hartree–Fock [3] and hybrid density functional theory [4] to correlated wavefunction [5], quantum-embedding [6], and lattice model [7] approaches. When the orbitals are localized, the number of non-negligible pair densities is $O(N_o)$ rather than $N_o^2$, and the ERI tensor is correspondingly sparse, with $O(N_o^2)$ non-negligible entries. It can be compressed further by low-rank factorizations such as Cholesky decomposition [8], density fitting [9, 10], tensor hypercontraction [11, 12, 13, 14, 15], or the use of structured/adaptive basis sets [16, 17, 18, 19, 20].

In standard settings, $W(x,x') = \frac{1}{4\pi\lvert x-x'\rvert}$ is the Coulomb kernel, corresponding to electrons interacting in vacuum or a homogeneous medium, and the evaluation of the resulting integrals is routine and highly optimized. For Gaussian basis sets, for example, they admit analytic recurrences [21, 2]. In this work, we consider the case in which $W$ is the Green’s function for the Poisson equation in a spatially varying dielectric environment, rather than the bare Coulomb kernel. Such a kernel is not usually available in closed form, and can be applied only by solving an electrostatic interface boundary-value problem (BVP).

Spatially varying dielectric environments arise in many experimental platforms of contemporary interest. For example, atomically thin crystals such as graphene [22] and transition metal dichalcogenide monolayers [23, 24] are regularly exposed to oxide substrates, encapsulated in hexagonal boron nitride [25], and gated by nearby conducting layers. Because a two-dimensional (2D) material screens the in-plane Coulomb interaction only weakly on its own [26, 27, 28], the effective interaction between its electrons is determined to a large degree by the screening properties of its polarizable dielectric surroundings. This environmental screening can change quasiparticle band gaps and exciton binding energies by hundreds of meV [29, 30], and can be patterned laterally to define junctions within a single homogeneous monolayer, an approach known as Coulomb engineering [31, 30]. In moiré superlattices such as magic-angle twisted bilayer graphene [32], tuning the screening environment alone can even suppress correlated insulating phases in favor of superconducting ones [33, 34].

Capturing such an environment fully from first principles is infeasible, since the surrounding device is far too large to treat quantum-mechanically alongside the active material. Practical descriptions therefore retain only a small set of localized orbitals within the “active” part of the system, namely the part treated quantum-mechanically, such as Wannier functions of its low-energy bands, giving an extended Hubbard-type model whose interaction coefficients are precisely the ERIs [35]. Downfolding schemes such as the constrained random phase approximation [36] and quantum embedding theories for correlated defect states [37, 38] account for the screening from the retained electronic degrees of freedom, while the surrounding device structure must instead be described approximately at the level of continuum electrostatics [39, 40].

We model the dielectric environment as a collection of (possibly-adjacent) dielectric cuboids $\Omega:= \{\Omega_k\}_{k=1}^{N_k}$, where $\Omega_k$ has constant relative permittivity $\epsilon_k$, embedded in a vacuum background $\Omega_0 = \mathbb{R}^3 \setminus\bigcup_{k=1}^{N_k}\Omega_k$ with $\epsilon_0 = 1$. Figure 1 shows an example with $N_k = 3$ materials. This setup covers many device geometries of interest, and requires addressing aspect ratios and edge/corner singularities, while allowing certain other simplifications in the boundary integral equation (BIE) implementation. While the BIE method also applies to more general geometries with curved interfaces, for which various accurate quadratures could be used [41 §9.3] [42, 43, 44], the focus of the present paper is on several other challenges that arise in the ERI application. We assume that all retained pair densities $\rho$ are supported inside a single “active” region $\Omega_m$ (for example, $m = 3$ in Figure 1). In the intended application, the basis functions $\varphi$ are localized orbitals of this active material, obtained from an outer electronic structure calculation such as a Wannier representation of the low-energy bands derived from density functional theory [45, 46, 47]. Thus each pair density $\rho_{ij}$ is smooth and compactly supported inside the region. Dropping the $ij$ indices, a given pair density function $\rho$ induces the potential function $\phi$ which solves the Poisson BVP

where $\Gamma_{kk'} := \partial\Omega_k \cap\partial\Omega_{k'}$ is the interface between two adjacent dielectric regions $\Omega_k$ and $\Omega_{k'}$, and $n$ denotes its normal vector (oriented to point from $\Omega_k$ into $\Omega_{k'}$ when $k < k'$). The first two equations are the (nondimensionalized) Poisson equation. The third and fourth equations are continuity and flux-matching conditions. The last line imposes zero potential at infinity, ensuring uniqueness. Once $\phi$ has been obtained, the corresponding ERI entries can be evaluated by the inner product (3) with a “target” pair density.

Existing continuum treatments of dielectric screening often exploit analytically solvable models, such as effective dielectric constants, image-charge constructions for a single planar interface, or Rytova–Keldysh kernels for an infinite uniform slab [27, 39]. These descriptions are effective for idealized planar or otherwise separable geometries, but the finite, multi-region, laterally-structured configurations considered here demand a numerical solution of (6).

One standard approach is to discretize the Poisson equation throughout the volume, using, for example, interface-fitted finite differences or finite elements [48, 49], multigrid continuum-solvation solvers [50], or particle–mesh techniques [51, 52, 53]. For the present problem, however, such a volumetric solve faces three difficulties:

1. The localized orbital densities vary on a much finer length scale than the surrounding dielectric structure of the device: this ratio is often $\gtrsim10^3$. Any volume mesh must therefore resolve this range in each of three dimensions.

2. The dielectric geometry contains edges and corners at which the electric field is singular. Resolving these local features would require additional 3D refinement.

3. The problem is posed in free space with a decay condition at infinity; volumetric discretization thus requires artificial truncation to a sufficiently large domain, or specialized boundary conditions or elements for its accurate handling.

These limitations motivate a BIE approach. BIE solvers are widely used in computational electrostatics for piecewise-homogeneous media [54, 55], with applications including Poisson–Boltzmann models of solvated biomolecules [56, 57], polarizable continuum models in quantum chemistry [58, 59], and dielectric inclusions with edges and corners [60]. Because the permittivity is piecewise constant, the BVP (6) can be reformulated via Green’s theorem as a BIE for an induced surface charge density defined only on interfaces. The reduction in dimensionality brings a large reduction in the number of unknowns—especially those induced by multi-scale refinement—and in the cost of meshing, compared to finite elements. The free-space condition at infinity is already built into the integral representation. Furthermore, the BIE may be arranged to be Fredholm equation of the second kind, so that the condition number of the resulting linear system is typically small and fixed as the discretization is refined. Thus iterative solution is effective and, although the system matrix is dense, it may be applied at each iteration in linear cost by the fast multipole method (FMM) [61].

Although a BIE approach removes many of the difficulties of volumetric discretization, some challenges remain, and motivate the techniques developed in this work. First, the induced surface charge is singular at dielectric edges and corners, while a localized source close to an interface produces rapidly-varying boundary data. We handle both issues with adaptive panel refinement. Second, evaluating layer potentials between nearby panels, such as parallel facing panels in high-aspect geometries, requires corrected near-singular quadratures. Third, each source pair density (supplied on a quadrature grid of millions of points) generates a free-space incident potential that must be evaluated accurately both near and far from its source box. For this, we combine an FFT-based truncated-kernel method (TKM) [62] for near targets with FMM-accelerated quadrature for far targets. Finally, computing the ERI tensor requires many thousands of interface solves and potential evaluations on the fixed geometry. We therefore batch the solves and reuse geometry-, source-, and target-dependent work throughout a distributed computation pipeline.

While most of these ingredients are established, a major methodological contribution is their integration into a scalable workflow for the automated, efficient evaluation of many dielectric-screened ERIs in a given geometrically complex environment, to a user-prescribed accuracy. In particular, the workflow combines adaptive interface discretization tailored to localized pair-density sources and edge singularities, FMM-compatible corrections for close-to-touching panels, hybrid near/far incident-field evaluation, reuse of fixed-geometry data across many right-hand sides, and batched distributed evaluation of the retained ERI tensor.

An outline of the paper is as follows. In § 2, we present the BIE formulation of the Poisson BVP (6), including the locally corrected Nyström discretization, the adaptive panel refinement strategy, and the near-panel quadrature correction of the discrete operator. § 3 presents the hybrid TKM and FMM quadrature methods for evaluating the incident potential and its normal gradient. § 4 explains how these tools are used in a scalable fashion: batched interface solves over a fixed geometry, repeated potential evaluation on many target densities, and distributed computation of the sparse four-index ERI tensor. § 5 presents numerical results: convergence and dielectric-contrast studies on a three-material benchmark with edges, corners, and a buried material junction, followed by computation of the full retained ERI tensor for graphene orbital-pair densities on a $\mathrm{Si}/\mathrm{SiO}_{2}$ heterojunction substrate.

Boundary integral formulation and discretization

In this section, we reformulate the BVP (6) as a second-kind BIE posed on the dielectric interfaces, and describe its discretization. § 2.1 derives the integral equation. § 2.2 overviews the locally corrected Nyström discretization and the resulting linear system, assuming a given panel mesh. § 2.3 describes the adaptive construction of this panel mesh, which is guided by the orbital structure as well as the edge and corner singularities of the interface density. § 2.4 details the local quadrature corrections for nearby panel pairs.

Second-kind boundary integral formulation

We begin by splitting the physical potential in (4) into two parts (using language common in wave scattering theory):

The first term is called the incident potential, and is generated by the density $\rho$ in a homogeneous space with permittivity $\epsilon_m$,

so that it satisfies

Note that the first equation holds on all of $\mathbb{R}^3$, since $\rho$ vanishes outside $\Omega_m$. Moreover, since $\rho$ is smooth and compactly supported, $u_{\mathrm{inc}}$ is smooth throughout $\mathbb{R}^3$. Let $\Gamma=\bigcup_{1\leq k<k'\leq N_k}\Gamma_{kk'}$ be the set of all interfaces. (Note that interfaces between regions with equal permittivity are not dielectric interfaces and may be removed from $\Gamma$.)

Subtracting the conditions (7) satisfied by $u_{\mathrm{inc}}$ from the BVP (4) for the total potential $\phi$, and using the splitting (5), we see that the scattered potential $u$ solves a Laplace BVP with inhomogeneous transmission data:

Here $u^-$ and $u^+$ denote the one-sided limits from $\Omega_k$ and $\Omega_{k'}$, respectively, consistent with the normal orientation on $\Gamma_{kk'}$.

We introduce a single layer potential representation of the scattered potential, a standard approach for solving the Laplace transmission problem with a continuous potential:

Here $\sigma$ is an unknown single layer density defined on $\Gamma$, and $\mathrm{d}S$ is surface area measure. Any such representation is harmonic off $\Gamma$, continuous across $\Gamma$, and decays at infinity, so it automatically satisfies the first, second, and fourth conditions of (8). For the flux-matching condition, the third line of (8), we note that the single layer potential satisfies the standard jump relation

where $D^T$ is the on-surface adjoint double layer operator, defined via the improper integral

For mathematical background in the Sobolev setting $\sigma\in H^{-1/2}(\Gamma)$, sufficient for the integrable singularities arising with edges/corners of Lipschitz surfaces, we refer the reader to, e.g., [63 Ch. 3]. Following [55], substituting this jump relation into the flux-matching condition reduces it to a second-kind BIE,

where $\gamma$ is an inverse contrast function defined on $\Gamma$ by

obeying $|\gamma| > 1$ for bounded real permittivities. Solving (14) yields the single layer density $\sigma$, from which $u$ is generated by (11), and finally the desired total potential $\phi$ via (7).

Locally corrected Nyström discretization

We discretize the integral equation (14) using a panelized Nyström method. We assume that $\Gamma$ has been partitioned into a disjoint collection $\mathcal{P}$ of $M = |\mathcal{P}|$ flat rectangular panels, each living on one of the interfaces $\Gamma_{kk'}$. The adaptive construction of this panelization will be described in Section 2.3. Each panel $P \in\mathcal{P}$ uses a $p \times p$ tensor product Gauss–Legendre quadrature rule, with nodes mapped onto the panel and weights scaled by its area. Let

denote the mapped quadrature nodes, surface weights, and unit normals at each node on $P$. We denote the union of all such nodes, weights, and normals as $\{x_i,w_i,n_i\}_{i=1}^{N}$, where $N = Mp^2$. We abbreviate

for the nodal values of the single-layer density and right-hand side data on panel $P$, and collect them into vectors $\boldsymbol{\sigma}_P,\boldsymbol{f}_P \in\mathbb{R}^{p^2}$. Stacking these vectors panel-by-panel gives

The Nyström discretization of the integral equation (14) takes the form of a block linear system

where $A \in\mathbb{R}^{N\times N}$ is composed of square $p^2 \times p^2$ blocks $A_{PQ}$, each indexed by a “target” panel $P$ (block row) and a “source” panel $Q$ (block column). We make the decomposition

where the three terms are, respectively, the diagonal contrast term, the standard panelwise Gauss–Legendre discretization of the adjoint double layer operator, and the sparse local quadrature correction. These terms are defined blockwise as follows:

1. $\boldsymbol{G}$ is block diagonal, and discretizes the first (scaled identity) term of (14). The $P$th diagonal block is $G_{PP} = \frac{1}{2}\gamma_P I$, with $\gamma_P$ the contrast of the interface containing $P$.

2. $D^{T}$ discretizes the second term of (14) using standard Gauss–Legendre quadrature:

Because all interfaces are flat in our geometry, $D^{T}_{PQ}\equiv0$ whenever the source and target panels are coplanar; this includes panels on the same interface, and in particular diagonal blocks $D^{T}_{PP}$.

3. $\mathbf{C}^{T}$ is a block-sparse local correction: its block $\mathbf{C}^{T}_{PQ}$ is nonzero only when the panels $P$ and $Q$ are near each other, or touching at an edge or corner, in which case it corrects the standard block $\mathbf{D}^{T}_{PQ}$ using a target-specific quadrature

where $\mathbf{D}^{T,\mathrm{acc}}_{PQ}$ contains the accurate local quadrature weights, whose construction is given in § 2.4.

The standard rule (21) is high-order accurate (with respect to panel refinement) whenever the integrand is $C^{\infty}$ smooth on the source panel, which holds when $P$ and $Q$ are well-separated relative to the size of $P$. This is true for all but a small fraction of panel pairs. However, for close and touching pairs, which occur at edges, and for the large faces of thin slabs (as in Figure 4), the fixed $p \times p$ rule becomes inaccurate, and a new scheme is needed to preserve accuracy. This is the reason for introducing the correction term $\mathbf{C}^{T}$. This correction, and the classification of close and touching pairs, will be described in § 2.4.

With the system matrix defined, we solve the square linear system (19) iteratively with GMRES [64], using $\mathbf{G}^{-1}$ as a diagonal right preconditioner. The resulting preconditioned system is

so the matrix-vector multiplication at each iteration is with the matrix

$\mathbf{G}^{-1}$ is block-diagonal and can therefore be applied with $\mathcal{O}(N)$ cost. $\mathbf{D}^{T}$ is dense but involves source (charge) weights that are independent of the target point, thus can be applied to a vector at $\mathcal{O}(N)$ cost using a single FMM call, with the target-normal derivative of the potential as the output. The correction term $\mathbf{C}^{T}$ is block sparse and can therefore be applied by a direct sparse matrix–vector product, with cost proportional to the number of close and touching pairs. The latter cost is usually also linear in $N$ (since each panel has $\mathcal{O}(1)$ near neighbors), and much cheaper than an FMM call.

If $\Gamma$ were smooth, $\mathbf{D}^{T}$ would be compact, so the eigenvalues of $\mathbf{A}\mathbf{G}^{-1}$ would cluster only at $1$, and the GMRES convergence rate would be expected to be fixed under mesh refinement. Corners and edges add a non-compact contribution to $\mathbf{D}^{T}$. However, its norm is known to be less than one [65 §6.5], so we still only expect at most weak growth in iteration count. We verify this numerically in § 5.1.

Adaptive panel refinement

In this section, we describe the adaptive panelization $\mathcal{P}$ of the set of interfaces $\Gamma$, needed to discretize the integral equation (14). We fix the quadrature order $p$ and construct the panelization $\mathcal{P}$ using an $h$-adaptive refinement strategy, recursively subdividing each selected panel into four child panels. Refinement is driven by two independent considerations. First, the right-hand side $f=-\partial_n u_{\mathrm{inc}}$ of Eq. (14) can vary rapidly when the source density $\rho$ is close to the interface, and the layer density $\sigma$ inherits this rapid variation. Second, $\sigma$ develops integrable singularities near edges and corners of the dielectric interfaces. Our corresponding refinement criteria will be data-driven and geometry-driven, respectively.

Right-hand-side-driven refinement

The first component of the refinement criterion is based on the local interpolation error of $f$. On each panel $P$, let $\Pi_P f$ be the tensor-product polynomial interpolant of degree at most $p-1$ in each local coordinate, such that

We refine each panel until

given a prescribed local refinement threshold $\tau_{\mathrm{rhs}}$. We estimate the norm in Eq. (26) using an oversampled $3p \times3p$ tensor-product Gauss–Legendre rule.

This adaptive refinement turns out to be crucial in our setting, as we now illustrate. We study the absolute global interpolation $L^2$-error of the right-hand side,

for the test geometry of Figure 2 with data from a narrow Gaussian source $\rho$ (a simple pair density model) a typical distance above one face. The left picture shows that the refinement concentrates near the source, as expected. The right plot confirms that although both uniform and adaptive refinement exhibit the expected high-order convergence, the adaptive mesh needs $N$ about 20 times smaller than for the uniform mesh to reach the same error. This holds essentially independently of error level and $p$.

Edge-driven refinement

For a rectangular panel $P$ with side lengths $s_1(P)$ and $s_2(P)$, we use

as a measure of its size. Our geometric refinement strategy is simply to recursively subdivide into four equal panels any panel adjacent to an edge, until every such panel $P$ has $h_P\le h_{\min}$, where $h_{\min}$ is a global convergence parameter. This results in a dyadic grading with respect to distance from the edge (see Figure 4 inset), giving $O(h_{\min}^{-2})$ panels in total.

To predict the resulting solution error convergence rate with respect to $h_{\min}$, we note that a dielectric edge generically induces a weak, integrable power-law singularity in the solution density $\sigma$, matching the electric field singularity. For simplicity, we consider the case of a single right-angle dielectric edge (noting that the case of multi-junctions is similar [66 (10)] [67 §5.3]). The following is derived in Appendix A by combining classical wedge analysis with the jump relation.

We now use this for a heuristic prediction of the far-field solution error induced by a finite $h_{\min}$. It is known that the discretized Nyström error is of the same order as the quadrature error applied to the true solution [65 p. 227]. Although panel-panel local corrections are used, the panels touching an edge cannot accurately resolve the edge singularity. Thus we expect the quadrature error to be dominated by the contribution from the unresolved innermost strip within distance $0 < r < h_{\min}$ from the edge. By the above Proposition, the $L^1$-norm of this contribution is bounded by

for some unknown constant $C$, giving our convergence rate prediction for the edge-related error. A rigorous analysis of this would be technical and is beyond the scope of this paper, but the power law (31) is useful in practice.

To verify this, we place a unit point charge inside a dielectric cube of relative permittivity $\epsilon\in\{1.87, 4.08, 21.0\}$, surrounded by vacuum, so that corresponding singularity exponents are $\beta\in\{-0.1, -0.2, -0.3\}$. We fix the quadrature order $p=4$ and vary $h_{\min}$. By Gauss’s law, the integral of surface charge equals $1-1/\epsilon$, so we choose

as an error indicator. Figure 3(a) shows that the density blows up with the predicted power-law $O(d^{\beta})$ over the resolved distance range, where $d$ is distance to the edge. The singularity is stronger for larger dielectric contrast. Figure 3(b) shows that the predicted $\mathcal{O}(h_{\min}^{\beta+1})$ error convergence order is remarkably accurate in practice. We have not attempted to analyze the corner singularities, but expect their error contribution to be of lower order than the edges, and the above results suggest this.

Locally corrected quadrature for near-panel pairs

The adaptive mesh of § 2.3 is designed to resolve the layer density on the $p \times p$ nodes of each panel, but the plain quadrature rule (21) is inaccurate for close source and target panels. We correct the affected blocks $\mathbf{D}_{PQ}^{T}$ as follows (note that this does not change the number of unknowns nor their node locations). For an ordered pair with source panel $Q$ and target panel $P$, the contribution of $Q$ at a target collocation node $x_{\mu}^{P}$ is

To express this contribution directly in terms of the unknown nodal values $\{\sigma_{\nu}^{Q}\}_{\nu=1}^{p^{2}}$, we interpolate

where $\ell_{\nu}^{Q}$ are the product Lagrange polynomial basis functions [68 Ch. 5] for the $p \times p$ panel nodes. Substituting this into the panel integral defines replacement matrix elements

These integrals must now be approximated, and touching and non-touching near ordered pairs require different constructions.

For a touching pair, sources on $Q$ and targets on $P$ can approach each other closely. We therefore compute each integral in (34) by adaptive 2D quadrature over $Q$ (a quad-panel version of the triangle-panel neighbor correction strategy in [44]). The source panel $Q$ is recursively subdivided, and the same $p \times p$ tensor-product Gauss–Legendre rules are applied on its sub-panels until the estimated local quadrature error falls below a preset parameter $\tau_{\mathrm{quad}}$. We vectorize over the $\nu$ index for speed.

This adaptive scheme is accurate for any near pair, but expensive: each target point in each block requires its own recursive subdivision of the source panel. Thus, for near but non-touching pairs we apply a cheaper upsampling strategy exploiting the fact that the kernel is now analytic on $Q$. Let $Q_{\mathrm{up}}$ denote a $p_{\mathrm{up}} \times p_{\mathrm{up}}$ Gauss–Legendre upsampled grid on the same panel, where $p_{\mathrm{up}} > p$, and let $\mathbf{U}_{Q_{\mathrm{up}}\leftarrow Q} \in\mathbb{R}^{p_{\mathrm{up}}^{2}\times p^{2}}$ be the interpolation matrix from the original $p \times p$ grid to this new grid. Using quadrature on the upsampled grid gives

where $\mathbf{D}_{P,Q_{\mathrm{up}}}^{T} \in\mathbb{R}^{p^{2}\times p_{\mathrm{up}}^{2}}$ has entries as in (21) but instead using the upsampled nodes and weights. It remains to identify the non-touching ordered pairs that require correction, then to select $p_{\mathrm{up}}$ for each. For this, we use convergence theory involving the complex plane singularities of the integrand (following [69] and other recent work), building on analytic 1D Gauss–Legendre error estimates. For a fixed target point $x$ and a parametrization $x'(\boldsymbol{\xi})$ of $Q$, with $\boldsymbol{\xi}=(\xi_{1},\xi_{2})\in[-1,1]^{2}$, define the squared-distance function $R^{2}(\boldsymbol{\xi},x):=\lvert x'(\boldsymbol{\xi})-x\rvert^{2}$. For fixed $\xi_{2}$, the nearest complex roots in $\xi_{1}$ form a conjugate pair; denote the conjugate root in the upper half-plane by $t_{0}^{(1)}(\xi_{2};x)$. Define $t_{0}^{(2)}(\xi_{1};x)$ analogously by fixing $\xi_{1}$ and solving in $\xi_{2}$. For a complex singularity $t_0 \notin[-1,1]$, its Bernstein parameter is [69 Eq. 64]

This gives the largest Bernstein ellipse with foci at $\pm1$ within which the 1D integrand remains analytic. The associated $p$-point Gauss–Legendre error has the exponential rate [68 Thm. 19.3]

For the tensor-product rule, the smaller radius in the two reference directions should control the asymptotic rate. This motivates a heuristic “2D Bernstein” convergence rate (noting that more elaborate variants exist [69]),

As the target approaches $Q$, then $\varrho(x,Q)\to1^+$, making the convergence arbitrarily slow. Conversely, a larger $\varrho$ corresponds to a larger analyticity neighborhood and faster Gauss–Legendre convergence. Taking the worse-case Bernstein parameter over the target points in panel $P$ leads to

with expected quadrature error $\mathcal{O}\left(\varrho_{\min}(P,Q)^{-2p}\right)$ for a $p\times p$ rule. Equating this error with $\tau_{\mathrm{quad}}$, ignoring the prefactor, motivates our classifying a non-touching ordered pair $(P,Q)$ as near when

Note that since $\varrho_{\min}(P,Q)$ need not equal $\varrho_{\min}(Q,P)$, nearness is not a symmetric relation. For each near ordered pair $(P,Q)$ we then choose the upsampled order by inverting the asymptotic rate formula:

With $p_{\mathrm{up}}$ chosen, $\mathbf{D}_{PQ}^{T,\mathrm{acc}}$ is filled via (33), then $\mathbf{C}_{PQ}^{T}$ is filled via (20). For an efficient implementation, a $k$-d tree on panel centroids is used to generate spatially-nearby candidate pairs. For each ordered panel pair $(P,Q)$ and each target node $x_P^\mu$, we approximate the minimization in Eq. (38) over the $p$ Gauss–Legendre nodes in each direction on $Q$. We then apply Eq. (39) and the criterion in Eq. (40). Figure 4 highlights the resulting set of target panels $P$ that require near corrections for one large source panel $Q$ in a thin cuboid geometry.

We now verify the above error estimates and the resulting near quadrature correction. Figure 5(a) provides evidence that the above 2D Bernstein radius gives a reasonable prediction of the error using the standard (not upsampled) $p\times p$ rule. Their agreement over the plotted range shows that the Bernstein radius captures the observed spatial error landscape. Figure 5(b) validates the non-touching-panel upsampling strategy, for parallel panels separated by a range of gaps $d$:

We plot the relative $L^2$ error

Here $\mathbf{D}_{PQ}^{T,\circ}$ is either the standard $p=8$ block or the upsampled replacement, choosing $\tau_{\mathrm{quad}}=10^{-12}$. The uncorrected rule loses accuracy as $d$ shrinks, whereas the upsampled rule maintains accuracy throughout the tested non-touching regime. Over the same range, $p_{\mathrm{up}}$ increases from 16 at $d=1$ to 139 at $d=0.1$.

Evaluation of incident potential and its normal gradient

This section describes the evaluation of the incident potential Eq. (8) and its gradient, required respectively for the evaluation of the physical potential (5) and the right-hand side of the BIE (12). This is done most efficiently by separation into targets near to the pair density support, handled by the truncated-kernel method (TKM), and the remaining far targets, handled by particle-FMM accelerated quadrature.

Far-field quadrature evaluation and near/far splitting

We first describe a trapezoidal quadrature approximation for the incident potential Eq. (8) and its gradient at far-field targets. Let $\rho$ be a smooth source density whose support is contained in a finite cuboidal box $\mathcal{B}$. Choosing the center of $\mathcal{B}$ as the origin, we write

We sample $\rho$ on a uniform Cartesian grid in $\mathcal{B}$ with spacing $h$ in all three directions. The numbers of grid points in each direction are $N_x$, $N_y$, and $N_z$, respectively. The grid nodes are

We enumerate these nodes as $\{y_j\}_{j=1}^{N_\rho}$, where $N_\rho=N_xN_yN_z$. Let $\rho_h$ denote the trigonometric interpolant [71 §9] constructed from the grid samples of $\rho$ with spacing $h$ (since $\rho$ decays smoothly to zero, it may be treated as periodic). Then we choose the grid spacing $h$ small enough so that

where $\tau_\rho$ is a prescribed tolerance. This is equivalent to insuring that the tails of the Fourier transform of $\rho$ are small beyond the Nyquist frequency magnitude $\pi/h$.

For far-field targets, we simply apply the tensor-product trapezoidal rule, discretizing the volume integral (6) on $\{y_j\}_{j=1}^{N_\rho}$. All weights are $\omega_\rho=h^3$, giving

which together are evaluated efficiently at all targets $x$ with a single FMM call.

For targets far from the support of $\rho$, the integrands are smooth, so this trapezoidal rule converges superalgebraically [72 §10]. However, this convergence is not uniform as the target approaches the support of $\rho$, because the integrands are nearly singular and become underresolved on the fixed $h$ grid. To distinguish between far and near targets, we define a near-field region $B_{\mathrm{pad}}$ by padding $\mathcal{B}$ by $h_{\mathrm{pad}}=c_{\mathrm{pad}}h$ in each of the six directions. Motivated by the analogous trapezoidal-rule error analysis for layer potentials [73], which proves that the error decays to machine precision as the target moves several grid spacings from the source, we take $c_{\mathrm{pad}}=5$ as a practical choice. For the set of targets inside $B_{\mathrm{pad}}$, we instead apply the TKM to evaluate the volume integral.

Truncated kernel method for near-field volume integral evaluation

We introduce a Fourier-based truncated kernel method, inspired by Vico et al. [62], to evaluate the volume integral in (8) for targets inside $B_{\mathrm{pad}}$. By the convolution theorem, this volume integral has the Fourier representation

where $|k|^{-2}$ is the Fourier transform of the Coulomb kernel, and the source Fourier transform $\widehat{\rho}$ is

This representation avoids direct integration against the singular real-space kernel at $x'=x$. Since $\rho$ is smooth and compactly supported, $\widehat{\rho}$ decays rapidly so that the Fourier integral can be truncated to a bounded frequency domain, which is a cube of side length $2k_{\max}$ centered at the origin. A uniform-grid trapezoidal discretization then permits FFT-based evaluation. Two issues need to be handled: (48) has a singularity at $k=0$, and the aliasing error induced by Fourier-space sampling must be controlled.

To remove the Fourier-space singularity without changing the potential in $\mathcal{B}_{\mathrm{pad}}$, we truncate the Coulomb kernel in real space beyond the largest distance from a source $x^{\prime} \in\mathcal{B}$ to a target $x \in\mathcal{B}_{\mathrm{pad}}$,

The radially truncated Coulomb kernel is then defined by

We then define

Since for every $x^{\prime}\in\mathcal{B}$ and $x\in\mathcal{B}_{\mathrm{pad}}$, $G_L(x-x^{\prime})=G(x-x^{\prime})$, $u_{\mathrm{inc}}(x)=u_L(x)$ for all $x\in\mathcal{B}_{\mathrm{pad}}$. The Fourier transform of $G_L$ is given by [62]

which is now smooth. We discretize the Fourier integral (52) using a truncated trapezoidal rule on a Cartesian grid with spacing $\Delta k_\alpha$ in each direction (again here $a,b,c\in\mathbb{Z}$)

where $M_\alpha=\left\lceil\frac{k_{\max}}{\Delta k_\alpha}\right\rceil$ for $\alpha\in\{x,y,z\}$. We explain the optimal selection of the Fourier grid spacings $\Delta k_\alpha$ in the next section. We enumerate this Fourier grid as $\{k_\mu\}_{\mu=1}^{M_k}$, where $M_k=(2M_x+1)(2M_y+1)(2M_z+1)$. At each Fourier grid point, we approximate the Fourier transform of the source density using its existing grid (43):

which can be evaluated via a single padded 3D FFT. Using this discretized $\widehat{\rho}$ grid, the approximation to (52) gives the incident potential

Analytic differentiation of (56) gives the corresponding gradient approximation:

We use a 3D type-2 non-uniform fast Fourier transform (NUFFT) [74, 75, 76] to evaluate Eqs. (56) and (57) at a large, arbitrary set of target points simultaneously, at a cost that is quasi-linear in $M_k$ and number of target points. This is typically only one order of magnitude more expensive than an FFT. The inputs to the transform for the potential are the kernel-weighted Fourier coefficients $\widehat{u}=\{\widehat{u}(k_\mu)\}_{\mu=1}^{M_k}$, where

For the gradient components, the inputs are multiplied by the components of $ik_\mu$. The stack of four type-2 NUFFTs can be done more efficiently in a single call.

Remark 3.1. In practice, orbitals are often provided on regular but skew grids associated with the crystal lattice, so using a 3D FFT for (55) would need resampling to a compatible Cartesian grid. Since there are typically $10^{6}$ such points, this is cumbersome. Instead we work directly with $\{y_j\}$ the input set of skew grid points, and $\omega_\rho$ their quadrature weight, evaluating (55) using a single 3D type-1 NUFFT call. The speed difference vs an FFT is outweighed by the flexibility that comes from using NUFFTs for both input and target point sets.

Fourier-grid selection and error control

If $\Delta k_\alpha$ are chosen to obey certain bounds, then the aliasing error introduced by the discreteness of the Fourier quadrature in (56) and (57) vanishes entirely, leaving only the high-frequency truncation error (due to the finite $k_{\max}$).

The proof relies on the finite spatial support of $G_L$, adapting arguments of [62] to the padded cuboidal case. The 3D version of the Poisson summation formula [77 §7, Cor. 2.6] gives

where $T_n=(2\pi n_x/\Delta k_x,2\pi n_y/\Delta k_y,2\pi n_z/\Delta k_z)$, and the second term on the right-hand side is the error. Since $\rho$ is supported in $\mathcal{B}$ and $G_L$ is supported in the ball of radius $L$, the support of $u_L$ lives within $\prod_{\alpha\in\{x,y,z\}}[-l_\alpha/2-L,l_\alpha/2+L]$. Under condition (59), for all targets $x\in\mathcal{B}_{\mathrm{pad}}$ and $n\ne(0,0,0)$, each $x+T_n$ lies outside $u_L$’s support. Thus the second term on the right-hand side of Eq. (62) vanishes, and so does the error. Recall that $u_{\mathrm{inc}}(x)=u_L(x)$ for all $x\in\mathcal{B}_{\mathrm{pad}}$, so that Eq. (60) follows. Similar considerations apply to the gradient of $u_L$, leading to Eq. (61). Figure 6 illustrates the geometric separation that excludes contributions from periodic source images.

In practice, for efficiency we take equality in Eq. (59), giving the largest $\Delta k_\alpha$ permitted by Proposition 3.2, hence the coarsest Fourier lattice satisfying the sufficient condition for eliminating aliasing error on $\mathcal{B}_{\mathrm{pad}}$.

Efficient, parallel evaluation of the ERI tensor

In this section, we describe the overall procedure for evaluating the ERI tensor in Eq. (1) for a set of orbitals, which combines the BIE solver in Section 2 and the TKM/FMM solvers in Section 3.

Batched BIE solves

Evaluating the ERI tensor requires computing the potential generated by the orbital pair densities defined in Eq. (4). For localized orbitals, many pair densities are small and can be neglected. For a given threshold $\tau_{\mathrm{pair}}$, we define the set of active orbital pairs as

where $\lVert\rho_{ij} \rVert_2$ is the $L^2$ norm of the pair density $\rho_{ij}$. Each $(i,j) \in\mathcal{I}_{\mathrm{pair}}$ defines an active density $\rho_{ij}$ and a corresponding BIE to solve. We enumerate the active densities as $\{\rho_a\}_{a=1}^{N_p}$, where $N_p = \lvert\mathcal{I}_{\mathrm{pair}}\rvert$.

All active densities lead to BIEs with the same continuous operator but different right-hand sides and interface discretizations. Instead of solving each of them separately, we group spatially nearby densities into batches so that their right-hand sides can be efficiently resolved on a common mesh. The discrete operator $\mathbf{A}$ and FMM tree are shared across all source densities in a batch. In practice, $\mathcal{I}_{\mathrm{pair}}$ is divided into $N_b$ batches $\{\mathcal{B}_q\}_{q=1}^{N_b}$ with size $K_q \le K_{\max}$, where $K_{\max}$ is a prescribed maximum batch size.

The geometry-driven criterion of § 2.3.2 is source-independent, but the right-hand-side criterion of § 2.3.1 is not, since each right-hand-side $f_a$ requires refinement near its own source. We therefore build the common mesh for a batch $\{\rho_a\}_{a=1}^{K}$ using the following joint refinement criterion:

Here we recall from (25) that $\Pi_P f_a$ is the interpolant of $f_a$. This plays the role of (26), i.e., we subdivide each panel until the criterion holds for all panels $P$. Then the $K$ linear systems can be stacked as columns to give

with $f_a = -\partial_n u_{\mathrm{inc},a}$ sampled at the interface nodes. The preconditioner (23) can then be applied.

Eq. (65) is solved using block GMRES [78, 79, 80], which constructs a common block Krylov subspace for all $K$ systems, and is stopped at the selected residual tolerance $\tau_{\mathrm{GMRES}}$, as summarized in Algorithm 1. The FMM evaluates all $K$ strength vectors on the same source and target nodes, sharing the tree, interaction lists, and M2L matrix-vector multiplies, etc., across the batch. The vectorized FMM has given several-fold acceleration in BEM applications [81]. The user-selected $K_{\max}$ controls the tradeoff between the FMM cost per right-hand side and the orthogonalization work and memory required by a wider block Krylov basis. We assess this tradeoff numerically in Section 5.2.

Algorithm 1 Interface Poisson BIE solve for one source batch.

Input: pair-density batch $B=\{\rho_a\}_{a=1}^{K}$ on uniform volume grids; interface $\Gamma$; permittivities $\{\epsilon_k\}$; order $p$; thresholds $\tau_{\mathrm{rhs}},\tau_{\mathrm{quad}},\tau_{\mathrm{GMRES}}$ and $\tau_{\rho}$; Fourier/NUFFT/FMM tolerance; smallest panel size $h_{\min}$

Output: layer densities $\Sigma$ and the coefficient arrays $\{\widehat{u}_a\}$ of (58)

1. procedure SolveBatch$(B,\Gamma,\{\epsilon_k\})$

2. for all $a=1,\ldots,K$ do

3. form $\widehat{u}_a$ of (58) from the input samples of $\rho_a$, using a type-1 NUFFT ▷ Sec. 3

4. end for

5. partition $\Gamma$ into flat panels carrying $p\times p$ Gauss–Legendre rules; refine until the batch criterion (64) holds, then recursively subdivide near edges and corners until those panels have size $h_{\min}$ ▷ Sec. 2.3

6. construct $A=G+D^{\mathsf{T}}+C^{\mathsf{T}}$ of (20) ▷ Secs. 2.2, 2.4

7. form RHS $F$ and solve $A\Sigma=F$ by block GMRES, applying $A$ via FMM plus the block-sparse correction

8. return $\Sigma,\{\widehat{u}_a\}$

9. end procedure

Post-evaluation of the potential

For each active source density $\rho_a$, once the corresponding BIE has been solved, we evaluate the physical potential $\phi_a=u_{\mathrm{inc},a}+u_a$ on the volume grids (43) of all active target densities $\{\rho_b\}_{b=1}^{N_p}$. For each target density $\rho_b$, we then compute the corresponding ERI by the trapezoidal rule:

where $\{y_m^{(b)}\}_{m=1}^{N_{\rho_b}}$ are the volume grid nodes of $\rho_b$.

More specifically, to compute $\{V_{ab}\}_{b=1}^{N_p}$ for a given $a$, we consider the union of all target grid points:

The incident potential $u_{\mathrm{inc},a}$ defined in (8) is then evaluated on $\mathcal{T}$ by the near/far procedure of Section 3. The scattered potential $u_a$ is evaluated from the discretization of the single-layer representation (9), using the density column vector $\boldsymbol{\sigma}_a$ from the block linear solve (63).

Since the target set $\mathcal{T}$ can be close to the interface, special quadrature for near interactions is also needed. For a panel $P$ and a target point $y$, the same criterion as (40) is used to determine whether $y$ is near to or far from $P$. Far-field contributions are evaluated by FMM-accelerated quadrature using the existing interface nodes and weights. Integrals over near panels are evaluated by adaptive quadrature using the known density $\boldsymbol{\sigma}_a$ to the desired accuracy, as a correction to the far-field contributions. Unlike the on-surface correction, these near-field contributions are computed on-the-fly without forming or storing a correction matrix. The resulting evaluation procedure is summarized in Algorithm 2.

Algorithm 2 Evaluation of one source batch on the merged target set.

Input: batch $B$ with global pair-index set $I_B$ and its solve panels and volume-lattice samples; layer densities $\Sigma$ and coefficient arrays $\{\widehat{u}_a\}$ from Algorithm 1; merged set of all targets $\mathcal{T}$.

Output: $\{V_{ab}\}_{b=1}^{N_p}$, $a \in I_B$

1: procedure EvalBatch$(B, \Sigma, \{\widehat{u}_a\}, \mathcal{T})$ 2: for all $a \in I_B$ do 3: evaluate $u_{\mathrm{inc},a}$ on $\mathcal{T}$ from $\widehat{u}_a$ and the samples of $\rho_a$ by the near–far split of Section 3 4: evaluate $u_a$ on $\mathcal{T}$ from the discretized Eq. (11), with near panels selected by Eq. (40) 5: $\phi_a \leftarrow u_{\mathrm{inc},a} + u_a$ on $\mathcal{T}$ 6: for $b = 1, \ldots, N_p$ do 7: gather $\phi_a$ at $\{y_m^{(b)}\}$ and form $V_{ab}$ by Eq. (66) 8: end for 9: end for 10: return $\{V_{ab}\}_{b=1}^{N_p}$, $a \in I_B$ 11: end procedure

Distributed pipeline and computational cost

The overall procedure for evaluating the ERI tensor is summarized as a three-stage pipeline:

1. Prepare. Sample and batch the pair densities, and construct target grids.

2. Solve. Construct the common discretization for each batch and solve Eq. (65).

3. Evaluate. Evaluate each source potential to form the entries $V_{ab}$ with $b \ge a$.

The independent solve and evaluation stages are distributed across workers; the full procedure is summarized in Algorithm 3.

Algorithm 3 Distributed computation of the four-index tensor $\mathbf{V}$.

Input: orbitals $\varphi$; dielectric regions $\Omega$; pair-retention tolerance $\tau_{\mathrm{pair}}$; numerical controls $\tau_{\mathrm{rhs}}, \tau_{\mathrm{quad}}, \tau_{\mathrm{GMRES}}, \tau_{\rho}$, and the Fourier/NUFFT/FMM tolerance; maximum batch size $K_{\max}$; smallest panel size $h_{\min}$; quadrature order $p$

Output: non-negligible elements of the four-index tensor $\mathbf{V}$ in sparse form

1: set $\Gamma$; sample the pair densities on volume grids Eq. (45); select the pair densities using Eq. (63); form the $N_b$ localized batches of size at most $K_{\max}$; set $\mathcal{T}$ by Eq. (67) ▷ serial 2: for all batches $B$ in parallel do ▷ distributed across workers 3: $\Sigma_B, \{\widehat{u}_a\} \leftarrow$ SolveBatch$(B, \Gamma, \{\epsilon_k\})$ ▷ Alg. 1 4: end for 5: for all batches $B$ in parallel do ▷ distributed across workers 6: $\{V_{ab}\} \leftarrow$ EvalBatch$(B, \Sigma_B, \{\widehat{u}_a\}, \mathcal{T})$ ▷ Alg. 2 7: end for 8: return $\mathbf{V}$

We estimate the main computational costs of Algorithm 3. Let $N_p$ denote the number of active pair densities, $N$ a typical number of surface Nyström nodes, $N_\rho$ a typical number of volume-grid points per pair density, and $N_{\mathcal{T}} = |\mathcal{T}|$ the number of distinct target points. For localized orbitals, $N_p = O(N_o)$. The estimates below describe total computational work at fixed $p$ and tolerances, assuming that the Fourier mode count $M_k = O(N_\rho)$ and that local quadrature corrections only require linear work. The solve and evaluation stages account for most of the measured runtime.

For the solve stage, evaluating the right-hand sides at the interface nodes costs

where the first term accounts for the FMM-accelerated far-field in Eq. (47), and the second term accounts for the NUFFT-based near-field evaluation in Eq. (57). This is quasilinear in the total number of grid points $N_{\mathcal{T}} = O(N_pN_\rho)$. The cost of BIE solves is dominated by applying $\mathbf{A}$ in GMRES iterations, given by

where $N_{\mathrm{iter}}$ is a representative number of block GMRES iterations, which is studied in § 5.1. For the evaluation stage, the scattered potentials are evaluated on $\mathcal{T}$, while the incident potentials are evaluated by the near/far procedure in § 3, giving the cost

The first term is $O(N_p^2N_\rho)$, quadratic in the number of active pair densities (but still quasilinear with respect to number of grid points per pair). This contrasts the previous stages, which were all $O(N_p)$.

Numerical Results

In this section, we demonstrate the convergence, accuracy, timing and memory usage of the scheme. We also show the method’s performance on a large-scale, realistic system with complex dielectric geometry. Our implementation is in Julia [82], and is available at [83].

Convergence and effect of dielectric contrast

We first study the convergence with respect to the quadrature order $p$ and the minimum panel size $h_{\min}$ for edge refinement at fixed dielectric contrast. We then hold $p$ and $h_{\min}$ fixed and study the effect of increasing dielectric contrast on accuracy and GMRES iteration count, also showcasing the effect of the diagonal right preconditioner of Eq. (23).

A simple three-material geometry is used, as shown in Figure 7(a): a $10 \times10 \times1$ dielectric slab ($\Omega_3$) rests on two $10 \times10 \times10$ cubes sharing a vertical face ($\Omega_1$ and $\Omega_2$), all embedded in vacuum. A Gaussian density of width 0.05 at the center of the slab is used as the source density $\rho$, which is represented by the red dot in Figure 7(a), We compute $\rho$’s self-interaction $U_0$ as the quantity of interest for convergence studies, given by

where $\mathcal{B}$ is the bounding box of $\rho$, and $W$ is the Green’s function of the three-material system. For the convergence study, we set $\epsilon_1 = 4$, $\epsilon_2 = 12$, and $\epsilon_3 = 10$, and decrease $h_{\min}$ from $1/2$ to $1/32$ for each of $p \in\{2, 4, 6\}$. Other numerical tolerances are all set to $10^{-4}$, including the RHS-refinement, near-quadrature, GMRES, and volume-solver tolerances. The relative error

is reported against the number of unknowns $N$. Here $U_{\mathrm{ref}}$ is from a reference calculation with $p = 8$, $h_{\min} = 1/64$ and other tolerances set to $10^{-6}$. The corresponding errors and GMRES iteration counts $N_{\mathrm{iter}}$ are shown as solid curves in Figure 7(b) and (c), respectively. At comparable $N$, the error curves for different $p$ are almost identical, with convergence rate approximately $O(N^{-0.7})$, suggesting that the error is dominated by resolving edge singularities. Indeed, the analysis of § 2.3.2 predicts $O(h_{\min}^{\beta+1}) = O(N^{-(\beta+1)})$, where $\beta+ 1 \approx0.72$ using the dominant 12:1 dielectric ratio. Meanwhile, $N_{\mathrm{iter}}$ remains between 9 and 13, growing only mildly under mesh refinement. Dashed curves show results obtained without the touching-panel correction for $p = 2$ and $p = 4$. At $p = 2$, the error and $N_{\mathrm{iter}}$ are similar to the corrected case, but at $p = 4$, the error is about an order of magnitude larger, and $N_{\mathrm{iter}}$ is larger. Although $p = 2$ yields comparable accuracy without the touching-panel correction, a small $p$ requires more extensive near-quadrature corrections, substantially increasing the computational cost. This motivates using a larger $p$ together with touching-panel corrections to reduce errors associated with edge singularities.

We next isolate the effect of dielectric contrast, fixing $p = 6$, $h_{\min} = 1/32$, $\epsilon_1 = 4$, and $\epsilon_3 = 10$, and varying $\epsilon_2 \in\{6, 20, 60, 200\}$. As shown in Table 1, the relative error $\mathcal{E}_r$ increases mildly with $\epsilon_2$, from about $10^{-4}$ at $\epsilon_2 = 6$ to $10^{-3}$ at $\epsilon_2 = 200$. The preconditioned GMRES iteration count also increases, from 11 to 19. Thus increasing dielectric contrast increases the error and iteration count mildly at this fixed resolution. We also see that the diagonal right preconditioner in Eq. (23) reduces the iteration count by 39%–56%, with the largest reduction at the highest contrast.

$\epsilon_1$	$\epsilon_2$	$\epsilon_3$	\multicolumn{2}{c}{N_{\mathrm{iter}}}	$\mathcal{E}_r$
			unpreconditioned	preconditioned
4	6	10	18	11	$1.1 \times10^{-4}$
4	20	10	28	13	$1.7 \times10^{-4}$
4	60	10	36	17	$3.5 \times10^{-4}$
4	200	10	43	19	$1.0 \times10^{-3}$

Performance of the batched interface solve

We benchmark the batched interface solve described in § 4.1, focusing on strong scaling across CPU cores and the reduction in cost per source density obtained by solving $K$ right-hand sides simultaneously. We report single-node timings for the surface discretization, BIE solve, and evaluate steps of Algorithms 1 and 2. The discretization step includes the mesh discretization and the construction of the discrete operator $\mathbf{A}$, and the BIE solve refers to the block GMRES solve of Eq. (65). Timings for the full distributed pipeline are presented in § 5.3.

We consider a three-material system modeling a trilayer of graphene placed above a Si/SiO$_2$ heterojunction, as shown in Figure 8. The heterojunction consists of two $270 \times270 \times270$ Å cubes with relative permittivities of 11.9 for Si and 3.9 for SiO$_2$. The graphene has thickness corresponding to three graphene layers, and is represented by a $43 \times44 \times9$ Å slab with effective relative permittivity $\epsilon_{\mathrm{slab}} = 2.4$, the r\mathrm{P}A value for graphene [35].

We focus orbitals in the central graphene layer, which contains $N_o = 198$ Wannier $p_z$ orbitals (see Appendix § 8) centered on carbon atoms, as shown in the right panel of Figure 8. Each orbital density is supplied on a Wannier90 [84] real-space skew grid covering a $5 \times5 \times1$ graphene supercell with dimensions $12.325 \times12.325 \times15$ Å. The grid contains $150 \times150 \times192 = 4.32 \times10^{6}$ points, with spacings of about $0.08$ Å. From those orbitals, we form the $N_p = 2613$ pair densities that are non-negligible according to the criterion of Eq. (63), with $\tau_{\mathrm{pair}} = 10^{-4}$. This threshold corresponds approximately to only keep the pairs whose centers lie within $5$ Å of one another. Unless otherwise stated, all numerical tolerances are set to $10^{-5}$. We use quadrature order $p = 6$ and refine the interface mesh, and the edge panels are refined to a minimum size of $h_{\min} = 1.14$ Å. All benchmarks are performed on a single compute node equipped with two 48-core AMD EPYC 9474F processors and 1.5 TB of RAM.

We first examine the strong scaling of the method for a single source density ($K = 1$). We select the orbital closest to the center of the slab (see Figure 8), denoted by $\varphi_1$, and use $\rho_{11} = \varphi_1^2$ as the source density. We solve the corresponding BIE, whose discretization requires $N \approx1.35 \times10^{6}$ interface unknowns. We then evaluate the resulting potential on the merged target set $\mathcal{T}$, which contains $7.17 \times10^{6}$ points, and compute the 2613 tensor elements $V_{ij11}$, where $(i,j) \in\mathcal{I}_{\mathrm{pair}}$ (Eq. (63)). Keeping all other parameters fixed, we vary the number of threads from 1 to 96.

The results are shown in Table 2. The total wall-clock time decreases from 1020 s with 1 thread to 62 s with 96 threads, corresponding to a speedup of $16.4\times$ and a parallel efficiency of 17%. The decrease in parallel efficiency is consistent with memory-bandwidth saturation during the FMM-dominated iterative solve and potential evaluation.

$N_{\mathrm{threads}}$	discretization (s)	BIE solve (s)	eval (s)	total (s)	speedup	parallel efficiency
1	133.6	759	127	1020	$1.0\times$	100%
4	40.2	381	47	469	$2.2\times$	54%
16	14.4	146	22	183	$5.6\times$	35%
64	8.7	50	15	74	$13.9\times$	22%
96	8.9	39	14	62	$16.4\times$	17%

We then measure the benefit of batching multiple right-hand sides in the BIE solve and potential evaluation. We fix the number of threads at 96 and vary the batch size $K$ from 1 to 79, by including more neighboring pair densities. The results are shown in Table 3. The discretization time only increases slightly with $K$, from 10 s to 16 s. Both the BIE solve and evaluation phases benefit from the batched FMM calls. The end-to-end time per right-hand side decreases from 76 s at $K = 1$ to approximately 25 s at $K = 37$, and remains nearly constant through $K = 58$, corresponding to a 2.5× increase in throughput. Larger batches provide no further improvement, with the speedup declining slightly to 2.8× at $K = 79$. Meanwhile, peak memory usage grows approximately linearly with $K$ from 24 GB to about 1 TB. Thus, intermediate batch sizes, in this case specifically $K$ from 37 to 58, achieve nearly the highest measured throughput, with substantially lower memory usage than larger batches.

$K$	$N_{\mathrm{iter}}$	discretization (s)	BIE solve (s)	eval (s)	total (s)	peak RAM (GB)	speedup
1	11	8.9	39	14	62	24	–
4	10	11.0	88	40	140	62	1.8×
19	10	10.8	317	180	507	263	2.3×
37	9	12.0	552	371	934	512	2.5×
58	9	13.5	846	593	1453	783	2.5×
79	9	15.8	1258	829	2104	1064	2.3×

Large-scale ERI tensor computation

We apply the complete ERI computation pipeline, Algorithm 3, to the model introduced in § 5.2. We first show the wall-clock time for the distributed computation of the four-index tensor $\mathbf{V}$, and then show some representative entries of $\mathbf{V}$ and their convergence properties. We use the same numerical parameters: $p = 6$, $h_{\min} = 1.14$ Å and the RHS-refinement, near-quadrature, and GMRES tolerances all set to $10^{-5}$. The $N_o = 198$ orbitals in the active graphene layer yield $N_p = 2613$ pair densities. From Table 3, we select the batch size to be $K = 46$ for the distributed computation, which gives $N_b = 57$ batches in total. Those batches are distributed across ten 96-core nodes of the Rusty cluster at the Flatiron Institute, and processed in parallel.

Table 4 shows the per-stage wall-clock time for the distributed computation of $\mathbf{V}$. The calculation takes 2 h 10 m, or approximately 22 node-hours (about 2100 core-hours) in total. The solve and evaluation stages account for 50% and 33% of the wall-clock time, respectively. We see that even with $N_p = 2613$ pairs, the evaluation stage is not yet dominant over the other $O(N_p)$-scaling stages despite its $O(N_p^2)$ scaling.

stage	wall-clock
prepare	18 m
solve	65 m
eval	47 m
total	130 m

Figure 9(a) shows the computed ERI tensor element magnitudes $\lvert V_{ijkl}\rvert$, where $(i,j)$ and $(k,l)$ enumerate the $N_p$ orbital pairs. Because the retained pairs are ordered according to spatial locality, the matrix shows broad concentrations of relatively large entries together with repeated horizontal and vertical structure. To estimate the effect of discretization error at the production resolution, we perform a self-convergence study of $U_1 := V_{1111}$. Holding all other parameters fixed, we decrease $h_{\min}$ to refine the interface mesh near edges and corners, increasing the number of interface unknowns from $6.5 \times10^5$ to $2.2 \times10^7$. Figure 9(b) shows that the relative self-convergence error decays nearly as $N^{-0.8}$ over this range. For $h_{\min} = 1.14$ Å and $N \approx1.35 \times10^6$, the estimated relative error is approximately $4 \times10^{-5}$ and the absolute error is about 0.3 meV.

We next examine how the dielectric environment modifies the on-site interactions $U_i := V_{iiii}$ and inter-site density–density interactions $U_{ij} := V_{iijj}$. Figure 10 plots $\Delta U := U - U^{\mathrm{bare}}/\epsilon_{\mathrm{slab}}$, where $U^{\mathrm{bare}}$ is the bare Coulomb integral for the same orbitals. In panel (a), the on-site and first- to third-nearest-neighbor interactions lie approximately 0.19–0.20 eV below this reference on the Si side and increase by approximately 0.22–0.25 eV across the junction, becoming slightly positive on the SiO$_2$ side. Panel (b) compares $\Delta U_{ij}$ as function of orbital distance $r_{ij}$ for source orbitals $\varphi_i$ located in the Si region, near the junction, and in the SiO$_2$ region. At short range, screening is strongest on the Si side, intermediate near the junction, and weakest on the SiO$_2$ side. The results show that the Coulomb interactions are modified by the dielectric environment in a position-dependent manner, which explains the need for a full ERI tensor calculation in realistic device simulations.

Conclusion

We have presented a boundary integral method and a scalable HPC implementation for computing the four-index electron repulsion integral tensor in piecewise-constant dielectric environments composed of cuboidal regions. The method reformulates the governing Poisson interface problem as a second-kind boundary integral equation, discretized by a quad-panel-based Nyström scheme, with $h$-adaptivity driven by both the incident data and the edge singularities. To handle high aspect ratio cuboids efficiently, a Bernstein-radius criterion is used to engage special quadrature corrections between relatively close or touching panels. The incident fields of the smooth orbital-pair densities are evaluated by a truncated-kernel NUFFT method inside or near source supports, and quadrature farther away. All far-field interactions are accelerated by a vectorized FMM, with block GMRES as the iterative solver, showing significant amortization of costs over multiple right-hand sides.

Our extensive numerical studies included: i) convergence studies, showing rates matching an analytic density edge-singularity analysis, ii) validation that GMRES iteration counts remain moderate under refinement or large dielectric contrasts, iii) a strong scaling study with respect to thread count, and iv) parallel timings for the full ERI tensor (6.9 million outputs) on $10^3$ cluster cores, given 198 orbitals embedded in a realistic three-dielectric structure. In this final set-up, the ratio of system length to orbital grid-spacing is about $10^4$.

Much of our focus was on the scalable many-right-hand-side workflow required for ERI construction. While experiments confirmed that edge singularities induce a quite low convergence order with respect to $N$ (lower than the expected $O(N^{p/2})$ for uniform panels on a smooth surface), an important conclusion is that 3–4 digit accuracy is achieved with $N \sim10^6$, reasonable solution times, and $h_{\min} \approx1$Å. This value of $h_{\min}$ has reached the atomic scale at which a sharp-edge model becomes physically unimportant. However, we plan to develop edge-singularity-capturing panel quadratures, which will allow us to increase the convergence order, and reduce $N$. We also plan to extend the scheme to more general polyhedral or curved interfaces using existing BIE surface discretization tools.

For more general physics, possible future directions include periodic boundary conditions for extended multilayered materials, and the integration of the solver into self-consistent many-body frameworks, for which the screened interaction tensor computed here enters as a fundamental building block.

Density singularity exponent for a right-angle dielectric wedge

Here we derive the result on the density singularity in Proposition 2.1, by combining the jump relation with the right-angle special case of classical dielectric-wedge analysis [66, 85] (see also [67 §5.3]). The generic singularity results from assuming invariance along the edge, leaving a transverse 2D sectoral problem that is separable in polar coordinates $(r,\theta)$, with the edge at the origin and the interfaces at $\theta= \pm\pi/4$. The physical potential $\phi$ satisfies the 2D Laplace equation in each material,

Substituting the separated-variable form $\phi(r,\theta) = R(r)\psi(\theta)$ with separation constant $g$ gives

For $g > 0$, the bounded, nonconstant radial branch is $R(r) = r^g$. Thus in each material sector $\phi(r,\theta) = Ar^g \cos(g\theta+ C)$ for some coefficients $A,C$ that may differ between sectors. To find eigenvalues $g > 0$, we note that across either interface, $\phi$ and $\epsilon_0\partial_n\phi$ are continuous. Solutions fall into even and odd symmetry about $\theta= 0$. An even mode can be written

Setting $a = \pi g/4$, the two interface conditions at $\theta= \pi/4$ imply

After using addition formulae, nontrivial solutions demand that

The sin factor vanishing gives only regular modes, while the other factor vanishing gives

recalling (15), where the arccos is taken on the standard branch $[0,\pi]$. Similarly, the smallest positive exponent for a odd mode is $g_o = \frac{2}{\pi}\arccos\left(-\frac{1}{2\gamma}\right)$. The smallest (leading) positive wedge exponent is therefore the smaller of the two,

as in [66].

Since $u_{\mathrm{inc}}$ is smooth, (5) means that the scattered potential $u$ has the same singularity as $\phi$. Since the normal derivative is $\partial_n = r^{-1}\partial_\theta$, then the radial dependence of the jump in $\partial_n u$ is $\mathcal{O}(r^{g_*-1})$. By the jump relation for the single-layer representation (9), this jump determines $\sigma$, giving Proposition 2.1 with $\beta= g_* - 1$. The condition $\beta> -1/3$ follows from $|\gamma| > 1$.

If smooth potential variation along the edge is allowed, one still gets the same singularity as the above 2D analysis, since a Fourier transform along the edge direction gives a 2D modified Helmholtz equation in each material. Then $R(r) = I_g(r)$, the regular modified Bessel function of real order $g$, whose asymptotic power law is the same.

Wannier function generation

The Wannier functions were obtained from density functional theory (DFT) calculations performed with Quantum ESPRESSO (v7.4) [86, 87] using the PBE exchange–correlation functional [88] and Grimme’s DFT-D3 dispersion correction [89]. Carbon was described by an optimized norm-conserving Vanderbilt pseudopotential from the SG15 library (C_ONCV_PBE-1.2) [90, 91], with a plane-wave kinetic-energy cutoff of 100 Ry. We modeled monolayer graphene in a hexagonal unit cell with lattice constant $a = 2.465$ Å, and periodic images were separated by 15 Å of vacuum. The self-consistent ground state was converged to $10^{-10}$ Ry on a $\Gamma$-centered $32 \times32 \times1$ Monkhorst–Pack $k$-mesh, with 0.05 eV (3.7 mRy) Gaussian smearing. From these we built two maximally localized Wannier functions [46, 47] with C $p_z$ initial projections using Wannier90 (v3.1.0) [84]. The inner (frozen) energy window was set to $E_F \pm3$ eV. The Wannier functions were plotted on a $5 \times5 \times1$ real-space supercell for visualization.