Sliced Basis Density Matrix Renormalization Group for Electronic Structure

TL;DR

This paper presents a hybrid DMRG method combining grid and Gaussian basis sets for efficient electronic structure calculations of chain-like molecules, achieving near-linear scaling for large systems.

Contribution

The paper introduces a novel hybrid approach to DMRG that combines grid and Gaussian bases, enabling scalable calculations for long chain molecules.

Findings

Linear scaling of computational time with chain length

Near-exact results for large hydrogen chains within the basis

Effective handling of long-range Coulomb interactions

Abstract

We introduce a hybrid approach to applying the density matrix renormalization group (DMRG) to continuous systems, combining a grid approximation along one direction with a finite Gaussian basis set along the remaining two directions. This approach is especially useful for chain-like molecules, where the grid is used in the long direction, and we demonstrate the approach with results for hydrogen chains. The computational time for this system scales approximately linearly with the length of the chain, as we demonstrate with minimal basis set calculations with up to 1000 atoms, which are near-exact within the basis. The linear scaling comes from the combination of localization of the basis and a compression method with controlled accuracy for the long-ranged Coulomb terms in the Hamiltonian.

Equations

$$V_{ijkl}={\int }_{{\mathbf{r}}_1}{\int }_{{\mathbf{r}}_2}\frac{{\varphi }_i({\mathbf{r}}_1){\varphi }_l({\mathbf{r}}_1){\varphi }_j({\mathbf{r}}_2){\varphi }_k({\mathbf{r}}_2)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}$$

$${\hat{H}}_{\text{el}}$$

$$\text{mbox}+\frac{1}{2}{\int }_{\mathbf{r},{\mathbf{r}}^{\prime }}\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\hat{\psi }}_{\sigma }^{†}(\mathbf{r}){\hat{\psi }}_{{\sigma }^{\prime }}^{†}({\mathbf{r}}^{\prime }){\hat{\psi }}_{{\sigma }^{\prime }}({\mathbf{r}}^{\prime }){\hat{\psi }}_{\sigma }(\mathbf{r}).$$

$$\hat{H}$$

$$+\frac{1}{2}\sum _{n{n}^{\prime }}\sum _{ijkl}{V}_{ijkl}^{n{n}^{\prime }}{\hat{c}}_{ni\sigma }^{†}{\hat{c}}_{n^{\prime }j{\sigma }^{\prime }}^{†}{\hat{c}}_{n^{\prime }k{\sigma }^{\prime }}{\hat{c}}_{nl\sigma }.$$

$$V_{ijkl}^{n{n}^{\prime }}$$

$$t_{ij}^{n{n}^{\prime }}$$

$$-{\delta }_{ij}\frac{1}{2a^2}{\Delta }_{n{n}^{\prime }}.$$

$${\hat{c}}_{nj\sigma }=\sqrt{a}{\int }_{x,y}{\varphi }_j(x,y){\hat{\psi }}_{n\sigma }(x,y,z_n).$$

$$(x^2-y^2)\exp [-\zeta (x^2+y^2)]$$

$$2xy\exp [-\zeta (x^2+y^2)].$$

$${\eta }_i^k={\int }_{x,y}{\varphi }^k(x,y,z_n){\xi }_i(x,y)$$

$${\rho }_{i{i}^{\prime }}=\sum _k{\eta }_i^k{\eta }_{i^{\prime }}^k.$$

$$\sum _{n\le n^{\prime }}{V}_{n{n}^{\prime }}{\hat{n}}_n{\hat{n}}_{n^{\prime }}.$$

$$V^{(k)}=U^{(k)}{S}^{(k)}{W}^{(k)}$$

$$U^{(k+1)}=P(U^{(k)})X^{(k+1)}.$$

$$V_{ijkl}^{n{n}^{\prime }}$$

$${\tilde{t}}_{ij}^{nn}$$

$$\text{mbox}+{\int }_{\mathbf{\rho }}{\varphi }_i(\mathbf{\rho })\left[v(\mathbf{\rho },z_n)\right]{\varphi }_j(\mathbf{\rho })$$

$$\frac{1}{r}\approx \sum _{i=1}^P{c}_i\exp (-a_i{r}^2).$$

$$\sum _{n\le n^{\prime }}{V}_{n{n}^{\prime }}{\hat{n}}_n{\hat{n}}_{n^{\prime }}$$

$$V_{n{n}^{\prime }}={\lambda }^{|n-n^{\prime }|}.$$

$$V_{n{n}^{\prime }}={\lambda }^{|n-n^{\prime }|}={\lambda }^{n^{\prime }-n}={\lambda }^{-n}{\lambda }^{n^{\prime }}(n^{\prime }\ge n).$$

$$V_{r,c}^{(p)}=V_{r,(c+p-1)}$$

$$V^{(p)}=U^{(p)}{S}^{(p)}{W}^{(p)}.$$

$$V^{(p+1)}$$

$$\displaystyle=\big{[}U^{(p)}S^{(p)}C^{-}(W^{(p)})\big{]}\oplus r^{(p+1)}$$

$$\displaystyle=P(U^{(p)})\big{[}S^{(p)}C^{-}(W^{(p)})\oplus r^{(p+1)}\big{]}.$$

$$U^{(p+1)}=P(U^{(p)})X^{(p+1)}$$

