Custom Hamiltonians

[ ]:

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

In this tutorial, we provide an example python scripts for performing DMRG using custom Hamiltonians, where the operators and states at local Hilbert space at every site can be redefined. It is also possible to use different local Hilbert space for different sites. New letters can be introduced for representing new operators (the operator name can only be a single lower or upper case character).

Note the following examples are only supposed to work in the Abelian symmetry modes (SZ, SZ|CPX, SGF, SGFCPX, SAny, or SAny|CPX) for Abelian symmetries, and non-Abelian symmetry modes (SAnySU2 or SAnySU2|CPX) for non-Abelian symmetries, respectively.

The Hubbard Model

In the following example, we implement a custom Hamiltonian for the Hubbard model. In the standard implementation, the on-site term was represented as cdCD. Here we instead introduce a single letter N for the cdCD term. For each letter in cdCDN (representing elementary operators), we define its matrix representation in the local basis in site_ops. The quantum number and number of states in each quantum number at each site (which defines the local Hilbert space) are set in site_basis.

[2]:

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

L = 8
U = 2
N_ELEC = 8

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

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []
Q = driver.bw.SX # quantum number wrapper (n_elec, 2 * spin, point group irrep)

for k in range(L):
    basis = [(Q(0, 0, 0), 1), (Q(1, 1, 0), 1), (Q(1, -1, 0), 1), (Q(2, 0, 0), 1)] # [0ab2]
    ops = {
        "": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]),   # identity
        "c": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]),  # alpha+
        "d": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]),  # alpha
        "C": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+
        "D": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta
        "N": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]]),  # cdCD
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

b.add_term("cd", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("CD", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("N", np.array([i for i in range(L)]), U)

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy = %20.15f" % energy)


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.167 | E =      -6.2256341447 | DW = 2.65116e-16

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.240 | E =      -6.2256341447 | DE = -8.88e-15 | DW = 4.93007e-16

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.305 | E =      -6.2256341447 | DE = 1.78e-15 | DW = 9.66767e-17

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.367 | E =      -6.2256341447 | DE = 1.78e-15 | DW = 1.20251e-16

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.431 | E =      -6.2256341447 | DE = -8.88e-16 | DW = 3.71514e-20

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.493 | E =      -6.2256341447 | DE = -3.55e-15 | DW = 6.52472e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.560 | E =      -6.2256341447 | DE = -1.78e-15 | DW = 4.11603e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.626 | E =      -6.2256341447 | DE = 0.00e+00 | DW = 4.46198e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      0.674 | E =      -6.2256341447 | DE = 7.11e-15 | DW = 4.96657e-20

DMRG energy =   -6.225634144657917

The Hubbard-Holstein Model

The above script can be easily extended to treat phonons.

[3]:

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

N_SITES_ELEC, N_SITES_PH, N_ELEC = 4, 4, 4
N_PH, U, OMEGA, G = 11, 2, 0.25, 0.5
L = N_SITES_ELEC + N_SITES_PH

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

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []
Q = driver.bw.SX # quantum number wrapper (n_elec, 2 * spin, point group irrep)

for k in range(L):
    if k < N_SITES_ELEC:
        # electron part
        basis = [(Q(0, 0, 0), 1), (Q(1, 1, 0), 1), (Q(1, -1, 0), 1), (Q(2, 0, 0), 1)] # [0ab2]
        ops = {
            "": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]),   # identity
            "c": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]),  # alpha+
            "d": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]),  # alpha
            "C": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+
            "D": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta
        }
    else:
        # phonon part
        basis = [(Q(0, 0, 0), N_PH)]
        ops = {
            "": np.identity(N_PH), # identity
            "E": np.diag(np.sqrt(np.arange(1, N_PH)), k=-1), # ph+
            "F": np.diag(np.sqrt(np.arange(1, N_PH)), k=1),  # ph
        }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms in Hubbard-Holstein model
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

