# Fast and spectrally accurate construction of
adaptive diagonal basis sets for electronic structure

Michael Lindsey<EMAIL_ADDRESS>Department of Mathematics,
University of California, Berkeley, Berkeley, CA 94720, USA  Sandeep Sharma
<EMAIL_ADDRESS>Department of Chemistry, University of Colorado, Boulder, CO
80302, USA

###### Abstract

In this article, we combine the periodic sinc basis set with a curvilinear
coordinate system for electronic structure calculations. This extension allows
for variable resolution across the computational domain, with higher
resolution close to the nuclei and lower resolution in the inter-atomic
regions. We address two key challenges that arise while using basis sets
obtained by such a coordinate transformation. First, we use pseudospectral
methods to evaluate the integrals needed to construct the Hamiltonian in this
basis. Second, we demonstrate how to construct an appropriate coordinate
transformation by solving the Monge-Ampère equation using a new approach that
we call the cyclic Knothe-Rosenblatt flow. The solution of both of these
challenges enables mean-field calculations at a cost that is log-linear in the
number of basis functions. We demonstrate that our method approaches the
complete basis set limit faster than basis sets with uniform resolution. We
also emphasize how these basis sets satisfy the diagonal approximation, which
is shown to be a consequence of the pseudospectral method. The diagonal
approximation is highly desirable for the solution of the electronic structure
problem in many frameworks, including mean field theories, tensor network
methods, quantum computing, and quantum Monte Carlo.

##  I Introduction

The first step in the solution of the electronic structure problem involves
the discretization of the Hamiltonian using a basis set. Two types of basis
sets are most commonly used: first, atom-centered basis functions (with
Gaussian type orbitals (GTO) being the most common in this category(molpro, ;
molcas, ; g16, ; orca, ; dalton, ; dirac, ; nwchem, ; gamess, ; psi4, ; pyscf,
; cp2k, ; turbomole, ; crystal, ), although numerical atomic orbitals are also
used(fhi-aims, )) and, second, plane wave basis sets(quantum_espresso, ; vasp,
). These basis sets have complementary strengths. With GTOs one is able to
perform all-electron calculations because they are well-suited to describing
the nuclear cusp in the wavefunction. They usually deliver impressive accuracy
even when very few basis functions are used, i.e., roughly $O(10)$ per atom,
although converging results to the complete basis set limit with error smaller
than $\mu$Eh can be challenging. Plane wave (PW) basis functions, on the other
hand, usually require the use of a pseudopotential (PP) as well as a much
larger number of basis functions, roughly $O(1000)$ per atom when PPs are
used. However, plane waves enjoy spectral convergence to the exact result
(within the PP approximation), and the complete basis set limit can be
achieved routinely. In addition, the various integrals in the Hamiltonian take
on a particularly simple form with PW. Moreover, the use of fast Fourier
transforms (FFTs) enables mean-field calculations with a computational cost
that scales better than mean-field calculations with GTOs by a factor
proportional to the system size, without further approximations such as tensor
hypercontraction (THC)THC1 ; THC2 . The disadvantage of the PW basis is that
it is not adapted to the molecular geometry. As a result, PW calculations
require uniformly high resolution across the entire computational domain,
unlike GTOs which provide high resolution close to nuclear centers and lower
resolution in the interatomic space.

More recently, other constructions including multiwavelet(madness, ; bigdft, ;
mrchem, ), finite differenceoctopus ; gpaw , and finite element
methods(GaviniEnriched1, ; rufus2021fast, ; dft-fe-1.0, ; lin_dg, ) have risen
to prominence, offering variable resolution over different parts of the
computational domain and enabling all-electron calculations with guaranteed
accuracy. Some of the most accurate, benchmark-quality all-electron
calculations for large system sizes have been performed using these basis
sets. However, to the best of the authors’ knowledge, such approaches have not
yet been used to write the full electronic Hamiltonian in second quantization.
This makes it difficult to use them using general correlated calculations,
although we note that perturbation theory calculations using multiwavelet
basis sets are feasible(kottmann2020, ). These approaches also do not deliver
spectral convergence with increasing basis set size, in contrast with the PW
basis.

In 1992, Gygi introduced a PW basis defined with respect to a curvilinear
coordinate system(Gygi1992, ; Gygi1993, ; Gygi1995, ) that deforms the domain
such that more basis functions are concentrated around the nucleus and fewer
basis functions are included in the interatomic space. Such a perspective aims
to build the adaptivity of an atom-centered basis into the PW framework,
maintaining its benefits such as spectral accuracy and fast implementation.
The curvilinear coordinates were optimized to minimize the energy of the mean-
field calculations using a double loop: in the inner loop the coefficients of
the molecular orbitals were optimized, and in the outer loop the map
describing the curvilinear coordinates was optimized. This work inspired other
followup calculations by various groups(Devenyi1994, ; JoseM.Perez-Jorda1995,
; Perez-Jorda1998, ; Perez-Jorda2011, ; Perez-Jorda2014, ; Hamann1995, ;
Hamann1996, ; Hamann2001, ; Hamman1995, ; Modine1997, ; Zumbach1996, ;
Rodriguez2008, ).

Our work here is inspired by the pioneering work of Gygi, and we modify the
approach by incorporating the so-called diagonal approximation, as well as a
new framework for computing the deformation. The concept of the diagonal
approximation was formalized in the work of White et al in the development of
the Gausslet basis functions(White2017, ) and various extensions (White2019, ;
Qiu2021, ). This work demonstrates that for basis functions which are
orthogonal and interpolating (i.e., behaving like delta functions, cf. Section
II for more details), the two-electron integrals can be effectively replaced
by a two-index, rather than four-index, quantity. The original Gausslets (like
PWs) offer uniform resolution across the entire computational domain and thus
suffer from the same problems as PW in terms of the efficiency of the basis
representation. We highlight that the Gausslets were designed in part for
deployment within tensor network approaches to electronic structure, where
they nonetheless enjoy advantages over PWs in maintaining locality of
entanglement and allowing for efficient representation of the many-body
Hamiltonian. It is worth pointing out that such ‘diagonal’ basis sets also
deliver significant advantage for quantum computing applicationsbabbush2018low
; mcclean2020discontinuous .

All the same, the original shortcoming of non-adaptivity was partially
overcome in the work of White and Stoudenmire(White2019, ) by introducing
curvilinear coordinates while maintaining the diagonal approximation. The
deformation used in this work is of separable or multi-sliced (cf. Section
IV.1 below) form, enabling essentially exact evaluation of all required
integrals. However, for a fully 3-D system, multi-sliced deformations are not
flexible enough to distribute basis functions ideally throughout the
computational domain. This shortcoming was addressed in a recent paper by
White and Lindsey(White2023, ), which is not simply based on a deformation of
an underlying grid of basis functions. This work introduces nested Gausslets
basis sets, which are carefully constructed using a procedure based on the
coarsening of 1-D diagonal basis sets constructed in terms of a separable
deformation near atom centers. In particular, the ‘COMX theorem’ demonstrates
that this coarsening can be achieved with a simple and tractable
diagonalization. Ultimately the nested Gausslet construction eliminates
unnecessary functions in regions where high resolution is not needed,
delivering controlled accuracy with a small number of basis functions and
maintaining the diagonal approximation. In addition, the construction can be
combined with GTOs while still maintaining diagonality. However, so far the
construction has only been realized for atomic and diatomic molecules,
although extensions to fully 3-D systems may be possible with suitable
extensions. A side effect of using nested Gausslets is that one needs to store
a dense matrix of two-electron integrals, with memory cost scaling
quadratically in the number of nested Gausslets. Although the basis set size
significantly improves on that of the original Gausslets, the relatively large
number of basis functions compared to a typical GTO basis set could make the
storage of two-index quantities a potential bottleneck for large systems.

In this work we combine the periodic sinc basis setboyd2001chebyshev  with a
curvilinear coordinate system. The sinc functions are identical in span to
plane waves but allow us to achieve the diagonal approximation, which we
explain as a consequence of pseudospectral approximation boyd2001chebyshev .
Our curvilinear deformation is fully 3-D, allowing for the treatment of
general molecular geometries in a straightforward and automatic fashion.
However, due to the use of general 3-D deformations, it is nontrivial to
calculate the integrals required to construct the Hamiltonian. Following the
work of GygiGygi1992 , we demonstrate that they can be computed efficiently
within a pseudospectral approximation framework that introduces little
additional error on top of the basis truncation error itself. In addition, we
eliminate the non-analyticity due to nuclear cusp by using all-electron
pseudopotentials of Gygi(Gygi2023, ) that are highly accurate and can be
systematically improved to arbitrary precision, enabling our pseudospectral
approach to enjoy spectral convergence as the discretization is refined. Using
our implementation, the cost of mean-field calculations scales log-linearly
with the number of basis functions, just like calculations in a traditional PW
basis. However, in our approach the number of basis functions is dramatically
reduced.

A novel aspect of our work is the computation of the map or deformation
defining the curvilinear coordinate system. In contrast to preceding work, we
take as given a prescribed density of grid points that the deformation should
achieve. Such a density can be supplied either from an approximate preliminary
calculation or from a formula inspired by multiresolution analysis. In
particular, the functional form that we have used was suggested in the nested
Gausslet paper(White2023, ), cf. (VI.1) below. We shall demonstrate that this
approach based on matching a prescribed density yields rapidly converging
energies as the discretization is refined. Moreover, the invertibility of the
deformation is guaranteed by the construction, even when the deformation is
quite sharp. This property is difficult to maintain in more heuristic
approaches when extreme deformations are required. Although the goal of
matching a prescribed density poses an apparently difficult computational
problem, namely, the problem of solving a Monge-Ampère equation Villani_2003 ,
we show how this problem can be solved efficiently with spectral accuracy
using a novel approach. We call our transformation solving this problem the
_cyclic Knothe-Rosenblatt flow_ , which can be viewed as a collocation method
for the Monge-Ampère equation in both forward and inverse formulations. In
particular, we demonstrate that computing the cyclic Knothe-Rosenblatt flow
achieving the map from uniform to curvilinear coordinates can be achieved with
Chebyshev or PW functions at a cost that scales linearly with the number of
function, so that the step of computing the deformation does not impose a
computational bottleneck in the broader context of our work.

The rest of the paper is organized as follows. In Section 2, we begin with
preliminaries and briefly describe the pseudospectral method and diagonal
approximation. In Section 3, we show how the use of the pseudospectral method
leads to a diagonal approximation for our deformed sinc basis set. In Section
4, we show how to solve the Monge-Ampère equation using the cyclic Knothe-
Rosenblatt flow. In Section 5, we show how the integrals needed for mean-field
calculations are evaluated, and we also describe our log-linear-scaling
algorithm for solving the Hartree-Fock equations. Then in Section 6 we present
numerical results demonstrating that rapid energy convergence is obtained for
the various systems we have studied. Finally, we close with conclusions and a
discussion of future work.

##  II Preliminaries

Consider an electronic stucture problem, discretized with the orthonormal
basis set $\\{\eta_{i}(\boldsymbol{x})\\}$. The quantum chemistry Hamiltonian

$\sum_{i,j}\mathbf{H}_{ij}\,a_{i}^{\dagger}a_{j}+\frac{1}{2}\sum_{ijkl}\mathbf{V}_{ijkl}\,a_{i}^{\dagger}a_{k}^{\dagger}a_{l}a_{j}$
(II.1)

is specified by one-electron and two-electron integrals
$\\{\mathbf{H}_{ij}\\}$ and $\\{\mathbf{V}_{ijkl}\\}$, respectively. In turn
$\mathbf{H}_{ij}=\mathbf{T}_{ij}+\mathbf{U}_{ij}$ can be written as a sum of
kinetic and potential contributions

$\mathbf{T}_{ij}=-\int\eta_{i}(\boldsymbol{x})\Delta\eta_{j}(\boldsymbol{x})\,d\boldsymbol{x},\quad\mathbf{U}_{ij}=\sum_{I}\int\eta_{i}(\boldsymbol{x})V_{I}(\boldsymbol{x})\eta_{j}(\boldsymbol{x})\,d\boldsymbol{x},$
(II.2)

where $\Delta$ denotes the Laplacian operator, $I$ indexes the atoms with
atomic numbers $Z_{I}$ and centers $R_{I}$, and $V_{I}$ is the potential due
to nucleus $I$. Note that $V_{I}$ can either be the Coulomb potential
$V_{I}(\boldsymbol{x})=\frac{-Z_{I}}{|\boldsymbol{x}-R_{I}|}$ or a
pseudopotential. We will explore both possibilities below.

Meanwhile, the two-electron integrals are given by

$\mathbf{V}_{ijkl}=\int\int\eta_{i}(\boldsymbol{x})\eta_{j}(\boldsymbol{x})\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}\eta_{k}(\boldsymbol{y})\eta_{l}(\boldsymbol{y})\,d\boldsymbol{x}\,d\boldsymbol{y}.$
(II.3)

The structure of the quantum chemistry Hamiltonian, especially the storage of
the two-electron integrals, simplifies considerably for certain ‘diagonal’
basis sets, which behave like smooth delta functions in the sense that

$\int
f(\boldsymbol{x})\eta_{i}(\boldsymbol{x})d\boldsymbol{x}=w_{i}f(\boldsymbol{x}_{i}),$
(II.4)

where $\boldsymbol{x}_{i}$ is the spatial center of the $i$-th basis function
and $w_{i}$ is its weight, independent of $f$. This implies that for a
sufficiently smooth function $f$ one has

$f(\boldsymbol{x})=\sum_{i}w_{i}f(\boldsymbol{x}_{i})\eta_{i}(\boldsymbol{x}),$
(II.5)

i.e., the coefficients of the expansion can simply be read off by evaluating
the function at grid points $\boldsymbol{x}_{i}$. Diagonal basis sets include
(periodic) sinc functionsboyd2001chebyshev  as well as more recently
introduced localized alternatives such as the gaussletsWhite2017  and their
various extensionsWhite2019 ; Qiu2021 ; White2023 .

With a diagonal basis one needs to retain only a one-index quantity to store
the nuclear integrals (instead of a two-index quantity) and two-index quantity
to store the 2-electron Coulomb integrals (instead of a four-index quantity).
In fact, as we shall explore below, it is possible to deal with the Coulomb
integrals in matrix-free fashion, avoiding the need to even form this two-
index quantity.

To see this, let us assume that the subscripts $p,q,r,s$ are used to represent
an arbitrary molecular orbital (MO). Then, abusing notation slightly by
overloading $\mathbf{U}$, the nuclear matrix in the MO basis is given by

$\displaystyle\mathbf{U}_{pq}$ $\displaystyle=\sum_{I}\int
f_{p}(\boldsymbol{x})f_{q}(\boldsymbol{x})V_{I}(\boldsymbol{x})\,d\boldsymbol{x}$
$\displaystyle\approx\sum_{I}\int\sum_{i}w_{i}f_{p}(\boldsymbol{x}_{i})f_{q}(\boldsymbol{x}_{i})\eta_{i}(\boldsymbol{x})V_{I}(\boldsymbol{x})\,d\boldsymbol{x}$
$\displaystyle=\sum_{i}w_{i}^{2}f_{p}(\boldsymbol{x}_{i})f_{q}(\boldsymbol{x}_{i})\left[\sum_{I}\frac{1}{w_{i}}\int\eta_{i}(\boldsymbol{x})V_{I}(\boldsymbol{x})\,d\boldsymbol{x}\right]$
$\displaystyle=\sum_{i}f_{pi}\,f_{qi}\,\mathbf{u}_{i},$ (II.6)

where $\mathbf{u}_{i}$ is defined as suggested by the last equality and
$f_{pi}:=w_{i}f_{p}(\boldsymbol{x}_{i})$ is the coefficient of $f_{p}$ in the
$\\{\eta_{i}\\}$ basis, following (II.5). The equation (II.6) emphasizes that
one is only required to evaluate and store the single-index quantity
$\mathbf{u}_{i}$, which is independent of the MOs.

Similarly, the two-electron integrals in the MO basis (once again abusing
notation slightly) are given by

