Time-reversal symmetry adaptation in relativistic density matrix renormalization group algorithm

We assume the one-electron basis has a Kramers paired structure $\{\psi_{p},\psi_{\bar{p}}\}$, which can either be spin-orbitals or spinors computed from spin-restricted or Kramers-restricted self-consistent field calculations, respectively. The action of the time-reversal symmetry operator $\mathcal{T}$ on spin-orbitals/spinors reads

$\displaystyle\mathcal{T}|\psi_{p}\rangle=|\psi_{\bar{p}}\rangle,\quad\mathcal{T}|\psi_{\bar{p}}\rangle=-|\psi_{p}\rangle.$

(23)

In the absence of magnetic field, $\mathcal{T}$ commutes with the Hamiltonian $\hat{H}$, which is written in a second quantized form as

$\displaystyle\hat{H}=\sum_{pq}h_{pq}a_{p}^{\dagger}a_{q}+\frac{1}{4}\sum_{pqrs}\langle pq\|sr\rangle a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}.$

(24)

To get the representation of $H$ in the direct product space $V_{l}\otimes V_{r}$, $\hat{H}$ is usually rewritten as

$$\begin{array}[]{rcl}\hat{H}&=&\hat{H}^{\mathrm{L}}+\hat{H}^{\mathrm{R}}\\ &+&\sum_{p_{\mathrm{L}}}(a_{p_{\mathrm{L}}}^{\mathrm{L}\dagger}S_{p_{\mathrm{L}}}^{\mathrm{R}}+\mathrm{h.c.})+\sum_{q_{\mathrm{R}}}(a_{q_{\mathrm{R}}}^{\mathrm{R}\dagger}S_{q_{\mathrm{R}}}^{\mathrm{L}}+\mathrm{h.c.})\\ &+&\sum_{p_{\mathrm{L}}<q_{\mathrm{L}}}(A_{p_{\mathrm{L}}q_{\mathrm{L}}}^{\mathrm{L}}P_{p_{\mathrm{L}}q_{\mathrm{L}}}^{\mathrm{R}}+\mathrm{h.c.})\\ &+&\sum_{p_{\mathrm{L}}\leq s_{\mathrm{L}}}w_{p_{\mathrm{L}}s_{\mathrm{L}}}(B_{p_{\mathrm{L}}s_{\mathrm{L}}}^{\mathrm{L}}Q_{p_{\mathrm{L}}s_{\mathrm{L}}}^{\mathrm{R}}+\mathrm{h.c.}),\end{array}$$

(25)

where $w_{pq}=1-\frac{1}{2}\delta_{pq}=w_{qp}$, $p_{\mathrm{L}}$ ($p_{\mathrm{R}}$) represents the index of one-electron basis for the subspace $V_{l}$ ($V_{r}$), and the introduced intermediates (normal and complementary operators^{56, 13}) are

$$\begin{array}[]{rcl}A_{pq}&\triangleq&a_{p}^{\dagger}a_{q}^{\dagger},\\ B_{ps}&\triangleq&a_{p}^{\dagger}a_{s},\\ P_{pq}&\triangleq&\sum_{s<r}\langle pq\|sr\rangle a_{r}a_{s},\\ Q_{ps}&\triangleq&\sum_{qr}\langle pq\|sr\rangle a_{q}^{\dagger}a_{r},\\ S_{p}&\triangleq&\sum_{q}\frac{1}{2}h_{pq}a_{q}+\sum_{q,s<r}\langle pq\|sr\rangle a_{q}^{\dagger}a_{r}a_{s}.\end{array}$$

(26)

The time-reversal invariance of the Hamiltonian $\hat{H}=\mathcal{T}\hat{H}\mathcal{T}^{-1}$ is not obvious in this form. Since the integrals satisfy time-reversal symmetry, viz., $h_{pq}^{*}=h_{\bar{p}\bar{q}}$ and $h_{p\bar{q}}^{*}=-h_{\bar{p}q}$, we can rewrite Eq. (25) as

