# Hubbard Model

```
[1]:
```

```
!pip install block2==0.5.2rc13 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
```

## Introduction

In this tutorial we introduce the python interface of `block2`

. This interface allows the user to design custom workflow for running DMRG. It is also more convenient and efficient for model Hamiltonians like Hubbard and Heisenberg models, and non-standard quantum chemistry models with non-Hermitian or complex number integrals.

In `block2`

we support a few different symmetry modes. For fermionic models, if the Hamiltonian conserves the total spin \(S\), the projected spin \(S_z\), and the number of electrons \(N\), we can label the states using quantum numbers \((N, S_z)\) or \((N, S)\) or simply \(N\).

The symmetry group for \((N, S_z)\) is \(\mathrm{U(1)\times U(1)}\), which is Abelian. This symmetry mode is activated using

`symm_type=SymmetryTypes.SZ`

. We will first introduce the DMRG calculation using this symmetry since it is the most convenient choice.If we label states using \((N, S)\), the symmetry group is \(\mathrm{U(1)\times SU(2)}\). SU(2) is a non-Abelian symmetry, which is slightly difficult to work with, but it can be two to four times faster than the

`SZ`

mode. This symmetry mode is activated using`symm_type=SymmetryTypes.SU2`

.If we label states using only \(N\), we can remove the

`alpha`

and`beta`

subscripts in the Hamiltonian, and the symmetry group is only \(\mathrm{U(1)}\). This symmetry mode is activated using`symm_type=SymmetryTypes.SGF`

(general spin fermionic). Note that the fermion statistics is considered in this mode.If we label states using only \(N\), we can also transform the Fermionic Hamiltonian to qubit Hamiltonian using JW transformation and then do DMRG on qubits without considering the fermion statistics. This symmetry mode is activated using

`symm_type=SymmetryTypes.SGB`

(general spin bosonic). This mode has the worst efficiency (most of the time the efficiency of`SGF`

and`SGB`

should be the same), but it may show a clearer connection to other JW-based methods.For 2D Hubbard model with only nearest-neighbor interactions, the particle number is also a pseudo-spin (particle-hole) \(\mathrm{SU(2)}\) symmetry. If we label states using \((N, S_z)\), the symmetry group is \(\mathrm{SU(2)\times U(1)}\). This symmetry mode is activated using

`symm_type=SymmetryTypes.SAnyPHSU2`

.For 2D Hubbard model with only nearest-neighbor interactions, we can use the spin and pseudo-psin symmetries simultaneously. If we label states using \((N, S)\), the symmetry group is \(\mathrm{SU(2)\times SU(2)~/~Z_2} = \mathrm{SO(4)}\). This symmetry mode is activated using

`symm_type=SymmetryTypes.SAnySO4`

.

Next, we will explain the settings of the Hubbard Hamiltonian in each of the modes (1) (2) (3) (4) (5) and (6).

```
[2]:
```

```
import numpy as np
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
```

## The `SZ`

Mode

### Initialization

The SZ mode should be the most convenient for the Hubbard model. We first set a driver using this mode. `scratch`

sets a folder for writing temparary scratch files. `symm_type`

sets the symmetry mode (and whether the complex number should be used. Here we allow only real numbers). `n_threads`

sets the number of threads for shared-memory parallelization.

Next, we use the `initialize_system`

method to set the property of the target wavefunction (namely, which symmetry sector is targeted).

`n_sites`

is the number of sites \(L\) (or length of the lattice). Here each site has four local states: \(|0\rangle, |\alpha\rangle, |\beta\rangle\), and \(|\alpha\beta\rangle\).`n_elec`

is the number of electrons \(N\) in the target wavefunction. For half-filling Hubbard model, we should have`n_elec == n_sites`

.`spin`

is two times the total spin \(2S\) (in`SU2`

mode) or two times the projected spin \(2S_z\) (in`SZ`

mode) quantum number of the target wavefunction. Since we want the final wavefunction to have a equal number of alpha and beta electrons (namely, 4 alpha and 4 beta electrons for \(N=8\)), this number is set to zero here. Note that \(2S_z = N_{\alpha} - N_{\beta}\).

```
[3]:
```

