Spin-projected matrix product states (SP-MPS): a versatile tool for strongly correlated systems

TL;DR

The paper introduces spin-projected matrix product states (SP-MPS), a flexible wavefunction ansatz that simplifies achieving spin-adaptation and enables efficient exploration of complex magnetic systems in strongly correlated electronic structures.

Contribution

It develops a new SP-MPS ansatz combining spin projection with MPS, offering a simpler alternative to non-Abelian symmetry incorporation and facilitating studies of magnetic and open-shell systems.

Findings

Provides a simpler route to spin-adapted wavefunctions.

Enables straightforward generation of broken symmetry states.

Supports parallel computation of expectation values.

Abstract

We present a new wavefunction ansatz that combines the strengths of spin projection with the language of matrix product states (MPS) and matrix product operators (MPO) as used in the density matrix renormalization group (DMRG). Specifically, spin-projected matrix product states (SP-MPS) are constructed as $|\Psi^{(N,S,M)}_{SP-MPS}\rangle=\mathcal{P}_S|\Psi_{MPS}^{(N,M)}\rangle$, where $\mathcal{P}_S$ is the spin projector for total spin $S$ and $|\Psi_{MPS}^{(N,M)}\rangle$ is an MPS wavefunction with a given particle number $N$ and spin projection $M$. This new ansatz possesses several attractive features: (1) It provides a much simpler route to achieve spin-adaptation (i.e. to create eigenfunctions of $\hat{S}^2$) compared to explicitly incorporating the non-Abelian SU(2) symmetry into the MPS. In particular, since the underlying state in the SP-MPS uses only Abelian symmetries, one does not need the singlet embedding scheme for non-singlet states, as normally employed in spin-adapted DMRG, to achieve a single consistent variationally optimized state. (2) Due to the use of $|\Psi_{MPS}^{(N,M)}\rangle$ as its underlying state, the SP-MPS can be closely connected to broken symmetry mean-field states. This allows to straightforwardly generate the large number of broken symmetry guesses needed to explore complex electronic landscapes in magnetic systems. Further, this connection can be exploited in the future development of quantum embedding theories for open-shell systems. (3) The sum of MPOs representation for the Hamiltonian and spin projector $\mathcal{P}_S$ naturally leads to an embarrassingly parallel algorithm for computing expectation values and optimizing SP-MPS. (4) Optimizing SP-MPS belongs to the variation-after-projection (VAP) class of spin projected theories.

Equations

$$|\Psi ⟩=\sum _{\{n_k\}}{\Psi }^{n_1{n}_2\cdots n_K}|n_1{n}_2\cdots n_K⟩,$$

$${\Psi }^{n_1{n}_2\cdots n_K}=\sum _{\{{\alpha }_k\}}{A}_{{\alpha }_1}^{n_1}[1]A_{{\alpha }_1{\alpha }_2}^{n_2}[2]\cdots A_{{\alpha }_{K-1}}^{n_K}[K],$$

$$\hat{O}=\sum _{\{n_k\}\{n_k^{\prime }\}}{O}_{n_1^{\prime }{n}_2^{\prime }\cdots n_K^{\prime }}^{n_1{n}_2\cdots n_K}|n_1{n}_2\cdots n_K⟩⟨n_1^{\prime }{n}_2^{\prime }\cdots n_K^{\prime }|,O_{n_1^{\prime }{n}_2^{\prime }\cdots n_K^{\prime }}^{n_1{n}_2\cdots n_K}\triangleq ⟨n_1{n}_2\cdots n_K|O|n_1^{\prime }{n}_2^{\prime }\cdots n_K^{\prime }⟩,$$

$$O_{n_1^{\prime }{n}_2^{\prime }\cdots n_K^{\prime }}^{n_1{n}_2\cdots n_K}=\sum _{\{{\beta }_k\}}{W}_{{\beta }_1}^{n_1{n}_1^{\prime }}[1]W_{{\beta }_1{\beta }_2}^{n_2{n}_2^{\prime }}[2]\cdots W_{{\beta }_{K-1}}^{n_K{n}_K^{\prime }}[K].$$

$$\hat{H}=\sum _{pq}{h}_{pq}{a}_p^{†}{a}_q+\frac{1}{2}\sum _{pqrs}{g}_{pq,sr}{a}_p^{†}{a}_q^{†}{a}_r{a}_s,g_{pq,sr}=⟨pq|sr⟩.$$

$$\hat{H}=\sum _p{\hat{H}}_p,{\hat{H}}_p=a_p^{†}\left(\sum _q{h}_{pq}{a}_q+\sum _{q,r<s}{g}_{pqrs}^A{a}_q^{†}{a}_r{a}_s\right),g_{pqrs}^A\triangleq \frac{1}{2}(g_{pq,sr}-g_{pq,rs}).$$

$$\hat{V}=\sum _q{h}_q{a}_q+\sum _{qrs}{V}_{qrs}{a}_q^{†}{a}_r{a}_s,V_{q,r\ge s}=0.$$

$${\hat{V}}^{[k,K]}$$

$${\hat{P}}_l^{[k,K]}$$

$${\hat{Q}}_l^{[k,K]}$$

$$\displaystyle\left[\begin{array}[]{c}\hat{V}^{[k,K]}\\ \hline\cr\hline\cr(a^{\dagger}_{[k]})^{[k,K]}\\ (a^{\dagger}_{[k+1,K]})^{[k,K]}\\ \hline\cr(a_{[k]})^{[k,K]}\\ (a_{[k+1,K]})^{[k,K]}\\ \hline\cr\hat{Q}^{[k,K]}_{[1,k-1]}\\ \hline\cr\hat{P}^{[k,K]}_{[1,k-1]}\\ \hline\cr\hline\cr I^{[k,K]}\\ \end{array}\right]=\left[\begin{array}[]{c||c|c|c|c|c|c||c}\widetilde{I}_{k}&0&V_{kkr}\widetilde{n}_{k}&0&-a_{k}&0&a_{k}^{\dagger}&h_{k}a_{k}\\ \hline\cr\hline\cr 0&0&0&0&0&0&0&a^{\dagger}_{k}\\ 0&\widetilde{I}_{k}&0&0&0&0&0&0\\ \hline\cr 0&0&0&0&0&0&0&a_{k}\\ 0&0&\widetilde{I}_{k}&0&0&0&0&0\\ \hline\cr 0&-V_{rlk}\widetilde{a}_{k}&V_{klr}\widetilde{a}^{\dagger}_{k}&I_{k}&0&0&0&V_{klk}n_{k}\\ \hline\cr 0&0&V_{lkr}\widetilde{a}_{k}&0&0&I_{k}&0&0\\ \hline\cr\hline\cr 0&0&0&0&0&0&0&I_{k}\\ \end{array}\right]\left[\begin{array}[]{c}\hat{V}^{[k+1,K]}\\ \hline\cr\hline\cr(a^{\dagger}_{[k+1,K]})^{[k+1,K]}\\ \hline\cr(a_{[k+1,K]})^{[k+1,K]}\\ \hline\cr\hat{Q}^{[k+1,K]}_{[1,k-1]}\\ \hat{Q}^{[k+1,K]}_{[k]}\\ \hline\cr\hat{P}^{[k+1,K]}_{[1,k-1]}\\ \hat{P}^{[k+1,K]}_{[k]}\\ \hline\cr\hline\cr I^{[k+1,K]}\\ \end{array}\right],$$

$$[I]$$

$$[a^{†}]$$

$$[a]$$

$$[n]$$

$$|{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}⟩={\mathcal{P}}_S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩,$$

$${\mathcal{P}}_S=\prod _{S^{\prime }\ne S}{\hat{P}}_{S^{\prime }},{\hat{P}}_{S^{\prime }}=\frac{{\hat{S}}^2-S^{\prime }(S^{\prime }+1)}{S(S+1)-S^{\prime }(S^{\prime }+1)},$$