$\displaystyle\mathbf{V}_{pqrs}$ $\displaystyle=\int\int
f_{p}(\boldsymbol{x})f_{q}(\boldsymbol{x})\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}f_{r}(\boldsymbol{x}^{\prime})f_{s}(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}\,d\boldsymbol{x}^{\prime}$
$\displaystyle\approx\sum_{ij}w_{i}^{2}f_{p}(\boldsymbol{x}_{i})f_{q}(\boldsymbol{x}_{i})\left[\frac{1}{w_{i}w_{j}}\int\int\eta_{i}(\boldsymbol{x})\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}\eta_{j}(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}\,d\boldsymbol{x}^{\prime}\right]w_{j}^{2}f_{r}(\boldsymbol{x}_{j})f_{s}(\boldsymbol{x}_{j})$
$\displaystyle=\sum_{ij}f_{pi}\,f_{qi}\,\mathbf{v}_{ij}\,f_{rj}\,f_{sj},$
(II.7)

i.e., the two-electron integrals in the MO basis are determined by the two-
index quantity $\mathbf{v}_{ij}$ which is independent of the MOs. Note the
similarity of (II.7) to the tensor hypercontraction (THC) formTHC1 ; THC2  of
the two-electron integrals integrals. However, unlike conventional THC, used
as a tool for compressing a given GTO basis, this approach will yield complete
basis set results when the number of diagonal basis functions is taken to be
sufficiently large.

Note that the approximations (II.6) and (II.7) for the nuclear and two-
electron integrals in the MO basis ($p,q,r,s$) would hold if it were the case
that the nuclear and two-electron integrals (II.2) and (II.3) in the diagonal
basis $\\{\eta_{i}\\}$ were diagonal i.e.
$\mathbf{U}_{ij}\approx\mathbf{u}_{i}\,\delta_{ij}$ and
$\mathbf{V}_{ijkl}\approx\mathbf{v}_{ik}\,\delta_{ij}\,\delta_{kl}$, as can be
verified by substitution of (II.5). Such a condition is however sufficient but
not necessary, and indeed even while such approximations for $\mathbf{U}_{ij}$
and $\mathbf{V}_{ijkl}$ do not hold elementwise to high accuracy, the
approximations (II.6) and (II.7) for $\mathbf{U}_{pq}$ and $\mathbf{V}_{pqrs}$
can still hold to high accuracy due to the smoothness of the MOs. Moreover,
any computational procedure for ground state calculations can treat the
substitutions $\mathbf{U}_{ij}\leftarrow\mathbf{u}_{i}\,\delta_{ij}$ and
$\mathbf{V}_{ijkl}\leftarrow\mathbf{v}_{ik}\,\delta_{ij}\,\delta_{kl}$ as if
they were exact, in the sense of pseudospectral methods.

Such a ‘diagonal approximation’ is not variational, since the effective one-
and two-electron integrals used in such a calculation are not obtained by
exact Galerkin truncation. Thus in particular the Hartree-Fock results
presented in this work do not converge _from above_ to the exact (i.e.,
complete basis set) Hartree-Fock results. However, the non-variational error
incurred due to the diagonal approximation is of similar magnitude to the
error due to truncation error of the underlying diagonal basis set (for
example, see Figures 8 and 9 in RefWhite2017 ). At the same time, the energy
obtained from these basis sets (treated either variationally or within the
diagonal approximation) converges to the exact result only algebraically,
rather than exponentially, due to the non-analyticity of the Coulomb
potential. To overcome this slow convergence, we will make use of the all-
electron pseudopotential (PP) proposed recently by GygiGygi2023 . This PP can
be made arbitrarily accurate by modifying just a single parameter (called $a$
in the original paper). For the calculations presented here we use $a=4$ which
is expected to give results that are within a few tens of $\mu$Eh of the exact
results for elements in the first two rows of the periodic table.

In this paper we show how to extend existing grid-based diagonal basis sets by
composing with an arbitrary invertible transformation. Specifically, we will
focus on the case where the original diagonal basis set is a set of periodic
sinc functions, or sinc functions for short, whose span is identical to that
of a suitable set of planewaves adapted to a periodic box. For periodic sinc
functions the weights satisfy $w_{i}=w=\sqrt{\frac{\mathcal{V}}{M}}$, where
$\mathcal{V}$ is the volume of the computational box and $M$ is the number of
basis functions or grid points in the domain, so $w_{i}^{2}$ can be viewed as
a quadrature weight for integrals performed over the computational grid. Note
that sinc functions have the property that
$\phi_{i}(\boldsymbol{x}_{j})=w^{-1}\delta_{ij}$. Moreover, for sinc functions
it is well-known that (II.4) is satisfied for suitably band-limited functions.
Further details on this choice of basis are given in Appendix A.1.

##  III Diagonal approximation with the deformed basis

One of the shortcomings of sinc and Gausslet basis is that it provides uniform
resolution in the entire domain of the problem. Thus even though high
resolution is needed only near the nucleus, one has to uniformly increase the
number of basis function over the entire computational domain to get accurate
results. This is wasteful, and it was noted is the early 1990s (by
Gygi(Gygi1992, ; Gygi1993, ; Gygi1995, ; Gygi1995a, ), followed by
others(Devenyi1994, ; JoseM.Perez-Jorda1995, ; Perez-Jorda1998, ; Perez-
Jorda2011, ; Perez-Jorda2014, ; Hamann1995, ; Hamann1996, ; Hamann2001, ;
Hamman1995, ; Modine1997, ; Zumbach1996, ; Rodriguez2008, )) that one can
overcome this shortcoming by introducing a deformation that preferentially
introduces a larger number of basis functions or grid points near the nucleus.

To see how this works, let us say we are given a collection of sinc functions
$\\{\phi_{i}\\}$ centered at uniform rectangular grid points
$\\{\boldsymbol{y}_{i}\\}$. Let us now introduce an invertible transformation
$T:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$. Then a deformed diagonal basis
set $\\{\eta_{i}\\}$ with centers
$\\{\boldsymbol{x}_{i}=T^{-1}(\boldsymbol{y}_{i})\\}$ can be constructed
formally as

$\displaystyle\eta_{i}(\boldsymbol{x})$
$\displaystyle:=\psi_{i}(\boldsymbol{x})\sqrt{J(\boldsymbol{x})},$
$\displaystyle\psi_{i}(\boldsymbol{x})$
$\displaystyle:=\phi_{i}(T(\boldsymbol{x}))$

where $J(\boldsymbol{x}):=\det DT(\boldsymbol{x})$ is the Jacobian determinant
of the deformation map and can be interpreted as the specifying the density of
point at position $\boldsymbol{x}$. In the following it is useful to set the
following notation:

  * •

$\phi_{i}$ denote the sinc basis functions.

  * •

$\boldsymbol{y}_{i}$ denote the uniform grid of sinc centers.

  * •

$\boldsymbol{x}_{i}=T^{-1}(\boldsymbol{y}_{i})$ denote the deformed grid
points.

  * •

$\psi_{i}(\boldsymbol{x}):=\phi_{i}(T(\boldsymbol{x}))$ denote the transformed
sinc functions, which inherit the interpolating property.

  * •

$\eta_{i}(\boldsymbol{x}):=\phi_{i}(T(\boldsymbol{x}))\,\sqrt{J(\boldsymbol{x})}=\psi_{i}(\boldsymbol{x})\,\sqrt{J(\boldsymbol{x})}$
denote the orthogonal diagonal basis functions.

  * •

$J_{i}:=J(\boldsymbol{x}_{i})$ denote the determinant of the Jacobian at the
deformed grid points.

It is important to note the distinction between $\psi_{i}$ and $\eta_{i}$. The
latter form the orthogonal diagonal basis that we actually use as the quantum
chemistry basis set. However, for a ‘suitably smooth’ $f$ we can calculate

$f(T^{-1}(\boldsymbol{y}))\approx
w\sum_{j}f(T^{-1}(\boldsymbol{y}_{j}))\,\phi_{j}(\boldsymbol{y})$

by evaluating (II.5) at $T^{-1}(\boldsymbol{y})$, which implies that for
arbitrary $\boldsymbol{x}$,

$\displaystyle f(\boldsymbol{x})$ $\displaystyle\approx
w\sum_{j}f(\boldsymbol{x}_{j})\,\psi_{j}(\boldsymbol{x}).$ (III.1)

Therefore the functions $\psi_{i}$ (not $\eta_{i}$) inherit the interpolating
property (II.5) from the sinc functions, as well as the property that
$\psi_{i}(\boldsymbol{x}_{j})=w^{-1}\delta_{ij}$.

In the above discussion, the precise meaning of ‘suitably smooth’ is that
$f\circ T^{-1}$ is sufficiently bandlimited. Therefore the $\psi_{i}$
functions are good interpolators even for functions with sharp features as
long as the the density $J(\boldsymbol{x})$ of points is also large in these
regions—and moreover with only dull features where $J(\boldsymbol{x})$ is
small.

Evidently we can also approximate

$f(\boldsymbol{x})\approx
w\sum_{j}\frac{f(\boldsymbol{x}_{j})}{\sqrt{J_{j}}}\eta_{j}(\boldsymbol{x}),$
(III.2)

as can be obtained by using the interpolating basis $\\{\psi_{j}\\}$ to expand
$f(\boldsymbol{x})/\sqrt{J(\boldsymbol{x})}$ instead of $f(\boldsymbol{x})$.
For electronic structure calculations we note that the basis set
$\\{\eta_{j}\\}$ is convenient because it is orthogonal. For mean-field
calculations we often go between representations in the $\\{\eta_{j}\\}$ and
$\\{\psi_{j}\\}$ basis representations given in (III.2) and (III.1),
respectively, by simply multiplying the coefficients in the $\\{\eta_{j}\\}$
expansion by $\sqrt{J_{j}}$ to get the coefficients in the $\\{\psi_{j}\\}$
expansion. The expansion in terms of $\\{\psi_{j}\\}$ is more convenient when
one wants to apply the kinetic and Coulomb operators (cf. Section V below),
while expansion in terms of $\\{\eta_{j}\\}$ is more convenient for
calculating overlaps between orbitals and maintaining orthogonality.

Using (III.2) it can be readily shown that the nuclear integrals of MOs will
be given (up to the approximation introduced above) by

$\mathbf{U}_{pq}=\sum_{i}w^{2}\frac{f_{p}(\boldsymbol{x}_{i})}{\sqrt{J_{i}}}\frac{f_{q}(\boldsymbol{x}_{i})}{\sqrt{J_{i}}}\left[\frac{\sqrt{J_{i}}}{w}\sum_{I}\int\eta_{i}(\boldsymbol{x})V_{I}(\boldsymbol{x})\,d\boldsymbol{x}\right]=\sum_{i}f_{pi}\,f_{qi}\,\mathbf{u}_{i},$
(III.3)

where $\mathbf{u}_{i}$ is the bracketed quantity in the second expression,
adapting the definition in (II.6) to the deformed basis. Meanwhile
$f_{pi}=f_{p}(\boldsymbol{x}_{i})/\sqrt{J_{i}}$ are once again the coordinates
of the MO $f_{p}$ in the orthonormal $\\{\eta_{i}\\}$ basis.

Similarly the 2-electron integrals can be written as

$\displaystyle\mathbf{V}_{pqrs}$
$\displaystyle=\sum_{ij}w^{2}\frac{f_{p}(\boldsymbol{x}_{i})}{\sqrt{J_{i}}}\frac{f_{q}(\boldsymbol{x}_{i})}{\sqrt{J_{i}}}\left[\frac{\sqrt{J_{i}\,J_{j}}}{w^{2}}\int\int\eta_{i}(\boldsymbol{x})\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}\eta_{j}(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}\,d\boldsymbol{x}^{\prime}\right]w^{2}\frac{f_{r}(\boldsymbol{x}_{j})}{\sqrt{J_{j}}}\frac{f_{s}(\boldsymbol{x}_{j})}{\sqrt{J_{j}}}$
$\displaystyle=\sum_{ij}f_{pi}\,f_{qi}\,\mathbf{v}_{ij}\,f_{rj}\,f_{sj}$
(III.4)

where the two-index quantity $\mathbf{v}_{ij}$, adapting (II.7), is given by
the bracketed quantity in the first line.

Alternatively, by instead using (III.1) we can similarly derive

$\mathbf{U}_{pq}=\sum_{i}f_{pi}\,f_{qi}\,\tilde{\mathbf{u}}_{i},\qquad\mathbf{V}_{pqrs}=\sum_{ij}f_{pi}\,f_{qi}\,\tilde{\mathbf{v}}_{ij}\,f_{rj}\,f_{sj},$

where

$\tilde{\mathbf{u}}_{i}:=\frac{{J_{i}}}{w}\sum_{I}\int\psi_{i}(\boldsymbol{x})V_{I}(\boldsymbol{x})\,d\boldsymbol{x},\qquad\tilde{\mathbf{v}}_{ij}:=\frac{{J_{i}\,J_{j}}}{w^{2}}\int\int\psi_{i}(\boldsymbol{x})\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}\psi_{j}(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}\,d\boldsymbol{x}^{\prime}.$
(III.5)

Similarly to the original sinc basis, the deformed sinc basis, together with
the suitable diagonal approximation, allows one to reduce the two-electron
integrals from a 4-index quantity to a 2-index quantity. Once again, to
construct an electronic structure Hamiltonian (II.1) with accurate ground-
state properties, we can perform the effective substitutions
$\mathbf{U}_{ij}\leftarrow\mathbf{u}_{i}\,\delta_{ij}$ and
$\mathbf{V}_{ijkl}\leftarrow\mathbf{v}_{ik}\,\delta_{ij}\,\delta_{kl}$, or
alternatively $\mathbf{U}_{ij}\leftarrow\tilde{\mathbf{u}}_{i}\,\delta_{ij}$
and
$\mathbf{V}_{ijkl}\leftarrow\tilde{\mathbf{v}}_{ik}\,\delta_{ij}\,\delta_{kl}$,
within the one- and two-electron integrals, in place of the exact Galerkin
formulations of (II.2) and (II.3). However, compared to the sinc basis (as we
shall verify in the numerical studies of Section VI), significantly fewer
deformed sinc basis functions are needed for high accuracy results, given an
appropriately chosen deformation $T$.

Unfortunately, the deformed sinc basis functions do not enjoy the same
separable form as the original sinc basis functions, if the deformation $T$ is
not of separable form. Thus the integrals $\mathbf{u}_{i},\mathbf{v}_{ij}$, as
well as the suitable matrix elements $\mathbf{T}_{ij}$ of the kinetic operator
defined in (II.2), are nontrivial to evaluate for a general deformation $T$.
Furthermore, although it is possible to specify a density $J(\boldsymbol{x})$
of deformed grid points, it is not trivial to specify a map $T$ that is
consistent with the prescribed density.

In the following two sections, we will explain how to approach these two key
computational difficulties of working with such a basis. First, in Section IV,
we will show how to efficiently construct a map $T(\boldsymbol{x})$ achieving
a prescribed Jacobian $J(\boldsymbol{x})$, i.e., a prescribed density of
deformed basis function centers. Second, in Section V, we will show how the
computation of the integrals $\mathbf{u}_{i},\mathbf{v}_{ij},\mathbf{T}_{ij}$
can be simplified with a pseudospectral approach. In fact we will show that
matrix-vector multiplication by suitable matrix quantities can be performed
efficiently without forming them as dense matrices, which will imply that for
mean-field calculations all operations can be performed at a cost that is log-
linear in the number of basis functions instead of quadratic. Thus we shall
achieve the same asymptotic scaling with respect to the number of basis
function as planewave calculations, but with a smaller basis set.

##  IV Constructing the coordinate transformation

Suppose we are given an arbitrary probability density $\rho$, which is not to
be conflated with the electron density, on a box
$\mathcal{B}:=\prod_{i=1}^{d}[a_{i},b_{i}]$ of volume $\mathcal{V}$ in
$\mathbb{R}^{d}$. We want to find a map $T:\mathcal{B}\rightarrow\mathcal{B}$
such that the Monge-Ampère equation Villani_2003

$\det(DT(\boldsymbol{x}))=\mathcal{V}\,\rho(\boldsymbol{x})$ (IV.1)

holds for $x\in\mathcal{B}$. The solution map $T$ pushes forward the density
$\rho$ on $\mathcal{B}$ to the uniform density $\mathcal{V}^{-1}$ on
$\mathcal{B}$. Note that the right-hand side is the density ratio of
$\rho(\boldsymbol{x})$ to the target density at $T(\boldsymbol{x})$, which is
simply the uniform density $\mathcal{V}^{-1}$ on the box. It is useful to
define $L_{i}:=b_{i}-a_{i}$ to be the length of the box in each dimension, so
$\mathcal{V}=\prod_{i=1}^{d}L_{i}$.