$\displaystyle\begin{array}[]{rcl}\hat{H}&=&\tilde{H}+\mathcal{T}\tilde{H}\mathcal{T}^{-1},\\ \tilde{H}&\triangleq&\frac{1}{2}(\hat{H}^{\mathrm{L}}+\hat{H}^{\mathrm{R}})+(H_{aS}+H_{AP}+H_{BQ}),\\ H_{aS}&\triangleq&\sum_{p_{\mathrm{L}}}a_{p_{\mathrm{L}}}^{\mathrm{L}\dagger}S_{p_{\mathrm{L}}}^{\mathrm{R}}+\sum_{p_{\mathrm{R}}}a_{p_{\mathrm{R}}}^{\mathrm{R}\dagger}S_{p_{\mathrm{R}}}^{\mathrm{L}}+\mathrm{h.c.},\\ H_{AP}&\triangleq&\sum_{p_{\mathrm{L}}<q_{\mathrm{L}}}A^{\mathrm{L}}_{p_{\mathrm{L}}q_{\mathrm{L}}}P^{\mathrm{R}}_{p_{\mathrm{L}}q_{\mathrm{L}}}\\ &&+\sum_{p_{\mathrm{L}}\leq q_{\mathrm{L}}}w_{p_{\mathrm{L}}q_{\mathrm{L}}}A^{\mathrm{L}}_{p_{\mathrm{L}}\bar{q}_{\mathrm{L}}}P^{\mathrm{R}}_{p_{\mathrm{L}}\bar{q}_{\mathrm{L}}}+\mathrm{h.c.},\\ H_{BQ}&\triangleq&\sum_{p_{\mathrm{L}}\leq s_{\mathrm{L}}}w_{p_{\mathrm{L}}s_{\mathrm{L}}}B^{\mathrm{L}}_{p_{\mathrm{L}}s_{\mathrm{L}}}Q^{\mathrm{R}}_{p_{\mathrm{L}}s_{\mathrm{L}}}\\ &&+\sum_{p_{\mathrm{L}}\leq s_{\mathrm{L}}}w_{p_{\mathrm{L}}s_{\mathrm{L}}}B^{\mathrm{L}}_{p_{\mathrm{L}}\bar{s}_{\mathrm{L}}}Q^{\mathrm{R}}_{p_{\mathrm{L}}\bar{s}_{\mathrm{L}}}+\mathrm{h.c.},\end{array}$

(34)

using the derived time-reversal symmetry properties of intermediates such as

	$\displaystyle\mathcal{T}A_{pq}\mathcal{T}^{-1}$	$\displaystyle=$	$\displaystyle\mathcal{T}(p^{\dagger}q^{\dagger})\mathcal{T}^{-1}=\bar{p}^{\dagger}\bar{q}^{\dagger}=A_{\bar{p}\bar{q}},$		(35)
	$\displaystyle\mathcal{T}A_{p\bar{q}}\mathcal{T}^{-1}$	$\displaystyle=$	$\displaystyle\mathcal{T}(p^{\dagger}\bar{q}^{\dagger})\mathcal{T}^{-1}=-\bar{p}^{\dagger}q^{\dagger}=-A_{\bar{p}q}.$		(36)

The obtained ’skeleton’ operator $\tilde{H}$ is Hermitian but not time-reversal invariant, but the number of operators in $\tilde{H}$ is roughly half of that in $\hat{H}$. This form will be used later to reduce the computational cost of R-DMRG.

To solve Eq. (22) in $V_{l}\otimes V_{r}$ with $\hat{H}$ in Eq. (34), we can use the iterative Davidson algorithm, where in each step the so-called $\sigma$ vector needs to be formed $|\sigma\rangle=\hat{H}|\Psi\rangle$. For the even-electron case, the coefficients of $\sigma$ are

	$\displaystyle\sigma_{lr}^{\mathrm{e}}$	$\displaystyle\triangleq$	$\displaystyle\langle lr|\tilde{H}+\mathcal{T}\tilde{H}\mathcal{T}^{-1}|\Psi_{\mathrm{e}}\rangle$		(37)
		$\displaystyle=$	$\displaystyle\langle lr|\tilde{H}|\Psi_{\mathrm{e}}\rangle+\langle\bar{l}\bar{r}|\tilde{H}|\Psi_{\bar{\mathrm{e}}}\rangle^{*}.$

Since we require $|\Psi_{\mathrm{e}}\rangle=|\Psi_{\bar{\mathrm{e}}}\rangle$, it simply becomes

$\displaystyle\sigma_{lr}^{\mathrm{e}}=\tilde{\sigma}_{lr}^{\mathrm{e}}+\tilde{\sigma}_{\bar{l}\bar{r}}^{\mathrm{e}*},$

(38)

