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:

\[\hat{H} = -t \sum_{i=1,\sigma}^{L-1} \left( a_{i\sigma}^\dagger a_{i+1\sigma} + a_{i+1\sigma}^\dagger a_{i\sigma}\right) + U \sum_{i=1}^L a^\dagger_{i\alpha}a_{i\alpha} a^\dagger_{i\beta}a_{i\beta}\]

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

\[n_{i\uparrow} n_{i\downarrow} = a_{i\alpha}^\dagger a_{i\alpha} a_{i\beta}^\dagger a_{i\beta} = a_{i\alpha}^\dagger a_{i\beta}^\dagger a_{i\beta} a_{i\alpha} = a_{i\beta}^\dagger a_{i\alpha}^\dagger a_{i\alpha} a_{i\beta} = \frac{1}{2} \sum_{\sigma} \sum_{\sigma' \neq \sigma} a_{i\sigma}^\dagger a_{i\sigma'}^\dagger a_{i\sigma'} a_{i\sigma} = \frac{1}{2} \sum_{\sigma\sigma'} a_{i\sigma}^\dagger a_{i\sigma'}^\dagger a_{i\sigma'} a_{i\sigma}\]

We can rewrite the Hubbard Hamiltonian as

\[\hat{H} = -t\sum_{\langle i j \rangle,\sigma} a^\dagger_{i\sigma} a_{j\sigma} + \frac{U}{2} \sum_{i,\sigma\sigma'} a_{i\sigma}^\dagger a_{i\sigma'}^\dagger a_{i\sigma'} a_{i\sigma}\]

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

\[\begin{split}\big(a_p^\dagger\big)^{[1/2]} := \begin{pmatrix} a_{p\alpha}^\dagger \\ a_{p\beta}^\dagger \end{pmatrix}^{[1/2]}\quad \big(a_p\big)^{[1/2]} := \begin{pmatrix} -a_{p\beta} \\ a_{p\alpha} \end{pmatrix}^{[1/2]}\end{split}\]

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:

\[\begin{split}\big(a_p^\dagger\big)^{[1/2]} \otimes_{[0]} \big(a_q\big)^{[1/2]} \quad \text{(singlet product)} \\ \big(a_p^\dagger\big)^{[1/2]} \otimes_{[1]} \big(a_q\big)^{[1/2]} \quad \text{(triplet product)}\end{split}\]

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)

\[\hat{H}^{[0]} = -\sqrt{2} \ t\sum_{\langle i j \rangle} \big(a_i^\dagger\big)^{[1/2]} \otimes_{[0]} \big(a_j\big)^{[1/2]} + U \sum_i \bigg( \big(a_i^\dagger\big)^{[1/2]} \otimes_{[1/2]} \Big( \big(a_i^\dagger\big)^{[1/2]} \otimes_{[0]} \big(a_i\big)^{[1/2]} \Big) \bigg) \otimes_{[0]} \big(a_i\big)^{[1/2]}\]

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

\[\hat{H} = -t \sum_{i=1}^{L-1} \left( a_{2i}^\dagger a_{2i+2} + a_{2i+2}^\dagger a_{2i} + a_{2i-1}^\dagger a_{2i+1} + a_{2i+1}^\dagger a_{2i-1}\right) + U \sum_{i=1}^L a^\dagger_{2i-1}a_{2i-1} a^\dagger_{2i}a_{2i}\]

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

\[\hat{H} = -t \sum_{i=1}^{L-1} \Big( S^+_{2i} S^z_{2i} S^z_{2i+1} S^-_{2i+2} - S^-_{2i} S^z_{2i} S^z_{2i+1} S^+_{2i+2} + S^+_{2i-1} S^z_{2i-1} S^z_{2i} S^-_{2i+1} - S^-_{2i-1} S^z_{2i-1} S^z_{2i} S^+_{2i+1} \Big) + U \sum_{i=1}^L S^+_{2i-1} S^-_{2i-1} S^+_{2i} S^-_{2i}\]

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

\[\begin{split}\big( E_p \big)^{[1/2]} := \begin{pmatrix} a^\dagger_{p,\uparrow} \\ (-1)^p a_{p,\downarrow} \end{pmatrix}^{[1/2]} \quad \big( F_p \big)^{[1/2]} := \begin{pmatrix} -(-1)^p a^\dagger_{p,\downarrow} \\ a_{p,\uparrow} \end{pmatrix}^{[1/2]}\end{split}\]

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

\[\begin{split}\big( E_p \big)^{[1/2]} \otimes_{[0]} \big( F_q \big)^{[1/2]} = \begin{pmatrix} a^\dagger_{p,\uparrow} \\ (-1)^p a_{p,\downarrow} \end{pmatrix}^{[1/2]} \otimes_{[0]} \begin{pmatrix} -(-1)^q a^\dagger_{q,\downarrow} \\ a_{q,\uparrow} \end{pmatrix}^{[1/2]} \\ = \frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{q,\uparrow} - a_{p,\downarrow} a^\dagger_{q,\downarrow} \right) = \frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{q,\uparrow} + a^\dagger_{q,\downarrow} a_{p,\downarrow} \right) \\\end{split}\]

So