In the context of the deformed basis construction outlined above, our goal is
to solve (IV.1) where $\rho(\boldsymbol{x})\propto J(\boldsymbol{x})$ and $J$
is a prescribed density of deformed basis function centers. However, we will
consider (IV.1) more abstractly because ultimately our approach will view the
solution of (IV.1), for arbitrary right-hand side $\rho(\boldsymbol{x})$, as a
computational primitive that we shall wrap in several layers of abstraction in
order to obtain a deformation of high quality, as efficiently as possible.

We will build our approach in several stages. In Section IV.1 we review an
essential ingredient used to build our map: the Knothe-Rosenblatt transport
Rosenblatt ; Knothe . In Section IV.2, we explain how the Knothe-Rosenblatt
transport extends to the setting of a periodic computational domain. In
Section IV.3, we show how the Knothe-Rosenblatt transport can be computed
exactly when $\rho$ is presented as a sum of separable functions. In Section
IV.4, we explain how several Knothe-Rosenblatt transports (with the roles of
the variables cyclically permuted) can be composed to produce solutions to the
Monge-Ampère equation (IV.1) which are preferable to the unmodified Knothe-
Rosenblatt transport. In particular we explain how to fit each transport
problem in the sequence in the form (IV.1) with separable right-hand side
$\rho$. We call this approach the cyclic Knothe-Rosenblatt flow.

Finally, in Section IV.5, we wrap the cyclic Knothe-Rosenblatt flow in one
more layer of abstraction. This final complication is motivated by the fact
that, as we shall explain, it is computationally preferable to compute the
inverse map $S(\boldsymbol{y})=T^{-1}(\boldsymbol{y})$, which necessarily
satisfies the ‘inverse’ Monge-Ampère equation

$\det(DS(\boldsymbol{y}))\propto\frac{1}{J(S(\boldsymbol{y}))},$ (IV.2)

where $J$ is the prescribed density of basis function centers. This
formulation is more difficult to solve directly due to the dependence of the
right-hand side on the unknown map $S$ itself. However, motivated by the
approach of Lindsey and Rubinstein LindseyRubinstein2017 , we will show how
(IV.2) can be solved by computing a convergence sequence of solutions of
problems of the form (IV.1), where $\rho$ is determined self-consistently.

###  IV.1 Review of Knothe-Rosenblatt transport

The Monge-Ampère equation is underdetermined. Many solutions can be
constructed as optimal transport maps Villani_2003 , but there is in fact a
semi-explicit solution that forms the building block of our approach. To wit,
the Knothe-Rosenblatt transport Rosenblatt ; Knothe  is the unique solution of
(IV.1) among ‘triangular maps’ i.e., maps of the form
$T=(T_{1},\ldots,T_{d})$, where each $T_{k}$ depends only on
$x_{1},\ldots,x_{k}$. Thus, e.g., for $d=3$,

$DT(\boldsymbol{x})=\left(\begin{array}[]{ccc}\frac{\partial T_{1}}{\partial
x_{1}}(x_{1})&0&0\\\ &\frac{\partial T_{2}}{\partial x_{2}}(x_{1},x_{2})&0\\\
&*&\frac{\partial T_{3}}{\partial
x_{3}}(x_{1},x_{2},x_{3})\end{array}\right),$

and $\det(DT)$ is the product of the diagonal terms. This solution was used
previous by Perez-JordaJoseM.Perez-Jorda1995  to solve the Monge-Ampère
equation in the context of solution of electronic structure with adaptive
coordinates.

Let $\rho_{1}(x_{1})$ be the first marginal distribution of $\rho$, and
$\rho_{2}(x_{1},x_{2})$, etc. denote the marginals for the first several
variables. Moreover, let
$\rho_{2}(x_{2}|x_{1})=\rho_{2}(x_{1},x_{2})/\rho_{1}(x_{1})$ and
$\rho_{3}(x_{3}|x_{1},x_{2})=\rho(x_{1},x_{2},x_{3})/\rho(x_{1},x_{2})$, etc.
denote the conditional densities.

It follows from the Monge-Ampère equation (IV.1) and our restriction on the
form of $T$ that the equations

$\frac{\partial T_{1}}{\partial
x_{1}}(x_{1})=L_{1}\,\rho_{1}(x_{1}),\quad\frac{\partial T_{2}}{\partial
x_{2}}(x_{1},x_{2})=L_{2}\,\rho_{2}(x_{2}|x_{1}),\quad\frac{\partial
T_{3}}{\partial x_{3}}(x_{1},x_{2},x_{3})=L_{3}\,\rho_{3}(x_{3}|x_{1},x_{2}),$

etc., must hold.

The first equation is solved simply with an antiderivative:

$T_{1}(x_{1})=a_{1}+L_{1}\int_{a_{1}}^{x_{1}}\rho_{1}(y_{1})\,dy_{1}.$ (IV.3)

Evidently, $T_{1}$ satisfies the boundary conditions $T_{1}(a_{1})=a_{1}$ and
$T_{1}(b_{1})=b_{1}$.

Next we solve the second equation:

$T_{2}(x_{1},x_{2})=a_{2}+L_{2}\int_{a_{2}}^{x_{2}}\rho_{2}(y_{2}|x_{1})\,dy_{2},$
(IV.4)

which also achieves the second boundary condition $T_{2}(x_{1},a_{2})=a_{2}$
and $T_{2}(x_{1},b_{2})=b_{2}$.

Then the third equation:

$T_{3}(x_{1},x_{2},x_{3})=a_{3}+L_{3}\int_{a_{3}}^{x_{3}}\rho_{2}(y_{3}|x_{1},x_{2})\,dy_{2},$
(IV.5)

and so on.

It follows from the construction that

$T_{i}(x)=\begin{cases}a_{i}&\text{if}\ x_{i}=a_{i}\\\ b_{i}&\text{if}\
x_{i}=b_{i},\end{cases}$

i.e., $T$ fixes the boundary of $\mathcal{B}$.

It may be useful to compute the inverse map $S=T^{-1}$. In fact
$S=(S_{1},\ldots,S_{d})$ has the same form in that $S_{1}=S_{1}(y_{1})$,
$S_{2}=S_{2}(y_{1},y_{2})$, etc. Evidently $S_{1}(y_{1})=T_{1}^{-1}(y_{1})$
can be obtained by univariate inversion of $T_{1}$. Once we have found
$x_{1}=S_{1}(y_{1})$, we want to find $x_{2}$ such that

$T_{2}(x_{1},x_{2})=y_{2},$

which is achieved by univariate inversion of the function
$T_{2}(x_{1},\,\cdot\,)$. The result $x_{2}$ is $S_{2}(y_{1},y_{2})$.
Similarly, once we have found $x_{1},x_{2}$, we can obtain
$x_{3}=S_{3}(y_{1},y_{2},y_{3})$ by inverting the function
$T_{2}(x_{1},x_{2},\,\cdot\,)$.

###  IV.2 Periodic extension

In the case that $\rho$ extends smoothly and periodically from the box
$\mathcal{B}$ to $\mathbb{R}^{d}$, in fact the map
$T:\mathcal{B}\rightarrow\mathcal{B}$ also extends smoothly and periodically,
in the sense that the coordinate functions of $T-\mathrm{Id}$ are all
periodic, and the extension preserves each unit cell. We prove this claim in
Appendix B.1.

###  IV.3 Explicit solution for sum of separable functions

Suppose $\rho$ has the functional form

$\rho(\boldsymbol{x})=\sum_{\alpha}c^{\alpha}g^{\alpha}(\boldsymbol{x}),$

where each term
$g^{\alpha}(\boldsymbol{x})=\prod_{i=1}^{d}g_{i}^{\alpha}(x_{i})$ is
separable. Let
$G_{i}^{\alpha}(x_{i})=\int_{a_{i}}^{x_{i}}g_{i}^{\alpha}(y_{i})\,dy_{i}$ be
the appropriate antiderivative, and assume that we can evaluate
$g_{i}^{\alpha},G_{i}^{\alpha}$ easily. We show in Appendix B.2 how the exact
Knothe-Rosenblatt transport can be computed using such evaluations.

It is convenient to think of the basis index $\alpha$ as a multi-index
$\alpha=(n_{1},\ldots,n_{d})$ and define

$g_{i}^{\alpha}=P_{n_{i}}\circ u_{i},$

where $P_{n}$ is the $n$-th Chebyshev polynomial, defined on the standard
reference domain $[-1,1]$, and $u_{i}$ is the linear map sending
$[a_{i},b_{i}]$ to $[-1,1]$. This choice defines an efficient basis for the
expansion of the density $\rho$, and in Appendix B.3 we explain how the exact
computation of Appendix B.2 can be carried out explicitly in this setting.
This discussion amount to the computation of $G_{i}^{\alpha}$ in the special
case at hand.

In the case of periodic $\rho$, we can likewise expand in a basis
$g_{k}^{\alpha}=e_{n_{k}}\circ u_{k}$, where $e_{n}(x)=e^{2\pi in_{k}x}$ is a
suitable Fourier mode and $u_{k}$ maps $[a_{k},b_{k}]$ linearly to $[0,1]$. In
this case the computation of $G_{k}^{\alpha}$ is even simpler, cf. Appendix
B.4.

###  IV.4 Cyclic Knothe-Rosenblatt flow

A single Knothe-Rosenblatt transport is precisely the type of category of map
considered in the context of the multi-sliced Gausslet basis White2019 , which
was applied to quasi-1D molecular geometries. However, the Knothe-Rosenblatt
transport treats the coordinates $x_{1},x_{2},x_{3}$ asymetrically and can
produce awkward deformations in more general cases where it is not possible to
define a geometrically natural ordering of the coordinates, as illustrated in
Figure IV.1.

Figure IV.1: Consider $\rho(\boldsymbol{x})\propto
0.01+\sum_{I=1}^{3}e^{-3\sqrt{0.3+|\boldsymbol{x}-R_{I}|^{2}}}$ on the non-
periodic 3-dimensional box $[-4,4]^{3}$, where $R_{1}=(0,1,0)$,
$R_{2}=(1,-1,0)$, and $R_{3}=(-1,-1,0)$. We plot the curvilinear coordinates
$\boldsymbol{x}=T^{-1}(\boldsymbol{y})$ furnished by naive Knothe-Rosenblatt
transport. Left: image of $y_{2}$-$y_{3}$ coordinate lines in $x_{2}$-$x_{3}$
plane for fixed $y_{1}=0$. Right: image of $y_{1}$-$y_{3}$ coordinate lines in
$x_{1}$-$x_{3}$ plane for fixed $y_{2}=0$. Observe that vertical lines are
preserved by the coordinate transformation.

As such it is desirable to restore the permutation symmetry among the
coordinates. We will achieve this goal by incrementally building a transport
map as a composition of small Knothe-Rosenblatt transports in which the roles
of the coordinates are cycled. We call the resulting map the _cyclic Knothe-
Rosenblatt flow_.

Concretely, we will build

$T=T^{(N)}\circ\cdots\circ T^{(1)},$

where we choose $T^{(1)}$ such that

$\det(DT^{(1)}(\boldsymbol{x}))\propto\rho^{1/N}(\boldsymbol{x})$

for all $\boldsymbol{x}$.

Then we choose $T^{(2)}$ such that

$\det(DT^{(2)}(\boldsymbol{x}))\propto\rho^{1/N}\left([T^{(1)}]^{-1}(\boldsymbol{x})\right),$

or equivalently

$\det(DT^{(2)}(T^{(1)}(\boldsymbol{x})))\propto\rho^{1/N}\left(\boldsymbol{x}\right),$

for all $\boldsymbol{x}$.

This ensures that

$\displaystyle\det\left(D\left[T^{(2)}\circ
T^{(1)}\right](\boldsymbol{x})\right)$
$\displaystyle=\det\left(DT^{(2)}(T^{(1)}(\boldsymbol{x}))\,DT^{(1)}(\boldsymbol{x})\right)$
$\displaystyle=\det\left(DT^{(2)}(T^{(1)}(\boldsymbol{x}))\right)\det\left(DT^{(1)}(\boldsymbol{x})\right)$
$\displaystyle\propto\rho^{1/N}(x)\,\rho^{1/N}(\boldsymbol{x})$
$\displaystyle=\rho^{2/N}(\boldsymbol{x}).$

More generally, we construct $T^{(n)}$ inductively in terms of $T^{(n-1)}$ to
ensure that

$\det(DT^{(n)})\propto\rho^{(n)},$

where

$\rho^{(n)}\propto\rho^{1/N}\circ[T^{(n-1)}\circ\cdots\circ T^{(1)}]^{-1}.$

Once $\rho^{(n)}$ is fit with a separable expansion (e.g., Chebyshev or
Fourier), the map $T^{(n)}$ can be computed according to the procedure in
Section IV.3 for computing the Knothe-Rosenblatt transport. However, as we
iterate $n=1,\ldots,N$, we cycle the roles of the coordinates before computing
the Knothe-Rosenblatt transport.

How do we fit each $\rho^{(n)}$ with a separable expansion? It is useful to
note that

$\rho^{(n)}\propto\rho^{(n-1)}\circ\left[T^{(n-1)}\right]^{-1}.$

Therefore we can evaluate $\rho^{(n)}$ on Chebyshev nodes provided we can
invert the incremental map $T^{(n-1)}$ and evaluate the previous right-hand
side $\rho^{(n-1)}$, which we assume inductively is already fit with a
Chebyshev expansion. Once we have evaluated $\rho^{(n)}$ on the Chebyshev
nodes, the fitting with Chebyshev polynomials follows from standard
techniques. The same approach applies _mutatis mutandi_ to the periodic case
where the right-hand sides may be fit with Fourier expansions by evaluation on
an equispaced grid.

Figure IV.2: Just as in Figure IV.1, we consider $\rho(\boldsymbol{x})\propto
0.01+\sum_{I=1}^{3}e^{-3\sqrt{0.3+|\boldsymbol{x}-R_{I}|^{2}}}$ on the non-
periodic 3-dimensional box $[-4,4]^{3}$, where $R_{1}=(0,1,0)$,
$R_{2}=(1,-1,0)$, and $R_{3}=(-1,-1,0)$. We plot the curvilinear coordinates
$\boldsymbol{x}=T^{-1}(\boldsymbol{y})$ furnished by cyclic Knothe-Rosenblatt
flow with $N=40$. Left: image of $y_{2}$-$y_{3}$ coordinate lines in
$x_{2}$-$x_{3}$ plane for fixed $y_{1}=0$. Right: image of $y_{1}$-$y_{3}$
coordinate lines in $x_{1}$-$x_{3}$ plane for fixed $y_{2}=0$.

The cyclic Knothe-Rosenblatt flow is illustrated in Figure IV.2 using the same
choice of $\rho$ as in our demonstration of the ordinary Knothe-Rosenblatt
transport in Figure IV.1. Compared to the ordinary Knothe-Rosenblatt
transport, our flow yields a smoother transformation with significantly
better-conditioned coordinate boxes.

Once the cyclic Knothe-Rosenblatt flow $T$ is constructed,
$T^{-1}(\boldsymbol{y})$ can be evaluated by successively inverting the
incremental maps $T^{(n)}$. Since these are all triangular maps, they can be
inverted in three space dimensions by solving $d=3$ successive univariate
equations. However, since each incremental map is a small perturbation of the
identity, it can also be inverted efficiently be a few iterations of fixed-
point iteration. Concretely, if we want to solve
$\boldsymbol{y}=T^{(n)}(\boldsymbol{x})$ for $\boldsymbol{x}$, we can look for
a fixed point of the map
$\boldsymbol{x}\mapsto\boldsymbol{x}-\left[T^{(n)}(\boldsymbol{x})-\boldsymbol{y}\right]$
using fixed-point iteration (possibly with mixing to aid convergence).

We comment that the evaluation of a single layer of Knothe-Rosenblatt
transport is easier to perform on a separable grid than an arbitrary
collection of points, as can be readily observed from the conclusion of
Appendix B.4. However, several steps in our approach require evaluation of
such a layer on a non-separable collection of points. Specifically, the
evaluation of the inverse maps $T^{(n)}$ by fixed-point iteration requires
evaluations on unstructured grids. Moreover, a full forward pass through the
transport as a composition of the incremental maps $T^{(n)}$ necessarily
requires evaluation on unstructured grids. In addition, the fitting step
wherein $\rho^{(n)}$ is evaluated on a separable grid requires us to evaluate
the previous $\rho^{(n-1)}$ on an unstructured grid, which is similarly
inefficient. We accelerate all of these steps using non-uniform fast Fourier
transforms (NUFFTs) BOYD1992243 ; DuttRokhlin  to perform fast and highly
accurate interpolation of such unstructured evaluations from a fast evaluation
on a separable grid, maintaining log-linear scaling in the size of our
computational grid throughout the implementation. We use the FiNUFFT library
barnett_parallel_2019 ; cufinufft ; barnett_aliasing_2021  to evaluate our
NUFFTs.