```
L = 8
N = 8
TWOSZ = 0
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=N, spin=TWOSZ)
```

### Build Hamiltonian

Now we set up the Hamiltonian. We first write the nearest-neighbor Hubbard model as the following:

Where we have considered the open boundary condition, and \(\sigma \in \{ \alpha, \beta\}\). In the code we use characters `c`

, `d`

, `C`

, `D`

to represent \(a^\dagger_\alpha, a_\alpha, a^\dagger_\beta, a_\beta\), respectively.

The method `expr_builder`

will return a `ExprBuilder`

object.

The method `add_term`

of the `ExprBuilder`

object has three parameters.

The first parameter is a string of characters in

`cdCD`

.The second parameter is a list of integers indicating the site indices associated with each charactor (operator) in the first parameter

The third parameter is a floating point number (or a list of floating point numbers), indicating the coefficient of the added term.

For example, `b.add_term("dCc", [1, 3, 0], 0.7)`

will add the term: \(0.7 \ a_{1\alpha} a^\dagger_{3\beta} a^\dagger_{0\alpha}\).

Most of the time, there can be multiple terms differ only in the operator indices and coefficients. For this case, One can invoke `add_term`

only once. So `b.add_term("CD", [1, 3, 3, 1, 2, 2, 2, 4], [0.7, 0.6, 0.5, 0.4])`

is equivalent to

```
b.add_term("CD", [1, 3], 0.7)
b.add_term("CD", [3, 1], 0.6)
b.add_term("CD", [2, 2], 0.5)
b.add_term("CD", [2, 4], 0.4)
```

`b.finalize`

will do some necessary operator reordering and simplification and `driver.get_mpo`

will create an optimal MPO for the given Hamiltonian. So we can set the MPO for the above Hubbard Hamiltonian using the following:

```
[4]:
```

```
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("cd", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(), -t)
b.add_term("CD", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(), -t)
# onsite term
b.add_term("cdCD", np.array([[i, ] * 4 for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
```

```
Build MPO | Nsites = 8 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 2 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 3 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 4 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 5 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 6 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 6 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.019 Tmvc-total = 0.000 MPO bond dimension = 6 MaxDW = 0.00e+00
NNZ = 78 SIZE = 228 SPT = 0.6579
Rank = 0 Ttotal = 0.047 MPO method = FastBipartite bond dimension = 6 NNZ = 78 SIZE = 228 SPT = 0.6579
```

Note that the above method should allow the user to define arbitrary second-quantized fermionic (particle-number conserving) Hamiltonians (with possibly long-range and high-order terms).

### Run DMRG

Finally, we can run DMRG to find the ground state energy.

We first use `driver.get_random_mps`

to create an initial guess for the wavefunction (MPS).

`tag`

is part of the scratch file name for this MPS. If there are multiple MPS objects, their tags must be all different.`bond_dim`

is the bond dimension of the initial guess MPS. The final optimized MPS may have a larger bond dimension.`nroots`

is the number of roots. If`nroots > 1`

, it will compute the ground state and`nroots - 1`

low-energy excited states.

The DMRG algorithm consists of many “sweeps”. The bond dimension of the MPS can increase after several sweeps. After a few tens of the sweeps, the MPS gradually converges to the ground state if `nroots == 1`

. To prevent stuck in local minima, we may add noise in most sweeps, but the last sweep will not have any noise (namely, `noises[-1] == 0`

) so that the accuracy of the final printed energy is not affected. In each sweep, the DMRG algorithm will optimize \(L\) tensors in the MPS, one
by one. For optimizing each tensor, an effective Hamiltonian will be created and the Davidson algorithm is used to sovle the eigenvalue problem of this effective Hamiltonian. The eigenvalue (energy) of this effective Hamiltonian is the same as the eigenvalue of the full many-body Hamiltonian (if the bond dimension of the MPS is big enough).

Therefore, to run a DMRG calculation, we need to set up a “sweep schedule”, namely, we need to set for each sweep: (a) the MPS bond dimension used in this sweep; (b) the noise used in this sweep; and (c) the Davidson algorithm convergence threshold for this sweep. Typically, these parameters should be changed every 4 to 5 sweeps (for Hubbard model with several hundred sites they may be changed every 30 to 50 sweeps), and the MPS bond dimension increases, while the noise and Davidson algorithm convergence decrease. The Davidson algorithm convergence cannot be set to zero.