$$\hat{V}=\sum _{n\le n^{\prime }}{V}_{n{n}^{\prime }}{\hat{n}}_n{\hat{n}}_{n^{\prime }}.$$

$$\hat{V}=\sum _{\{\alpha \}}{\hat{M}}_{{\alpha }_1}^{(1)}{\hat{M}}_{{\alpha }_1{\alpha }_2}^{(2)}{\hat{M}}_{{\alpha }_2{\alpha }_3}^{(3)}\cdots {\hat{M}}_{{\alpha }_{N-1}}^{(N)}$$

$${\hat{M}}^{(1)}=\left[\begin{array}{ccc}{V}_{11}({\hat{n}}_1{)}^2& X_{11}^{(1)}{\hat{n}}_1& {\hat{I}}_1\end{array}\right],$$

$${\hat{M}}^{(2)}=\left[\begin{array}{cccc}{\hat{I}}_2& 0& 0& 0\\ {\Omega }_{11}^{(2)}{\hat{n}}_2& X_{11}^{(2)}{\hat{n}}_2& X_{12}^{(2)}{\hat{n}}_2& 0\\ V_{22}({\hat{n}}_2{)}^2& X_{21}^{(2)}{\hat{n}}_2& X_{22}^{(2)}{\hat{n}}_2& {\hat{I}}_2\end{array}\right].$$

$${\hat{M}}^{(3)}=\left[\begin{array}{ccccc}{\hat{I}}_3& 0& 0& 0& 0\\ {\Omega }_{11}^{(3)}{\hat{n}}_3& X_{11}^{(3)}{\hat{n}}_3& X_{12}^{(3)}{\hat{n}}_3& X_{13}^{(3)}{\hat{n}}_3& 0\\ {\Omega }_{21}^{(3)}{\hat{n}}_3& X_{21}^{(3)}{\hat{n}}_3& X_{22}^{(3)}{\hat{n}}_3& X_{23}^{(3)}{\hat{n}}_3& 0\\ V_{33}({\hat{n}}_3{)}^2& X_{31}^{(3)}{\hat{n}}_3& X_{32}^{(3)}{\hat{n}}_3& X_{33}^{(3)}{\hat{n}}_3& {\hat{I}}_3\end{array}\right]$$

Full text

In the last decade the density matrix renormalization group (DMRG) has become a powerful method for computing the electronic structure of molecules Chan and Sharma (2011). The now standard quantum chemistry DMRG approach (QCDMRG) works with a discrete Hamiltonian defined by an orthogonalized, contracted Gaussian basis set White and Martin (1999). For systems with strong correlation, problems of inaccuracy and poor convergence plaguing other approaches are not a serious problem for DMRG. But QCDMRG has major limitations associated with basis set size and dimensionality. Calculation times grow rapidly with the number of active basis functions, and the current practical limit is about 100-200 basis functions. In addition, there are fundamental limitations for DMRG when the transverse size of the system becomes large, which we do not try to address here.

The Hilbert space used in QCDMRG is the same as that of the Hubbard model, equating a Hubbard site with a single basis function. However, the rapid scaling of computation time with the number of basis functions in QCDMRG does not occur for a one-dimensional Hubbard model, for which the calculation time is approximately linear (when keeping a fixed number of states in DMRG). The main reason for the poor scaling of QCDMRG is the complexity of the Hamiltonian in the basis, particularly the two-electron terms. The electron-electron Coulomb interaction terms are defined by two-electron integrals

[TABLE]

where the $\phi_{i}(\mathbf{r})$ are orthonormal basis functions. If the basis functions are delocalized, as they are when using molecular orbitals from a Hartree Fock calculation, the number of significant $V_{ijkl}$ terms scales as $N_{b}^{4}$, where $N_{b}$ is the number of basis functions. This leads to a computation time for QCDMRG which scales as $N_{b}^{4}m^{2}+N_{b}^{3}m^{3}$, where $m$ is the number of many-body states kept. 111In long molecules, with minimal basis sets, truncation of the interactions can improve the scaling of QCDMRG to $O(N_{b}^{2})$.

The nonlocality of the orthogonal basis functions also increases the $m$ needed for a given accuracy. DMRG is a low-entanglement approximation, and the entanglement of ground states is governed by the area law Evenbly and Vidal (2011); Hastings (2007). The area law is a property that holds for ground states described in terms of local, “real space” degrees of freedom. In a delocalized basis, a volume law of entanglement holds instead (except for non-interacting systems, a special point where the entanglement is zero in the eigenstate basis). To capture volume-law states, $m$ must grow exponentially with the system size, even in one dimension. For this reason, some effort should be made to localize the basis before applying standard QCDMRG, except on very small molecules. The localization is always imperfect—the basis functions have oscillating tails which tend to be slowly decaying.

Hypothetically, one could get rid of both the $N_{b}^{4}$ scaling and the increase in entanglement from extended basis functions by going to a real-space grid defined by finite differences. In such a grid the interactions are defined as $V_{ij}\hat{n}_{i}\hat{n}_{j}$, where $\hat{n}_{i}$ is the density operator on site $i$. For model one-dimensional continuum systems, this is currently the most powerful approach, and we have used it to simulate systems of 100 pseudo-hydrogen atoms with about 20 grid points per atom Stoudenmire et al. (2012). A key part of using a one-dimensional grid is compressing the interactions by approximating long-range interactions as a sum of exponentials Crosswhite et al. (2008); Pirvu et al. (2010). With this compression, the calculation time grows only linearly with the number of atoms. The problem with such a grid approach for three dimensions is that the number of grid points would be very high, for example of order $10^{6}$ for a system of modest size.