###  IV.5 Self-consistent inverse transport

As mentioned in the opening of Section IV, we wrap the cyclic Knothe-
Rosenblatt transport in one final computational layer.

Our motivation is as follows. Suppose we try to solve the Monge-Ampère
equation (IV.1), where the right-hand side is given by
$\rho(\boldsymbol{x})\propto J(\boldsymbol{x})$. Within each layer of the
computation of the cyclic Knothe-Rosenblatt transport, we must fit the right-
hand side $\rho^{(n)}$ with a sum of separable functions, by pointwise
evaluation on a separable computational grid. However, $\rho(\boldsymbol{x})$,
and therefore the intermediate functions $\rho^{(n)}(\boldsymbol{x})$ as well,
are very sharp precisely in the case where $J(\boldsymbol{x})$ is sharp. Just
as we want to avoid resolving the wavefunction with uniform resolution across
the entire computational domain, so also we want to avoid resolving the
nonuniform density $J(\boldsymbol{x})$ of basis functions with uniform
accuracy across the entire computational domain.

Figure IV.3: Image of the subset $[-4,4]^{2}\times\\{0\\}$ of the
computational domain $[-5,5]^{3}$ under $T^{-1}$, where $T$ is computed via
naive cyclic Knothe-Rosenblatt flow and $\rho$ is specified in the main text.
Left: result using 21 computational grid points per dimension. Right: result
using 65 computational grid points per dimension. Observe that significant
aliasing error is still present even for the finer computational grid.

The difficulty is illustrated in Figure IV.3, which illustrates results for
$\rho(\boldsymbol{x})$ chosen following the functional form (VI.1) that we use
in our experiments below, where we take $a=0.1$, $b=4$, $c=0.01$, and we
consider three atom centers $I=1,2,3$ with $Z_{I}=1$, positioned at the
vertices of an equilateral triangle on the unit circle in the $\\{x_{3}=0\\}$
plane. The computational box is $[-5,5]^{3}$ with periodic boundary
conditions. The prescribed $\rho$ is very sharp, with over 3 orders of
magnitude of variation in the density over the computational domain, and as a
result it is too expensive to overcome aliasing error using the naive cyclic
Knothe-Rosenblatt transport of Section IV.4. In other words, the computational
grid is too coarse in some regions to capture all qualitative features of the
intermediate functions $\rho^{(n)}$, and as a result the transformations
furnished by the method display visually apparent artifacts.

By comparison, suppose for the sake of argument that we already knew the
right-hand side
$\tilde{\rho}(\boldsymbol{y})\propto\frac{1}{\rho(S(\boldsymbol{y}))}$ of
(IV.2). In fact $\tilde{\rho}(\boldsymbol{y})$ is significantly smoother than
$\rho(\boldsymbol{x})$ because the composition of $\rho$ with $S$ precisely
smooths out the nonuniformity in $\rho$.

Thus if we apply the cyclic Knothe-Rosenblatt flow to the computation of $S$
solving

$\det DS(\boldsymbol{y})\propto\tilde{\rho}(\boldsymbol{y}),$ (IV.6)

it is reasonable to use a separable grid with uniform resolution to fit each
successive right-hand side $\rho^{(n)}$ in the flow.

The only catch is that $\tilde{\rho}$ is not known _a priori_ , so it must be
determined self-consistently. Therefore we propose the following procedure,
alternating the following steps until convergence is reached:

  1. 1.

Define $\tilde{\rho}(\boldsymbol{y})\propto\frac{1}{\rho(S(\boldsymbol{y}))}$,
fitting with a sum of separable functions by evaluation of $S$ on a suitable
grid.

  2. 2.

Construct $S$ solving (IV.2) by cyclic Knothe-Rosenblatt transport.

In practice the convergence is rapid and requires less than 20 iterations. We
comment that the same self-consistent loop was used in the work of Lindsey and
Rubinstein LindseyRubinstein2017  to solve general Monge-Ampère equations of
the form $\det DT(\boldsymbol{x})=f(\boldsymbol{x})/g(T(x))$, by likewise
reducing to the case where the right-hand side is independent of $T$.

In fact, the cyclic Knothe-Rosenblatt flow of Section IV.4, in which the
intermediate right-hand sides $\rho^{(n)}$ are fit by interpolating on the
uniform computational grid (via Chebyshev or Fourier interpolation), can be
viewed as a _collocation_ method for the Monge-Ampère equation (IV.1) in the
sense that the solution $T$ solves the PDE (IV.1) _exactly_ at the
computational grid points. Likewise, if the self-consistent procedure for
computing $S=T^{-1}$ is solved to convergence, this method can be viewed as a
collocation method for (IV.2), in the sense that $\det
DS(\boldsymbol{y})\propto 1/\rho(S(\boldsymbol{y}))$ holds _exactly_ at the
computational grid points in the $\boldsymbol{y}$ domain. Since
$\rho(S(\boldsymbol{y}))$ has more uniform smoothness, this approach is
preferable as a collocation method.

In turn $T$ can be recovered if desired from the inversion procedure outlined
at the end of Section IV.4, if desired. However, for downstream purposes all
that we shall require will be the evaluation on our uniform grid
$\\{\boldsymbol{y}_{n}\\}$ of $S$, i.e.,
$S(\boldsymbol{y}_{n})=T^{-1}(\boldsymbol{y}_{n})=\boldsymbol{x}_{n}$, as well
as of the full Jacobian matrix
$DS(\boldsymbol{y}_{n})=\left[DT(\boldsymbol{x}_{n})\right]^{-1}$, which we
obtain by analytical differentiation of the layers of Knothe-Rosenblatt as
computed in Section IV.3.

Figure IV.4: Image of the subset $[-4,4]^{2}\times\\{0\\}$ of the
computational domain $[-5,5]^{3}$ under $S=T^{-1}$, where $S$ is computed via
the self-consistent procedure outlined in Section IV.5 and $\rho$ is the same
as in Figure IV.3 Left: result using 21 computational grid points per
dimension. Right: result using 65 computational grid points per dimension.
Observe that, in contrast with Figure IV.3, there is no noticeable aliasing
error for the result with the coarser grid, which is almost converged to the
exact limit.

##  V Pseudospectral method for integral evaluation

Before describing how the integrals are evaluated, we provide a brief outline
of the Hartree-Fock (HF) algorithm used in this work, in order to clarify our
exact computational requirements. We refer the reader to several excellent
textbooksszabo1996modern ; martin2020electronic ; lin2019mathematical .

###  V.1 Review of SCF for HF equations

The HF equations are solved using the self-consistent field (SCF) iteration.
Within each iteration of the most basic approach, we perform the following two
steps, which are then repeated over many iterations until convergence of the
occupied molecular orbitals $f_{p}$ for $p=1,\ldots,N_{e}/2$ is achieved.

  1. 1.

Given current guesses for the occupied molecular orbitals $f_{p}$ for
$p=1,\ldots,N_{e}/2$, we construct a subroutine that can compute the action of
the Fock operator on an arbitrary smooth function $f(\boldsymbol{x})$. The
Fock operator is equal to the sum of kinetic ($\hat{T}$), nuclear ($\hat{N})$,
Coulomb ($\hat{J}$), and exchange ($\hat{K}$) operators, whose actions are
defined by:

$\displaystyle\hat{T}f(\boldsymbol{x})$ $\displaystyle=-\frac{1}{2}\Delta
f(\boldsymbol{x})$ $\displaystyle\hat{N}f(\boldsymbol{x})$
$\displaystyle=\sum_{I}V_{I}(\boldsymbol{x})f(\boldsymbol{x})$
$\displaystyle\hat{J}f(\boldsymbol{x})$
$\displaystyle=\left[\int\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}\,\rho(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}^{\prime}\right]f(\boldsymbol{x})$
$\displaystyle\hat{K}f(\boldsymbol{x})$
$\displaystyle=-\sum_{p=1}^{N_{e}/2}f_{p}(\boldsymbol{x})\int\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}f_{p}(\boldsymbol{x}^{\prime})f(\boldsymbol{x}^{\prime})\,d\boldsymbol{x}^{\prime}.$
(V.1)

Here $V_{I}(\boldsymbol{x})$ is once again the nuclear potential due to the
$I$-th nucleus, $N_{e}$ is the number of electrons, the index $p$ runs over
occupied orbitals (each of which corresponds to two spin-orbitals), and
$\rho(\boldsymbol{x})=2\sum_{p=1}^{N_{e}/2}\left|f_{p}(\boldsymbol{x})\right|^{2}$
is the electron density. We will represent the function $f$ as well as the
molecular orbitals $f_{p}$ via their evaluations $f(\boldsymbol{x}_{i})$ and
$f_{p}(\boldsymbol{x}_{i})$ on the deformed grid points, which can be related
to coefficients in the deformed interpolating basis $\psi_{i}$ via (III.1) or
in the deformed orthonormal basis $\eta_{i}$ via (III.2). By ‘computing the
action of the Fock operator,’ we mean that given these representations, we can
construct (within our pseudospectral approximation framework) the left-hand
sides of (V.1) evaluated at the same points
$\boldsymbol{x}=\boldsymbol{x}_{i}$.

  2. 2.

Given a subroutine computing the action of the Fock operator, we use the
Davidson method(davidson1975iterative, ) to obtain the $N_{e}/2$ lowest
eigenpairs of the Fock operator. The eigenvectors define the next guesses
$f_{p}$, $p=1,\ldots,N_{e}/2$, for the occupied molecular orbitals.

The basic 2-step algorithm outlined above often fails to converge, and we use
Pulay’s DIIS(pulay1980diis, ; pulay1982diis, ) to stabilize and accelerate the
convergence. The computational bottleneck in these steps is the application of
the exchange operator. To minimize the number of exchange operator actions
required, we use the adaptively compressed exchange (ACE) approach due to
Lin(lin2016adaptively, ; hu2017adaptively, ), which partitions the algorithm
into two nested loops. In the inner loop, the ACE exchange is kept fixed, and
it is only updated in the outer loop once the inner loops converges. Typically
one needs to perform $O(10)$ iterations of the inner loop between for each
outer loop iteration. In each iteration of the outer loop, we must solve
$N_{e}^{2}$ Poisson equations to construct the ACE operator. The overall cost
of the calculation will scale as $\tilde{O}(MN_{e}^{2})$, where the tilde
indicates the omission of log factors and $M$ is the number of basis
functions.

In this publication we will not focus on improving the convergence rate of the
SCF cycles or the Davidson solver. Instead we will only describe how the
various operator actions in (V.1) can be performed.

###  V.2 General operator matrix elements

As we will see below, the operator of central importance to us is the kinetic
energy operator, because it is used not only to evaluate the kinetic energy
but also to solve the Poisson equation for evaluating the Coulomb, exchange,
and nuclear terms. It will also become apparent in Section V.3 that the matrix
describing the action of the kinetic energy on the deformed sinc basis
functions $\psi_{i}(\boldsymbol{x})$ can be evaluated efficiently using FFTs.
However, because the $\psi_{i}(\boldsymbol{x})$ functions are not orthogonal,
care is needed while using this matrix as we describe below.

For a general operator $\hat{O}$ (e.g., $-\Delta$ for the kinetic energy or
$\Delta^{-1}$ for the solution operator of the Poisson equation), consider the
matrix elements $O_{ij}=\langle\eta_{i},\hat{O}\eta_{j}\rangle$ in our
orthonormal basis. Here $\langle\,\cdot\,,\,\cdot\,\rangle$ denotes the
ordinary $L^{2}$ inner product over the computational domain. The matrix
elements satisfy the defining property that for a basis expansion
$f(\boldsymbol{x})=\sum_{j}f_{j}\,\eta_{j}(\boldsymbol{x})$ of some function
$f$, we can expand $(\hat{O}f)(\boldsymbol{x})$ in the basis as
$(\hat{O}f)(\boldsymbol{x})=\sum_{i}\left(\sum_{j}O_{ij}f_{j}\right)\eta_{i}(\boldsymbol{x})$.
Equivalently, it holds that

$\langle\eta_{i},\hat{O}f\rangle=\sum_{j}O_{ij}f_{j}$ (V.2)

We claim that the defining property (V.2) still holds for the alternative
matrix elements

$O_{ij}=\sqrt{J_{i}}\,\langle\psi_{i},\hat{O}\,\psi_{j}\rangle\sqrt{J_{j}},$
(V.3)

provided that $f$ is suitably smooth relative to our deformed coordinates, as
if we were effectively identifying
$\eta_{i}(\boldsymbol{x})=\psi_{i}(\boldsymbol{x})\,\sqrt{J(\boldsymbol{x})}$
with $\psi_{i}(\boldsymbol{x})\,\sqrt{J_{i}}$. The claim is justified in
Appendix C.

Thus, it is easy to show that we can compute the action of $\hat{O}$ on a
smooth function $f$ by simply performing matrix-vector multiplications by the
matrix $\langle\psi_{i},\hat{O}\,\psi_{j}\rangle$. To see this, recall that
the coefficients $f_{i}$ in the $\\{\eta_{i}\\}$ basis correspond to deformed
grid point evaluations $f(\boldsymbol{x}_{i})$ via
$f_{i}=wf(\boldsymbol{x}_{i})/\sqrt{J_{i}}$, and likewise
$\langle\eta_{i},\mathcal{O}f\rangle=w\,(\hat{O}f)(\boldsymbol{x}_{i})/\sqrt{J_{i}}$.
Therefore we can alternatively express the action of $\hat{O}$ in terms of
evaluations on the deformed grid points via the identity

$(\hat{O}f)(\boldsymbol{x}_{i})={J_{i}}\sum_{j}\langle\psi_{i},\hat{O}\,\psi_{j}\rangle\,f(\boldsymbol{x}_{j})$
(V.4)

from which the matrix elements $\langle\eta_{i},\mathcal{O}f\rangle$ can be
readily obtained from $(\hat{O}f)(\boldsymbol{x}_{i})$.

###  V.3 Kinetic matrix

Following the discussion in the preceding Section V.2, we want to show how to
multiply by the Laplacian matrix elements
$\mathbf{L}_{ij}:=-\langle\psi_{i},\Delta\,\psi_{j}\rangle$. In Appendix A.2,
we explain how this operation can be performed with pseudospectral accuracy,
using 4 FFTs and 4 inverse FFTs.

The matrix-vector multiplication by $\mathbf{L}$ not only allows us to perform
the action of the kinetic operator $\hat{T}$, but also will serve as a key
computational primitive that we leverage in the calculation of both nuclear
and two-electron integrals. Both of these types of integral can be understood
in terms of the Poisson equation, which we shall discuss first in an abstract
context.

###  V.4 Solving the Poisson equation

In this section we explain how to use the Laplacian matrix $\mathbf{L}$ to
solve a general Poisson equation

$-\Delta u(\boldsymbol{x})=v(\boldsymbol{x})$ (V.5)

within the pseudospectral approximation, for general right-hand side
$v(\boldsymbol{x})$. The method is motivated by stipulating that the Poisson
equation (V.5) holds on the deformed grid points, i.e.,

$-\Delta u(\boldsymbol{x}_{i})=v(\boldsymbol{x}_{i})$ (V.6)

for all $i$.

To solve (V.6), take $\hat{O}=-\Delta$ in (V.4) to obtain equations

${J_{i}}\sum_{j}\mathbf{L}_{ij}\,u(\boldsymbol{x}_{j})=v(\boldsymbol{x}_{i}),$
(V.7)

which are equivalent to (V.6) up to the pseudospectral approximation.
Therefore we can solve the linear system

$\mathbf{L}\tilde{u}=b,$ (V.8)

where the right-hand side is defined by $b_{i}=v(\boldsymbol{x}_{i})/J_{i}$,
and the solution vector $\tilde{u}_{j}=(u(\boldsymbol{x}_{j}))$ is the vector
of evaluations of the solution $u(\boldsymbol{x})$ of the Poisson equation
(V.5) on the deformed grid points. In turn $u(\boldsymbol{x})$ can be
recovered in either the $\\{\psi_{i}\\}$ or $\\{\eta_{i}\\}$ basis via
interpolation, cf. (III.1) and (III.2).