# electron part
b.add_term("cd", np.array([[i, i + 1, i + 1, i] for i in range(N_SITES_ELEC - 1)]).ravel(), -1)
b.add_term("CD", np.array([[i, i + 1, i + 1, i] for i in range(N_SITES_ELEC - 1)]).ravel(), -1)
b.add_term("cdCD", np.array([[i, i, i, i] for i in range(N_SITES_ELEC)]).ravel(), U)

# phonon part
b.add_term("EF", np.array([[i + N_SITES_ELEC, ] * 2 for i in range(N_SITES_PH)]).ravel(), OMEGA)

# interaction part
b.add_term("cdE", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)
b.add_term("cdF", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)
b.add_term("CDE", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)
b.add_term("CDF", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy = %20.15f" % energy)


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     94.146 | E =      -6.9568929542 | DW = 3.61640e-09

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =    126.358 | E =      -6.9568932112 | DE = -2.57e-07 | DW = 3.01342e-19

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =    139.180 | E =      -6.9568932112 | DE = 4.44e-15 | DW = 1.23451e-19

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =    149.110 | E =      -6.9568932112 | DE = 1.78e-15 | DW = 7.41629e-20

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    158.521 | E =      -6.9568932112 | DE = 3.55e-15 | DW = 7.71387e-20

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    166.802 | E =      -6.9568932112 | DE = -5.33e-15 | DW = 6.44981e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    175.980 | E =      -6.9568932112 | DE = 8.88e-16 | DW = 8.11423e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    184.501 | E =      -6.9568932112 | DE = -2.66e-15 | DW = 6.41555e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =    192.172 | E =      -6.9568932112 | DE = 0.00e+00 | DW = 7.78301e-20

DMRG energy =   -6.956893211180857

Custom Symmetry Groups

In the following we show how to set the custom symmetry groups. The symmetry mode SymmetryTypes.SAny (or SymmetryTypes.SAny | SymmetryTypes.CPX if complex number is required) should be used for this purpose. Currently we support the definition of symmetry group as an arbitrary direct product of up to six Abelian symmetry sub-groups. Possible sub-group names are “U1”, “Z1”, “Z2”, “Z3”, …, “Z2055”, “U1Fermi”, “Z1Fermi”, “Z2Fermi”, “Z3Fermi”, …, “Z2055Fermi”, “LZ”, and “AbelianPG”. The names with the suffix “Fermi” should be used for Fermion symmetries. The names without the suffix “Fermi” should be used for spin or Boson symmetries. The non-Abelian sub-group names “SU2” and “SU2Fermi” can be used with some restrictions, see the example in later sections. The DMRGDriver.set_symmetry_groups(sub_group_name_1: str, sub_group_name_2: str, ...) method can be used to set the symmetry sub-groups. The number of arguments in the quantum number wrapper Q should then match the number of sub-group names given in DMRGDriver.set_symmetry_groups.

As a first example, we use the custom symmetry group syntax to recompute the Hubbard model. We first use \(U(1) \times U(1)\) symmetry, which should be equivalent to the previous SZ mode example.

[4]:

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

L = 8
U = 2
N_ELEC = 8
TWO_SZ = 0

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAny, n_threads=4)

# quantum number wrapper (U1 / n_elec, U1 / 2*Sz)
driver.set_symmetry_groups("U1Fermi", "U1")
Q = driver.bw.SX

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []

for k in range(L):
    basis = [(Q(0, 0), 1), (Q(1, 1), 1), (Q(1, -1), 1), (Q(2, 0), 1)] # [0ab2]
    ops = {
        "": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]),   # identity
        "c": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]),  # alpha+
        "d": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]),  # alpha
        "C": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+
        "D": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(N_ELEC, TWO_SZ), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