\[\begin{split}\big( E_p \big)^{[1/2]} \otimes_{[0]} \big( F_{p+1} \big)^{[1/2]} + \big( E_{p+1} \big)^{[1/2]} \otimes_{[0]} \big( F_p \big)^{[1/2]} \\ =\frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{p+1,\uparrow} + a^\dagger_{p+1,\downarrow} a_{p,\downarrow} \right) + \frac{1}{\sqrt{2}} \left( a^\dagger_{p+1,\uparrow} a_{p,\uparrow} + a^\dagger_{p,\downarrow} a_{p+1,\downarrow} \right) \\ =\frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{p+1,\uparrow} + a^\dagger_{p,\downarrow} a_{p+1,\downarrow} + a^\dagger_{p+1,\uparrow} a_{p,\uparrow} + a^\dagger_{p+1,\downarrow} a_{p,\downarrow} \right) \\ =\frac{1}{\sqrt{2}} \sum_\sigma \left( a^\dagger_{p,\sigma} a_{p+1,\sigma} + \mathrm{H.c.} \right)\end{split}\]

\[-t \sum_{\langle i,j \rangle, \sigma} \left( c_{i,\sigma}^\dagger c_{j,\sigma} + \mathrm{H.c.} \right) = -\sqrt{2} t \sum_{\langle i,j \rangle} \left[ \big( E_i \big)^{[1/2]} \otimes_{[0]} \big( F_{j} \big)^{[1/2]} + \big( E_{j} \big)^{[1/2]} \otimes_{[0]} \big( F_i \big)^{[1/2]} \right]\]

Now consider the onsite term

\[\begin{split}\big( E_p \big)^{[1/2]} \otimes_{[0]} \big( F_p \big)^{[1/2]} = \begin{pmatrix} a^\dagger_{p,\uparrow} \\ (-1)^p a_{p,\downarrow} \end{pmatrix}^{[1/2]} \otimes_{[0]} \begin{pmatrix} -(-1)^p a^\dagger_{p,\downarrow} \\ a_{p,\uparrow} \end{pmatrix}^{[1/2]} \\ = \frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{p,\uparrow} + a_{p,\downarrow} a^\dagger_{p,\downarrow} \right) = \frac{1}{\sqrt{2}} \left( a^\dagger_{p,\uparrow} a_{p,\uparrow} - a^\dagger_{p,\downarrow} a_{p,\downarrow} + 1\right) \\ \big( E_p \big)^{[1/2]} \otimes_{[0]} \big( E_p \big)^{[1/2]} = \begin{pmatrix} a^\dagger_{p,\uparrow} \\ (-1)^p a_{p,\downarrow} \end{pmatrix}^{[1/2]} \otimes_{[0]} \begin{pmatrix} a^\dagger_{p,\uparrow} \\ (-1)^p a_{p,\downarrow} \end{pmatrix}^{[1/2]} \\ = \frac{1}{\sqrt{2}} \left( (-1)^p a^\dagger_{p,\uparrow} a_{p,\downarrow} -(-1)^p a_{p,\downarrow} a^\dagger_{p,\uparrow} \right) = (-1)^p \sqrt{2} a^\dagger_{p,\uparrow} a_{p,\downarrow} \\ \big( F_p \big)^{[1/2]} \otimes_{[0]} \big( F_p \big)^{[1/2]} = \begin{pmatrix} -(-1)^p a^\dagger_{p,\downarrow} \\ a_{p,\uparrow} \end{pmatrix}^{[1/2]} \otimes_{[0]} \begin{pmatrix} -(-1)^p a^\dagger_{p,\downarrow} \\ a_{p,\uparrow} \end{pmatrix}^{[1/2]} \\ = \frac{1}{\sqrt{2}} \left( -(-1)^p a^\dagger_{p,\downarrow} a_{p,\uparrow} +(-1)^p a_{p,\uparrow} \right) = -(-1)^p \sqrt{2} a^\dagger_{p,\downarrow} a_{p,\uparrow}\end{split}\]

So the product of above terms

\[\begin{split}\left( \big( E_p \big)^{[1/2]} \otimes_{[0]} \big( E_p \big)^{[1/2]} \right) \otimes_{[0]} \left( \big( F_p \big)^{[1/2]} \otimes_{[0]} \big( F_p \big)^{[1/2]} \right)\\ = -2 a^\dagger_{p,\uparrow} a_{p,\downarrow} a^\dagger_{p,\downarrow} a_{p,\uparrow} = -2 a^\dagger_{p,\uparrow} a_{p,\uparrow} a_{p,\downarrow} a^\dagger_{p,\downarrow} = 2 a^\dagger_{p,\uparrow} a_{p,\uparrow} a^\dagger_{p,\downarrow} a_{p,\downarrow} -2 a^\dagger_{p,\uparrow} a_{p,\uparrow} = 2 n_{p\uparrow} n_{p\downarrow} - 2 n_{p\uparrow}\end{split}\]

[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.

\[\begin{split}\big( G_p \big)^{[1/2,1/2]} := \begin{pmatrix} a^\dagger_{p,\uparrow} & a^\dagger_{p,\downarrow} \\ (-1)^p a_{p,\downarrow} & -(-1)^p a_{p,\uparrow} \end{pmatrix}^{[1/2,1/2]}\end{split}\]

So the components can be represented as

\[\begin{split}\begin{pmatrix} [1/2,1/2] & [1/2,-1/2] \\ [-1/2,1/2] & [-1/2,-1/2] \end{pmatrix}^{[1/2,1/2]}\end{split}\]

[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

Hubbard Model

Introduction

The SZ Mode

Initialization

Build Hamiltonian

Run DMRG

The SU2 Mode

Initialization

Build Hamiltonian

Run DMRG

The SGF Mode

Build Hamiltonian

Run DMRG

The SGB Mode

Initialization

Build Hamiltonian

Run DMRG

The PHSU2 Mode

The SO4 Mode

The `SZ` Mode

The `SU2` Mode

The `SGF` Mode

The `SGB` Mode

The `PHSU2` Mode

The `SO4` Mode