Since $\mathbf{L}$ is positive definite by construction, we can solve the
linear system (V.8) using the conjugate gradient (CG) method, relying on the
subroutine discussed in Section V.3 to perform matrix-vector multiplications
by $\mathbf{L}$ without every forming the matrix fully. For our preconditioner
we use the exact inverse $\tilde{\mathbf{L}}^{-1}$ of the ‘undeformed’
Laplacian matrix
$\tilde{\mathbf{L}}_{ij}=\langle\phi_{i},-\Delta\phi_{j}\rangle$, which is
diagonalized by the discrete Fourier transform. Thus the predconditioner
$\tilde{\mathbf{L}}^{-1}$ can be applied to a vector using a single FFT and a
single inverse FFT. In our experiments we require on roughly 20 to 30 CG
iterations to solve (V.8). Since each iteration requires 5 forward and inverse
FFTs, the solve is about 100 times more expensive than the analogous solve for
a uniformly spaced sinc basis.

###  V.5 Coulomb and exchange operators

The action of the Coulomb and exchange operators $\hat{J}$ and $\hat{K}$
defined in in (V.1) can be computed by leveraging our subroutine for solving
the Poisson equation outline in the preceding Section V.4. Indeed, the action
of the integral Coulomb kernel
$\frac{1}{|\boldsymbol{x}-\boldsymbol{x}^{\prime}|}$ appearing in the
definitions of $\hat{J}$ and $\hat{K}$ can be viewed (up to a proportionality
constant) as a Poisson solve. In fact, the formulas of (V.1) involving the
Coulomb kernel are valid only for open boundary conditions, but more generally
the understanding of these operations as Poisson solves holds even when
periodic boundary conditions are considered, provided that the Laplacian is
considered with appropriate boundary conditions.

To apply $\hat{J}$, we must form the electron density values
$\rho(\boldsymbol{x}_{j})$ on our deformed grid and then performs the Poisson
solve

$\Delta u(\boldsymbol{x})=\rho(\boldsymbol{x})$

following the method of Section V.4, which recovers the values of the solution
$u(\boldsymbol{x}_{i})$ on the deformed grid points. Then
$\hat{J}f(\boldsymbol{x}_{i})$ can be constructed from $u(\boldsymbol{x}_{i})$
and $f(\boldsymbol{x}_{i})$ by pointwise multiplication on the deformed grid.

To form $\hat{K}f(\boldsymbol{x}_{i})$, the logic is much the same, except
that we must perform $N_{e}/2$ Poisson solves:

$\Delta u_{p}(\boldsymbol{x})=f_{p}(\boldsymbol{x})f(\boldsymbol{x})$

to furnish $u_{p}(\boldsymbol{x}_{i})$ for $p=1,\ldots,N_{e}/2$.

###  V.6 Pseudopotential-based nuclear operator

The action of the nuclear operator $\hat{N}$ in (V.1) is quite simple to
perform under the assumption that the nuclear potentials $V_{I}$ are smooth
relative to our deformed coordinates. Note that this is not exactly the case
for the bare nuclear potential
$V_{I}(\boldsymbol{x})=-Z_{I}/|\boldsymbol{x}-R_{I}|$, due to the singularity
at the atom center. By considering more strongly deformed coordinates near the
origin, we improve the quantitative smoothness near the center but cannot
remove the singularity. Thus we are motivated to introduce an all-electron
pseudopotential (PP) following the recent work of Gygi(Gygi2023, ). For the
purposes of this section, we view $V_{I}$ as smooth, having been constructed
as a PP following the discussion in Appendix D, and in this construction we
only need to assume that $V_{I}(\boldsymbol{x}_{i})$ can be evaluated on our
deformed grid. The details of an alternate approach using the bare nuclear
potential will be provided in the following Section V.7.

Under the PP assumptions outlined above, the action of the nuclear operator
$\hat{N}$ defined in (V.1) can be performed simply as

$\hat{N}f(\boldsymbol{x}_{i})=\sum_{I}V_{I}(\boldsymbol{x}_{i})f(\boldsymbol{x}_{i}).$

###  V.7 Bare nuclear operator

Consider a general diagonal operator $\hat{V}$ which operates as
$\hat{V}f(\boldsymbol{x})=V(\boldsymbol{x})f(\boldsymbol{x})$. Compared to the
discussion of Section V.2, we will consider an alternative defining property
for the matrix elements $V_{ij}$. Namely, given expansions
$f(\boldsymbol{x})=\sum_{i}f_{i}\,\eta_{i}(\boldsymbol{x})$ and
$g(\boldsymbol{x})=\sum_{i}g_{i}\,\eta_{i}(\boldsymbol{x})$ for suitably
smoooth functions $f$ and $g$ in the orthonormal $\\{\eta_{i}\\}$ basis, we
insist that

$\langle f,\hat{V}g\rangle=\int
f(\boldsymbol{x})g(\boldsymbol{x})V(\boldsymbol{x})\,d\boldsymbol{x}=\sum_{ij}f_{i}g_{j}V_{ij}.$
(V.9)