with $\tilde{\sigma}_{lr}^{\mathrm{e}}\triangleq\langle lr|\tilde{H}|\Psi_{\mathrm{e}}\rangle$ This shows that only $\tilde{\sigma}_{lr}^{\mathrm{e}}$ needs to be constructed, whose computational cost is roughly half of that for constructing $\sigma_{lr}^{\mathrm{e}}$. Similarly, for odd-electron systems, we can find

	$\displaystyle\sigma_{lr}^{\mathrm{o}}$	$\displaystyle\triangleq$	$\displaystyle\langle lr|\tilde{H}+\mathcal{T}\tilde{H}\mathcal{T}^{-1}|\Psi_{\mathrm{o}}\rangle$		(39)
		$\displaystyle=$	$\displaystyle\langle lr|\tilde{H}|\Psi^{\mathrm{o}}\rangle+\langle\bar{l}\bar{r}|\tilde{H}|\Psi_{\bar{\mathrm{o}}}\rangle^{*}=\tilde{\sigma}_{lr}^{\mathrm{o}}+\tilde{\sigma}_{\bar{l}\bar{r}}^{\bar{\mathrm{o}}*},$
	$\displaystyle\sigma_{lr}^{\bar{\mathrm{o}}}$	$\displaystyle\triangleq$	$\displaystyle\langle lr|\tilde{H}+\mathcal{T}\tilde{H}\mathcal{T}^{-1}|\Psi_{\bar{\mathrm{o}}}\rangle$		(40)
		$\displaystyle=$	$\displaystyle\langle lr|\tilde{H}|\Psi_{\bar{\mathrm{o}}}\rangle-\langle\bar{l}\bar{r}|\tilde{H}|\Psi_{\mathrm{o}}\rangle^{*}=\tilde{\sigma}_{lr}^{\bar{\mathrm{o}}}-\tilde{\sigma}_{\bar{l}\bar{r}}^{\mathrm{o}*}.$

Thus, it suffices to construct $\tilde{\sigma}_{lr}^{\mathrm{o}}$ and $\tilde{\sigma}_{lr}^{\bar{\mathrm{o}}}$, which reduces the computational cost for constructing $\sigma_{lr}^{\mathrm{o}}$ and $\sigma_{lr}^{\bar{\mathrm{o}}}$ roughly by half.

In summary, the full $\sigma$ can be recovered from the skeleton one $\tilde{\sigma}$ by a ’time-reversal symmetrization’ in both even- and odd-electron cases. This is in a similar spirit to the construction of Fock matrix using the skeleton-matrix algorithm^{57, 58, 59}. Expressions for $\tilde{H}$ in Eq. (34) immediately show that only half of the intermediate operators are necessary for constructing $\tilde{H}$, which reduces the memory and computational cost for renormalized operators by half compared with an implementation without using time-reversal symmetry.

To use such reduction, we need to maintain the structure (1) for basis vectors of the subspace in Davidson algorithm. For the even-electron case, suppose the current subspace $V=\mathrm{span}\{|b_{k}\rangle\}$ in Davidson algorithm is spanned by time-reversal invariant basis, then the representation of $\hat{H}$ is a real matrix,

$\displaystyle\langle b_{k}|H|b_{l}\rangle=(\mathcal{T}\langle b_{k}|H|b_{l}\rangle)^{*}=\langle\bar{b}_{k}|H|\bar{b}_{l}\rangle=\langle b_{k}|H|b_{l}\rangle.$

(41)

Consequently, the eigenvectors $\mathbf{X}$ of $\mathbf{H}$ are real and the states $|x_{i}\rangle=\sum_{k}|b_{k}\rangle X_{ki}$ are time-reversal invariant. It can be verified that the residual $|r_{i}\rangle=\hat{H}|x_{i}\rangle-|x_{i}\rangle E_{i}$ is also time-reversal invariant, so does the precondition residual as $\langle lr|H|lr\rangle=\langle\bar{l}\bar{r}|H|\bar{l}\bar{r}\rangle$.