Here we introduce a hybrid approach, which we call sliced basis DMRG (SBDMRG). Along one particular “z” direction we use a grid. This grid direction is chosen to be the direction over which the molecule extends furthest. At each grid point, the remaining transverse dimensions, $x$ and $y$, are captured by a small number of basis functions derived from standard Gaussian basis sets, making what we call a “slice”—see Fig. 1. The total number of DMRG “sites” is therefore $N_{b}=N_{z}N_{o}$, where $N_{z}$ is the number of grid points, and $N_{o}$ is the number of transverse functions (“orbitals”) per grid point. The DMRG path progresses through all orbitals on a slice, then moves to the next. This approach has several major advantages. First, all interaction terms $V_{ijkl}$ where $i$ and $l$ are not on the same slice are zero, and similarly for $j$ and $k$. Thus the number of terms scales as $N_{z}^{2}$. Second, the remaining interactions can be compressed very efficiently, making the dominant part of the calculation time linear in $N_{z}$. Third, since there is no spatial extent of the basis functions in the $z$ direction, there is no extra entanglement due to nonlocality, potentially reducing the number of states $m$ needed for a given accuracy.

We demonstrate our method by simulating linear chains of hydrogen atoms. Although these are three-dimensional systems, their linear nature makes them especially well suited for both SBDMRG and QCDMRG. They also exhibit strong correlation, and can be quite challenging for electronic structure methods. The electronic density in a plane through the nuclei for a typical calculation is presented in Fig. 2.

To define the sliced basis approach in detail, consider the electronic structure Hamiltonian for fixed nuclei in atomic units

[TABLE]

Summation over spin labels $\sigma$ is implied above and in what follows, and $v(\mathbf{r})$ is the single particle potential generated by the nuclei.

Along the $z$ direction, we make a grid approximation by taking $z_{n}=n\!\cdot\!a$ with $n$ an integer and $a$ a small grid spacing. Then on each slice $n$, we introduce a finite, orthonormal basis of functions $\{\phi_{j}(x,y)\}$ where $j=1,2,\ldots,N_{o}$. For simplicity, we use the same $N_{o}$ and functions $\{\phi_{j}(x,y)\}$ on every slice $n$. At a later stage one can perform a change of basis to adapt the basis for each slice, possibly reducing the number of functions. We introduce discrete operators $\hat{c}^{\dagger}_{nj\sigma}$ and $\hat{c}_{nj\sigma}$ which create and destroy electrons in a slice orbital. In terms of these operators, the discretized Hamiltonian takes the form

[TABLE]

Introducing the notation $\boldsymbol{\rho}=(x,y)$ for convenience, the interaction integrals are defined as

[TABLE]

Note that the $i,j,k,l$ indices only run over the small number of functions $N_{o}$ on each slice. Thus, the Hamiltonian is defined by just $N_{z}^{2}N_{o}^{4}$ interaction integrals. The single-particle couplings are defined to be

[TABLE]

Our discrete Hamiltonian treats the $z$-direction kinetic energy terms Eq. (7) on a different footing than the “integral” terms. For the $z$-direction kinetic energy, we treat the basis functions as being smooth functions of $z$, and think of the slices as sampling those functions. Thus we use standard finite difference formulas, defined via $\Delta_{nn^{\prime}}$. One could take a second order approximation for $\Delta$, with nonzero terms $\Delta_{nn}=-2$ and $\Delta_{n,n+1}=\Delta_{n+1,n}=1$. However, to reduce the grid error to $a^{4}$ we use a fourth-order approximation. For the “integral” terms, we think of the basis functions as being completely localized and nonoverlapping between slices, i.e. $\varphi_{nj}(\mathbf{r})=\delta^{\frac{1}{2}}(z-z_{n})\phi_{j}(x,y)$. This corresponds to taking

[TABLE]

and then transforming Eq. (2) accordingly. The distinct treatments of the terms means that the results are not strictly variational at finite $a$. However, we find finite-$a$ errors for hydrogen chains of only about 0.1 mH per atom for $a=0.1$, and in the limit of $a\to 0$, the results are variational.

In what follows, we construct the transverse basis functions on a slice $\{\phi_{j}(x,y)\}$ out of standard atom-centered Gaussian basis sets. We assume all the atoms are identical. In going from the spherical symmetry used in standard Gaussians to slices, we switch to cylindrical symmetry. Thus, an $S$-function becomes a $\sigma$ function, $P$-functions become $\pi$ functions, etc. Whereas there are $2\ell+1$ functions in a spherical set with angular momentum $\ell$, there are only two cylindrical functions for any $\ell>0$. For example, a set of $D$ functions, with coefficient $\zeta$, becomes the two slice basis functions

[TABLE]

We leave out functions like $P_{z}$, which looks like a $\sigma$ function on a slice, or any other function looking like a function of smaller $\ell$. (In principle, $P_{z}$ could be kept as an additional $S$ function.) The slice basis functions are only orthogonal between different slices. This means the functions within each slice must be orthogonalized.