By the same manipulations as in Section III (cf. (III.3), (III.4), and
(III.5), this property is satisfied for either the choice
$V_{ij}=\mathbf{v}_{i}\,\delta_{ij}$ or
$V_{ij}=\tilde{\mathbf{v}}_{i}\,\delta_{ij}$, where
$\mathbb{v}_{i}=\frac{\sqrt{J_{i}}}{w}\int\eta_{i}(\boldsymbol{x})\,V(\boldsymbol{x})\,d\boldsymbol{x}$
or
$\tilde{\mathbb{v}}_{i}=\frac{J_{i}}{w}\int\psi_{i}(\boldsymbol{x})\,V(\boldsymbol{x})\,d\boldsymbol{x}$.
We choose the latter as the more convenient, leading to formulation

$\hat{N}f(\boldsymbol{x}_{i})=\sum_{I}V_{I,i}\,f(\boldsymbol{x}_{i}),$ (V.10)

for the action of the nuclear operator (V.1), where

$V_{I,i}:=\frac{J_{i}}{w}\int\psi_{i}(\boldsymbol{x})\,V_{I}(\boldsymbol{x})\,d\boldsymbol{x}$
(V.11)

can be evaluated even in the case of the bare nuclear potential
$V_{I}=-Z_{I}/|\boldsymbol{x}-R_{I}|$, which has an integrable singularity.
Note that the formulation (V.10) for the nuclear action defined in (V.1)
holds, via our defining property, in the sense of being accurate up to
projection onto the space of suitably smooth functions.

In Appendix E, we demonstrate how the calculation of $\sum_{I}V_{I,i}$ as in
(V.11) can be reduced to Poisson solves and in turn a single linear solve
involving the matrix $\mathbf{L}$.

##  VI Results

Here we present the results for calculations on He, He2, Be, Be2, and CH4
systems. All calculations are performed on a cubic unit cell of size $10$ Bohr
with periodic boundary conditions.

For each system we compute a single transformation $T$ by constructing the
self-consistent inverse cyclic Knothe-Rosenblatt flow following the discussion
of Section IV.5. The choice of $\rho$ that determines the transformation will
be outlined below for each system. Then the same transformation $T$ is used
with various choices $M$ of deformed sinc basis set size. For He, He2, Be, and
Be2, we compute $T$ using a $20\times 20\times 20$ computational grid. We use
a $30\times 30\times 30$ computational grid to compute the transformation for
CH4.

###  VI.1 He and He2

Potential | Energy
---|---
PP, $a=4$ | -2.58304 _2_
PP, $a=5$ | -2.58302 _7_
PP, $a=6$ | -2.58302 _3_
Bare Nuclear | -2.5830 _15_

Figure VI.1: The left figure shows the convergence of the energy of He atom as
the number of basis functions is increased. The four curves with the label
“distorted” show the results obtained with distorted sinc functions using both
the PP (with parameter values $a=4,5,6$) and the bare nuclear potential. The
difference between the energies is not visible on the scale of this graph. The
right graph is a zoomed-in version of the left graph where the convergence of
the various curves is more clearly visible. For additional clarity the
numerical values of the converged energies obtained with $M=70^{3}$ basis
functions are given in the table on the right. The italicized number
corresponds to the digit which is uncertain in our calculations due to the
finiteness of the basis set. The purple curve at left plots the convergence of
energy when plane waves (PWs) are used together with the bare nuclear
potential. Points on the purple curve are sufficiently far away from
convergence that they would fall outside of the range of values displayed on
the right curve.

For He and He2 we prescribe $J(\boldsymbol{x})\propto\rho(\boldsymbol{x})$
where

$\rho(\boldsymbol{x})=\sum_{I}\frac{\mathrm{erf}(Z_{I}|\boldsymbol{x}-R_{I}|/a_{I})-\mathrm{erf}(Z_{I}|\boldsymbol{x}-R_{I}|/b_{I})}{|\boldsymbol{x}-R_{I}|}+c.$
(VI.1)

Here $R_{I}$ and $Z_{I}$ are again the position and charge of the $I$-th
nucleus, and $a_{I}=0.1$, $b_{I}=4$, and $c=0.01$ are parameters that
determine the spatial extent of the deformation. This deformation ensures that
the density of grid points increases approximately as $1/r$ between the
distances of $b/Z_{I}$ Bohr and $a/Z_{I}$ Bohr from the center of the nucleus.
Outside of this region, the density is roughly a non-zero constant. The total
number of distorted sinc functions $M$ is progressively increased until
convergence is achieved. The Monge-Ampère equation (IV.1) is solved by the
method of Section IV.5 using $N=15$ steps and $20^{3}$ computational grid
points in each of the $N=15$ intermediate fitting steps.

Figure VI.1 concerns the He atom exclusively and is used to justify the
accuracy of our practical choice of all-electron pseudopotential (PP) for the
remaining calculations. The numbers in the table in Figure VI.1 show that the
energies of the pseudopotential (PP) calculations converge to that of the bare
nuclear calculations as the PP parameter $a$ is increased. In all the
calculations shown in the table we have used a large number ($M=70^{3}$) of
basis functions, and our estimated uncertainty due to basis set incompleteness
with PP are on the order of a few $\mu$Eh. The table shows that when $a=4$ is
used in the PP we can expect to get an energy which is within a few tens of
$\mu$Eh of the converged result compared to the bare nuclear potential. The
left graph in Figure (VI.1) shows that the convergence of the energy with
distorted grids is vastly improved compared the convergence with undistorted
grids. One would expect that the difference in performance would increase for
larger atoms, as all-electron calculations with plane waves for heavier atoms
quickly become impractical without the use of large supercomputers. The right
graphs shows the zoomed in version of the left graph where we have not plotted
the energies with the undistorted grids. This graphs shows that the converge
of energy is much more smoother when PP is used and the bare nuclear energies
are still not fully converged to $\mu$Eh even with $M=70^{3}$ number of basis
functions.

Figure VI.2: The figure on the left (resp., right) shows the convergence of
the energy with respect to the distorted basis set size for the He atom
(resp., He2 molecule), with comparisons to GTO energies plotted as horizontal
lines.

Figure (VI.2) shows that we can converge the energy of He atom and He dimer
(with a bond length of 2 Å) to accuracy surpassing that of the 5Z basis by
using $30^{3}$ and $40^{3}$ grid points, respectively, and the PP with $a=4$.
Interestingly, the He dimer requires less than double the number of basis
functions compared to the He atom, likely because the density contributions in
$\rho(\boldsymbol{x}$), cf. (VI.1), of the two atoms overlap. It is also worth
noting that the convergence of the energy with respect to the number of grid
points is significantly less smooth for the bare nuclear potential than for
the PP.

###  VI.2 Be and Be2

Next we use the procedure outlined above to evaluate the energy of the
beryllium atom and beryllium dimer (with a bond length of 2.44 Å). For these
systems we have used the same deformation as (VI.1) with $a_{I}=0.1$,
$b_{I}=8$, and $c=0.03$. Figure (VI.3) shows that the energy surpasses the
accuracy of the 5Z basis set with grid sizes of $30^{3}$ and $42^{3}$,
respectively, for the atom and the dimer. These numbers are not significantly
different than the ones we found for He atom, suggesting that with an
appropriately chosen deformation, heavy nuclei can be treated without
significant overhead. Further tests will be needed to confirm this claim for a
range of heavier atoms. We point out that, once again, the number of grid
points needed for the dimer is only slighly higher than the number needed for
the atom because the density contributions of the two nuclei in the
specification of $\rho(\boldsymbol{x})$ have significant overlap.

Figure VI.3: The figure on the left (resp., right) shows the convergence of
the energy with respect to the distorted basis set size for the Be atom
(resp., Be2 molecule), with comparisons to GTO energies plotted as horizontal
lines.

###  VI.3 CH4

For our final test we perform a calculation on methane, which does not have a
linear geometry, suggesting that simpler approaches based on 1-D deformations
will not be highly efficient. Again we use the deformation (VI.1), choosing
$a_{I}=0.1,\,b_{I}=18$ for the carbon atom and $a_{I}=0.1,\,b_{I}=1.5$ for the
hydrogen atoms, as well as $c=0.01$. These parameters are chosen to ensure
that a higher density of basis functions surrounds the carbon atom, which is
bonded to all the hydrogen atoms surrounding it. Results are plotted in Figure
VI.4. For this system we observe that for $M=50^{3}$ grid points, the energy
is converged to within a few tenths of mEh. We can further converge the
results within a few tens of $\mu$Eh by taking $M=80^{3}$ grid points.

Figure VI.4: The figure shows the convergence of the energy with respect to
the distorted basis set size for the methane molecule.

##  VII Discussion and conclusion

In this article we have shown how the two major challenges associated with the
use of the distorted sinc basis set can be overcome. These challenges are,
namely: (1) computing the transformation maps that yield a prescribed density
of basis functions and (2) calculating and operating with integrals for the
Hamiltonian. We have demonstrated that these techniques can be used to
evaluate mean-field energies at a cost that is log-linear in the number of
basis functions.

Several questions motivate further investigation in this direction. First, in
this work we have designed the deformation by prescribing a density of basis
functions with an analytical functional form. A more systematic approach with
less hand-tuning would be desirable. Such an approach might be achieved by
minimizing the variance of energy for a fixed number of grid points. (One
cannot directly minimize the energy because it is not variational due to the
use of diagonal approximation.) Another approach might be to minimize the
error in the one- and two-electron integrals for GTO basis functions incurred
by representing the GTOs in the interpolating basis of deformed sinc
functions.

Second, we expect to be able to get away with a smaller number of basis
functions if we can combine the deformed sinc basis set with sharp GTO basis
functions to describe the nuclear cusp. This perspective has been applied
fruitfully in several works White2023 ; GaviniEnriched1 ; dft-fe-1.0 . It is
an open question as to whether such a mixed basis set can be used while
keeping the cost of the calculation linear in the number of basis functions as
we have done here.

Third, currently the most significant cost of the calculation is the solution
of the Poisson equation. The use of the preconditioner as described in Section
V.4 significantly speeds up the calculations, but we still require on the
order of a few 10’s of CG iterations per solve. We plan to investigate more
sophisticated proconditioners in future work.

Finally, we plan to investigate applications of our basis set in frameworks
for correlated calculations, including tensor networks methods, auxiliary
field quantum Monte Carlo, variational Monte Carlo, and quantum algorithms.

##  VIII Acknowledgements

The authors would like to thank Francois Gygi and Steve White for helpful
discussions. S.S. was supported through the National Science Foundation grant
CHE-2145209. M.L. was supported by the Applied Mathematics Program of the US
Department of Energy (DOE) Office of Advanced Scientific Computing Research
under contract number DE-AC02-05CH11231.

## References

  * [1] H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al.  Molpro, version 2019.2, a package of ab initio programs.  J. Chem. Phys., 152:144107, 2020.
  * [2] Frank Aquilante, Leifeng De Vico, Nicholas Ferré, Gunnar Delcey, Ilias M. Angelis, Thomas B. Pedersen, Riccardo Lindh, and Roland Lindh.  Openmolcas: From source code to insight.  J. Chem. Theory Comput., 16(2):694–709, 2020.
  * [3] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox.  Gaussian 16, Revision C.01.  Gaussian, Inc., Wallingford, CT, 2019.
  * [4] Frank Neese.  Orca - an ab initio, density functional and semiempirical program package.  J. Comput. Chem., 30(4):167–178, 2012.
  * [5] A. S. Christensen, K. Kristensen, M. A. A. Jensen, and P. Jørgensen.  Dalton, a molecular electronic structure program, release dalton2018.2 (2018).  J. Chem. Theory Comput., 16(1):1944–1951, 2020.
  * [6] L. Visscher, H. J. Aa. Jensen, J. Autschbach, K. G. Dyall, I. S. P. Mancini, et al.  Dirac, a relativistic ab initio electronic structure program, release dirac18 (2018).  J. Chem. Phys., 152:144108, 2020.
  * [7] Marat Valiev, Dmytro Bykov, Doenoe Kwak, Justin M. Johnson, Zachary R. Kann, and Wyatt Slack.  Nwchem: Past, present, and future.  J. Chem. Phys., 152(20):204110, 2020.
  * [8] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery.  Advances in electronic structure theory: Gamess a decade later.  J. Comput. Chem., 14(11):1347–1363, 1993.
  * [9] Daniel G. A. Smith, Daniel R. Burns, Robert N. Parrish, Matthew C. Simmonett, Lori A. Schieber, Roseann M. Galvin, Edward D. Turner, Matthew S. Dellapiazza, Edward G. O’Brien, James S. Sirianni, Jarrod N. Nascimento, Taylor D. Kraus, Alexei Sokolov, Roberto Di Remigio, Antonio Alenaizan, Xavier Sergeev, John F. F. Mathewes, Hrant P. Lotrich, Adrian E. DePrince, Utiu F. Y., Rollin A. King, Edward F. Valeev, Francesco A. Evangelista, Justin M. Turney, Thom Dunning, and C. David Sherrill.  Psi4 1.1: An open-source electronic structure program emphasizing automation, advanced libraries, and interoperability.  Journal of Chemical Theory and Computation, 16(10):5612–5628, 2020\.
  * [10] Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan.  Pyscf: the python-based simulations of chemistry framework.  Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1):e1340, 2018.
  * [11] Jürg Hutter, Michele Iannuzzi, Florian Schiffmann, and Joost VandeVondele.  Cp2k: An electronic structure and molecular dynamics software package \- quickstep: Efficient and accurate electronic structure calculations.  Wiley Interdisciplinary Reviews: Computational Molecular Science, 4(1):15–25, 2014.
  * [12] R. Ahlrichs, M. Bär, H.-P. Baron, R. Bauernschmitt, S. Böcker, L. von Arnim, R. Farman, M. Häser, P. Horn, C. Huber, H. Stoll, W. König, C. Körzdörfer, T. Kußmann, K. Kümmel, U. Huniar, O. Treutler, A. Vöhringer, M. Deglmann, A. Grüneis, C. van Wüllen, and M. Ehrig.  Turbomole: Modular program suite for ab initio quantum-chemical and condensed-matter simulations.  J. Chem. Phys., 152(12):144107, 2020.
  * [13] R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, S. Casassa, L. Maschio, M. Ferrabone, M. De La Pierre, P. D’Arco, Y. Noel, M. Causa, M. Rerat, and B. Kirtman.  Crystal17: A program for the ab initio investigation of crystalline solids.  J. Chem. Phys., 152(20):204111, 2020.
  * [14] Volker Blum, Ralf Gehrke, Felix Hanke, Peter Havu, Ville Havu, Xinguo Ren, Karsten Reuter, and Matthias Scheffler.  Fhi-aims: an all-electron full-potential code from atoms to solids.  Comput. Phys. Commun., 180(11):2175–2196, 2009.
  * [15] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch.  Quantum espresso: a modular and open-source software project for quantum simulations of materials.  J. Phys.: Condens. Matter, 21(39):395502, 2009.
  * [16] G. Kresse and D. Joubert.  Projector augmented-wave method.  Phys. Rev. B, 59(3):1758–1775, 1999.
  * [17] Edward G. Hohenstein, Robert M. Parrish, and Todd J. MartÃnez.  Tensor hypercontraction density fitting. I. Quartic scaling second- and third-order MÃžller-Plesset perturbation theory.  The Journal of Chemical Physics, 137(4):044103, 07 2012.
  * [18] Robert M. Parrish, Edward G. Hohenstein, Todd J. Martínez, and C. David Sherrill.  Tensor hypercontraction. II. Least-squares renormalization.  Journal of Chemical Physics, 137(22):1–12, 2012.
  * [19] Robert J. Harrison, George I. Fann, Takeshi Yanai, Zhihua Gan, and Gregory Beylkin.  Madness: A multiresolution, adaptive numerical environment for scientific simulation.  SIAM Journal on Scientific Computing, 28(6):2461–2486, 2006.
  * [20] L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, and R. Schneider.  Bigdft: a flexible wavelet-based first-principles code for massively parallel computations.  Journal of Chemical Physics, 129(1):014109, 2008.
  * [21] Simen Kvaal, Christian Angeli, Damien Caliste, Morten P. Gauss, Klaus W. Jacobsen, Michel H. Linscheid, and Bálint Aradi.  Mrchem: A numerical multi-resolution approach to density functional theory.  Journal of Chemical Physics, 156(15):154103, 2022.
  * [22] Xavier Andrade, David A. Strubbe, Umberto De Giovannini, Ask H. Larsen, Micael J. T. Oliveira, Jose Alberdi-Rodriguez, Alejandro Varas, Iris Theophilou, Nicole Helbig, Matthieu J. Verstraete, Lorenzo Stella, Fernando Nogueira, AlÃ¡n Aspuru-Guzik, Alberto Castro, Miguel A. L. Marques, and Angel Rubio.  Real-space grids and the octopus code as tools for the development of new simulation approaches for electronic systems.  Journal of Physics: Condensed Matter, 24(23):233202, 2012.
  * [23] Jussi Enkovaara, Carsten Rostgaard, Jens Jørgen Mortensen, Jingzhe Chen, Marcin Dulak, Lara Ferrighi, Jeppe Gavnholt, Christian Glinsvad, Ville Haikola, Heine A. Hansen, Henrik H. Kristoffersen, Mikael Kuisma, Ask H. Larsen, Lauri Lehtovaara, Mathias Ljungberg, Poul G. Moses, Janne Ojanen, Thomas Olsen, Victor Petzold, Nichols A. Romero, Jesper Stausholm-Møller, Mikkel Strange, Georgios A. Tritsaris, Marco Vanin, Michael Walter, Bjørk Hammer, Hannu Häkkinen, Georg K. H. Madsen, Risto M. Nieminen, Jens K. Nørskov, Martti Puska, Tapio T. Rantala, Jakob Schiøtz, Kristian S. Thygesen, and Karsten W. Jacobsen.  Electronic structure calculations with gpaw: a real-space implementation of the projector augmented-wave method.  Journal of Physics: Condensed Matter, 22(25):253202, 2010.
  * [24] Bikash Kanungo and Vikram Gavini.  Large-scale all-electron density functional theory calculations using an enriched finite-element basis.  Phys. Rev. B, 95:035112, Jan 2017.
  * [25] Nelson D. Rufus, Bikash Kanungo, and Vikram Gavini.  Fast and robust all-electron density functional theory calculations in solids using orthogonalized enriched finite elements.  Physical Review B, 104(8):085112, 2021.
  * [26] Sambit Das, Phani Motamarri, Vishal Subramanian, David M. Rogers, and Vikram Gavini.  Dft-fe 1.0: A massively parallel hybrid cpu-gpu density functional theory code using finite-element discretization.  Comput. Phys. Commun., 264:107801, 2021.
  * [27] Lin Lin, Wei Hu, Chao Yang, and John E. Pask.  Chebyshev polynomial filtered subspace iteration in the discontinuous galerkin method for large-scale electronic structure calculations.  Journal of Chemical Physics, 145:154101, 2016.
  * [28] Jakob S. Kottmann, Florian A. Bischoff, and Edward F. Valeev.  Direct determination of optimal pair-natural orbitals in a real-space representation: The second-order Moller-Plesset energy.  The Journal of Chemical Physics, 152(7):074105, 02 2020.
  * [29] F. Gygi.  Adaptive riemannian metric for plane-wave electronic-structure calculations.  Epl, 19(7):617–622, 1992.
  * [30] Frangois Gygi.  Electronic-structure calculations in adaptive coordinates.  Phys. Rev. B, 48(16):11692–11700, 1993.
  * [31] Franqois Gygi and Giulia Galli.  Real-space adaptive-coordinate electronic-structure calculations.  Physical Review B, 52:2229, 1995.
  * [32] A. Devenyi, K. Cho, T. A. Arias, and J. D. Joannopoulos.  Adaptive Riemannian metric for all-electron calculations A.  Phys. Rev. B, 49(19):13373, 1994.
  * [33] Jose M. Perez-Jorda.  Change of variables for fully numerical electronic-structure calculations that concentrates the grid points in the atomic-core regions and is adequate for fast Fourier transforms.  Phys. Rev. A1, 52(4):2778, 1995.
  * [34] José M. Pérez-Jordá.  Variational plane-wave calculations in adaptive coordinates.  Physical Review B - Condensed Matter and Materials Physics, 58(3):1230–1235, 1998.
  * [35] José M. Pérez-Jordá.  Variational solution of the three-dimensional Schrdinger equation usingplane waves in adaptive coordinates.  Journal of Chemical Physics, 135(20), 2011.
  * [36] José M. Pérez-Jordá.  Variational solution of Poisson’s equation using plane waves in adaptive coordinates.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 90(5):1–6, 2014.
  * [37] D. R. Hamann.  Application of adaptive curvilinear coordinates to the electronic structure of solids.  Physical Review B, 51(11):7337–7340, 1995.
  * [38] D. Hamann.  Generalized-gradient functionals in adaptive curvilinear coordinates.  Physical Review B - Condensed Matter and Materials Physics, 54(3):1568–1574, 1996.
  * [39] D. R. Hamann.  Comparison of global and local adaptive coordinates for density-functional calculations.  Physical Review B - Condensed Matter and Materials Physics, 63(7):1–6, 2001.
  * [40] D. R. Hamman.  Band structure in adaptive curvilinear coordinates.  Physical Review B, 51(15):9508–9514, 1995.
  * [41] N A Modine, Gil Zumbach, and Efthimios Kaxiras.  Adaptive-coordinate real-space electronic-structure calculations for atoms, molecules, and solids.  Technical report, 1997.
  * [42] Gil Zumbach, N A Modine, and Efthimios Kaxiras.  Adaptive coordinate, real-space electronic structure calculations on parallel computers.  Solid State Communications, 99(2):57–61, 1996.
  * [43] Juan I. Rodríguez, David C. Thompson, Paul W. Ayers, and Andreas M. Köster.  Numerical integration of exchange-correlation energies and potentials using transformed sparse grids.  Journal of Chemical Physics, 128(22):0–10, 2008.
  * [44] Steven R. White.  Hybrid grid/basis set discretizations of the Schrödinger equation.  Journal of Chemical Physics, 147(24), 2017.
  * [45] Steven R. White and E. Miles Stoudenmire.  Multisliced gausslet basis sets for electronic structure.  Physical Review B, 99(8), 2019.
  * [46] Yiheng Qiu and Steven R. White.  Hybrid gausslet/Gaussian basis sets.  The Journal of Chemical Physics, 155(18):184107, 11 2021.
  * [47] Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven, and Garnet Kin-Lic Chan.  Low depth quantum simulation of electronic structure.  Physical Review X, 8(1):011044, 2018.
  * [48] Jarrod R. McClean, Fabian M. Faulstich, Qinyi Zhu, Bryan O’Gorman, Yiheng Qiu, Steven R. White, Ryan Babbush, and Lin Lin.  Discontinuous galerkin discretization for quantum simulation of chemistry.  New Journal of Physics, 22(9):093015, 2020.
  * [49] Steven R. White and Michael J. Lindsey.  Nested gausslet basis sets.  The Journal of Chemical Physics, 159(23):234112, 12 2023.
  * [50] John P. Boyd.  Chebyshev and Fourier Spectral Methods.  Dover Publications, Mineola, New York, second revised edition edition, 2001.
  * [51] François Gygi.  All-Electron Plane-Wave Electronic Structure Calculations.  Journal of Chemical Theory and Computation, 19(4):1300–1309, 2023\.
  * [52] Cédric Villani.  Topics in Optimal Transportation, volume 58 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, Rhode Island, March 2003.
  * [53] Francois Gygi.  Ab initio molecular dynamics in adaptive coordinates.  Technical report, 1995.
  * [54] Murray Rosenblatt.  Remarks on a Multivariate Transformation.  The Annals of Mathematical Statistics, 23(3):470 – 472, 1952.
  * [55] Herbert Knothe.  Contributions to the theory of convex bodies.  Michigan Mathematical Journal, 4(1):39 – 52, 1957.
  * [56] Michael Lindsey and Yanir A. Rubinstein.  Optimal transport via a monge–ampÃšre optimization problem.  SIAM Journal on Mathematical Analysis, 49(4):3073–3124, 2017.
  * [57] John P Boyd.  A fast algorithm for chebyshev, fourier, and sinc interpolation onto an irregular grid.  Journal of Computational Physics, 103(2):243–257, 1992.
  * [58] A. Dutt and V. Rokhlin.  Fast fourier transforms for nonequispaced data.  SIAM Journal on Scientific Computing, 14(6):1368–1393, 1993.
  * [59] Alexander H. Barnett, Jeremy Magland, and Ludvig Af Klinteberg.  A Parallel Nonuniform Fast Fourier Transform Library Based on an “Exponential of Semicircle” Kernel.  SIAM Journal on Scientific Computing, 41(5):C479–C504, January 2019\.
  * [60] Y. Shih, G. Wright, J. Anden, J. Blaschke, and A. H. Barnett.  cuFINUFFT: a load-balanced GPU library for general-purpose nonuniform FFTs.  In 2021 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 688–697, Los Alamitos, CA, USA, jun 2021\. IEEE Computer Society.
  * [61] Alex H. Barnett.  Aliasing error of the $\exp(\beta\sqrt{1-z^{2}})$ kernel in the nonuniform fast Fourier transform.  Applied and Computational Harmonic Analysis, 51:1–16, March 2021\.
  * [62] Attila Szabo and Neil S. Ostlund.  Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory.  Dover Publications, Inc., Mineola, first edition, 1996.
  * [63] Richard M. Martin.  Electronic Structure: Basic Theory and Practical Methods.  Cambridge University Press, Cambridge, 2nd edition, 2020.
  * [64] Lin Lin and Jianfeng Lu.  A Mathematical Introduction to Electronic Structure Theory.  SIAM, 2019.
  * [65] Ernest R. Davidson.  The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices.  Journal of Computational Physics, 17(1):87–94, 1975.
  * [66] P. Pulay.  Convergence acceleration of iterative sequences. the case of scf iteration.  Chemical Physics Letters, 73(2):393–398, 1980.
  * [67] P. Pulay.  Improved scf convergence acceleration.  Journal of Computational Chemistry, 3(4):556–560, 1982.
  * [68] Lin Lin.  Adaptively compressed exchange operator.  Journal of Chemical Theory and Computation, 12(5):2242–2249, 2016\.
  * [69] Wei Hu, Lin Lin, Amartya S. Banerjee, Eugene Vecharynski, and Chao Yang.  Adaptively compressed exchange operator for large-scale hybrid density functional calculations with applications to the adsorption of water on silicene.  Journal of Chemical Physics, 147(14):144107, 2017.

##  Appendix A Periodic sinc basis

###  A.1 Definitions

Let $\boldsymbol{L}=(L_{1},L_{2},L_{3})$ be the unit cell side lengths, and
let $\mathcal{V}=L_{1}L_{2}L_{3}$ denote the volume of the unit cell. Let
$\boldsymbol{m}=(m_{1},m_{2},m_{3})$ denote the numbers of grid points per
dimension. We assume that these numbers are odd in order to construct the
periodic sinc basis. To index the Fourier basis, it is useful then to define
the integers $\ell_{a}=(m_{a}-1)/2$ for $a=1,2,3$. We also let
$M=m_{1}m_{2}m_{3}$ denote the total number of basis functions and
$\Delta\mathcal{V}=\mathcal{V}/M$ the discrete volume element.

For the purposes of this discussion, let us identify the indices $i,j$ for the
grid points / basis functions with the multi-indices
$\boldsymbol{i},\boldsymbol{j}\in\prod_{a=1}^{3}\\{0,\ldots,m_{a}-1\\}$.