$$\displaystyle\left[\begin{array}[]{c}(\hat{S}^{2})^{[k,K]}\\ \hline\cr\hline\cr(\hat{S}_{+})^{[k,K]}\\ \hline\cr(\hat{S}_{-})^{[k,K]}\\ \hline\cr(\hat{S}_{z})^{[k,K]}\\ \hline\cr\hline\cr I^{[k,K]}\\ \end{array}\right]=\left[\begin{array}[]{c||c|c|c||c}I_{k}&\hat{S}_{-,k}&\hat{S}_{+,k}&2\hat{S}_{z,k}&\hat{S}^{2}_{k}\\ \hline\cr\hline\cr 0&I_{k}&0&0&\hat{S}_{+,k}\\ \hline\cr 0&0&I_{k}&0&\hat{S}_{-,k}\\ \hline\cr 0&0&0&I_{k}&\hat{S}_{z,k}\\ \hline\cr\hline\cr 0&0&0&0&I_{k}\\ \end{array}\right]\left[\begin{array}[]{c}(\hat{S}^{2})^{[k+1,K]}\\ \hline\cr\hline\cr(\hat{S}_{+})^{[k+1,K]}\\ \hline\cr(\hat{S}_{-})^{[k+1,K]}\\ \hline\cr(\hat{S}_{z})^{[k+1,K]}\\ \hline\cr\hline\cr I^{[k+1,K]}\\ \end{array}\right],$$

$$[I]$$

$$[{\hat{S}}_z]$$

$${\mathcal{P}}_{M,M^{\prime }}^S=\frac{2S+1}{8{\pi }^2}\int d\Omega D_{M,M^{\prime }}^{S\ast }(\Omega )\hat{R}(\Omega ),\hat{R}(\Omega )=e^{-\text{mathbbm}i\alpha {\hat{S}}_z}{e}^{-\text{mathbbm}i\beta {\hat{S}}_y}{e}^{-\text{mathbbm}i\gamma {\hat{S}}_z},$$

$${\mathcal{P}}_{M,M^{\prime }}^S$$

$${\mathcal{P}}_M$$

$${\hat{P}}_{M,M^{\prime }}^S$$