In the parent 3D Gaussian bases, usually some of the functions (particularly $S$-type) are contracted, meaning out of $N_{g}$ original Gaussians, one uses a smaller number $N_{o}$ of linear combination of functions for each atom: $\phi^{j}=\sum_{m=1}^{N_{g}}c^{j}_{m}\exp[-\zeta_{m}(\vec{r}-\vec{r}_{A})^{2}]$ where $j=1\ldots N_{o}$, and $N_{o}<N_{g}$. In this case, to define the transverse basis on a slice, we follow an approach that is useful very generally: we form a local orbital density matrix for each slice. Let $i$ and $i^{\prime}$ run over an orthonormal uncontracted basis for the slice at $z_{n}$, defined by functions $\xi_{i}(x,y)$. Let $\phi^{k}(x,y,z)$ be a particular 3D contracted basis function attached to one of the atoms, and let

[TABLE]

Then let

[TABLE]

The leading eigenvectors of $\rho$ form optimal local functions for representing the contracted 3D basis. More generally, $\rho$ could come from the interacting ground state, as a block of the single particle reduced density matrix $\langle c^{\dagger}_{i}c_{i^{\prime}}\rangle$, and we would call the eigenvectors of $\rho$ “slice natural orbitals” (SNOs). A subset with only $N_{j}$ of these SNOs would be an ideal reduced local basis. Our procedure for contractions is conceptually similar to this, but with equal weighting for all 3D contracted basis functions. In this case, for example, the sharp Gaussians used to represent the nuclear cusps only appear significantly in the slices close to nuclei. In our hydrogen chain calculations, if the basis has $N_{S}$ contracted $S$ functions per atom, we keep $N_{S}$ contracted functions per slice.

We perform DMRG with the Hamiltonian represented as a sum of matrix product operators (MPOs), one of which represents the long-ranged two-electron interactions. For this MPO we use a compression technique giving an MPO with matrix dimension $D$ which is nearly independent of system length, leading to a linear scaling of the computation time. (The other MPOs, say for $v(\mathbf{r})$, are naturally of constant dimension.) Consider the simplest case of a single basis function per slice such that the interaction part of the Hamiltonian Eq. (4) simplifies to

[TABLE]

Here we will focus on the compression of the upper triangle of the matrix $V_{nn^{\prime}}$, giving just an outline; more details are given in Appendix B. Note that since the local basis varies from slice to slice, $V$ is not translationally invariant; if it was, an MPO could be constructed based on fitting $V(n-n^{\prime})$ to a sum of exponentials Crosswhite et al. (2008). We use a more general method based on a sequence of singular value decompositions (SVDs). This is a simplification of more general SVD approaches for potentially more complicated Hamiltonians Zaletel et al. (2015); Chan et al. (2016).

For a particular diagonal index $k$, let $V^{(k)}$ be the rectangular block of $V$ with the lower left corner at $V_{kk}$, and extending to the upper right corner of $V$. An SVD gives

[TABLE]

where $S^{(k)}$ is the diagonal matrix of singular values. The smoothness of $V(n-n^{\prime})$ away from the diagonal makes this SVD have a small number $D$ of significant singular values, allowing us to approximate $S^{(k)}$ as a $D\times D$ matrix, with appropriate reductions in the number of columns of $U^{(k)}$ and rows of $W^{(k)}$.

This factorized representation at index $k$ can be related to a similar representation at $k+1$. Define $P(U^{(k)})$ to be the direct sum of $U^{(k)}$ and a $1\times 1$ identity matrix, that is add an extra column and row of zeros to the bottom and right of $U^{(k)}$ and set the new diagonal element to 1. Then a matrix $X^{(k+1)}$ can be computed such that

[TABLE]

The matrix $X^{(k+1)}$ is of dimension $(D+1)\times D$. We see that we can recover all the $U^{(k)}$ if we know all the $X^{(k)}$ and $U^{(1)}$. Similarly, all of the $W^{(k)}$ can be generated in terms of a reverse recursion involving $D\times(D+1)$ matrices $Y^{(k)}$. This means we can reconstruct every $V^{(k)}$, and thus the entire $N_{z}\times N_{z}$ matrix $V$ out of the $O(N_{z}D^{2})$ parameters in $X^{(k)}$, $S^{(k)}$, and $Y^{(k)}$. In Appendix B, we detail how to compute the $X^{(k)}$ and $Y^{(k)}$ matrices, and show how they lead to an MPO representation of the interactions with MPO matrix dimension $D+2$.

In Fig. 3 we show results for chains of 10 equally-spaced hydrogen atoms as a function of separation $R$, for several different basis sets with $a=0.1$ and for comparison, standard QCDMRG results for parent 3D basis sets Zheng et al. . The STO-6G basis is a minimal basis, contracting 6 Gaussians to one function per atom; the sliced version also has one function per slice. One can see that the completeness of the standard and sliced bases are similar; which basis gives a lower energy varies with $R$. The double $\zeta$ basis (cc-pVDZ) has five functions per atom Jr. (1989), and the sliced version has four per slice (no $P_{z}$). Here the energies are even closer, but the sliced version is consistently slightly lower. The triple $\zeta$ basis (cc-pVTZ) has 14 functions per atom, or 140 functions total, making this a somewhat challenging QCDMRG calculation. The sliced version has 9 functions per slice, with up to 561 slices. To get the SBDMRG total energy errors to within 1 mH took from 4-10 days (depending on $R$), with bond dimensions $m\sim 300-1000$, running on a 2013 quad core Mac mini with 16Gb. For triple $\zeta$ the sliced and non-sliced energies are also very close, but with the sliced version slightly lower. All DMRG calculations were performed using the ITensor library ITe .