b.add_term("cd", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("CD", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("cdCD", np.array([i for i in range(L) for _ in range(4)]), U)

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy = %20.15f" % energy)


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.247 | E =      -6.2256341447 | DW = 5.28457e-16

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.344 | E =      -6.2256341447 | DE = -8.88e-15 | DW = 6.68285e-16

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.446 | E =      -6.2256341447 | DE = -4.44e-15 | DW = 1.63432e-16

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.543 | E =      -6.2256341447 | DE = 1.78e-15 | DW = 2.16706e-16

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.623 | E =      -6.2256341447 | DE = 3.55e-15 | DW = 2.88426e-20

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.692 | E =      -6.2256341447 | DE = -5.33e-15 | DW = 2.28873e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.764 | E =      -6.2256341447 | DE = 8.88e-16 | DW = 3.56876e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      0.831 | E =      -6.2256341447 | DE = 1.78e-15 | DW = 3.86969e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      0.880 | E =      -6.2256341447 | DE = -2.66e-15 | DW = 5.43709e-20

DMRG energy =   -6.225634144658391

As a second example, we recompute the Hubbard model using \(Z_2 \times Z_2\) symmetry. This time we cannot easily target the \(N_{\mathrm{elec}} = 8\) symmetry sector. Instead, we compute a few excited states and compute the \(\langle N\rangle\) to identify the correct state.

[5]:

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

L = 8
U = 2
N_ELEC = 8
TWO_SZ = 0

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAny, n_threads=4)

# quantum number wrapper (Z2 / n_elec, Z2 / 2*Sz)
driver.set_symmetry_groups("Z2Fermi", "Z2")
Q = driver.bw.SX

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []

for k in range(L):
    basis = [(Q(0, 0), 2), (Q(1, 1), 2)] # [02ab]
    ops = {
        # note the order of row and column is different from the U1xU1 case
        "": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]),   # identity
        "c": np.array([[0, 0, 0, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 0, 0]]),  # alpha+
        "d": np.array([[0, 0, 1, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 1, 0, 0]]),  # alpha
        "C": np.array([[0, 0, 0, 0], [0, 0, -1, 0], [0, 0, 0, 0], [1, 0, 0, 0]]), # beta+
        "D": np.array([[0, 0, 0, 1], [0, 0, 0, 0], [0, -1, 0, 0], [0, 0, 0, 0]]), # beta
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(N_ELEC % 2, TWO_SZ % 2), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

b.add_term("cd", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("CD", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)
b.add_term("cdCD", np.array([i for i in range(L) for _ in range(4)]), U)

# [Part C] Perform state-averaged DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=10)
energies = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=200, iprint=1)

# [Part D] Check particle number expectations
b = driver.expr_builder()
b.add_term("cd", np.array([[i, i] for i in range(L)]).ravel(), 1)
b.add_term("CD", np.array([[i, i] for i in range(L)]).ravel(), 1)
partile_n_mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)

kets = [driver.split_mps(mps, ir, tag="KET-%d" % ir) for ir in range(mps.nroots)]
for ir in range(mps.nroots):
    n_expt = driver.expectation(kets[ir], partile_n_mpo, kets[ir])
    print("Root = %d <E> = %20.15f <N> = %10.3f" % (ir, energies[ir], n_expt))


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     36.478 | E[ 10] =      -6.9958183038     -6.5705972881     -6.5705972845     -6.5705972809     -6.2773135425     -6.2256341359     -6.1122535707     -6.1122534910     -6.1122534358     -6.0320900884 | DW = 4.44265e-09

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     40.236 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -1.62e-07 | DW = 1.12165e-09

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     43.608 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = 3.43e-12 | DW = 1.18398e-09

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     46.245 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -1.75e-12 | DW = 1.12063e-09

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     49.844 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -1.88e-12 | DW = 7.51709e-17

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     52.386 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -8.88e-15 | DW = 1.55303e-17

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     55.928 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = 0.00e+00 | DW = 1.59309e-17

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     58.510 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = 1.24e-14 | DW = 6.42978e-18

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     60.677 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = 0.00e+00 | DW = 9.28475e-19

