Custom Hamiltonians

[1]:

!pip install block2==0.5.3rc6 -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).

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) is 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.379 | E =      -6.2256341447 | DW = 2.65e-16

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.487 | E =      -6.2256341447 | DE = -1.15e-14 | DW = 4.93e-16

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

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

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

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

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

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

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

DMRG energy =   -6.225634144657922

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 =     90.266 | E =      -6.9568929201 | DW = 3.62e-09

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

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

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

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

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    158.074 | E =      -6.9568932112 | DE = 8.88e-16 | DW = 7.14e-20

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

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

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

DMRG energy =   -6.956893211182315

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 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", "U1Fermi")
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.176 | E =      -6.2256341447 | DW = 5.28e-16

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

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.313 | E =      -6.2256341447 | DE = -8.88e-16 | DW = 1.63e-16

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

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

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

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

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

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

DMRG energy =   -6.225634144658390

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", "Z2Fermi")
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 =     35.522 | E[ 10] =      -6.9958183038     -6.5705972882     -6.5705972844     -6.5705972810     -6.2773135424     -6.2256341358     -6.1122535707     -6.1122534909     -6.1122534352     -6.0320900882 | DW = 4.44e-09

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     39.172 | 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.21e-09

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

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

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     49.388 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -1.55e-13 | DW = 8.14e-17

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

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     54.521 | E[ 10] =      -6.9958183059     -6.5705972895     -6.5705972895     -6.5705972895     -6.2773135514     -6.2256341447     -6.1122535739     -6.1122535739     -6.1122535739     -6.0320902499 | DE = -1.78e-15 | DW = 2.12e-17

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     57.156 | 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 = 5.80e-18

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     59.782 | 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-16 | DW = 1.09e-18

Root = 0 <E> =   -6.995818305900327 <N> =      6.000
Root = 1 <E> =   -6.570597289542966 <N> =      6.000
Root = 2 <E> =   -6.570597289542286 <N> =      6.000
Root = 3 <E> =   -6.570597289541241 <N> =      6.000
Root = 4 <E> =   -6.277313551398623 <N> =      6.000
Root = 5 <E> =   -6.225634144677048 <N> =      8.000
Root = 6 <E> =   -6.112253573865427 <N> =      6.000
Root = 7 <E> =   -6.112253573864907 <N> =      6.000
Root = 8 <E> =   -6.112253573863808 <N> =      6.000
Root = 9 <E> =   -6.032090249941684 <N> =      4.000

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

[6]:

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.422 | E =      -0.5128539883 | DW = 1.32e-05

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     18.071 | E =      -0.5145284723 | DE = -1.67e-03 | DW = 3.73e-08

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     23.571 | E =      -0.5145507220 | DE = -2.22e-05 | DW = 3.84e-07

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     25.998 | E =      -0.5145508502 | DE = -1.28e-07 | DW = 4.73e-07

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     33.865 | E =      -0.5145512426 | DE = -3.92e-07 | DW = 8.02e-09

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     47.131 | E =      -0.5145512486 | DE = -5.93e-09 | DW = 1.24e-08

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     60.265 | E =      -0.5145512511 | DE = -2.51e-09 | DW = 5.33e-09

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     73.415 | E =      -0.5145512506 | DE = 4.86e-10 | DW = 1.74e-09

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     83.304 | E =      -0.5145512503 | DE = 3.24e-10 | DW = 2.58e-18

DMRG energy (per site) =   -0.514551250253768