In Fig. 4, we present results for very long chains, demonstrating the linear scaling of SBDMRG. These calculations were at the stretched distance $R=3.6$, using a sliced STO-6G basis with one basis function per slice, and grid spacing $a=0.2$. The inset shows the calculation time per sweep on a single core of a 2013 3.5GHz Mac Pro, for a sweep keeping $m=100$ states. The calculation time not only grows very close to linearly in the number of atoms, it is also quite modest. The largest system, with 1000 atoms, had over 18,000 sliced basis functions, and an $m=100$ sweep took a little more than an hour. The number of states kept was slowly ramped up, with 30 smaller-$m$, faster sweeps occuring before three $m=100$ sweeps. Subsequent sweeps up to $m=400$ showed that at $m=100$, the energy per atom was in error by only 0.06 mH (DMRG error only, excluding the finite basis and finite $a$ errors). The main part of the figure shows the energy per site, in comparison with QCDMRG STO-6G. The energy results show the modest difference in completeness of STO-6G and sliced STO-6G, and also demonstrate that the sliced DMRG is converged to high accuracy.

The sliced basis set approach we have introduced here can be seen to be very well suited to DMRG calculations. Coupled with a compression method for the interactions, this approach gives linear scaling of computation time with the length of the system, allowing very long systems to be treated. This formulation brings DMRG for electronic structure closer to DMRG for models, and new approaches introduced for models (such as working directly with an infinite chain) can probably be adapted to SBDMRG with little difficulty. We also anticipate that extending SBDMRG to more complicated molecules will be reasonably straightforward.

We acknowledge support from the Simons Foundation through the Many-Electron Collaboration, and from the U.S. Department of Energy, Office of Science, Basic Energy Sciences under award #DE-SC008696.

Appendix A Interaction Integrals for Sliced Basis Sets

Recall that to construct a Hamiltonian in a sliced basis set, one must compute the integrals

[TABLE]

for the interaction terms and

[TABLE]

for the single-particle terms. (Recall the full expression for $t^{nn}_{ij}$ includes the grid kinetic energy $t^{nn}_{ij}=\tilde{t}^{nn}_{ij}-\frac{1}{2a^{2}}\Delta_{nn}$.)

To use a sliced basis, we need to evaluate integrals between basis function representing:

the overlap of two nonorthogonal function on a slice 2. 2.

kinetic energy matrix elements on a slice, Eq. (17) 3. 3.

single particle potential matrix elements from the Coulomb potential of the nuclei, Eq. (18) 4. 4.

the two particle terms $V^{nn^{\prime}}_{ikjl}$, Eq. (16)

The integrals for (1) and (2) for Gaussian functions have simple analytic formulas. The matrix elements (3) can be considered a limiting case of (4), where we consider a nucleus as an $S$-type Gaussian of vanishing width on one slice, and then the terms from the second coordinate define the one particle potential; thus we need only consider case (4).

A.1 Gaussian Fitting For $\ell>0$ Integrals

For $S$ functions, the $V^{nn^{\prime}}_{ijkl}$ have analytic formulas. However, for other types of orbitals, the formulas get both tedious to derive and very time consuming to evaluate. Instead, we implemented another approach: fit the function $1/r$ to a sum of Gaussians

[TABLE]

The widths of these Gaussians $a_{i}^{-1/2}$ were taken to be equally spaced on a logarithmic scale, except for the ten largest widths, which were optimized over both $a_{i}$ and $c_{i}$. Using $P=220$, we obtained a fit good to $O(10^{-10})$ over the range $10^{-8}$ to $10000$. The integrals were evaluate by taking the sum over $i$ outside the integrals, turning them into simple analytic Gaussian integrals which also separated by dimension $x,y,z$. The separation meant that the integral formulas for each single dimension could be calculated and stored quickly, and then each $V^{nn^{\prime}}_{ijkl}$ evaluation could be done as a loop of length $P$ involving only multiplications and additions, making it very fast.

A.2 Smoothing Procedure for Integrals with Cusps

If a continuous function does not have any frequency components above $\pi/a$, sampling it with grid spacing $a$ is exact. In contrast, sampling a function with a slope discontinuity leads to errors in the function of order $a$. The divergence of the $1/r$ interaction at short distances makes some of the two electron interaction integrals have slope discontinuities at $z=z^{\prime}$.

To accelerate the convergence with $a$, we adopt a pre-filtering technique, which is done before any contractions, when the integrals are still a function of $z-z^{\prime}$. The interaction is first computed at a finer grid spacing of $2^{-r}a$ for a small integer $r$. Then the interactions are put through a low-pass filter and factor-of-two decimation $r$ separate times, giving a final spacing of $a$. The low pass filter is designed to reproduce exactly all frequencies up to half the maximum frequency. Thus this smoothing procedure does not alter any low frequency parts of the interaction, but smoothly removes components at frequencies higher than $\pi/a$. The same smoothing procedure is also used for the nucleus-electron interaction integrals. We tested the accuracy of this procedure on H2, and found $r=3$ nicely accelerates convergence with $a$ while not increasing the computation time too much. The errors in the resulting total energies shown in Fig. (5) scale approximately as $a^{2.5}$ to $a^{3.1}$, and are approximately 0.1 mH per atom at $a=0.1$.