For the odd-electron case, suppose the current subspace $V=\mathrm{span}\{|b_{k}\rangle\}\oplus\mathrm{span}\{|\bar{b}_{k}\rangle\}$ is spanned by a Kramers paired basis, i.e., $|\bar{b}_{k}\rangle=\mathcal{T}|b_{k}\rangle$, then the representation of Hamiltonian has a quaternion structure (6). Diagonalizing it with a structure-preserving algorithm will produced Kramers paired eigenvectors $\{|x_{i}\rangle,|\bar{x}_{i}\rangle\}$. Similarly, one can show that the residuals also form a Kramers pair, $|\bar{r}_{i}\rangle=\mathcal{T}|r_{i}\rangle$. An important point for constructing new Kramers paired orthonormal basis is that if $|r_{i}\rangle$ is already orthonormalized against $V=\mathrm{span}\{|b_{k}\rangle\}\oplus\mathrm{span}\{|\bar{b}_{k}\rangle\}$, then $|\bar{r}_{i}\rangle$ will be automatically orthogonal to the basis vectors in $V$, such that the pair $(|r_{i}\rangle,|\bar{r}_{i}\rangle)$ can be added simultaneously. This property can be seen from

	$\displaystyle\langle\bar{r}_{i}|b_{k}\rangle$	$\displaystyle=$	$\displaystyle(\mathcal{T}\langle\bar{r}_{i}|b_{k}\rangle)^{*}=-\langle r_{i}|\bar{b}_{k}\rangle=0,$		(42)
	$\displaystyle\langle\bar{r}_{i}|\bar{b}_{k}\rangle$	$\displaystyle=$	$\displaystyle(\mathcal{T}\langle\bar{r}_{i}|\bar{b}_{k}\rangle)^{*}=\langle r_{i}|b_{k}\rangle=0,$		(43)

and $\langle\bar{r}_{i}|r_{i}\rangle=0$ due to Kramers’ theorem. In this way, we can maintain the structure (1) for the basis vectors of the subspace in the Davidson diagonalization algorithm.

Once the eigenvectors $|\Psi\rangle$ (18) of $\hat{H}$ have been found in the direct product space $V_{l}\otimes V_{r}$, we need to perform decimation to obtain an optimized basis for $V_{l}$ or $V_{r}$. This is done by diagonalizing the reduced density matrix $\rho_{l}$ or $\rho_{r}$. In the sweep optimization of DMRG, $V_{l}$ (or $V_{r}$) is also a direct product space denoted by $V_{D}$. It is formed by a direct product between the left environment and the left dot, referred as superblock. Simply diagonalizing the reduced density matrix will not yield a basis with the structure (1). We will show how to perform decimation in such space to produce time-reversal symmetry-adapted renormalized states.

For the even-electron subspace $V_{D}^{\mathrm{e}}$, we assume it is spanned by the following direct product basis

$\displaystyle V_{D}^{\mathrm{e}}=\mathrm{span}\{|D_{\mathrm{e}}\rangle,|D_{\bar{\mathrm{e}}}\rangle,|D_{0}\rangle\},$

(44)

with $|D_{\bar{\mathrm{e}}}\rangle=\mathcal{T}|D_{\mathrm{e}}\rangle$ and $|D_{0}\rangle=\mathcal{T}|D_{0}\rangle$. Here, $|D_{0}\rangle$ represents the part of direct product basis which is already time-reversal invariant such as $|l_{\mathrm{e}}r_{\mathrm{e}}\rangle$ in Eq. (9). The pair $|D_{\mathrm{e}}\rangle$ and $|D_{\bar{\mathrm{e}}}\rangle$ represent those parts which can be related by $\mathcal{T}$ such as $|l_{\mathrm{o}}r_{\mathrm{o}}\rangle$ and $|l_{\bar{\mathrm{o}}}\bar{r}_{\mathrm{o}}\rangle$. Suppose the reduced density matrix obtained from $|\Psi\rangle$ in $V_{D}$ is

$\displaystyle\bm{\rho}^{\Psi}=\left[\begin{array}[]{ccc}\rho_{\mathrm{e}\mathrm{e}}&\rho_{e\bar{\mathrm{e}}}&\rho_{e0}\\ \rho_{\bar{\mathrm{e}}e}&\rho_{\bar{\mathrm{e}}\bar{\mathrm{e}}}&\rho_{\bar{\mathrm{e}}0}\\ \rho_{0e}&\rho_{0\bar{\mathrm{e}}}&\rho_{00}\\ \end{array}\right],$

(48)

then that obtained from $|\bar{\Psi}\rangle=\mathcal{T}|\Psi\rangle$ is