Then the orthonormal periodic sinc basis functions can be expressed

$\phi_{\boldsymbol{j}}(\boldsymbol{y})=\frac{1}{\sqrt{M}}\sum_{\boldsymbol{k}}e^{-2\pi
i\boldsymbol{k}\cdot(\boldsymbol{j}/\boldsymbol{m})}e_{\boldsymbol{k}}(\boldsymbol{y}),$
(A.1)

where vector division is performed entrywise and $e_{\boldsymbol{k}}$ are the
orthonormal Fourier modes on our unit cell, defined by

$e_{\boldsymbol{k}}(\boldsymbol{y})=\frac{1}{\sqrt{\mathcal{V}}}\,e^{2\pi
i\boldsymbol{k}\cdot(\boldsymbol{y}/\boldsymbol{L})},\quad\boldsymbol{k}\in\prod_{a=1}^{3}\\{-\ell_{a},\ldots,\ell_{a}\\}.$

Note the indexing convention for $\boldsymbol{k}$.

###  A.2 Kinetic matrix multiplication

Here we retain the multi-indexing convention of Appendix A.1, to which we
refer the reader for the detailed specification of the sinc basis set
$\phi_{\boldsymbol{j}}$.

Now we concern ourselves with demonstrating how to perform matrix-vector
multiplications by the kinetic matrix
$\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}}=-\int\psi_{\boldsymbol{j}}(\boldsymbol{x})\Delta\psi_{\boldsymbol{j}^{\prime}}(\boldsymbol{x})\,dx$.
Substituting $\psi_{\boldsymbol{j}}:=\phi_{\boldsymbol{j}}\circ T$, we
compute:

$\displaystyle\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}}$
$\displaystyle=-\int\psi_{\boldsymbol{j}}(\boldsymbol{x})\Delta\psi_{\boldsymbol{j}^{\prime}}(\boldsymbol{x})\,d\boldsymbol{x}$
$\displaystyle=\int\nabla\psi_{\boldsymbol{j}}(\boldsymbol{x})\cdot\nabla\psi_{\boldsymbol{j}^{\prime}}(\boldsymbol{x})\,d\boldsymbol{x}$
$\displaystyle=\int\left[\nabla\phi_{\boldsymbol{j}}\circ
T\right](\boldsymbol{x})\,\left[DT(\boldsymbol{x})\,DT(\boldsymbol{x})^{\top}\right]\,\left[\nabla\phi_{\boldsymbol{j}^{\prime}}\circ
T\right](\boldsymbol{x})\,d\boldsymbol{x}$
$\displaystyle=\int\nabla\phi_{\boldsymbol{j}}(\boldsymbol{y})\,\left[\frac{K(\boldsymbol{y})K(\boldsymbol{y})^{\top}}{J(\boldsymbol{y})}\right]\,\nabla\phi_{\boldsymbol{j}^{\prime}}(\boldsymbol{y})\,d\boldsymbol{y}$
$\displaystyle=\int\nabla\phi_{\boldsymbol{j}}(\boldsymbol{y})\,\Sigma(\boldsymbol{y})\,\nabla\phi_{\boldsymbol{j}^{\prime}}(\boldsymbol{y})\,d\boldsymbol{y},$
(A.2)

where we have defined $K(\boldsymbol{y}):=[DT^{-1}(\boldsymbol{y})]^{-1}$ and
in turn
$\Sigma(\boldsymbol{y}):=\frac{K(\boldsymbol{y})K(\boldsymbol{y})^{\top}}{J(\boldsymbol{y})}$.

Next substitute the Fourier expansions (A.1) for
$\phi_{\boldsymbol{j}}=\overline{\phi_{\boldsymbol{j}}}$ and
$\phi_{\boldsymbol{j}^{\prime}}$ into (A.2), using the identity $\nabla
e_{\boldsymbol{k}}(\boldsymbol{y})=2\pi
i\frac{\boldsymbol{k}}{\boldsymbol{L}}e_{\boldsymbol{k}}(\boldsymbol{y})$, to
obtain:

$\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}}=\frac{(2\pi)^{2}}{M}\sum_{\boldsymbol{k}}e^{2\pi
i\boldsymbol{k}\cdot(\boldsymbol{j}/\boldsymbol{m})}\left[\frac{\boldsymbol{k}}{\boldsymbol{L}}\right]^{\top}\sum_{\boldsymbol{k}^{\prime}}\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}\left[\frac{\boldsymbol{k}^{\prime}}{\boldsymbol{L}}\right]e^{-2\pi
i\boldsymbol{k}^{\prime}\cdot(\boldsymbol{j}^{\prime}/\boldsymbol{m})},$ (A.3)

where

$\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}:=\int\overline{e_{\boldsymbol{k}}(\boldsymbol{y})}\Sigma(\boldsymbol{y})e_{\boldsymbol{k}^{\prime}}(\boldsymbol{y})\,d\boldsymbol{y}.$

Now we will rewrite the integral $d\boldsymbol{y}$ as a Riemann sum over the
real space grid points
$\boldsymbol{y}_{\boldsymbol{j}}=\left(\frac{j_{a}}{m_{a}}L_{a}\right)$:

$\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}\approx\mathcal{V}\,\sum_{\boldsymbol{i}}e^{-2\pi
i\boldsymbol{k}\cdot(\boldsymbol{i}/\boldsymbol{m})}\,\Sigma(y_{\boldsymbol{i}})\,M^{-1}e^{2\pi
i\boldsymbol{k}^{\prime}\cdot(\boldsymbol{i}/\boldsymbol{m})}.$ (A.4)

Note that $\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}$ is $(3\times
3)$-matrix-valued. We have split the volume element
$\Delta\mathcal{V}=\mathcal{V}/M$ into pieces to more naturally suggest
forward and inverse discrete Fourier transforms.

Now we are ready to describe how to perform a matrix vector product of
$(\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}})$ by
$(u_{\boldsymbol{j}^{\prime}})$. Observe first that from (A.3), we have

$\sum_{\boldsymbol{j}^{\prime}}\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}}v_{\boldsymbol{j}^{\prime}}=\frac{(2\pi)^{2}}{M}\sum_{\boldsymbol{k}}e^{2\pi
i\boldsymbol{k}\cdot(\boldsymbol{j}/\boldsymbol{m})}\left[\frac{\boldsymbol{k}}{\boldsymbol{L}}\right]^{\top}\sum_{\boldsymbol{k}^{\prime}}\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}\left[\frac{\boldsymbol{k}^{\prime}}{\boldsymbol{L}}\right]\hat{u}_{\boldsymbol{k}^{\prime}},$

where (subject to suitable indexing) $\hat{u}$ can be obtained precisely as an
FFT of $u$. (We have summed over $\boldsymbol{j}^{\prime}$.)

Then we form the 3-vector-valued Fourier grid function

$\left(\hat{v}_{\boldsymbol{k}^{\prime}}^{b}\right)=\left(\left[k^{\prime}_{b}/L_{b}\right]\hat{u}_{\boldsymbol{k}^{\prime}}\right)$

by pointwise operations over the Fourier grid.

In turn, to compute
$\hat{w}_{\mathbf{k}}:=\sum_{\boldsymbol{k}^{\prime}}\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}\hat{v}_{\boldsymbol{k}^{\prime}}$,
defined elementwise as

$\hat{w}_{\boldsymbol{k}}^{a}=\sum_{b,\boldsymbol{k}^{\prime}}\hat{\Sigma}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{ab}\hat{v}_{\boldsymbol{k}^{\prime}}^{b}$

for $a=1,2,3$, we insert into (A.4). Evidently this operation can be achieved
with one inverse FFT per coordinate $b=1,2,3$, followed by a pointwise
$(3\times 3)\times(3\times 1)$ matrix-vector multiplication at each real space
grid point, followed by one forward FFT per coordinate $a=1,2,3$.

Finally, we can recover

$\sum_{\boldsymbol{j}^{\prime}}\mathbf{L}_{\boldsymbol{j}\boldsymbol{j}^{\prime}}v_{\boldsymbol{j}^{\prime}}=\frac{(2\pi)^{2}}{M}\sum_{\boldsymbol{k}}e^{2\pi
i\boldsymbol{k}\cdot(\boldsymbol{j}/\boldsymbol{m})}\left[(\boldsymbol{k}/\boldsymbol{L})\cdot\hat{w}_{\boldsymbol{k}^{\prime}}\right]$

with more pointwise operations on the Fourier grid, followed by one final
inverse FFT.

##  Appendix B Knothe-Rosenblatt computations

###  B.1 Periodic extension

We assume that $\rho$ extends smoothly and periodically from the box
$\mathcal{B}$ to $\mathbb{R}^{d}$, and we want to show that the Knothe-
Rosenblatt transport $T:\mathcal{B}\rightarrow\mathcal{B}$ extends smoothly
and periodically—in the sense that the coordinate functions of $T-\mathrm{Id}$
are all periodic—to a map $\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that
preserves each unit cell.

To see this, simply adopt the same definitions (IV.3), (IV.4), (IV.5) for
$T_{1},T_{2},T_{3}$, but extend their domains to $\mathbb{R}^{d}$. Then
evidently the map $T$ is a smooth map
$\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. It remains to demonstrate that the
extension preserves the boundary of all unit cells.