Appendix B SVD Compression of Long-Range Interactions

In this section we give a more detailed discussion of the compression algorithm for long-range interactions described in the main body of the paper. The simplest case is compressing the interaction part of the sliced basis set Hamiltonian for the case of one transverse function per slice.

[TABLE]

Later below we discuss how to generalize the compression for the case of multiple transverse functions.

The basic idea of the compression algorithm is to use the singular value decomposition (SVD) to compress each of the rectangular blocks $V^{(k)}$ of the matrix $V$ extending from the element $V_{kk}$ to the upper-right corner of $V$. As a motivation, consider the case where the interactions decay exponentially:

[TABLE]

Restricting $V$ to an upper-right block constrains $n^{\prime}\geq n$, in which case $V$ factorizes as

[TABLE]

This factorization into the outer product of two vectors implies that each upper-right block $V^{(k)}$ has only one non-zero singular value (is rank 1) and will be maximally compressed by an SVD. The interaction matrix for a real system will be more complicated, but if one can approximate it as a sum of exponentials, then the number of significant singular values of $V^{(k)}$ should remain small. In practice, the SVD can uncover better compression strategies than just a sum of real exponentials.

B.1 Algorithm for an $N\times N$ Matrix

First we will detail the compression algorithm for the case where $V_{nn^{\prime}}$ is just an $N\times N$ matrix, and later generalize to the case where $V$ is a tensor (the latter corresponding in SBDMRG to having multiple functions on each slice). The compression deals with the upper-right blocks $V^{(p)}$ of $V$, defined such that

[TABLE]

where $r=1,2,...,p$ and $c=1,2,\ldots,(N-p+1)$, see Fig. 6.

For each of these blocks we define the matrices $U^{(p)}$, $S^{(p)}$, and $W^{(p)}$ by an SVD of $V^{(p)}$:

[TABLE]

The matrix $S^{(p)}$ is diagonal and contains the singular values. Assuming the smoothness of $V_{nn^{\prime}}$ away from the diagonal makes the $V^{(p)}$ have only $D$ significant singular values, the compression is achieved by truncating $S^{(p)}$ to be only a $D\times D$ matrix, reducing the corresponding columns of $U^{(p)}$ and rows of $W^{(p)}$.

We next seek a way to relate the SVD of any one of the blocks $V^{(p)}$ to another block $V^{(p+1)}$. It is helpful to define the following additional notation:

Define $C^{-}(M)$ to be the matrix $M$ with the first column removed (making a smaller matrix). 2. 2.

Define $M\oplus r$ to be $M$ with an extra row $r$ added at the bottom ($r$ is a vector). 3. 3.

For an $n\times m$ matrix $M$, define $P(M)$ to be $M$ with an extra row and column added at the bottom and right. The extra matrix elements are zero, except for a 1 on the diagonal (at position $(n+1),(m+1)$). 4. 4.

Define $r^{(p)}$ to be the bottom row of $V^{(p)}$.

Then it follows that

[TABLE]

Writing the SVD of the matrix in square brackets in Eq. (27) as $X^{(p+1)}S^{(p+1)}W^{(p+1)}$, we find that we have obtained the SVD of $V^{(p+1)}$, with

[TABLE]

Each matrix $X^{(p)}$ is of dimension $(D+1)\times D$. We see that we can recover all the $U^{(p)}$ if we know all the $X^{(p)}$ plus $U^{(1)}$. A similar calculation gives all the $W^{(p)}$ in terms of a reverse recursion involving $D\times(D+1)$ matrices $Y^{(p)}$. This means we can reconstruct the entire $N\times N$ matrix $V(n-n^{\prime})$ out of the $O(ND^{2})$ parameters in $X^{(p)}$, $S^{(p)}$, and $Y^{(p)}$.

In practice, to obtain the fully compressed representation of $V$, it is useful to start by computing the SVD of $V^{(2)}$ (the SVD of $V^{(1)}$ is trivial). The initial SVD has a cost only linear in $N$ since $V^{(2)}$ is a $2\times N$ matrix. To compute the $X^{(p)}$, one computes SVDs of the matrices $[S^{(p)}C^{-}(W^{(p)})\oplus r^{(p+1)}]$ which are of dimension $(D+1)\times(N-p+1)$. Thus the cost for each of these SVDs scales as $D^{2}N$ (assuming the entries of the matrix $V$ have already been computed). For a non-translationally invariant system, one must perform $N$ such SVDs, making the total cost $D^{2}N^{2}$. But the compression algorithm only has to be performed once, and thus does not dominate the scaling of a SBDMRG calculation. To achieve a linear scaling of the compression algorithm, one could start with a translationally invariant basis such that the SVD Eq. (24) is the same for every block of $V$. Following the compression, the basis can be contracted to a smaller number of functions in a non-translationally-invariant manner on each slice.

B.2 MPO Form of Compressed Interactions

A matrix product operator (MPO) is a compact rewriting of a sum of operators as a tensor network. An MPO resembles a matrix product state (MPS), but in an MPO each tensor has two physical indices. Thus each MPO tensor can be viewed as an operator valued matrix, which will be the notation we use below. Representing the Hamiltonian as an MPO, or as a sum of MPOs, not only makes a code more generic and flexible, but can also make calculations more efficient.