$\displaystyle\bm{\rho}^{\bar{\Psi}}=\left[\begin{array}[]{ccc}\rho^{*}_{\bar{\mathrm{e}}\bar{\mathrm{e}}}&\rho^{*}_{\bar{\mathrm{e}}e}&\rho^{*}_{\bar{\mathrm{e}}0}\\ \rho^{*}_{e\bar{\mathrm{e}}}&\rho^{*}_{\mathrm{e}\mathrm{e}}&\rho^{*}_{e0}\\ \rho^{*}_{0\bar{\mathrm{e}}}&\rho^{*}_{0e}&\rho^{*}_{00}\\ \end{array}\right].$

(52)

The average $\bm{\rho}=\frac{1}{2}(\bm{\rho}^{\Psi}+\bm{\rho}^{\bar{\Psi}})$ is time-reversal invariant

$\displaystyle\bm{\rho}=\frac{1}{2}\left[\begin{array}[]{ccc}\rho_{\mathrm{e}\mathrm{e}}+\rho^{*}_{\bar{\mathrm{e}}\bar{\mathrm{e}}}&\rho_{e\bar{\mathrm{e}}}+\rho^{*}_{\bar{\mathrm{e}}e}&\rho_{e0}+\rho^{*}_{\bar{\mathrm{e}}0}\\ \rho_{\bar{\mathrm{e}}e}+\rho^{*}_{e\bar{\mathrm{e}}}&\rho_{\bar{\mathrm{e}}\bar{\mathrm{e}}}+\rho^{*}_{\mathrm{e}\mathrm{e}}&\rho_{\bar{\mathrm{e}}0}+\rho^{*}_{e0}\\ \rho_{0e}+\rho^{*}_{0\bar{\mathrm{e}}}&\rho_{0\bar{\mathrm{e}}}+\rho^{*}_{0e}&\rho_{00}+\rho^{*}_{00}\\ \end{array}\right]\triangleq\left[\begin{array}[]{ccc}\mathcal{A}&\mathcal{B}&\mathcal{C}\\ \mathcal{B}^{*}&\mathcal{A}^{*}&\mathcal{C}^{*}\\ \mathcal{C}^{\dagger}&\mathcal{C}^{T}&\mathcal{E}\\ \end{array}\right],$

(59)

where $\mathcal{A}=\mathcal{A}^{\dagger}$, $\mathcal{B}=\mathcal{B}^{T}$, and $\mathcal{E}$ is real symmetric. However, simply diagonalizing it with a complex eigensolver will not produce time-reversal invariant basis function (1) due to the arbitrariness of the phase factor. To fix this problem, we can introduce a time-reversal invariant basis similar to Eq. (11),

$\displaystyle(|R_{-}\rangle,|R_{+}\rangle,|R_{0}\rangle)\triangleq(|R_{\mathrm{e}}\rangle,|R_{\bar{\mathrm{e}}}\rangle,|R_{0}\rangle)\mathbf{U}$

(60)

with $\mathbf{U}$ being

$\displaystyle\mathbf{U}=\left[\begin{array}[]{ccc}\frac{\mathbbm{i}}{\sqrt{2}}I&\frac{1}{\sqrt{2}}I&0\\ -\frac{\mathbbm{i}}{\sqrt{2}}I&\frac{1}{\sqrt{2}}I&0\\ 0&0&I\\ \end{array}\right],$

(64)

and transform $\bm{\rho}$ into this basis, which leads to a real symmetric reduced density matrix

$\displaystyle\tilde{\bm{\rho}}=\mathbf{U}^{\dagger}\bm{\rho}\mathbf{U}=\left[\begin{array}[]{ccc}(\mathcal{A}-\mathcal{B})_{\mathrm{R}}&(\mathcal{A}+\mathcal{B})_{\mathrm{I}}&\sqrt{2}\mathcal{C}_{\mathrm{I}}\\ -(\mathcal{A}-\mathcal{B})_{\mathrm{I}}&(\mathcal{A}+\mathcal{B})_{\mathrm{R}}&\sqrt{2}\mathcal{C}_{\mathrm{R}}\\ \sqrt{2}\mathcal{C}_{\mathrm{I}}^{T}&\sqrt{2}\mathcal{C}_{\mathrm{R}}^{T}&\mathcal{E}\\ \end{array}\right],$

(68)