Root = 0 <E> =   -6.995818305899173 <N> =      6.000
Root = 1 <E> =   -6.570597289543093 <N> =      6.000
Root = 2 <E> =   -6.570597289542542 <N> =      6.000
Root = 3 <E> =   -6.570597289541823 <N> =      6.000
Root = 4 <E> =   -6.277313551397241 <N> =      6.000
Root = 5 <E> =   -6.225634144677118 <N> =      8.000
Root = 6 <E> =   -6.112253573866418 <N> =      6.000
Root = 7 <E> =   -6.112253573864956 <N> =      6.000
Root = 8 <E> =   -6.112253573862747 <N> =      6.000
Root = 9 <E> =   -6.032090249939520 <N> =      4.000

Bose-Hubbard Model

In this example we will find the ground state of 1D Bose-Hubbard model. Unlike the Holstein model, in this model the boson number is conserved as an Abelian U(1) symmetry. The definition of the Hamiltonian can be found in https://en.wikipedia.org/wiki/Bose%E2%80%93Hubbard_model, which is given by

\[\hat{H} = -t \sum_{\langle ij\rangle} \Big( b^\dagger_i b_j + b^\dagger_j b_i \Big) + \frac{U}{2} \sum_i n_i (n_i - 1) - \mu \sum_i n_i\]

In the following script, we consider a 1D chain including \(L= 10\) sites with open boundary condition, \(t = 1, U = 0.1\) and \(\mu = 0\). The total number of boson in the many-body state is set to 10, and the maximal number of boson per site is 5.

[6]:

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

L = 10
T = 1
U = 0.1
MU = 0
NB_MAX = 5 # max n_boson per site
N_BOSON = 10

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAny, n_threads=4)

driver.set_symmetry_groups("U1")
Q = driver.bw.SX

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []

for k in range(L):
    basis = [(Q(i), 1) for i in range(NB_MAX + 1)] # [012..NB_MAX]
    ops = {
        "": np.identity(NB_MAX + 1),                           # identity
        "C": np.diag(np.sqrt(np.arange(1, NB_MAX + 1)), k=-1), # b+
        "D": np.diag(np.sqrt(np.arange(1, NB_MAX + 1)), k=1),  # b
        "N": np.diag(np.arange(0, NB_MAX + 1), k=0),           # particle number
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.initialize_system(n_sites=L, vacuum=Q(0), target=Q(N_BOSON), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

b.add_term("CD", np.array([j for i in range(L - 1) for j in [i, i + 1, i + 1, i]]), -T)
b.add_term("N", np.array([i for i in range(L)]), -(MU + U / 2))
b.add_term("NN", np.array([j for i in range(L) for j in [i, i]]), U / 2)

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=""), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy = %20.15f (per site = %10.6f)" % (energy, energy / L))


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.026 | E =     -18.5867775673 | DW = 4.33791e-12

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      3.090 | E =     -18.5873694219 | DE = -5.92e-04 | DW = 3.77544e-16

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      3.436 | E =     -18.5873694219 | DE = -7.11e-15 | DW = 3.97569e-17

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      3.813 | E =     -18.5873694219 | DE = -8.30e-12 | DW = 1.83165e-17

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.397 | E =     -18.5873694219 | DE = 3.55e-15 | DW = 1.35857e-19

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.703 | E =     -18.5873694219 | DE = -3.55e-15 | DW = 1.43599e-19

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.971 | E =     -18.5873694219 | DE = -1.07e-14 | DW = 8.50577e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      5.226 | E =     -18.5873694219 | DE = 3.55e-15 | DW = 1.44991e-19

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      5.415 | E =     -18.5873694219 | DE = 0.00e+00 | DW = 8.57382e-20

DMRG energy =  -18.587369421896540 (per site =  -1.858737)

SU(3) Heisenberg Model

In this example we will find the ground state of the 1D SU(3) Heisenberg model using the custom symemtry group syntax. We will only use Abelian symmetry groups (for quantum numbers \(S_z\) and \(Q_z\)) for this problem. The model used here can be found in Eq. (2) in Phys. Rev. B 79, 012408 (2009).

[7]:

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

L = 72

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAny, n_threads=4)