`bond_dims`

, `noises`

, and `thrds`

are three lists. The first number in each list is the parameter used in the first sweep, the second number is for the second sweep, etc.

We call `driver.dmrg`

to execute the DMRG algorithm.

`mpo`

is the Hamiltonian represented as an MPO.`ket`

is the initial guess of the MPS. After DMRG finishes, the content of the MPS will be changed and it will be the optimized state.`n_sweeps`

is the maximal number of sweeps.`bond_dims`

is a list of integers indicating the MPS bond dimension for each sweep.`noises`

is a list of real numbers indicating the noise for each sweep.`thrds`

is a list of real numbers indicating the Davidson algorithm convergence threshold for each sweep.`cutoff`

is a real number. Singular values below this number will be discarded during the optimization. Setting this to zero will not discard any small singular value in the MPS unless the MPS bond dimension becomes higher than the set value.`iprint`

can be 0 (quiet), 1 (print energy after each sweep), or 2 (print energy after each iteration in each sweep).

```
[5]:
```

```
def run_dmrg(driver, mpo):
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
return driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, cutoff=0, iprint=1)
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.288 | E = -6.2256341447 | DW = 5.89e-16
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.403 | E = -6.2256341447 | DE = 2.22e-14 | DW = 7.51e-16
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.519 | E = -6.2256341447 | DE = -7.99e-15 | DW = 1.19e-16
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.636 | E = -6.2256341447 | DE = 8.88e-16 | DW = 1.18e-16
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.868 | E = -6.2256341447 | DE = 5.33e-15 | DW = 9.04e-29
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.130 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.10e-28
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.382 | E = -6.2256341447 | DE = -8.88e-16 | DW = 1.27e-28
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.619 | E = -6.2256341447 | DE = 8.88e-16 | DW = 1.32e-28
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.844 | E = -6.2256341447 | DE = 2.66e-15 | DW = 2.76e-51
DMRG energy = -6.225634144662398
```

Note that in the outputs, `DE`

represents the energy difference between two sweeps. `DW`

is the maximal discarded weight (sum of discarded eigenvalues of the density matrix) during the current sweep.

The remaining part of this docuementation introduces some other alternative symmetry to solve the same problem.

## The `SU2`

Mode

### Initialization

In this section, we try to solve the same problem using the SU(2) symmetry, which has a lower computational cost. Now `spin`

represents the total spin (instead of the projected spin).

```
[6]:
```

```
TWOS = 0
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=N, spin=TWOS)
```

### Build Hamiltonian

We need to first do some rewriting of the Hamiltonian to use it with the SU(2) symmetry. This is mainly because operators like \(a_\alpha^\dagger\) do not conserve the total spin quantum number. We have to write the Hamiltonian in terms of only the total spin \(S\) preserving elementary operators.

First, note that

We can rewrite the Hubbard Hamiltonian as

Then we introduce the spin tensor operators which are spin-adaptation of elementary operators for alpha and beta spins. The superscript or subscript \([S]\) indicates the total spin quantum number of each spin tensor operator. Define

Because the total spin is a non-Abelian symmetry, we have the following notes:

The tensor product (coupling) of two spins may generate more than one possible resulting spins. For example, the tensor product between \(S=1/2\) (doublet) and \(S=1/2\) (doublet) can generate \(S=0\) (singlet) and \(S=1\) (triplet). Therefore, when we combine two spin operators \(\big(a_p^\dagger\big)^{[1/2]}\) and \(\big(a_q\big)^{[1/2]}\), we have to indicate which resulting total spin is desired. So we need the following notation to diffentiate the two cases:

There is no associative law. Namely, in general, \((S_1\otimes S_2) \otimes S_3 \neq S_1\otimes (S_2 \otimes S_3)\). The two products are connected by the spin recoupling, and the coefficients between the two products can be computed using Clebsch-Gordan factors. So when writing a string of spin tensor operators, we have to use parenthesis to indicate the coupling order of the operators.

