arXiv:2610.26745v1 Announce Type: new Abstract: We present a boundary integral method for evaluating electron repulsion integrals (ERIs) in dielectric environments of piecewise constant permittivity. This generalizes conventional bare-Coulomb ERIs by incorporating screening from structured dielectric surroundings, such as the substrates and/or encapsulating layers present in modern two-dimensional-material devices. We formulate the problem as a second-kind boundary integral equation for the induced surface charge on the dielectric interfaces, and discretize it with a Nystr\"om method on adaptively refined panels. The incident potential generated by each smooth orbital source density is evaluated by a Fourier truncated-kernel method in or near to its source box, and by fast multipole method (FMM) accelerated quadrature farther away. The linear system is solved by a Krylov iteration also accelerated by the FMM. We focus on cuboidal device geometries, and retain efficiency at high aspect ratios via close-to-touching panel quadratures. We further develop a distributed pipeline that computes the full ERI tensor, comprising thousands of orbital pairs, embedded in a single dielectric geometry. We demonstrate the method on an effective three-layer graphene slab over a Si/SiO$_2$ heterojunction, focusing on the central graphene layer represented by $198$ $p_z$-like orbitals, retaining the $2613^2$ non-negligible entries of the $198^4$-sized ERI tensor. The model is resolved on $1.3$ million surface unknowns, achieving a relative error of $0.004\%$, at a cost of $22$ node-hours.