# quantum number wrapper (2Sz / X, 2Qz / Y)
driver.set_symmetry_groups("U1", "U1")
Q = driver.bw.SX

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []

# Gell Mann operators
lambda_ops = {
    "L1": np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]),              # lambda_1
    "L2": np.array([[0, -1j, 0], [1j, 0, 0], [0, 0, 0]]),           # lambda_2
    "L3": np.array([[1, 0, 0], [0, -1, 0], [0, 0, 0]]),             # lambda_3
    "L4": np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]),              # lambda_4
    "L5": np.array([[0, 0, -1j], [0, 0, 0], [1j, 0, 0]]),           # lambda_5
    "L6": np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]),              # lambda_6
    "L7": np.array([[0, 0, 0], [0, 0, -1j], [0, 1j, 0]]),           # lambda_7
    "L8": np.array([[1, 0, 0], [0, 1, 0], [0, 0, -2]]) / 3 ** 0.5,  # lambda_8
}

for k in range(L):
    basis = [(Q(1, 1), 1), (Q(-1, 1), 1), (Q(0, -2), 1)]
    ops = {
        "": np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),        # identity
        "T": (lambda_ops["L1"] + 1j * lambda_ops["L2"]).real,   # T+
        "t": (lambda_ops["L1"] - 1j * lambda_ops["L2"]).real,   # T-
        "V": (lambda_ops["L4"] + 1j * lambda_ops["L5"]).real,   # V+
        "v": (lambda_ops["L4"] - 1j * lambda_ops["L5"]).real,   # V-
        "U": (lambda_ops["L6"] + 1j * lambda_ops["L7"]).real,   # U+
        "u": (lambda_ops["L6"] - 1j * lambda_ops["L7"]).real,   # U-
        "L": lambda_ops["L3"],                                  # L3
        "l": lambda_ops["L8"],                                  # L8
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(0, 0), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

b.add_term("Tt", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)
b.add_term("Vv", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)
b.add_term("Uu", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)
b.add_term("LL", np.array([[i, i + 1] for i in range(L - 1)]).ravel(), 0.25)
b.add_term("ll", np.array([[i, i + 1] for i in range(L - 1)]).ravel(), 0.25)
b.iscale(1 / L) # compute energy per site instead of total energy

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=""), algo_type=MPOAlgorithmTypes.FastBipartite)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)

energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy (per site) = %20.15f" % energy)


Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     11.778 | E =      -0.5120663491 | DW = 1.30017e-05

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     18.660 | E =      -0.5145229771 | DE = -2.46e-03 | DW = 4.16418e-08

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     23.883 | E =      -0.5145506979 | DE = -2.77e-05 | DW = 3.51827e-07

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     26.629 | E =      -0.5145508464 | DE = -1.48e-07 | DW = 4.79795e-07

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     34.588 | E =      -0.5145512453 | DE = -3.99e-07 | DW = 4.81635e-09

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     47.894 | E =      -0.5145512513 | DE = -6.07e-09 | DW = 4.36734e-09

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     61.178 | E =      -0.5145512510 | DE = 2.96e-10 | DW = 1.36107e-09

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     74.667 | E =      -0.5145512506 | DE = 4.30e-10 | DW = 8.13159e-10

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     84.889 | E =      -0.5145512503 | DE = 2.84e-10 | DW = 9.28021e-19

DMRG energy (per site) =   -0.514551250318804

SU(2) \(t-J\) Model

As an example of non-Abelian custom symmetry groups, in the following we will find the ground state of the 2D SU(2) \(t-J\) model using the custom symemtry group syntax. We will use the \(\mathrm{U(1) \times SU(2)}\) symmetry group (for particle number \(N\) and total spin \(S\)) for this problem.