Using spin tensor operators, we can rewrite the Hubbard Hamiltonian as (some derivation details can be found in https://block2.readthedocs.io/en/latest/theory/su2.html)

Next, we have to introduce a string notation to represent the above formula. We can use `C`

and `D`

for \(a^\dagger\) and \(a\), repsectively. Then we use `+`

to connect operators (namely, `+`

represents the tensor product). We use `()`

to indicate the coupling (associative) order. We add a number (2 times the total spin) after each `)`

to indicate the resulting spin number of the tensor product. So, for example, `(C+D)0`

represents the singlet product between
\(a^\dagger\) and \(a\), and `(C+D)2`

represents the triplet product.

As a result, we can set the MPO for the above SU(2) Hubbard Hamiltonian using the following:

```
[7]:
```

```
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("(C+D)0", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(),
np.sqrt(2) * -t)
# onsite term
b.add_term("((C+(C+D)0)1+D)0", np.array([[i, ] * 4 for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
```

```
Build MPO | Nsites = 8 | Nterms = 30 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Site = 1 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 2 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.001 T = 0.004
Site = 4 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 5 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 6 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.020 Tmvc-total = 0.001 MPO bond dimension = 4 MaxDW = 0.00e+00
NNZ = 50 SIZE = 104 SPT = 0.5192
Rank = 0 Ttotal = 0.045 MPO method = FastBipartite bond dimension = 4 NNZ = 50 SIZE = 104 SPT = 0.5192
```

### Run DMRG

Finally, we can run DMRG as normal:

```
[8]:
```

```
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.122 | E = -6.2256341447 | DW = 1.62e-31
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.193 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.34e-31
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.259 | E = -6.2256341447 | DE = 8.88e-16 | DW = 2.12e-32
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.325 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.46e-31
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.396 | E = -6.2256341447 | DE = -7.99e-15 | DW = 0.00e+00
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.466 | E = -6.2256341447 | DE = 1.78e-15 | DW = 0.00e+00
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.587 | E = -6.2256341447 | DE = 4.44e-15 | DW = 0.00e+00
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.703 | E = -6.2256341447 | DE = -5.33e-15 | DW = 0.00e+00
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.796 | E = -6.2256341447 | DE = 1.78e-15 | DW = 0.00e+00
DMRG energy = -6.225634144676156
```

## The `SGF`

Mode

In this section, we try to solve the same Hubbard problem using only particle number symmetry. The number of sites will now be the number of spin orbitals (which is equal to the number of qubits used in the next section), which is two times the number of the original (spatial) sites.

```
[9]:
```

```
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=L * 2, n_elec=N, heis_twos=1)
```

### Build Hamiltonian

Now we first split every Hubbard site into the alpha (odd) and beta (even) sites (next to each other). We have

We use characters `C`

and `D`

to represent \(a^\dagger\) and \(a\), respectively.

So we can set the MPO for the above Hubbard Hamiltonian using the following:

```
[10]:
```

```
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("CD", np.array([[[i, i + 2], [i + 2, i]] for i in range(2 * (L - 1))]).flatten(), -t)
# onsite term
b.add_term("CDCD", np.array([[2 * i, 2 * i, 2 * i + 1, 2 * i + 1] for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
```

```
Build MPO | Nsites = 16 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.004
Site = 2 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 3 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.004
Site = 4 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 5 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 6 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 7 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.006
Site = 8 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 9 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 10 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 11 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 12 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 13 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 14 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.002
Site = 15 / 16 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.044 Tmvc-total = 0.000 MPO bond dimension = 7 MaxDW = 0.00e+00
NNZ = 128 SIZE = 560 SPT = 0.7714
Rank = 0 Ttotal = 0.095 MPO method = FastBipartite bond dimension = 7 NNZ = 128 SIZE = 560 SPT = 0.7714
```

### Run DMRG

Finally, we can run DMRG as normal:

```
[11]:
```

```
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.921 | E = -6.2256341447 | DW = 2.51e-15
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.099 | E = -6.2256341447 | DE = -8.88e-15 | DW = 2.67e-15
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.272 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.33e-15
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.444 | E = -6.2256341447 | DE = 1.07e-14 | DW = 2.38e-15
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.853 | E = -6.2256341447 | DE = -8.88e-16 | DW = 9.93e-33
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.250 | E = -6.2256341447 | DE = -1.78e-15 | DW = 4.03e-32
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.697 | E = -6.2256341447 | DE = -8.88e-16 | DW = 1.85e-32
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.094 | E = -6.2256341447 | DE = -1.78e-15 | DW = 5.42e-32
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.330 | E = -6.2256341447 | DE = 3.55e-15 | DW = 0.00e+00
DMRG energy = -6.225634144669716
```

## The `SGB`

Mode

### Initialization

In this section, we try to solve the same Hubbard problem by first transfroming the model into a qubit (spin) model. The number of sites will now be the number of qubits (spins), which is two times the number of the original sites.

The parameter `heis_twos`

represents two times the spin in each site. Namely, for \(S = 1/2\) Heisenberg model, `heis_twos = 1`

, for \(S = 1\) Heisenberg model, `heis_twos = 2`

, etc.

```
[12]:
```

```
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGB, n_threads=4)
driver.initialize_system(n_sites=L * 2, n_elec=N, heis_twos=1)
```

### Build Hamiltonian

Now we first need to do the Jordan-wigner transform. We split every Hubbard site into the alpha (odd) and beta (even) sites (next to each other). We have

We use charactors `P`

(plus), `M`

(minus), and `Z`

to represent operators \(S^+\), \(S^-\) and \(\frac{1}{2}S^z\), respectively. Note that we assume \(S^z|\alpha\rangle = 1\cdot |\alpha\rangle\) and \(S^z|\beta\rangle = -1\cdot |\beta\rangle\) here.

So we can set the MPO for the above Hubbard Hamiltonian using the following:

```
[13]:
```

```
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("PZZM",
np.array([[i, i, i + 1, i + 2] for i in range(2 * (L - 1))]).flatten(),
[-t * 4] * 2 * (L - 1))
b.add_term("MZZP",
np.array([[i, i, i + 1, i + 2] for i in range(2 * (L - 1))]).flatten(),
[+t * 4] * 2 * (L - 1))
# onsite term
b.add_term("PMPM",
np.array([[i, i] for i in range(2 * L)]).flatten(),
[U] * 2 * L)
mpo = driver.get_mpo(b.finalize(), iprint=2)
```

```
Build MPO | Nsites = 16 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Site = 1 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.003
Site = 2 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 3 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 4 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 5 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 6 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 7 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 8 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.004
Site = 9 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 10 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 11 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 12 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 13 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 14 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.002
Site = 15 / 16 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.040 Tmvc-total = 0.000 MPO bond dimension = 7 MaxDW = 0.00e+00
NNZ = 128 SIZE = 560 SPT = 0.7714
Rank = 0 Ttotal = 0.087 MPO method = FastBipartite bond dimension = 7 NNZ = 128 SIZE = 560 SPT = 0.7714
```

### Run DMRG

Finally, we can run DMRG as normal:

```
[14]:
```

```
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.925 | E = -6.2256341447 | DW = 4.38e-15
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.078 | E = -6.2256341447 | DE = 2.13e-14 | DW = 4.26e-15
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.236 | E = -6.2256341447 | DE = 0.00e+00 | DW = 4.25e-15
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.432 | E = -6.2256341447 | DE = 5.33e-15 | DW = 3.95e-15
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.785 | E = -6.2256341447 | DE = 0.00e+00 | DW = 2.54e-32
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.172 | E = -6.2256341447 | DE = 4.44e-15 | DW = 7.02e-32
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.570 | E = -6.2256341447 | DE = 0.00e+00 | DW = 1.86e-32
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.929 | E = -6.2256341447 | DE = 1.78e-15 | DW = 5.74e-32
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.157 | E = -6.2256341447 | DE = -3.55e-15 | DW = 0.00e+00
DMRG energy = -6.225634144672112
```

## The `PHSU2`

Mode

Let

We have (assuming \(p\) and \(q\) differ by 1)

So

Now consider the onsite term

So the product of above terms

```
[15]:
```

```
t = 1
U = 2
n_elec = L
driver = DMRGDriver(scratch='./tmp', symm_type=SymmetryTypes.SAnyPHSU2, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=n_elec, spin=0)
b = driver.expr_builder()
b.add_term("(E+F)0", [g for i in range(L - 1) for g in [i, i + 1]], 2 ** 0.5)
b.add_term("(F+E)0", [g for i in range(L - 1) for g in [i, i + 1]], 2 ** 0.5)
b.add_term("((E+E)0+(F+F)0)0", np.array([i for i in range(L) for _ in range(4)]), 0.5 * U)
b.add_const(n_elec * U / 2)
mpo = driver.get_mpo(b.finalize(adjust_order=False), iprint=2)
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
energies = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energies)
```

```
Build MPO | Nsites = 8 | Nterms = 22 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Site = 1 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 2 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.005
Site = 4 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.005
Site = 5 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 6 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.028 Tmvc-total = 0.000 MPO bond dimension = 4 MaxDW = 0.00e+00
NNZ = 50 SIZE = 104 SPT = 0.5192
Rank = 0 Ttotal = 0.069 MPO method = FastBipartite bond dimension = 4 NNZ = 50 SIZE = 104 SPT = 0.5192
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.144 | E = -6.2256341447 | DW = 2.13e-20
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.187 | E = -6.2256341447 | DE = 4.26e-14 | DW = 2.06e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.232 | E = -6.2256341447 | DE = -8.88e-15 | DW = 1.49e-20
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.274 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.86e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.317 | E = -6.2256341447 | DE = -3.55e-15 | DW = 1.51e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.359 | E = -6.2256341447 | DE = 0.00e+00 | DW = 1.15e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.402 | E = -6.2256341447 | DE = 1.78e-15 | DW = 1.96e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.444 | E = -6.2256341447 | DE = 7.11e-15 | DW = 1.37e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.478 | E = -6.2256341447 | DE = -2.84e-14 | DW = 1.94e-20
DMRG energy = -6.225634144668852
```

## The `SO4`

Mode

Similarly, we can consider the \(\mathrm{SU(2) \times SU(2)}\) operator \(G\), where the row represents the pseudo spin symmetry and the column represents the spin symmetry.

So the components can be represented as

```
[16]:
```

```
t = 1
U = 2
n_elec = L
driver = DMRGDriver(scratch='./tmp', symm_type=SymmetryTypes.SAnySO4, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=n_elec, spin=0)
b = driver.expr_builder()
b.add_term("(G[1,1]+G[1,1])[0,0]",
[g for i in range(L - 1) for g in [i, i + 1]], [2.0 if i % 2 == 0 else -2.0 for i in range(L - 1)])
b.add_term("((G[1,1]+G[1,1])[0,2]+(G[1,1]+G[1,1])[0,2])[0,0]",
np.array([i for i in range(L) for _ in range(4)]), U / (2 * 3 ** 0.5))
b.add_const(n_elec * U / 2)
mpo = driver.get_mpo(b.finalize(adjust_order=False), iprint=2)
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
energies = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energies)
```

```
Build MPO | Nsites = 8 | Nterms = 15 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 3 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 2 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 4 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 5 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 6 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 3 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.022 Tmvc-total = 0.000 MPO bond dimension = 3 MaxDW = 0.00e+00
NNZ = 36 SIZE = 60 SPT = 0.4000
Rank = 0 Ttotal = 0.047 MPO method = FastBipartite bond dimension = 3 NNZ = 36 SIZE = 60 SPT = 0.4000
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.066 | E = -6.2256341447 | DW = 1.22e-20
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.113 | E = -6.2256341447 | DE = 3.02e-14 | DW = 2.80e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.159 | E = -6.2256341447 | DE = -3.55e-15 | DW = 7.31e-21
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.192 | E = -6.2256341447 | DE = 8.88e-15 | DW = 9.99e-21
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.229 | E = -6.2256341447 | DE = -1.95e-14 | DW = 1.22e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.261 | E = -6.2256341447 | DE = -1.78e-15 | DW = 5.33e-21
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.295 | E = -6.2256341447 | DE = 3.55e-15 | DW = 8.67e-21
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.333 | E = -6.2256341447 | DE = -1.78e-15 | DW = 7.17e-21
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.367 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.25e-20
DMRG energy = -6.225634144673016
```