$$|{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}⟩={\mathcal{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩={\mathcal{P}}_M{\hat{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩.$$

$$E[|{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}⟩]$$

$${\hat{P}}_{M,M}^S=\sum _{g=1}^{N_g}{w}_g{d}_{M,M}^S({\beta }_g)e^{-\text{mathbbm}i{\beta }_g{\hat{S}}_y},$$

$$\displaystyle[e^{-\mathbbm{i}\beta\hat{S}_{y,k}}]=\left[\begin{array}[]{cccc}1&0&0&0\\ 0&c&s&0\\ 0&-s&c&0\\ 0&0&0&1\\ \end{array}\right],\quad c=\cos(\beta/2),\quad s=\sin(\beta/2).$$

$$N_g$$

$${\Omega }_{max}$$

$$\frac{\partial ⟨{\Psi }_{\mathrm{MPS}}^{(N,M)}|\hat{H}{\hat{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩}{\partial A^{\ast }[k]}=E\frac{\partial ⟨{\Psi }_{\mathrm{MPS}}^{(N,M)}|{\hat{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩}{\partial A^{\ast }[k]},$$

$$H_{\mathrm{eff}}A[k]=N_{\mathrm{eff}}A[k]E.$$

$$0$$

$$B[k]$$

$$⟨{\hat{O}}_{sf}⟩=\frac{⟨{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}|{\hat{O}}_{sf}|{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}⟩}{⟨{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}|{\Psi }_{\mathrm{SP}\text{-}\mathrm{MPS}}^{(N,S,M)}⟩}=\frac{⟨{\Psi }_{\mathrm{MPS}}^{(N,M)}|{\hat{O}}_{sf}{\hat{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩}{⟨{\Psi }_{\mathrm{MPS}}^{(N,M)}|{\hat{P}}_{M,M}^S|{\Psi }_{\mathrm{MPS}}^{(N,M)}⟩}.$$

Full text

To illustrate the versatility of SP-MPS, we formulate algorithms for the optimization of ground and excited states, develop perturbation theory based on SP-MPS, and describe how to evaluate spin-independent and spin-dependent properties such as the reduced density matrices. We demonstrate the numerical performance of SP-MPS with applications to several models typical of strong correlation, including the Hubbard model, and [2Fe-2S] and [4Fe-4S] model complexes.

1 Introduction

Since its invention by White1, 2, the density matrix renormalization group (DMRG) has become the computational method of choice in strongly correlated one-dimensional systems. In quantum chemistry, despite the more complex structure of the entanglement between the orbitals in a general molecule, the DMRG has been applied successfully to a large class of chemical systems with strong multi-configurational effects and many open shells3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36. Part of its popularity in this context stems from its systematic improvability by changing a single parameter (the bond dimension), its lower computational scaling than other tensor network state (TNS) methods as a function of bond dimension, as well as various techniques that have significantly enhanced its efficiency in chemical applications, such as the use of orbital localization, reordering, point group and spin symmetry, and parallelization10. The identification of matrix product states (MPS) as the underlying variational wavefunction ansatz in DMRG37, 38 has further expanded the applications of DMRG and related techniques in quantum chemistry. For example, the similarities between MPS (DMRG) and Slater determinants (self-consistent field), has led to the identification of the Thouless parameterization for MPS 39 and the formulation of linear response theory and tangent space methods for excited states and for time evolution40, 41, 42, 43. Another example is the use of MPS to develop efficient multi-reference perturbation theories44, which avoid the computation of high-order ($>2$) reduced density matrices (RDMs) by directly minimizing the Hylleraas functional with the first-order wavefunction represented as a MPS. Finally, MPS have recently been reformulated in Hilbert space, providing a deep connection to configuration interaction approximations and their graphical representation 45. We emphasize that the MPS and DMRG languages are entirely complementary. The MPS language provides a precise mathematical description of the wavefunctions used in DMRG algorithms as a product of site tensors $A[k]$. The DMRG language, on the other hand, is the appropriate way to describe many of the efficient algorithms to evaluate and manipulate MPS.

In this paper, using MPS and related matrix product operator (MPO) techniques 46, 47, 48, 49, 50, we revisit the problem of constructing spin eigenfunctions from MPS, i.e. spin adaptation. Previously, spin adapted MPS were formulated by imposing non-Abelian SU(2) symmetry structure on the site tensors51, 52, 53, 54, 24, 31, 55. This amounts to generating spin-adapted reduced bases from direct products of two spin-adapted reduced bases during the DMRG blocking procedure. However, although it provides a practical advantage in terms of reducing the computational bond dimension by approximately half24, there are some formal and practical drawbacks to the explicitly spin-adapted MPS (SA-MPS) approach. First, for non-singlet states, the one-site DMRG optimization leads to different MPS wavefunctions, even at the same site of the sweep, when moving in the left and right directions. A simple example illustrates this24: consider a reduced triplet wavefunction written as $\|\Psi_{S=1}\rangle=\frac{1}{\sqrt{2}}\|S_{L}=1\rangle\times(\|S_{R}=0\rangle+\|S^{\prime}_{R}=2\rangle)$, where $\|S\rangle$ represents a reduced spin state with total spin quantum number $S$. In this case, the left reduced density matrix is $\rho^{L}=\left[1\right]$, while the right reduced density matrix is \rho^{R}=\frac{1}{2}\left[\begin{array}[]{cc}1&1\\ 1&1\\ \end{array}\right]. The density matrix $\rho^{L}$ is of rank one, which yields a discarded weight of zero in the DMRG optimization if the bond dimension $D$ (the number of states to be kept) is $\geq 1$. However, although the density matrix $\rho^{R}$ is also of rank one, in order to ensure that the renormalized states have well-defined spin quantum numbers, only the block diagonal part of $\rho^{R}$ (the so-called quasi-density matrix52, 53, 54, 56) enters the spin-adapted renormalization procedure, generating two renormalized states. This leads to a truncation of the wavefunction if $D=1$ and thus, even at the same site, the DMRG wavefunctions generated by a sweep to the left or to the right are different and yield different energies. While this is not a severe problem when $D$ is large as the difference is very small, the DMRG sweep no longer strictly corresponds to an energy minimization. This complicates the computation of properties, such as the nuclear gradients, since $\partial E/\partial A[k]$ may not be exactly zero. To avoid the above problem, one usually uses the singlet embedding scheme57, 24, 31 to create a singlet total wavefunction by coupling the non-singlet physical state to a set of noninteracting fictitious spins. Within this singlet total state, the spin quantum numbers and dimensions of the left and right renormalized Hilbert spaces must always match, ensuring that the one-site DMRG optimization algorithm converges to a single consistent MPS. However, the singlet embedding scheme introduces its own complications, for example in the evaluation of transition density matrices between two states with different spins, as required in the state interaction treatment of spin-orbit coupling with SA-MPS58.

The second, and perhaps more important in practice, deficiency of the spin-adapted MPS formulation is that it does not provide a simple connection to the intuitive spin structures and charge configurations of broken symmetry Slater determinants. This connection is important to retain for two reasons. First, it helps in interpreting the wavefunction when understanding the electronic structure of systems with many open shells, as is found in transition metal clusters59. Second, the connection to broken symmetry determinants can be used to easily prepare many different types of initial guess for the DMRG optimization. In complex systems where there are many competing low-energy spin states, the ability to systematically prepare many such initial guesses ensures that the subsequent optimization avoids physically irrelevant local minima.

To address these issues, in this work we propose an alternative way to construct spin eigenfunctions from MPS using spin projection. Spin projection techniques have a long history in quantum chemistry, dating back to the work by Löwdin60 using a spin projector with broken symmetry Slater determinants. The general idea has more recently been revived61, 62, 63, 64, 65, 66 using group theoretical projectors rather than Löwdin’s original projector. Here, we will use the wavefunction ansatz $|\Psi^{(N,S,M)}_{\mathrm{SP\text{-}MPS}}\rangle=\mathcal{P}_{S}|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ to obtain spin eigenstates, where $\mathcal{P}_{S}$ is the spin projector for total spin $S$ and $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ is an MPS wavefunction with given particle number $N$ and spin projection $M$. Such spin-projected matrix product states (SP-MPS) display the following interesting features:

• •

They allow for a consistent energy minimization procedure for non-singlet states without using the singlet embedding scheme as only Abelian symmetries are used in the underlying state $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$. As mentioned above, this is mainly a formal advantage over SA-MPS.

• •

The underlying MPS can reduce at $D=1$ to a “broken symmetry” determinant. Although in principle, different spatial orbitals for different spins (DODS) can be used for the underlying MPS, throughout this paper the same spatial orbitals for different spins (SODS) will be used, to simplify the representation of the spin projector $\mathcal{P}_{S}$ (see Sec. 2.2). Thus, the term “broken symmetry” used here refers to the breaking of spin symmetry at the level of configurations instead of orbitals. Compared with SA-MPS, this connection to a “broken symmetry” determinant helps in the interpretation of complex electronic wavefunctions, as well as in generating initial guesses for $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$, to systematically explore electronic landscapes with many competing states. Further, this connection may be used to extend the density matrix embedding theory (DMET)67, 68 to overall open-shell systems, using a mean-field state and an SP-MPS state for the low- and high-level descriptions of the open-shell system, respectively.

• •

It is a variation-after-projection (VAP) ansatz, where the underlying broken-symmetry MPS can be efficiently optimized with a DMRG sweep algorithm. Unlike other spin projected theories in quantum chemistry, it is easy to make the variational state more accurate, as one can simply increase the bond dimension in $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$. The resulting spin-projected MPS retains all the advantages of conventional MPS. For example, computing excited states is straightforward, as imposing orthogonality constraints is simple with MPS.

• •

A final advantage is that the computational implementation of SP-MPS algorithms is in many respects simpler than the implementation of algorithms involving spin-adapted MPS. This simplicity is potentially even more advantageous when working with more complex tensor network states, such as projected entangled pair states (PEPS)69.

The remainder of the paper is organized as follows. In Sec. 2.1, we recapitulate the MPS and MPO representations for many-body states and operators, which are essential for describing how to work with the SP-MPS ansatz. In particular, we provide an explicit recursive algorithm to derive an MPO representation for the second quantized ab initio Hamiltonian. The wavefunction ansatz for SP-MPS is then introduced and its properties discussed in Sec. 2.2. The DMRG-like sweep algorithm for optimizing the SP-MPS representation of ground and excited states, as well as a parallelization scheme based on the sum of MPOs representation for the Hamiltonian and spin projector $\mathcal{P}_{S}$, are also formulated in Sec. 2.3. The evaluation of properties such as spin-free and spin-dependent reduced density matrices is considered in Sec. 2.4. Pilot applications to some prototypical open shell strongly correlated systems are described in Sec. 3. Conclusions and outlines for future directions are presented in Sec. 4.

2 Theory

2.1 Matrix product states (MPS) and matrix product operators (MPO)

In this section, we recapitulate essential aspects of the MPS and MPO representations for many-body states and operators. This will be necessary to formulate the theory of SP-MPS. A more detailed description of this language can be found in Refs. 70, 50, 71, 55. As an important example of how to work with the MPO representation, we will give an explicit algorithm to derive the MPO representation for the second quantized ab initio Hamiltonian.

A simple way to introduce matrix product states is from the perspective of a repeated singular value decompositions (SVDs) of a general wavefunction defined in Fock space,

[TABLE]

where $\ket{n_{1}n_{2}\cdots n_{K}}$ ($n_{k}\in\{0,1\}$) is an occupancy basis state with $K$ spin-orbitals. A successive SVD for the coefficient tensor $\Psi^{n_{1}n_{2}\cdots n_{K}}$ leads to an MPS representation with open boundary conditions (OBC)70,

[TABLE]

where the site tensor $A[k]$ away from the two boundaries is a three-way tensor of dimension $2\times D\times D$, and the two tensors at the boundary are matrices of dimension $2\times D$ or $D\times 2$. The so-called bond dimension $D$ controls the number of retained (renormalized) states, and therefore controls the accuracy of the representation when used as a variational ansatz. Similarly to Eq. (1), within the occupation number representation, with basis vectors $\ket{n_{1}n_{2}\cdots n_{K}}$, an operator $\hat{O}$ is represented by its coefficient tensors $O^{n_{1}n_{2}\cdots n_{K}}_{n_{1}^{\prime}n_{2}^{\prime}\cdots n_{K}^{\prime}}$, viz.,

[TABLE]

and similarly to Eq. (2), an MPO representation for the coefficient can be obtained through successive SVDs as

[TABLE]

To simplify the discussion of MPS and MPO, it is convenient to use a graphical representation. In Figure 1, the graphical representations for a MPS (2) and MPO (4) are shown. A key simplification comes from the convention that the sum over auxiliary (dummy) indices is replaced by a bond between two tensors graphically represented by dots. This notation eliminates the need to write out nested algebraic summations to express contractions between MPS and MPO.

An important point that deserves emphasis is that in both the MPS and MPO representations, the site tensors such as $A[k]$ in Eq. (2) and $W[k]$ in Eq. (4) are not uniquely determined. Given a set of occupancies, Eq. (2) tells us that the coefficient is expressed as a product of matrices with vectors at the boundaries. There is thus a gauge degree of freedom arising from the trivial fact that the product of two matrices $AB$ is not affected by inserting a pair of matrices $GG^{-1}$, i.e. $AB=(AG)(G^{-1}B)$, if $G$ is nonsingular. In practice, we use certain types of gauges to simplify computations and improve numerically stability.

As an important example of the MPO formalism in practice, we now derive an explicit MPO representation for the generic second quantized Hamiltonian written as

[TABLE]

The MPO representation we obtain does not depend on the symmetry of the integrals $h_{pq}$ and $g_{pq,rs}$, and therefore also applies to other two-body operators, such as the doubles excitation operator $\hat{T}_{2}$ used in correlation methods. To begin, following Ref. 50, we rewrite $\hat{H}$ in Eq. (5) as a sum of $K$ operators,

[TABLE]

such that it suffices to examine the construction of the MPO for an operator $\hat{V}$ of the following form,

[TABLE]

Once the MPO form for $\hat{V}$ is obtained, we can substitute it in the bracket in Eq. (6), giving $\hat{H}$ as a sum of $K$ MPOs. As will be shown below, the bond dimension for representing $\hat{V}$ exactly as an MPO is of $O(K)$, and since $a_{p}^{\dagger}$ is a simple MPO with bond dimension 1, each $\hat{H}_{p}$ will also be an MPO with bond dimension of $O(K)$. Consequently, adding together the $\hat{H}_{p}$, we will have an exact MPO representation for $\hat{H}$ with bond dimension $O(K^{2})$, which is the correct scaling of the bond dimension for generic two-body operators50.

Assuming one spin orbital per site, there are only four basis operators $\{I_{k},a^{\dagger}_{k},a_{k},n_{k}=a_{k}^{\dagger}a_{k}\}$ per site. We will use the notation $\hat{V}^{[k,K]}$ to denote the operator $\hat{V}$ defined with spin orbitals only in the range from $k$ to $K$. To derive the formula for the MPO representations of $\hat{V}=\hat{V}^{[1,k]}$ (7), we start by examining the recurrence relation for $\hat{V}^{[k,K]}$,

[TABLE]

where the symbol $\otimes$ represents the Kronecker product, the tilde in $\widetilde{I}_{k}$ is used to indicate the presence of the fermionic sign factor that needs to be taken into account when going to the matrix representation of the operator (vide post), and the intermediates (complementary operators) $\hat{P}$ and $\hat{Q}$ are defined by

[TABLE]

Putting Eqs. (8), (9), and (10) together, we recast these recurrence relations in a compact matrix-vector product form,

[TABLE]

where the subscript in operators such as $\hat{Q}^{[k+1,K]}_{[1,k-1]}$ indicates the range for the index $l$ in Eq. (10). Since $\hat{V}^{[1,K]}$ is the operator to be written as an MPO, one can immediately recognize that the $2(K+1)$-by-$2(K+1)$ coefficient matrix in Eq. (35) with operator entries is just the operator counterpart of the site tensor $W[k]$ (2), while at the left and right boundaries, the first row and the last column in Eq. (35) can be read off for $W[1]$ and $W[K]$, respectively. To obtain the site tensor $W[k]$, we use the matrix representations for operators such as $a_{k}$ in the space of $\{|0\rangle,|1\rangle\}$. The necessary matrix representations for all the operators in Eq. (35) are as follows:

[TABLE]

The Pauli matrices in Eq. (55) remind us that the above matrix representations are simply an expression of the Jordan-Wigner transformation72 that maps fermions to spins and vice versa. A final remark is that to obtain the MPO in the spatial orbital basis, in which the local dimension is 4 instead of 2, one need only to merge two adjacent MPO site tensors defined in Eq. (35), placing the $\alpha$ and $\beta$ spin orbitals together. The resulting spatial MPO is the one we will use in the energy minimization algorithm with SP-MPS. In the following, we will assume the use of spatial orbitals as sites, the index $k$ for spatial orbitals, and $K$ for the total number of spatial orbitals.

2.2 Spin-projected matrix product states (SP-MPS)

Based on the above MPS and MPO formalism, we can now introduce the SP-MPS wavefunction $|\Psi^{(N,S,M)}_{\mathrm{SP\text{-}MPS}}\rangle$ as resulting from acting a spin projector $\mathcal{P}_{S}$ on an MPS $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ with given particle number $N$ and spin projection $M$, viz.,

[TABLE]

where its graphical representation is shown in Figure 2(a). It is possible to develop an even more general family of MPS via additional symmetry projections, such as particle number projection, which may allow to use smaller bond dimensions in the underlying MPS to achieve the same accuracy. However, since our purpose here is mainly to circumvent the problems arising from SA-MPS as mentioned in Sec. 1, we will only consider the simplest case, where the underlying MPS uses the physical Abelian symmetries of particle number and spin projection. We will also assume the SODS scheme which helps to significantly reduce the bond dimension to represent the projector $\mathcal{P}_{S}$. Physically, this means that the task of describing fluctuations around a classical broken symmetry determinant is largely performed by the underlying MPS, similarly to as in the normal DMRG case. This is conceptually rather different from other spin projected methods61, 62, 63, 64, 65, 66, where the spin projection itself is essential for restoring fluctuations and correlation via a deliberate symmetry breaking and restoration mechanism.

There are various choices for the spin projector $\mathcal{P}_{S}$ in Eq. (56). Löwdin’s spin projector60 takes the form,

[TABLE]

where for given numbers of $\alpha$ and $\beta$ electrons, assuming $N_{\alpha}\geq N_{\beta}$, the allowed values for $S^{\prime}$ range from $\frac{N_{\alpha}-N_{\beta}}{2}$ to $\frac{N_{\alpha}+N_{\beta}}{2}$, such that the projector in Eq. (57) is a product of $N_{\beta}$ terms. The use of this operator is a formidable task in conjunction with standard quantum chemistry methods due to its complicated operator form. However, within the MPS and MPO formalism, as long as the operators involved are local, the computation is tractable regardless of the rank of the operators. This is the case for the operator $\hat{S}^{2}$ in the projector (57). Using $\hat{S}^{2}=\frac{1}{2}(\hat{S}_{+}\hat{S}_{-}+\hat{S}_{-}\hat{S}_{+})+\hat{S}_{z}^{2}$ and following the same recursive method for $\hat{H}$, one can find a compact MPO representation for $\hat{S}^{2}$ with a bond dimension of only 5, viz.,

[TABLE]

where both $k$ and $K$ represent indices for spatial orbitals. It must be emphasized that this simplicity is associated with the SODS scheme, while using the DODS scheme would lead to a more complicated representation for $\hat{S}^{2}$. The matrix representations for the local operators appearing in Eq. (73) in the space $\mathrm{span}\{|0\rangle,|k_{\alpha}\rangle\}\otimes\mathrm{span}\{|0\rangle,|k_{\beta}\rangle\}=\mathrm{span}\{|0\rangle,|k_{\beta}\rangle,|k_{\alpha}\rangle,|k_{\alpha}k_{\beta}\rangle\}$ are the same for all $k$, and are

[TABLE]

Since adding a constant and multiplying by a factor to go from $\hat{S}^{2}$ to $\hat{P}_{S^{\prime}}$ in Eq. (57) does not change the MPO bond dimension in Eq. (73), we can conclude that the bond dimension for the Löwdin projector is at most $5^{N_{\beta}}$, with its graphical representation shown in Figure 2(b). Although formally the bond dimension is exponential in $N_{\beta}$, in practice one can expect that the wavefunction $\mathcal{P}_{S}|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ will be highly compressible at least for ground states, in the sense that when one projector $\hat{P}_{S^{\prime}}$ is applied to $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ with bond dimension $D$, the resulting wavefunction should be compressible back to an MPS without increasing the bond dimension too much in order to achieve a good accuracy, such as ca. 1mH in the ground state energies. Thus, unlike with other quantum chemistry methods, the combination of the Löwdin projector with MPS is in principle possible. However, this ansatz does not naturally fit into the DMRG sweep optimization algorithm, and hence needs to be optimized by other techniques.

In this paper, we consider another form of projector that is more compatible with DMRG sweep optimization, viz., the group theoretical projector73,

[TABLE]

where $\Omega=(\alpha,\beta,\gamma)$ are the Euler angles, $\hat{R}(\Omega)$ is the rotation operator, $D_{M,M^{\prime}}^{S}(\Omega)=\langle SM|\hat{R}(\Omega)|SM^{\prime}\rangle=e^{-\mathbbm{i}M\alpha}d_{M,M^{\prime}}^{S}(\beta)e^{-\mathbbm{i}M^{\prime}\gamma}$ is the Wigner $D$-matrix, and $d_{M,M^{\prime}}^{S}(\beta)$ is an element of Wigner’s small $d$-matrix, which can be chosen to be real for simplicity. The projector $\mathcal{P}^{S}_{M,M^{\prime}}$ (96) can then be rewritten as

[TABLE]

The calligraphic symbols $\mathcal{P}^{S}_{M,M^{\prime}}$ and $\mathcal{P}_{M}$ are used to indicate that these operators are Hermitian idempotent projectors (if $M=M^{\prime}$) in the $N$-electron Hilbert space, whereas $\hat{P}^{S}_{M,M^{\prime}}$, which induces local mixtures of $\alpha$ and $\beta$ orbitals, is not a projector. To construct the SP-MPS using $\mathcal{P}^{S}_{M,M^{\prime}}$ (96), we require $M^{\prime}=M$ in $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ such that Eq. (56) becomes,

[TABLE]

The energy to be variationally optimized can then be considered as a functional of the underlying MPS $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$, and its explicit functional form reads as,

[TABLE]

where the fact that $\hat{H}$ is spin-free ($[\mathcal{P}^{S}_{M,M},\hat{H}]=0$) has been used, and $\mathcal{P}_{M}$ has been dropped in the last two identities due to the left projection to the bra state with a good quantum number $M$. Before proceeding to the sweep algorithm for optimizing $E[|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle]$, we should point out that similar to other spin projected implementations61, 62, 63, 64, 65, 66, the integration in $\hat{P}^{S}_{M,M}$ is carried out in practice by numerical quadrature. Specifically, we employ Gauss-Legendre quadrature via the transformation $x=\cos\beta\in[-1,1]$,

[TABLE]

where $N_{g}$ is the number of quadrature points. Eq. (102) in the MPO language becomes a sum of $N_{g}$ simple MPOs with bond dimension 1, since the exponential $e^{-\mathbbm{i}\beta\hat{S}_{y}}$ of a sum of local operators $\hat{S}_{y,k}$ is just a product of local operators (as they commute with each other) $e^{-\mathbbm{i}\beta\hat{S}_{y}}=e^{-\mathbbm{i}\beta\sum_{k=1}\hat{S}_{y,k}}=\prod_{k=1}e^{-\mathbbm{i}\beta\hat{S}_{y,k}}$, where the matrix representation of $e^{-\mathbbm{i}\beta\hat{S}_{y,k}}$ is

[TABLE]

From Eq. (107), we note that only real algebra is needed to implement Eq. (102), even though the imaginary unit appears in the formulation. For a system of $N$ electrons in $K$ spatial orbitals, the quadrature (102) is exact if $N_{g}$ is chosen to be at least

[TABLE]

where $\Omega_{max}$ is the maximal seniority number (number of singly occupied orbitals) and $\lceil x\rceil$ is the ceiling function. Eq. (110) follows from the Gauss quadrature rule, which is constructed to be exact for polynomials of degree $2n-1$ or lower74 for an $n$-point quadrature, and the observation that the integrand in either $\langle\Psi_{\mathrm{MPS}}^{(N,M)}|\hat{H}\hat{P}^{S}_{M,M}|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ or $\langle\Psi_{\mathrm{MPS}}^{(N,M)}|\hat{P}^{S}_{M,M}|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$ is a polynomial in $x=\cos\beta$, see Eq. (107), whose maximal degree is $\Omega_{max}/2+S$ determined by the maximally singly occupied configurations.

Combined with the sum of MPOs representation for $\hat{H}$ (6), the operator $\hat{H}\hat{P}^{S}_{M,M}$ in the numerator of Eq. (101) becomes a sum of MPOs with bond dimension $2KN_{g}$. Distributing groups of MPOs to different processors leads to an embarrassingly parallel scheme to compute expectation values over $|\Psi_{\mathrm{MPS}}^{(N,M)}\rangle$, which will be exploited in the following DMRG sweep optimizations. Based on these MPO representations for both $\hat{H}$ and $\hat{P}^{S}_{M,M}$, the energy functional (101) possesses a very nice graphical representation, as shown in Figure 3(a), as a quotient of two fully contracted tensor networks.

2.3 Sweep algorithms for ground and excited states

2.3.1 Ground-state DMRG optimization

The DMRG sweep algorithm can be used to minimize the spin-projected energy functional (101). Specifically, assuming the one-site formalism for simplicity in this discussion, the optimization of the set of site tensors $A[k]$ is carried out one at a time, and at site $k$, the local minimization problem corresponds to solving the following stationary condition,

[TABLE]

which, when $A^{n_{k}}_{l_{k-1}r_{k}}[k]$ is viewed as a vector, leads to a generalized eigenvalue problem,

[TABLE]

While the explicit algebraic form of the effective Hamiltonian $H_{\mathrm{eff}}$ and the metric $N_{\mathrm{eff}}$, involving numerous sums of products, is very lengthy, the graphical representation introduced earlier enables a very compact expression. This is depicted in Figure 3, where the yellow dot represents $A[k]$, and deleting it from the Figures 3(b) and 3(c) leads to $H_{\mathrm{eff}}$ and $N_{\mathrm{eff}}$, respectively.

We mention here that the metric $N_{\mathrm{eff}}$ does not arise from the choice of gauge for the MPS, as it comes from the introduction of the projector, as clearly shown by Figure 3(c). This means that during the optimization sweeps, we are free to use the mixed canonical form for the underlying MPS as usual, which ensures that the renormalized configuration basis is orthonormal, leading to more stable numerical algorithms. This differs from the situation arising in the optimization of MPS with periodic boundary conditions (PBC)75, 76, where the metric arises from the impossibility of choosing an orthonormal gauge due to the cyclic structure of the MPS with PBC. Thus, the attractive feature of SP-MPS is that almost all the usual DMRG machinery can be reused without any modification. In particular, after solving the generalized eigenvalue problem (112), the site tensor can be chosen in left or right canonical form by using a SVD or density matrix renormalization in exactly the same way as in the usual DMRG. More importantly, since the underlying MPS $|\Psi^{(N,M)}_{\mathrm{MPS}}\rangle$ possesses Abelian symmetries only, there is no need to use the singlet embedding scheme for non-singlet states, and the one-dot algorithm naturally leads to a consistent MPS at convergence that fully minimizes the energy functional with respect to all site tensors(101).

The computational cost for optimizing the SP-MPS scales as $O(N_{g}(D^{3}K^{3}+D^{2}K^{4}))$, which is a factor of $N_{g}$ higher than the usual DMRG. By using the sum of MPO representations for $\hat{H}$ (6) and $\hat{P}^{S}_{M,M}$ (102), the calculations can be parallelized easily over up to $2KN_{p}$ processors, where $K$ represents the number of spatial orbitals here. In the present pilot implementation, the actual cost is higher than $N_{g}$ times that of a normal DMRG calculation. This is because due to the presence of $\exp(-\mathbbm{i}\beta\hat{S}_{y})$, the trial vector for $|\Psi^{(N,M)}_{\mathrm{MPS}}\rangle$ at a given site needs to first be subducted to a lower symmetry with particle number symmetry only, and only after the application of $\hat{H}\exp(-\mathbbm{i}\beta\hat{S}_{y})$ to the trial vector, is the resulting vector projected back to the space with symmetry $(N,M)$. Therefore, in the matrix-vector product step, there is less symmetry to use compared with a normal DMRG calculation with symmetries $(N,M)$. However, this can in principle be alleviated by exploiting block data sparsity, since it is reasonable to imagine that as the bond dimension becomes large in a finite system, the underlying state $|\Psi^{(N,M)}_{\mathrm{MPS}}\rangle$ should be close to the target state, such that the action of $\exp(-\mathbbm{i}\beta\hat{S}_{y})$ will not introduce too many states in the other $M$ sectors. This strategy for reducing the computational cost will be explored in our future studies.

2.3.2 State-specific excited-state optimizations

One important feature of the combination of MPS with spin projection is that compared with other spin projected methods based on Slater determinants, it is more straightforward to compute excited states, due to the simplicity in imposing orthogonality constraints using MPS. Specifically, to compute the first excited state using SP-MPS, aside from using the state-averaged algorithm, one can directly target the state by imposing the constraint $\langle\Psi_{\mathrm{SP\text{-}MPS},0}^{(N,S,M)}|\Psi_{\mathrm{SP\text{-}MPS},1}^{(N,S,M)}\rangle=0$ between the excited state $|\Psi_{\mathrm{SP\text{-}MPS},1}^{(N,S,M)}\rangle$ to be optimized and the ground state $|\Psi_{\mathrm{SP\text{-}MPS},0}^{(N,S,M)}\rangle$, which is assumed optimized by the method introduced in the previous section. This simply requires that at each local optimization step, a vector $B^{n_{k}}_{l_{k-1}r_{k}}[k]$ is constructed in the following way,

[TABLE]

such that the local optimization based on Eq. (112) for $A[k]$ is subject to the constraint (113). The graphical representation for such a condition is shown in Figure 4. It is seen that all the environmental tensors connected to $A[k]$ (in yellow in Figure 4(a)) can be contracted into a three-way tensor $B[k]$ (open diamond in Figure 4(b)). This constraint allows to tackle excited states within the same symmetry as the ground state, and the generalization to multiple constraints is straightforward, viz., only a set of $B[k]$ vectors (113) corresponding to each constraint needs to be constructed, then a projector $Q=1-VV^{T}$ to implement the set of orthogonality constraints can be defined within the set of orthonormal basis vectors $V$, which can be obtained from the QR decomposition of the set of $B[k]$ vectors, which are in general not orthonormal.

Finally, it should be mentioned that to reduce the cost of DMRG optimizations based on SP-MPS for both ground and excited states, instead of solving the optimization problem (112), perturbation theory can be used to approximately decompose the original optimization problem into a variational step plus a perturbation correction step. The first step can be carried out with a small bond dimension $D_{0}$, while the latter step can be performed with a large bond dimension $D_{1}$ but using a much simpler zeroth order Hamiltonian to reduce the computational cost. Such an idea has been used in the MPSPT44. Although not fully explored in this paper, we present a possible generalization based on SP-MPS in Appendix Appendix 1: SP-MPS perturbation theory.

2.4 Properties

While we have shown that the SP-MPS can be easily used in every situation where the standard SA-MPS is currently used, we emphasize that the SP-MPS is not intended to be a replacement for SA-MPS. Rather, these two classes of MPS have completely different sets of merits and demerits, making them quite complementary in their applicability. In particular, the unique feature of SP-MPS is its closer connection with mean-field states, which could be exploited, for example, in the future development of DMET67, 68 for open-shell systems. To this end, we discuss here the evaluation of two essential ingredients in DMET using SP-MPS, namely, the one-body reduced density matrix (1RDM) $\langle a_{p_{\sigma}}^{\dagger}a_{q_{\tau}}\rangle$ and the energy component for a fragment $X$ of the whole molecule defined as a sum of expectation values for $\hat{H}_{p}$, viz., $e_{X}=\sum_{p\in X}\langle\hat{H}_{p}\rangle$.

For spin-independent (spin-free) operators $\hat{O}_{sf}$, due to the commutation relation $[\mathcal{P}^{S}_{M,M},\hat{O}_{sf}]=0$, similarly to the energy functional (101) the expectation value for $\hat{O}_{sf}$ can be simplified as

[TABLE]

Operators belonging to this case include $E_{pq}=\sum_{\sigma}a_{p\sigma}^{\dagger}a_{q\sigma}$, whose expectation value gives rise to the spin-free 1RDM, and spin-spin correlation functions $\vec{S}_{X}\cdot\vec{S}_{Y}$. For the energy component, if both the $\alpha$ and $\beta$ orbitals for the same spatial orbital are selected in the same group $X$, then $\sum_{p\in X}\hat{H}_{p}$ is also spin-free, such that the total energy (101) can be rewritten as

[TABLE]

This expression also fits into the parallelization scheme using the sum of MPOs representation, meaning that each $e_{p}$ can be evaluated independently and in parallel.

The evaluation of expectation values of spin-dependent operators is in general more complicated than for spin-free operators. One of the most important spin-dependent properties is the spin density matrix for nonsinglet states $\langle\Psi_{\mathrm{SP\text{-}MPS}}^{(N,S,M)}|T_{pq}(1,0)|\Psi_{\mathrm{SP\text{-}MPS}}^{(N,S,M)}\rangle$. Another important property is the spin-orbit coupling matrix element $\langle\Psi_{\mathrm{SP\text{-}MPS},I}^{(N,S,M)}|H_{SOC}|\Psi_{\mathrm{SP\text{-}MPS},J}^{(N,S^{\prime},M^{\prime})}\rangle$ between two spin states $I$ and $J$, which at the one-electron level only requires the transition density matrices of form $\langle\Psi_{\mathrm{SP\text{-}MPS},I}^{(N,S,M)}|T_{pq}(1,\mu)|\Psi_{\mathrm{SP\text{-}MPS},J}^{(N,S^{\prime},M^{\prime})}\rangle$ ($\mu=1,0,-1$). More specifically, assuming that only the high-spin reference MPS $|\Psi_{\mathrm{MPS}}^{(N,M=S)}\rangle$ is used for each SP-MPS, then only transition density matrices of the form $\langle\Psi_{\mathrm{SP\text{-}MPS},I}^{(N,S,M)}|T_{pq}(1,\mu)|\Psi_{\mathrm{SP\text{-}MPS},J}^{(N,S^{\prime}=S-\mu,M^{\prime}=M-\mu)}\rangle$ need be computed by virtue of the Wigner-Eckart theorem77. These properties can be obtained in a straightforward approach using double integrations61 to discretize both the bra and ket spin projectors, or by using a single integration based on (114), in conjunction with higher-order spin-free density matrices78, 79. However, in both approaches, the computational scaling is higher than that for evaluating the spin-free 1RDM (114). Fortunately, by using the transformation properties of $T_{pq}(1,\mu)$ under the action of spin rotations $\hat{R}(\Omega)$, the following expressions for the spin-dependent one-body (transition) density matrices can be derived (for details, see Appendix Appendix 2: Derivations for Eqs. (116), (117), and (118)),

[TABLE]

which are sufficient for the state interaction treatment of spin-orbit coupling77, 58 as well as the computation of spin density matrices (with Eq. (117)) for nonsinglet states. These formulae show that rather than changing the computational scaling, the evaluation of spin-dependent properties within SP-MPS only changes the prefactor by a small factor compared with the evaluation of spin-free properties.

3 Numerical examples

3.1 Two-dimensional Hubbard model

As a proof-of-principle calculation, we consider the ground state of the two-dimensional Hubbard model on a 4$\times$4 cluster with the Hamiltonian $\hat{H}_{\mathrm{Hubbard}}\triangleq-t\sum_{\langle ij\rangle}\sum_{\sigma}(a_{i\sigma}^{\dagger}a_{j\sigma}+H.c.)+U\sum_{i}n_{i\alpha}n_{i\beta}$, where $\langle ij\rangle$ represents the sum over nearest neighbor pairs, for various values of $U$ ($U=1,2,4,8,16$, $t=1$) and at half-filling. This model can be solved by exact diagonalization80. Here we use it to compare the performance of MPS without spin-adaptation (denoted MPS for brevity), spin-adapted MPS (SA-MPS), and spin-projected MPS (SP-MPS) for the singlet ground state. For SP-MPS, to focus on the representational properties of SP-MPS, we eliminate the angular integration error in (102) by using $N_{g}=5$ in the numerical quadrature, which is exact in this model according to Eq. (110). The convergence of the ground state energies using these three kinds of MPS in the site basis (ordered in a zigzag ordering for rows) as the bond dimension $D$ is increased is shown in Figure 5. It is clear that without spin adaptation, the convergence is painfully slow, see Figure 5(a). For instance, it is not possible to converge to 10*-5* ($t$) even for $D$ greater than 7000. As shown in Figure 5(b), spin adaptation significantly improves the convergence with respect to the bond dimension due to the use of reduced renormalized states. It can be seen from Figure 5(c) that the SP-MPS also accelerates the convergence as a result of the spin projection. To obtain a better comparison of SA-MPS and SP-MPS, the convergence for $U$=1, 8, and 16 is compared in Figure 6. In general, we find that SP-MPS achieves an accuracy of 10*-5* in the energy with a bond dimension 1.3-1.4 times that needed with SA-MPS. For small $U$(=1), the SP-MPS yields lower energies than SA-MPS for bond dimensions in the range 4000-5000, whereas for large $U$, while both SA-MPS and SP-MPS converge faster than in the corresponding $U=1$ case and significantly better than MPS without spin adaptation, the SA-MPS generally tends to perform better than SP-MPS.

To understand such differences, we can compute various properties of the converged SP-MPS wavefunctions. Figure 7(a) plots the expectation value of the seniority number operator $\langle\hat{\Omega}\rangle$ for SP-MPS (solid lines) and its underlying MPS (dashed lines). The operator $\hat{\Omega}$ is used to measure the number of singly occupied (open-shell) orbitals, which is also a sum of local operators,

[TABLE]

and hence can be written as an MPO with bond dimension 2. Figure 7(b) displays the von Neumann entropy defined by $S_{\mathrm{von\;Neumann}}=-\mathrm{tr}(\rho\ln\rho)=-\sum_{i}\lambda_{i}\ln\lambda_{i}$ of the underlying MPS with bond dimension 7600 at each bond, where $\rho$ is the reduced density matrix of the system or environment, and $\lambda_{i}$ is its associated eigenvalue. In Figure 7(c), the weight of the singlet components in the underlying MPS $\langle\mathcal{P}_{S=0}\rangle\triangleq\langle\Psi^{(N,M=0)}_{\mathrm{MPS}}|\mathcal{P}^{S=0}_{M=0,M=0}|\langle\Psi^{(N,M=0)}_{\mathrm{MPS}}\rangle$, which is also the overlap between SP-MPS, that is, the denominator of the energy functional (101), is compared for different values of $U$. According to Figure 7(a), the number of open-shell orbitals decreases as $U$ decreases due to the enhanced hopping to other sites, and the underlying state in the site basis becomes more entangled for small $U$, as demonstrated by the increased $S_{\mathrm{von\;Neumann}}$ shown in Figure 7(b). [NB: We emphasize that the definition of entanglement depends on the one-particle basis used to define the partitioning when computing $S_{\mathrm{von\;Neumann}}$. Thus, while the large $U$ case is more entangled in the momentum or mean-field basis, the small $U$ case is more entangled with the site basis employed here.] In such a situation, the underlying MPS tends to break spin symmetry, as demonstrated by the small values of the singlet component (0.3-0.4) in Figure 7(c) at small $U$, in order to describe the increased entanglement in the ground state wavefunctions and to recover more correlation energy. Therefore, we can expect the SP-MPS to perform better than SA-MPS in highly entangled situations, and to become less superior in less entangled systems.

3.2 Iron-sulfur clusters

Having established the basic convergence properties and representational power of SP-MPS, we consider two practical applications of SP-MPS to iron-sulfur clusters. These have been discussed as typical of challenging strongly correlated systems, because they have competing low-energy states corresponding to different kinds of spin and charge fluctuations44. The first application is to \ce[Fe2S2(SCH3)4]^2-, where two different initial guesses and two different active spaces are constructed to illustrate the ability of SP-MPS to describe the correct spin states. The second is to a larger cluster with four iron atoms \ce[Fe4S4(SCH3)4]^2-. Here we examine the possibility of using different broken symmetry initial guesses in conjunction with SP-MPS with very small bond dimensions, to construct a map of the physically relevant low energy states in the Hilbert space.

First, a small active space, CAS(10e,10o) with only 3$d_{\mathrm{Fe}}$ orbitals constructed from localized DFT (density functional theory) orbitals using the BP86 functional81, 82, the TZP-DKH basis83, and the sf-X2C (spin-free exact two-component) Hamiltonian84, 85 to include scalar relativistic effects, is considered with two different initial configurations, viz., Fe(III)-Fe(III) and Fe(II)-Fe(IV) shown in Figure 8(a). The bond dimension $D=20$ is used in all calculations, and the absolute errors on a logarithmic scale are shown in Figure 8(b). While both the MPS without spin projection and SP-MPS for the singlet state converge to an accuracy of 10*-5* Hartrees, the final converged states are actually qualitatively different. This is because while the exact energy within the active space of the singlet state is -27.887643 Hartrees, the state with $S=5$ representing the ferromagnetically coupled iron centers is of lower energy (-27.890357 Hartrees). Consequently, without the spin projection, the normal MPS calculation without spin adaptation converges to the lowest energy (high spin) state. Including sulfur 3$p_{\mathrm{S}}$ orbitals is important to recover the correct spin-state ordering, as the superexchange effect stabilizes the singlet state29. Despite such deficiencies of the small active space, it is interesting to see that while the energy of the second initial configuration is higher as expected, after several sweeps the SP-MPS (blue line) starts to relax to the same state as when starting from the first initial configuration (red line), showing that even with this small bond dimension it is possible to recover from the poor initial guess.

Including sulfur 3$p_{\mathrm{S}}$ orbitals gives rise to an enlarged active space CAS(30e,20o), and the energy convergence is shown in Figure 9 as a function of the bond dimension $D$. In this case, both MPS and SP-MPS converge to the same singlet ground state, while the latter converges faster with $D$. For comparison, the convergence of spin-adapted MPS (SA-MPS) is also depicted in the same figure (black lines). We see that for this case the SA-MPS converges to an accuracy of 10*-5* Hartrees much more quickly, mainly because the ground state is not very highly entangled. However, if the accuracy required is 10*-3* Hartrees, then it may be advantageous to use SP-MPS with different initial guesses to sample the space spanned by low-lying states, especially for larger iron-sulfur clusters with many competing local minima in the energy landscape of MPS with a fixed small bond dimension. Here, the connection of SP-MPS to the underlying broken symmetry determinants is crucial as it makes it easy to setup and enumerate all the possible starting low-energy broken symmetry configurations.

We next consider the complex \ce[Fe4S4(SCH3)4]^2- shown in Figure 10(a) using an active space CAS(54e,36o) with all 3$d_{\mathrm{Fe}}$ orbitals and 3$p_{\mathrm{S}}$ orbitals. For this complex, 24(=4!) different physically meaningful initial guesses can be constructed by distributing four different kinds of iron oxidation and spin states, viz., spin-up/down Fe(II) and spin-up/down Fe(III), in four different positions. SP-MPS calculations for singlet states with a small bond dimension $D=200$ are carried out starting from these broken symmetry initial guesses for $|\Psi^{(N,M)}_{\mathrm{MPS}}\rangle$. After convergence, the spin-spin correlation functions $\langle\vec{S}_{\mathrm{Fe}_{i}}\cdot\vec{S}_{\mathrm{Fe}_{j}}\rangle$ ($i,j\in\{1,2,3,4\}$) between the four irons are computed by Eq. (114), and the results are depicted in Figure 10(a), where the odd rows contain the spin-spin correlation patterns for the SP-MPS state with initial product state configurations, and the even rows contain the corresponding converged SP-MPS results with $D=200$. Clearly, there are three distinct patterns in the final results, and in all cases the charges on the irons delocalize. In Figure 10(b), the energies relative to the lowest SP-MPS energy (the 5th state) for all SP-MPS states are compared for the three patterns. We observe that the blue bars are associated with lower energies. Thus, the corresponding spin-spin correlation pattern has a higher chance to be the true pattern for the ground state. In fact, this is indeed the case as demonstrated in an earlier study29 as well as in a state-averaged SA-MPS calculation that we have performed for the lowest two states with $D=2500$, with the results shown in the right panel of Figure 10(a). Both calculations demonstrate that the ground state spin-spin correlation pattern shown in Figure 10(a) is the same as the pattern with the lowest energy discovered using the SP-MPS, where two pairs of parallel spins are anti-ferromagetically coupled to form a global singlet state. We thus see that SP-MPS with small bond dimension, starting from different physically motivated broken symmetry initial guesses, can help in identifying the low energy electronic structure when there are many competing low-lying states. Extending this approach to study even larger iron-sulfur clusters will be pursued in the future.

4 Conclusions and outlook

In this work, we have developed a versatile tool for strongly correlated systems by combining the idea of spin projection in quantum chemistry with the simplest TNS, the MPS, which has an underlying one dimensional connectivity. Such an approach could be even more advantageous within two possible generalizations. The first generalization is to extend to non-Abelian point group symmetries, which would otherwise involve the use of the generalized 6$j$ and 9$j$ symbols for the non-Abelian point groups if the symmetry adaptation is carried out in a similar fashion as SA-MPS. However, using symmetry projection, only the representation matrix $D^{\Gamma}(\hat{R})$ is needed, where $\Gamma$ represents the target irreducible representation (irrep). In view of the fact that for most of the non-Abelian point groups, the dimension of the degenerate irrep is two (actually the largest dimension is 5 for the irrep $H_{g}$ of the $I_{h}$ group86), the bond dimension to represent each symmetry operation is small so long as orbitals belonging to the same irrep are placed in adjacent sites. The second generalization is to construct spin eigenfunctions for higher-dimensional tensor network states such as the PEPS69. This is also quite straightforward due to the simple structure of the “projector” $\hat{P}^{S}_{M,M}$ (102) that is a sum of products of local spin rotations. For instance, the overlap between two PEPS with projectors (102), $\langle\Psi_{\mathrm{PEPS}}|e^{-\mathbbm{i}\beta_{g}\hat{S}_{y}}|\Psi_{\mathrm{PEPS}}\rangle$ shown in Figure 11, reveals that the introduction of spin projectors adds no complications to the contraction of PEPS.

Regarding SP-MPS themselves, we have shown that they possess several distinct features that are not shared by the spin-adapted MPS with non-Abelian symmetry, such as a simple formulation and implementation which only requires the use of Abelian symmetries for the underlying MPS, and avoids the use of the singlet embedding scheme. Perhaps the most important feature is the close connection to traditional “broken symmetry” determinants. This gives the ability to seed SP-MPS from initial “broken symmetry” determinants, providing a route to connect chemical intuition about broken symmetry configurations to realistic calculations that properly incorporate fluctuations and correlations. This further opens up the possibility to fully map out the low energy landscape of competing states in finite chemical systems, in particular the polymetallic transition metal compounds, similar to what is already done in condensed phase problems 87. In such applications, the computed energies of SP-MPS can be improved by using the SP-MPS based perturbation theory, see Appendix Appendix 1: SP-MPS perturbation theory. Further, a combination of SP-MPS and SA-MPS is also possible by using the optimized SP-MPS with small bond dimension to initialize spin-adapted DMRG calculations with larger bond dimensions, to improve the computed energies efficiently. The connection to broken symmetry mean-field states may also help in the future development of quantum embedding methods for open-shell systems. Specifically, the SP-MPS can be used as an impurity solver for the embedded system in DMET, and such a combination may help overcome the difficulties in using MPS to describe very high-dimensional entanglement, as can be found in large transition metal clusters. These directions are being explored in our laboratory.

Appendix 1: SP-MPS perturbation theory

We here describe how to incorporate SP-MPS into the MPSPT framework44, and in particular, how to deal with the complications arising from the spin projectors. In the MPSPT, once a partition of $\hat{H}=\hat{H}_{0}+\hat{V}$ is given, the first order wavefunction represented in the MPS form is to be obtained by minimizing the Hylleraas functional44 (assuming real MPS for simplicity),

[TABLE]

In the SP-MPS case, both $|\Psi_{0}\rangle$ and $|\Psi_{1}\rangle$ are to be represented by the ansatz (56), and the only slight modification is that the reference $|\Psi_{0}\rangle$ is explicitly normalized via,

[TABLE]

The orthogonality constraint in Eq. (124) can be implemented as in the excited-state calculations discussed above. The key to avoid complications due to double integrations for $L[\Psi_{1}]$ is to choose a good spin-free $\hat{H}_{0}$. Then the perturbation $\hat{V}=\hat{H}-\hat{H}_{0}$ will also be spin-free, such that in both terms of $L[\Psi_{1}]$ (124), the commutator $[\mathcal{P}^{S}_{M,M},\hat{O}_{sf}]=0$ can be used to bring $L[\Psi_{1}]$ into a form that only involves a single integration in $\hat{P}^{S}_{M,M}$ similar to that in the energy functional (101). Then, the stationary condition for Eq. (124) will lead to a simple linear equation to be solved in each local optimization problem, i.e.,

[TABLE]

subject to the orthogonality constraint similar to Eq. (113). Here, all MPS in Eq. (126) refer to the underlying MPS such that the equation can be solved similarly to in MPSPT.

For spin-free $\hat{H}_{0}$, we propose to use a simple analog of the Epstein-Nesbet partition in the space of determinants,

[TABLE]

where $h_{pq}$ and $[pq|rs]$ are the one- and two-electron integrals over spatial orbitals. However, while the first two terms are singlet operators, which can be reexpressed by using $E_{pq}$, the last term containing $\sum_{\sigma}n_{p_{\sigma}}n_{q_{\sigma}}$ is not a pure singlet operator. To show this, it is rewritten into a combination of spin tensor operators,

[TABLE]

where $S_{pq}$ and $T_{pq}$ are singlet and triplet operators, respectively. Extracting the singlet component of $T_{pp}(1,0)T_{qq}(1,0)$, viz., $[T_{pp}(1)\times T_{qq}(1)]^{0}_{0}C^{00}_{10,10}$ with $C^{S_{3}M_{3}}_{S_{1}M_{1},S_{2}M_{2}}$ being the Clebsch-Gordan coefficient, gives the following form for the singlet component of $\hat{H}_{d}$,

[TABLE]

which can be rewritten as a sum of $2K$ MPO by the same splitting as in Eq. (6) for $\hat{H}$, each of which is of bond dimension $5$. Compared with the bond dimension $O(K^{2})$ for $\hat{H}$, it is seen that $\hat{H}_{d}^{\mathrm{singlet}}$ is a significant simplification. For $\hat{H}_{0}$ used in the Hylleraas function (124), taking into account the fact that $|\Psi_{0}\rangle$ is not an eigenfunction of $\hat{H}_{d}^{\mathrm{singlet}}$, the following form could be chosen,

[TABLE]

where $E_{0}$ could be chosen as the DMRG energy for $|\Psi_{0}\rangle$ for simplicity. The performance of such perturbation theory will be studied in future.

Appendix 2: Derivations for Eqs. (116), (117), and (118)

To derive a general expression for $\langle\Psi_{\mathrm{SP\text{-}MPS},I}^{(N,S,M)}|T_{pq}(1,\mu)|\Psi_{\mathrm{SP\text{-}MPS},J}^{(N,S^{\prime},M^{\prime})}\rangle$, we first examine the action of $\mathcal{P}^{S}_{M,M}$ on $T_{pq}(1,\mu)$,

[TABLE]

The term in the bracket can be recast into a linear combination of projectors, by using the following Clebsch-Gordan series88,

[TABLE]

where $C_{S_{1}M_{1},S_{2}M_{2}}^{S_{3}M_{3}}$ represents the (real) Clebsch-Gordan coefficients. With the decomposition (136), Eq. (135) becomes,

[TABLE]

Using this result, the matrix elements $\langle\Psi_{\mathrm{SP\text{-}MPS},I}^{(N,S,M)}|T_{pq}(1,\mu)|\Psi_{\mathrm{SP\text{-}MPS},J}^{(N,S^{\prime},M^{\prime})}\rangle$ are derived as

[TABLE]

where the projector property $\mathcal{P}^{S}_{M_{1},M_{2}}\mathcal{P}^{S^{\prime}}_{M_{1}^{\prime},M_{2}^{\prime}}=\delta_{SS^{\prime}}\delta_{M_{2}M_{1}^{\prime}}\mathcal{P}^{S}_{M_{1},M_{2}^{\prime}}$ has been used. Note that the nonvanishing condition for the Clebsch-Gordan coefficient imposes the triangular condition for angular momentum coupling in Eq. (138), $|S-1|\leq S^{\prime}\leq S+1$, otherwise the matrix elements are zero. If $M$ and $M^{\prime}$ are further limited to the high spin case, Eq. (138) simplifies to

[TABLE]

which gives rise to Eqs. (116)-(118) by substituting in the value of $\mu=1,0,-1$, respectively.

{acknowledgement}

This work was supported by the NSF through award SI2-SSI:1657286, with additional support from award 1650436.