The \(t-J\) model Hamiltonian is given by

\[\hat{H} = -t \sum_{\langle ij\rangle, \sigma} \Big( a^\dagger_{i\sigma} a_{j\sigma} + a^\dagger_{j\sigma} a_{i\sigma} \Big) + J \sum_{\langle ij\rangle} \Big( \mathbf{\hat{S}}_i \cdot \mathbf{\hat{S}}_j - \frac{1}{4} n_i n_j \Big)\]

We consider \(4\times 4\) square 2D lattice with open boundary condition, \(t = 1, J = 0.4\), and hole doping \(=1/8\). \(\langle ij\rangle\) only conatins one of \(i=1, j=2\) and \(i=2, j=1\), for example. The reference energy per site can be found in npj Quantum Materials 5, 28 (2020).

Written in SU(2) notation, the Hamiltonian is

\[\begin{split}\hat{H} = -\sqrt{2} t \sum_{\langle ij\rangle} \Big( \big(a^\dagger_i\big)^{[1/2]} \otimes_{[0]} \big(a_j\big)^{[1/2]} + \big(a^\dagger_j\big)^{[1/2]} \otimes_{[0]} \big(a_i\big)^{[1/2]} \Big) \\ + J \sum_{\langle ij\rangle} \Big[ -\frac{\sqrt{3}}{2} \Big( \big(a^\dagger_i\big)^{[1/2]} \otimes_{[1]} \big(a_i\big)^{[1/2]} \Big)\otimes_{[0]} \Big( \big(a^\dagger_j\big)^{[1/2]} \otimes_{[1]} \big(a_j\big)^{[1/2]} \Big) \\ - \frac{(\sqrt{2})^2}{4} \Big( \big(a^\dagger_i\big)^{[1/2]} \otimes_{[0]} \big(a_i\big)^{[1/2]} \Big)\otimes_{[0]} \Big( \big(a^\dagger_j\big)^{[1/2]} \otimes_{[0]} \big(a_j\big)^{[1/2]} \Big) \Big]\end{split}\]

[8]:

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

LX, LY = 4, 4
L = LX * LY
J = 0.4
N_ELEC = 14 # 1/8 doping
TWO_S = 0

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySU2, n_threads=4)

driver.set_symmetry_groups("U1Fermi", "SU2", "SU2")
Q = driver.bw.SX

# [Part A] Set states and matrix representation of operators in local Hilbert space
site_basis, site_ops = [], []

for k in range(L):
    basis = [(Q(0, 0, 0), 1), (Q(1, 1, 1), 1)] # [01]
    ops = {
        "": np.array([[1, 0], [0, 1]]),        # identity
        "C": np.array([[0, 0], [1, 0]]),       # a+
        "D": np.array([[0, 2**0.5], [0, 0]]),  # a
    }
    site_basis.append(basis)
    site_ops.append(ops)

# [Part B] Set Hamiltonian terms
driver.initialize_system(n_sites=L, vacuum=Q(0, 0, 0), target=Q(N_ELEC, TWO_S, TWO_S), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)
b = driver.expr_builder()

f = lambda i, j: i * LY + j if i % 2 == 0 else i * LY + LY - 1 - j

for i in range(0, LX):
    for j in range(0, LY):
        if i + 1 < LX:
            b.add_term("(C+D)0", [f(i, j), f(i + 1, j), f(i + 1, j), f(i, j)], -(2 ** 0.5))
            b.add_term("((C+D)2+(C+D)2)0", [f(i, j), f(i, j), f(i + 1, j), f(i + 1, j)], J * -(3 ** 0.5) / 2)
            b.add_term("((C+D)0+(C+D)0)0", [f(i, j), f(i, j), f(i + 1, j), f(i + 1, j)], J * -1 / 2)
        if j + 1 < LY:
            b.add_term("(C+D)0", [f(i, j), f(i, j + 1), f(i, j + 1), f(i, j)], -(2 ** 0.5))
            b.add_term("((C+D)2+(C+D)2)0", [f(i, j), f(i, j), f(i, j + 1), f(i, j + 1)], J * -(3 ** 0.5) / 2)
            b.add_term("((C+D)0+(C+D)0)0", [f(i, j), f(i, j), f(i, j + 1), f(i, j + 1)], J * -1 / 2)