Any sum of finite-range operators can be written exactly as an MPO using well-known conventions, which results in internal MPO indices whose sizes depend linearly on the range of the operators McCulloch (2007); Crosswhite and Bacon (2008). However, such an approach fails to be efficient when Hamiltonian terms do not have strictly finite support.

An interesting extension of the finite-range MPO construction allows MPOs to exactly capture sums of operators whose coefficients decay as pure exponentials McCulloch (2008); Crosswhite et al. (2008). By fitting other kinds of long-range terms, such as power-law decaying terms, to a sum of exponentials Pirvu et al. (2010), they can be approximated by MPOs in an efficient way.

But the exponential fitting approach leaves much to be desired. The best quality fits involve complex exponents, yet working with complex numbers incurs significant computational costs. Using a two-dimensional real matrix representation of the complex numbers avoids this issue, but complicates the method. Setting up the fits and the logic of the exponential decays for ladders and other quasi-one-dimensional systems with unit cells is also quite difficult.

Here we present an alternate approach to approximating sums of long-range operators as MPOs based on the SVD based compression algorithm discussed above. The approach here is closely related to the one proposed in Ref. Chan et al., 2016, especially in terms of the final MPO produced. But the present approach has some extra efficiencies arising from step that computes each $X^{(p)}$ from a matrix with only $(D+1)$ rows defined in square brackets Eq. (27). The cost of each associated SVD is at most linear in $N$, whereas the proposal in Ref. Chan et al., 2016 requires SVDs scaling as $N^{3}$. Both approaches are also related to a very general proposal for compressing MPOs in Ref. Zaletel et al., 2015.

In this section, we want to use the compression algorithm to produce an MPO for the sum of operators

[TABLE]

where $n,n^{\prime}=1,2,\ldots,N$. An MPO representation of $\hat{V}$ can be written as

[TABLE]

where each $M^{(n)}$ is an operator-valued matrix.

To make the following expressions more compact, it is convenient to define $\Omega^{(p)}=X^{(p)}S^{(p)}W^{(p)}$. Define the first MPO tensor to be:

[TABLE]

noting that $X^{(1)}\stackrel{{\scriptstyle\text{def}}}{{=}}U^{(1)}=1$. Define the second MPO tensor to be:

[TABLE]

And define the third MPO tensor to be:

[TABLE]

The general pattern for site $n$ is:

[TABLE]

From which we see the MPO has a matrix dimension of $(D+2)$.

Expanding this MPO, we can see that it represents $V$ as a sum of terms of the form

[TABLE]

where the notation $\tilde{M}$ means the first $D$ rows of a matrix $M$ (either $X$ or $\Omega$). To see how the expression in Eq. (35) recovers the matrix $V_{nn^{\prime}}$, note that row $(D+1)$ of each matrix $X^{(n)}$ is identical to row $n$ of $U^{(n)}$. Also note that

[TABLE]

for any $r\leq n$; the above equation can be seen to hold by omitting the last row of each of the matrices in Eq. (28). It follows that

[TABLE]

B.3 Generalization to Multiple Transverse Functions

For a sliced basis set with multiple transverse functions $\phi_{j}(x,y)$ on each slice, the interaction terms have the general form

[TABLE]

where $n,n^{\prime}=1,2,\ldots,N_{z}$ and $i,j,k,l=1,2,\ldots,N_{o}$. Thus to compress these interactions one must compress the tensor $V^{nn^{\prime}}_{ijkl}$. Because the indices $i,l$ label functions on slice $n$ and $j,k$ functions on slice $n^{\prime}$, reshape the tensor to an $(N_{z}N_{o}^{2})\times(N_{z}N_{o}^{2})$ matrix

[TABLE]

Then we can use a similar compression algorithm as that described above, with the key difference that one defines blocks of $V$ according to the $n,n^{\prime}$ indices, treating the $i,l$ or $j,k$ indices as a “unit cell” for each value of $n$ or $n^{\prime}$. So in contrast to the previous algorithm, where one would add a single row of $V$ in Eq. (25), for example, in the more general algorithm one adds $N_{o}^{2}$ rows of $V_{(nil)(n^{\prime}jk)}$.

The SVD one wants to obtain for each block $V^{(p)}$ of $V$ is of the form

[TABLE]

where $r=1,2,...,p$ and $c=1,2,\ldots,(N-p+1)$, and $i,l$ label the functions on slice $r$ while $j,k$ label the functions on slice $c$.

To compute matrices $X^{(p)}$ relating the SVD at one slice to that at another, make the following definitions:

Define $C^{-}_{N^{2}_{o}}(M)$ to be the matrix $M$ with the first $N_{o}^{2}$ columns removed. 2. 2.

Define $A\oplus B$ for an $a\times m$ matrix $A$ and an $b\times m$ matrix $B$ to be the $(a+b)\times m$ matrix whose first $a$ rows are those of $A$ and last $b$ rows are those of $B$. 3. 3.

Define $P_{N^{2}_{o}}$ for an $n\times m$ matrix $M$ to be the direct sum of $M$ and an $N_{o}^{2}\times N_{o}^{2}$ identity matrix. That is, append $N^{2}_{o}$ rows and columns to $M$ that are zero except for the diagonal elements which equal 1. 4. 4.

Define $v^{(p)}$ to be the last $N^{2}_{o}$ rows of $V^{(p)}$.