where $\mathcal{A}_{\mathrm{R}}$ (or $\mathcal{A}_{\mathrm{I}}$) represents the real (or imaginary) part. Diagonalizing $\tilde{\bm{\rho}}\mathbf{X}=\mathbf{X}\mathbf{\Lambda}$ yields a set of real vectors $\mathbf{X}$ in the time-reversal invariant basis, which can be back transformed to the original direct product basis (44) by

$\displaystyle\mathbf{UX}=\left[\begin{array}[]{ccc}\frac{\mathbbm{i}}{\sqrt{2}}I&\frac{1}{\sqrt{2}}I&0\\ -\frac{\mathbbm{i}}{\sqrt{2}}I&\frac{1}{\sqrt{2}}I&0\\ 0&0&1\\ \end{array}\right]\left[\begin{array}[]{c}\mathbf{X}_{-}\\ \mathbf{X}_{+}\\ \mathbf{X}_{0}\\ \end{array}\right]\equiv\left[\begin{array}[]{c}\mathbf{X}_{\mathrm{e}}\\ \mathbf{X}^{*}_{\mathrm{e}}\\ \mathbf{X}_{0}\\ \end{array}\right],$

(78)

with $\mathbf{X}_{\mathrm{e}}=(\mathbf{X}_{+}+\mathbbm{i}\mathbf{X}_{-})/\sqrt{2}$.

For the odd-electron subspace $V_{D}^{\mathrm{o}}$, we assume it is spanned by the following direct product basis

$\displaystyle V_{D}^{\mathrm{o}}=\mathrm{span}\{|D_{\mathrm{o}}\rangle,|D_{\bar{\mathrm{o}}}\rangle\},$

(79)

which is already Kramers paired, see Eq. (9). The reduced density matrices are

$\displaystyle\bm{\rho}^{\Psi}=\left[\begin{array}[]{cc}\rho_{\mathrm{o}\mathrm{o}}&\rho_{\mathrm{o}\bar{\mathrm{o}}}\\ \rho_{\bar{\mathrm{o}}\mathrm{o}}&\rho_{\bar{\mathrm{o}}\bar{\mathrm{o}}}\\ \end{array}\right],\quad\bm{\rho}^{\bar{\Psi}}=\left[\begin{array}[]{cc}\rho^{*}_{\bar{\mathrm{o}}\bar{\mathrm{o}}}&-\rho^{*}_{\bar{\mathrm{o}}\mathrm{o}}\\ -\rho^{*}_{\mathrm{o}\bar{\mathrm{o}}}&\rho^{*}_{\mathrm{o}\mathrm{o}}\\ \end{array}\right],$

(84)

and the average $\bm{\rho}$ has a quaternion structure

$\displaystyle\bm{\rho}=\frac{1}{2}\left[\begin{array}[]{cc}\rho_{\mathrm{o}\mathrm{o}}+\rho^{*}_{\bar{\mathrm{o}}\bar{\mathrm{o}}}&\rho_{\mathrm{o}\bar{\mathrm{o}}}-\rho^{*}_{\bar{\mathrm{o}}\mathrm{o}}\\ \rho_{\bar{\mathrm{o}}\mathrm{o}}-\rho^{*}_{\mathrm{o}\bar{\mathrm{o}}}&\rho_{\bar{\mathrm{o}}\bar{\mathrm{o}}}+\rho^{*}_{\mathrm{o}\mathrm{o}}\\ \end{array}\right]\triangleq\left[\begin{array}[]{cc}\mathcal{A}&\mathcal{B}\\ -\mathcal{B}^{*}&\mathcal{A}^{*}\\ \end{array}\right].$

(89)

Thus, diagonalizing it with a structural preserving algorithm^{47, 48, 49, 50, 51, 52, 53} will lead to a new set of Kramers paired renormalized states.

The above decimation procedure is quite general. We can even apply the decimation procedure for cases where $|\Psi\rangle$ is Kramers symmetry contaminated to produce a time-reversal symmetry-adapted basis. For instance, it can be applied to convert a Kramers symmetry contaminated selected configuration interaction wavefunctions to time-reversal symmetry-adapted MPS as the initial guess for R-DMRG⁴⁶. Then, by iterating the above diagonalization and decimation procedure, the time-reversal symmetry structure of basis functions (1) can be maintained recursively during the sweep optimization in DMRG. For other tensor network states without loops^{41, 42, 43, 44, 45, 46}, it is clear that the construction of time-reversal symmetry-adapted basis can be directly applied.