# [Part C] Perform DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True), algo_type=MPOAlgorithmTypes.FastBipartite, iprint=1)
mps = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)
print("DMRG energy = %20.15f (per site = %10.6f)" % (energy, energy / L))


Build MPO | Nsites =    16 | Nterms =         96 | Algorithm = FastBIP | Cutoff = 1.00e-14
 Site =     0 /    16 .. Mmpo =     5 DW = 0.00e+00 NNZ =        5 SPT = 0.0000 Tmvc = 0.000 T = 0.032
 Site =     1 /    16 .. Mmpo =    10 DW = 0.00e+00 NNZ =       13 SPT = 0.7400 Tmvc = 0.000 T = 0.021
 Site =     2 /    16 .. Mmpo =    14 DW = 0.00e+00 NNZ =       18 SPT = 0.8714 Tmvc = 0.000 T = 0.041
 Site =     3 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       22 SPT = 0.9127 Tmvc = 0.000 T = 0.012
 Site =     4 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       22 SPT = 0.9321 Tmvc = 0.000 T = 0.009
 Site =     5 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.012
 Site =     6 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.011
 Site =     7 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.012
 Site =     8 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       22 SPT = 0.9321 Tmvc = 0.000 T = 0.015
 Site =     9 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.013
 Site =    10 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.013
 Site =    11 /    16 .. Mmpo =    18 DW = 0.00e+00 NNZ =       26 SPT = 0.9198 Tmvc = 0.000 T = 0.012
 Site =    12 /    16 .. Mmpo =    14 DW = 0.00e+00 NNZ =       22 SPT = 0.9127 Tmvc = 0.000 T = 0.021
 Site =    13 /    16 .. Mmpo =    10 DW = 0.00e+00 NNZ =       18 SPT = 0.8714 Tmvc = 0.000 T = 0.018
 Site =    14 /    16 .. Mmpo =     5 DW = 0.00e+00 NNZ =       13 SPT = 0.7400 Tmvc = 0.000 T = 0.011
 Site =    15 /    16 .. Mmpo =     1 DW = 0.00e+00 NNZ =        5 SPT = 0.0000 Tmvc = 0.000 T = 0.008
Ttotal =      0.263 Tmvc-total = 0.001 MPO bond dimension =    18 MaxDW = 0.00e+00
NNZ =          316 SIZE =         3486 SPT = 0.9094

Rank =     0 Ttotal =      0.480 MPO method = FastBipartite bond dimension =      18 NNZ =          316 SIZE =         3486 SPT = 0.9094

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      3.850 | E =      -8.9796008520 | DW = 5.86154e-07

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      7.778 | E =      -9.0284250184 | DE = -4.88e-02 | DW = 1.36125e-06

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      9.742 | E =      -9.0298619274 | DE = -1.44e-03 | DW = 1.23409e-06

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     10.604 | E =      -9.0298670256 | DE = -5.10e-06 | DW = 1.13415e-06

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     11.898 | E =      -9.0298686867 | DE = -1.66e-06 | DW = 2.27179e-10

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     12.869 | E =      -9.0298686872 | DE = -4.35e-10 | DW = 1.12668e-10

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     13.782 | E =      -9.0298686872 | DE = -8.19e-12 | DW = 2.27110e-10

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     14.710 | E =      -9.0298686872 | DE = -1.02e-11 | DW = 1.12643e-10

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     15.473 | E =      -9.0298686872 | DE = 2.91e-11 | DW = 1.78156e-19

DMRG energy =   -9.029868687160421 (per site =  -0.564367)