Then the block $V^{(p+1)}$ is given by

[TABLE]

By computing an SVD of the matrix in square brackets in Eq. (49) above, and writing this SVD as $X^{(p+1)}S^{(p+1)}W^{(p+1)}$, it follows that

[TABLE]

similar to the algorithm for the $N_{o}=1$ case in Section B.1. It is helpful to note that the last $N^{2}_{o}$ rows of each matrix $X^{(p)}$ correspond to the last $N^{2}_{o}$ rows of $U^{(p)}$, which correspond to the indices $i,l$ labeling functions on slice $p$.

Finally, for the next section on constructing an MPO, it will be convenient to define

[TABLE]

B.4 MPO For $N_{o}>1$ Orbitals on Each Slice

Consider the case $N_{o}=3$ and $D=2$. To lighten the notation, we consider a particular slice $p$, suppressing the label $(p)$ and assuming all MPO matrices and matrices $X_{\alpha_{p-1}\alpha_{p}}$,$U^{il}_{r\alpha_{p}}$,$\Omega^{jk}_{\alpha_{p-1}c}$ are all associated with the same slice $p$. Subscripts on MPO matrices and on operators indicate the orbital number within the slice $p$. The entire MPO is formed by repeating these $N_{o}$ matrices for all $N_{z}$ slices, together with the boundary conditions given later below.

The MPO matrix for the first orbital on a slice is:

[TABLE]

where the “$\cdots$” indicate that the last eight columns are repeated, replacing $c^{\dagger}_{1\uparrow}\rightarrow c^{\dagger}_{1\downarrow}$ and $c_{1\uparrow}\rightarrow c_{1\downarrow}$.

The second MPO matrix is:

[TABLE]

where $F_{2}=(-1)^{n_{2}}$ is a fermion string operator. We include this detail to note that, at least in our ITensor implementation, the operators we denote here as $c$ and $c^{\dagger}$ only anticommute when acting on the same site, so the additional $F$ operators must be included between $c^{\dagger}$ and $c$ pairs acting on different sites. If the anticommutation bookkeeping is done in a more automatic way, where $c$ and $c^{\dagger}$ really do anticommute across different sites, one would replace these $F$ operators with identity operators.

The third, and last MPO matrix on this $N_{o}=3$ slice is:

[TABLE]

To make a well-defined MPO for a finite system, the first and last MPO tensors are contracted with the a boundary vector to the left of the first site:

[TABLE]

and a boundary vector to the right of the last site:

[TABLE]

To explain the design of the MPO above in words (recalling that it is a concrete example for the case $N_{o}=3$ and $D=2$, so that the row and column numbers are specific to that case):

Row 1 of each MPO matrix holds operators which begin an operator “string” on that site (these are the “starting” operators in a finite-state automaton picture of an MPO, Ref. Crosswhite et al., 2008). 2. 2.

The identity operator at element (2,2) of each matrix trails a completed string of operators (the “done” state in an automaton picture). 3. 3.

Rows and columns 3 and 4 correspond to the $\alpha$ indices formed from the SVDs in the compression algorithm (for general $D$ this would be rows and columns $3,4,\ldots,D+2$). For sites 1 and 2, operators from previous slices either connect with elements of $\Omega$ to form a completed operator string or are passed through to the next site. On site 3 (more generally site $N_{o}$), the $X$ matrix appears, transforming incomplete operator strings from the $\alpha_{p-1}$ basis into the $\alpha_{p}$ basis. 4. 4.

Columns 5–8 of $\hat{M}_{1}$, and rows and columns 5–8 of $\hat{M}_{2}$ collect pairs of operators on sites 1 and 2. In rows 5–8 of $\hat{M}_{3}$, these operator pairs get multiplied by elements of $U$ on site 3 to begin a new operator string connecting to a different slice. 5. 5.

Rows and columns 9 and 10 of $\hat{M}_{1}$ and $\hat{M}_{2}$ multiply nearly-complete operator strings from a previous slice by an element of $\Omega$ and a $c^{\dagger}_{\uparrow}$ operator. However these operator strings must be carried on to one of the remaining sites in the slice (sites 2 or 3) to be matched with the $c_{\uparrow}$ operators in column 2 of $\hat{M}_{2}$ and $\hat{M}_{3}$. 6. 6.

Rows and columns 11 and 12 begin operator strings starting with a $c^{\dagger}_{\uparrow}$ on either site 1 or 2, which will be paired with an element of $U$ and a $c_{3\uparrow}$ operator on site 3 to begin a new operator string.

The pattern of columns 9–12 of $\hat{M}_{1}$ and $\hat{M}_{2}$ repeats three more times, replacing $c^{\dagger}_{\uparrow}$ with $c_{\uparrow}$, $c^{\dagger}_{\downarrow}$, and $c_{\downarrow}$.

Note that in the third MPO matrix (more generally, the matrix on site number $N_{o}$ within a slice) we weighted new operator strings with elements of $U^{il}$ instead of elements of $X$ as in the $N_{o}=1$ MPO Eq. (34). This was for convenience as the elements $U^{(p)\,il}_{p\alpha_{p}}$ for $i,l=1,2,\ldots,N_{o}$ and $\alpha_{p}=1,2,\ldots,D$ correspond to the last $N_{o}^{2}$ rows of $X^{(p)}$, and listing these rows of $X^{(p)}$ would be unwieldy in the current notation.