First it is useful to note that the functions $\rho_{1}(\,\cdot\,)$,
$\rho_{2}(\,\cdot\,|x_{1})$, etc. are respectively $L_{1}$-periodic,
$L_{2}$-periodic, etc. To see this first observe that the marginals
$\rho_{1}(x$

For $k_{1}\in\mathbb{Z}$,

$\displaystyle T_{1}(a_{1}+k_{1}L_{1})$
$\displaystyle=a_{1}+L_{1}\int_{a_{1}}^{a_{1}+k_{1}L_{1}}\rho_{1}(y_{1})\,dy_{1}$
$\displaystyle=a_{1}+L_{1}\sum_{j_{1}=0}^{k_{1}-1}\int_{a_{1}+j_{1}L_{1}}^{a_{1}+(j_{1}+1)L_{1}}\rho_{1}(y_{1})\,dy_{1}$
$\displaystyle\overset{(\star)}{=}a_{1}+L_{1}\sum_{j_{1}=0}^{k_{1}-1}\int_{a_{1}}^{a_{1}+L_{1}}\rho_{1}(y_{1})\,dy_{1}$
$\displaystyle=a_{1}+k_{1}L_{1}\int_{a_{1}}^{a_{1}+L_{1}}\rho_{1}(y_{1})\,dy_{1}$
$\displaystyle=a_{1}+k_{1}L_{1},$

where the $(\star)$ step follows from the fact that $\rho_{1}$ is
$L_{1}$-periodic.

Similarly, for $k_{2}\in\mathbb{Z}$,

$\displaystyle T_{2}(x_{1},a_{2}+k_{2}L_{2})$
$\displaystyle=a_{2}+L_{2}\int_{a_{2}}^{a_{2}+k_{2}L_{2}}\rho_{2}(y_{2}|x_{1})\,dy_{2}$
$\displaystyle=a_{2}+L_{2}\sum_{j_{2}=0}^{k_{2}-1}\int_{a_{2}+j_{2}L_{2}}^{a_{2}+(j_{2}+1)L_{2}}\rho_{2}(y_{2}|x_{1})\,dy_{2}$
$\displaystyle=a_{2}+L_{2}\sum_{j_{2}=0}^{k_{2}-1}\int_{a_{2}}^{a_{2}+L_{2}}\rho_{2}(y_{2}|x_{1})\,dy_{2}$
$\displaystyle=a_{2}+k_{2}L_{2}.$

Note that we have used the fact that $\rho_{2}(\,\cdot\,|x_{1})$ is
$L_{2}$-periodic.

Extrapolating the same argument to $T_{3}$, etc., establishes that $T$
preserves the boundary of all unit cells.

###  B.2 Explicit solution for a sum of separable functions

Suppose $\rho$ has the functional form

$\rho(\boldsymbol{x})=\sum_{\alpha}c^{\alpha}g^{\alpha}(\boldsymbol{x}),$

where each term
$g^{\alpha}(\boldsymbol{x})=\prod_{i=1}^{d}g_{i}^{\alpha}(x_{i})$ is
separable. Let
$G_{i}^{\alpha}(x_{i})=\int_{a_{i}}^{x_{i}}g_{i}^{\alpha}(y_{i})\,dy_{i}$ be
the appropriate antiderivative, and assume that we can evaluate
$g_{i}^{\alpha},G_{i}^{\alpha}$ easily. We will show how to compute the
Knothe-Rosenblatt transport in terms of such evaluations.

Since $\rho$ must be a density on $\mathcal{B}$, we must have that

$1=\int_{\mathcal{B}}\rho\,d\boldsymbol{x}=\sum_{\alpha}c^{\alpha}\prod_{i=1}^{d}G_{i}^{\alpha}(b_{i})=1.$

Accordingly define the univariate masses

$m_{i}^{\alpha}=G_{i}^{\alpha}(b_{i})=\int_{a_{i}}^{b_{i}}g_{i}^{\alpha}\,dx_{i}$

and

$m^{\alpha}:=\prod_{i=1}^{d}m_{i}^{\alpha}=\int_{\mathcal{B}}g^{\alpha}\,d\boldsymbol{x},$

so we need

$\sum_{\alpha}m^{\alpha}c^{\alpha}=1.$

Therefore, given unnormalized $\tilde{c}^{\alpha}$, we can define normalized
weights

$c^{\alpha}=\frac{\tilde{c}^{\alpha}}{M},\quad\text{where
}M:=\sum_{\alpha}m^{\alpha}\tilde{c}^{\alpha},$

ensuring that $\rho$ is a valid probability density.

It is further useful to define

$\omega_{i}^{\alpha}:=\prod_{j=i+1}^{d}m_{i}^{\alpha}.$

We want to evaluate the expressions above for $T_{k}$. To begin, the marginal
$\rho_{1}$ can be computed as

$\rho_{1}(x)=\sum_{\alpha}c^{\alpha}\omega_{1}^{\alpha}g_{1}^{\alpha}(x).$

Then compute:

$\displaystyle T_{1}(x_{1})$
$\displaystyle=a_{1}+L_{1}\sum_{\alpha}c^{\alpha}\omega_{1}^{\alpha}G_{1}^{\alpha}(x_{1}).$

Next, simplify:

$\displaystyle T_{2}(x_{1},x_{2})$
$\displaystyle=a_{2}+L_{2}\int_{a_{2}}^{x_{2}}\rho_{2}(y_{2}|x_{1})\,dy_{2}$
$\displaystyle=a_{2}+\frac{L_{2}}{\rho_{1}(x_{1})}\int_{a_{2}}^{x_{2}}\rho_{2}(x_{1},y_{2})\,dy_{2},$

and compute

$\rho_{2}(x_{1},x_{2})=\sum_{\alpha}c^{\alpha}\omega_{2}^{\alpha}g_{1}^{\alpha}(x_{1})g_{2}^{\alpha}(x_{2}).$

Therefore

$T_{2}(x_{1},x_{2})=a_{2}+\frac{L_{2}}{\rho_{1}(x_{1})}\left[\sum_{\alpha}c^{\alpha}\omega_{2}^{\alpha}g_{1}^{\alpha}(x_{1})G_{2}^{\alpha}(x_{2})\right],$

or

$T_{2}(x_{1},x_{2})=a_{2}+L_{2}\frac{\sum_{\alpha}c^{\alpha}\omega_{2}^{\alpha}g_{1}^{\alpha}(x_{1})G_{2}^{\alpha}(x_{2})}{\rho_{1}(x_{1})}.$

It is clear by inspection that in fact $T_{2}(x_{1},a_{2})=a_{2}$ and
$T_{2}(x_{1},b_{2})=b_{2}$, so the boundary condition is satisfied.

Similarly

$T_{3}(x_{1},x_{2},x_{3})=a_{3}+L_{3}\frac{\sum_{\alpha}c^{\alpha}\omega_{3}^{\alpha}g_{1}^{\alpha}(x_{1})g_{2}^{\alpha}(x_{2})G_{3}^{\alpha}(x_{3})}{\rho_{2}(x_{1},x_{2})},$

etc.

###  B.3 Chebyshev polynomial case

Consider $\alpha$ as a multi-index $\alpha=(n_{1},\ldots,n_{d})$ and define

$g_{i}^{\alpha}=P_{n_{i}}\circ u_{i},$

where $P_{n}$ is the $n$-th Chebyshev polynomial, and $u_{i}$ is the
increasing linear map that sends $a_{i}\mapsto-1$ and $b_{i}\mapsto 1$, i.e.,

$u_{i}(t)=\frac{2}{b_{i}-a_{i}}t-\frac{a_{i}+b_{i}}{b_{i}-a_{i}}=\frac{t-\frac{a_{i}+b_{i}}{2}}{\frac{b_{i}-a_{i}}{2}}.$

We will show in this case how to explicitly evaluate the Knothe-Rosenblatt
transport computed in Section B.2. By the results of that section, it suffices
to derive an expression for the antiderivative $G_{i}^{\alpha}$.

Now

$\int P_{n}(u)\,du=\begin{cases}u,&n=0,\\\ \frac{u^{2}}{2},&n=1,\\\
\frac{1}{2(n+1)}P_{n+1}(u)-\frac{1}{2(n-1)}P_{n-1}(u),&n\geq 2,\end{cases}$

and

$\int P_{n}(u(t))\,dt=\int P_{n}(u)\,\frac{dt}{du}\,du.$

Moreover $\frac{du_{i}}{dt}\equiv\frac{2}{b_{i}-a_{i}}$, so we can define an
antiderivative

$\tilde{G}_{i}^{\alpha}(t)=\frac{b_{i}-a_{i}}{2}\begin{cases}u_{i}(t),&n_{i}=0,\\\
\frac{u_{i}(t)^{2}}{2},&n_{i}=1,\\\
\frac{1}{2(n_{i}+1)}P_{n_{i}+1}(u_{i}(t))-\frac{1}{2(n_{i}-1)}P_{n_{i}-1}(u_{i}(t)),&n_{i}\geq
2.\end{cases}$

in terms of which our desired antiderivative is

$G_{i}^{\alpha}(t)=\tilde{G}_{i}^{\alpha}(t)-\tilde{G}_{i}^{\alpha}(a_{i}).$

###  B.4 Fourier mode case

Consider $\alpha$ as a multi-index $\alpha=(n_{1},\ldots,n_{d})$ and define

$g_{i}^{\alpha}=e_{n_{i}}\circ u_{i},$

where $e_{n}(x)=e^{-2\pi inx}$ is the $n$-th Fourier mode polynomial, and
$u_{i}$ is the increasing linear map that sends $a_{i}\mapsto 0$ and
$b_{i}\mapsto 1$, i.e.,

$u_{i}(t)=\frac{t-a_{i}}{b_{i}-a_{i}}.$

We will repeat the computation of Section B.3 in this case.

Now

$\int e_{n}(u)\,du=\begin{cases}u,&n=0,\\\ \frac{i}{2\pi n}e_{n}(u)&n\neq
0.\end{cases}$

and

$\int P_{n}(u(t))\,dt=\int P_{n}(u)\,\frac{dt}{du}\,du.$

Moreover $\frac{du_{i}}{dt}\equiv\frac{1}{b_{i}-a_{i}}$, so we can define an
antiderivative

$\tilde{G}_{i}^{\alpha}(t)=(b_{i}-a_{i})\begin{cases}u_{i}(t),&n_{i}=0,\\\
\frac{i}{2\pi n_{i}}e_{n}(u_{i}(t))&n_{i}\neq 0.\end{cases}$

in terms of which our desired antiderivative is

$G_{i}^{\alpha}(t)=\tilde{G}_{i}^{\alpha}(t)-\tilde{G}_{i}^{\alpha}(a_{i})=(b_{i}-a_{i})\begin{cases}u_{i}(t),&n_{i}=0,\\\
\frac{i}{2\pi n_{i}}\left[e_{n}(u_{i}(t))-1\right]&n_{i}\neq 0.\end{cases}$

##  Appendix C Pseudospectral computations

In this section we prove the claim made in Section V.2 that the defining
property (V.2) for the matrix elements of an operator $\hat{O}$ holds for
$O_{ij}$ as in (V.3).

First observe that for arbitrary suitably smooth $g$, we have via (III.2) that
$\langle\eta_{i},g\rangle=w\,{g(\boldsymbol{x}_{j})}/{\sqrt{J_{i}}}$. Then by
replacing $g$ with $g/\sqrt{J}$ in this identity and shifting a factor of
$1/\sqrt{J}$ within the inner product, it follows that

$\langle\psi_{i},g\rangle=w\,{g(\boldsymbol{x}_{i})}/{{J_{i}}}$ (C.1)

for all suitably smooth $g$. Hence for such $g$, we can relate the two
relevant inner products via

$\langle\eta_{i},g\rangle=\sqrt{J_{i}}\,\langle\psi_{i},g\rangle,$ (C.2)

as if we were effectively identifying
$\eta_{i}(\boldsymbol{x})=\psi_{i}(\boldsymbol{x})\,\sqrt{J(\boldsymbol{x})}$
with $\psi_{i}(\boldsymbol{x})\,\sqrt{J_{i}}$.

Then for an expansion
$f(\boldsymbol{x})=\sum_{j}f_{j}\eta_{j}(\boldsymbol{x})$ of suitably smooth
$f$, we can alternatively write, following (III.1):
$f(\boldsymbol{x})=\sum_{j}f_{j}\sqrt{J_{j}}\,\psi_{j}(\boldsymbol{x})$.
Therefore
$\hat{O}f(\boldsymbol{x})=\sum_{j}f_{j}\sqrt{J_{j}}\,\hat{O}\psi_{j}(\boldsymbol{x})$.
Then under the assumption that $\hat{O}f$ is suitably smooth, by applying
(C.2) with $g=\hat{O}f$ we obtain
$\langle\eta_{i},\hat{O}f\rangle=\sqrt{J_{i}}\,\langle\psi_{i},\hat{O}f\rangle=\sum_{j}f_{j}\sqrt{J_{i}J_{j}}\,\langle\psi_{i},\hat{O}\psi_{j}\rangle$,
i.e., the defining property (V.2) holds for $O_{ij}$ as in (V.3).

##  Appendix D Pseudopotential computations

In this section we discuss how to evaluate $V_{I}(\boldsymbol{x})$ using the
all-electron pseudopotential of Gygi[51] in a periodic context.

The form of the potential is

$\displaystyle V_{\mathrm{PP}}(r)=$
$\displaystyle-\frac{1}{2}-\frac{\mathrm{erf}(ar)-2\left(a^{2}b+a/\sqrt{\pi}\right)re^{-a^{2}r^{2}}}{r}+\frac{\left[\mathrm{erf}(ar)-2\left(a^{2}b+a/\sqrt{\pi}\right)re^{-a^{2}r^{2}}\right]^{2}}{2}$
$\displaystyle\quad\quad+\
\frac{\left[-2a^{2}b-4a/\sqrt{\pi}+\left(4a^{4}b+4a^{3}/\sqrt{\pi}\right)r^{2}\right]e^{-a^{2}r^{2}}}{2},$
(D.1)

where $r$ is the distance from the nucleus and $a,b$ are parameters which can
be varied continuously to increase the accuracy of the potential. The PP (D.1)
is directly applicable for the hydrogen atom. To obtain the PP for more
general nuclear charge $Z$, one defines the expression

$V_{\mathrm{PP},Z}(r)=Z^{2}V_{\mathrm{PP}}(Zr).$

In this work we primarily consider $a=4$ and $b=-0.10200558466$, though we
consider the impact of varying the parameter $a$ in some of our experiments.
These choices are expected to yield energy error less than a few tens of $\mu
E_{h}$ for molecules containing elements in the first two rows of periodic
table. In terms of radially dependent functional form $V_{\mathrm{PP},Z}(r)$,
we define the PP as

$V_{I}(\boldsymbol{x})=V_{\mathrm{PP},Z_{I}}(|\boldsymbol{x}-R_{I}|).$ (D.2)

Since we perform periodic calculations in this work, we actually need to
consider the periodized pseudopotential. However simple summation in the real
space diverges and in practice one does a summation in the reciprocal space
and the divergent $G=0$ term is dropped. In other words the periodized
potential is

$U_{I}(\boldsymbol{x})=\sum_{\boldsymbol{G}\neq
0}\hat{V}_{I}(\boldsymbol{G})e^{-i\boldsymbol{G}\cdot\boldsymbol{x}},$ (D.3)

where $\boldsymbol{G}$ have the form
$n_{1}\boldsymbol{G}_{1}+n_{2}\boldsymbol{G}_{2}+n_{3}\boldsymbol{G}_{3}$ of
integer multiples of reciprocal lattice basis vectors, satisfying the
condition that $\boldsymbol{G}_{i}\cdot\boldsymbol{L}_{j}=2\pi\delta_{ij}$
where $\boldsymbol{L}_{1}$, $\boldsymbol{L}_{2}$, $\boldsymbol{L}_{3}$ are the
lattice basis vectors in real space.

In principle one can evaluate the Fourier transform of the PP and evaluate the
sum in (D.3) using NUFFTs to calculate the periodized potential
$U(\boldsymbol{x})$ at the grid points on the deformed coordinates. However,
this is expensive because despite the fact that the PP is band-limited the
number of $\boldsymbol{G}$ vectors needed will be extremely large. To avoid
this we use the Ewald summation trick whereby we split the summation into
short real and reciprocal space summations. This can be done as follows

$\displaystyle U_{I}(\boldsymbol{x})$ $\displaystyle=\sum_{\boldsymbol{G}\neq
0}\hat{V}_{I}(\boldsymbol{G})e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}$
$\displaystyle=\sum_{\boldsymbol{G}\neq
0}\left(\hat{V}_{I}(\boldsymbol{G})-\hat{V}_{1}(\boldsymbol{G})\right)e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}+\sum_{\boldsymbol{G}\neq
0}\hat{V}_{1}(\boldsymbol{G})e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}$
$\displaystyle=\sum_{\boldsymbol{G}}\left(\hat{V}_{I}(\boldsymbol{G})-\hat{V}_{1}(\boldsymbol{G})\right)e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}+\left(\hat{V}_{I}(0)-\hat{V}_{1}(0)\right)+\sum_{\boldsymbol{G}\neq
0}\hat{V}_{1}(\boldsymbol{G})e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}$
$\displaystyle=\sum_{\boldsymbol{L}}\left(V_{I}(\boldsymbol{x}+\boldsymbol{L})-V_{1}(\boldsymbol{x}+\boldsymbol{L})\right)+\sum_{\boldsymbol{G}\neq
0}\hat{V}_{1}(\boldsymbol{G})e^{-i\boldsymbol{G}\cdot\boldsymbol{x}}-\left(\hat{V}_{I}(0)-\hat{V}_{1}(0)\right)$

where ${V}_{1}(\boldsymbol{x})$ is defined by replacing $V_{\mathrm{PP},Z}(r)$
with the softened Coulomb potential $Z\frac{\mathrm{erf}(\omega r)}{r}$. The
parameter $\omega$ is chosen small enough such that the softened potential is
fairly smooth and the summation in the second term on the last line converges
quickly. Meanwhile, $\hat{V}_{I}$ and $\hat{V}_{1}$ are respectively the
Fourier transform of the PP and softened Coulomb potential. In going from the
third to the fourth line we have made use of the Poisson summation formula.

Both $V_{1}(\boldsymbol{r})$ and $V_{I}(\boldsymbol{r})$ converge to the usual
Coulomb potential with increasing $\boldsymbol{r}$ and thus their difference
goes to zero rapidly with increasing $\boldsymbol{r}$, which implies that the
first summation in the last line of the equation also rapidly converges.

In the last term both
$\hat{V}_{I}(0)=4\pi\int_{0}^{\infty}V_{\mathrm{PP},Z}(r)r^{2}dr$ and
$\hat{V}_{1}(0)=4\pi\int_{0}^{\infty}V_{1}(r)r^{2}dr$ individually diverge,
but their difference is bounded. We evaluate it by taking the difference
$4\pi(\int_{0}^{R}V_{\mathrm{PP},Z}(r)r^{2}dr-\int_{0}^{R}V_{1}(r)r^{2}dr)$
for a sufficiently large $R$ such that beyond this $R$ the two potentials are
equal and cancel each other. The two terms are

$\displaystyle\int_{0}^{R}V_{\mathrm{PP}}(r)r^{2}dr=$
$\displaystyle\frac{1}{96\pi a^{3}}\left[\sqrt{\pi}\sqrt{2}\left(9\pi
a^{2}b^{2}+78\sqrt{\pi}ab+49\right)\text{erf}\left(\sqrt{2}aR\right)+16\pi
a^{3}R^{3}\left(\text{erf}(aR)^{2}-1\right)\right.$ $\displaystyle-48\pi
a^{3}R^{2}\text{erf}(aR)-4aRe^{-2a^{2}R^{2}}\left(9\pi
a^{2}b^{2}+30\sqrt{\pi}ab+17\right)$
$\displaystyle-48a^{3}R^{3}e^{-2a^{2}R^{2}}\left(\sqrt{\pi}ab+1\right)^{2}$
$\displaystyle\left.-32\sqrt{\pi}e^{-a^{2}R^{2}}\left(\left(3\sqrt{\pi}ab+2\right)\left(a^{2}R^{2}+1\right)\text{erf}(aR)+3a^{4}R^{3}\left(\sqrt{\pi}ab+1\right)\right)\right]$
$\displaystyle\int_{0}^{R}V_{1}(r)r^{2}dr=$
$\displaystyle-\frac{Z}{4}\left[2e^{-R^{2}\omega^{2}}\frac{R}{\omega\sqrt{\pi}}+\left(2R^{2}-\frac{1}{\omega^{2}}\right)\mathrm{erf}(\omega
R)\right].$

The first term is only directly applicable for hydrogen atom, and the general
case can be recoverd via
$\int_{0}^{R}V_{\mathrm{PP},Z}(r)r^{2}dr=\frac{1}{Z}\int_{0}^{RZ}V_{\mathrm{PP}}(r)r^{2}dr$.

These expressions can be used to evaluate the PP $V_{I}(\boldsymbol{x}_{i})$
on our deformed grid for any given parameters $a$ and $b$ defining the PP. We
comment that these evaluations only have to be performed once at the beginning
of the calculation.

##  Appendix E Bare nuclear computations

Our goal in this section is to compute $V_{I,i}$ as defined in (V.11). Expand
this definition to obtain

$V_{I,i}=-\frac{4\pi
Z_{I}J_{i}}{w}\int\frac{1}{4\pi|\boldsymbol{x}-R_{I}|}\,\psi_{i}(\boldsymbol{x})\,d\boldsymbol{x}=4\pi
Z_{I}\,u_{i}(\boldsymbol{x}_{I}),$ (E.1)

where $u_{i}$ is defined to be the solution of the Poisson equation

$-\Delta
u_{i}(\boldsymbol{x})=f_{i}(\boldsymbol{x}):=-\frac{J_{i}}{w}\,\psi_{i}(\boldsymbol{x}).$

Recall from (V.7) that the solution $u_{i}$ satisfies

${J_{k}}\sum_{j}\mathbf{L}_{kj}\,u_{i}(\boldsymbol{x}_{j})=f_{i}(\boldsymbol{x}_{k})=-\frac{J_{i}}{w^{2}}\delta_{ik},$

where we have used the property that
$\psi_{i}(\boldsymbol{x}_{k})=w^{-1}\delta_{ik}$ in the last equality. It
follows that

$u_{i}(\boldsymbol{x}_{j})=-\frac{1}{w^{2}}\,\mathbf{L}^{-1}_{ij}.$

Then the solution $u_{i}(\boldsymbol{x})$, evaluated at $\boldsymbol{x}_{I}$,
can be expanded in the $\\{\psi_{j}\\}$ basis:

$u_{i}(\boldsymbol{x}_{I})=w\sum_{j}u_{i}(\boldsymbol{x}_{j})\,\psi_{j}(\boldsymbol{x}_{I})=-\frac{1}{w}\sum_{j}\mathbf{L}^{-1}_{ij}\psi_{j}(\boldsymbol{x}_{I}).$
(E.2)

Then it follows by plugging our expression (E.2) for
$u_{i}(\boldsymbol{x}_{I})$ into (E.1) and summing over $I$ that

$\sum_{I}V_{I,i}=\sum_{j}\mathbf{L}_{ij}^{-1}b_{j},$

where

$b_{i}=-\frac{4\pi}{w}\sum_{I}Z_{I}\,\psi_{j}(\boldsymbol{x}_{I}).$

In conclusion $\sum_{I}V_{I,i}$ can be computed for all $i$ with a single
linear solve of the form $\mathbf{L}x=b$.