Measure Expectation values

[cfengno1]

Dear Huanchen,

Thank you so much for the wonderful package!

I am trying to measure the electron occupation for a custom 40x6 system after DMRG calculations with bond dimension=4000. The following is my code:

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

LX, LY = 40, 6
L = LX * LY
J = 1
V = 0
t = 3
N_ELEC = L
N_db = 12
TWO_S = 0

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySU2, n_threads=12, stack_mem=164748364800)

driver.set_symmetry_groups("U1Fermi", "U1", "SU2", "SU2") # electron number, doublon number, total spin, total spin
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, 0), 1), (Q(1, 0, 1, 1), 1), (Q(2, 1, 0, 0) ,1)] # [01]
    ops = {
        "": np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),        # identity
        "B": np.array([[0, 0, 0], [0, 0, 0], [0, -2**0.5, 0]]),       # a+ create doublon
        "H": np.array([[0, 0, 0], [1, 0, 0], [0, 0, 0]]),  # a+ annhilate holon
        "b": np.array([[0, 0, 0], [0, 0, 1], [0, 0, 0]]),       # a annhilate doublon
        "h": np.array([[0, 2**0.5, 0], [0, 0, 0], [0, 0, 0]]),  # a create holon
        "n": np.array([[0, 0, 0], [0, 1, 0], [0, 0, 2]]),  # electron number
        "d": np.array([[0, 0, 0], [0, 0, 0], [0, 0, 1]]),  # doublon number
    }
    site_basis.append(basis)
    site_ops.append(ops)

driver.initialize_system(n_sites=L, vacuum=Q(0, 0, 0, 0), target=Q(N_ELEC, N_db, TWO_S, TWO_S), hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)

mps = driver.load_mps(tag="KET", nroots=1)
N = driver.get_npdm(mps, npdm_expr='n')[0]
np.savetxt('Ne.txt',N)

1. when I add fermionic_ops='' into N = driver.get_npdm(mps, npdm_expr='n', fermionic_ops='')[0], I get the error message:

File "/scratch/feng/tmp/Read.py", line 42, in
N = driver.get_npdm(mps, npdm_expr='n',fermionic_ops='')[0]
File "/home/cfengno1/local/block2/block2-CPX/block2-preview/pyblock2/driver/core.py", line 5539, in get_npdm
assert fermionic_ops is None
AssertionError

2. If I drop fermionic_ops='', and run the above mode, it takes huge memeory~300G and is killed eventually. I stopped my DMRG calculations by ctrl+c, so the MPS center is checked to be 86. Can this be the reason? How to fix it?

THank you so much!

Best,
Feng

[hczhai]

Hi Feng,

The current support for custom SU2 operators in block2 is unfortunately incomplete so certain setups may lead to incorrect results (with or without an error message). You may need to carefully check the correctness of the state and energy for small size systems.

1. We cannot do driver.get_npdm or driver.expectation for MPSs with center site in the middle. You may first follow the steps in https://block2.readthedocs.io/en/latest/tutorial/restarting-dmrg.html#Restarting-from-Interruption to finish a partial DMRG sweep, and then compute the expectation values. The restart_dir argument can be used to automatically backup a "normal" MPS after each sweep.
2. driver.get_npdm is not implemented for custom SU2 operators. You may need to use driver.expectation to compute the expectation values one at a time.
3. When you build the MPO for expectation or DMRG with custom SU2 operators, b.finalize(adjust_order=False) is required, and the site indices given in the second argument of b.add_term should be in non-descending order. If the original order is not non-descending, you may need to do the SU2 algebra manually to adjust the operator order, figuring out the correct SU2 factors, and then call b.add_term.
4. With the above setups you still need external code to check the correctness of the number. If you get unexpected results then the calculation is possibly beyond the current capability of block2.

[cfengno1]

Dear Huanchen,

Thank you so much for your reply!

Regarding point 3, I see that in you custom t-J model example, the second argument of b.add_term is not exactly nondescending. See

            b.add_term("(C+D)0", [f(i, j), f(i + 1, j), f(i + 1, j), f(i, j)], -(2 ** 0.5))

Can you why in this case we can call adjust_order=true?

Best,
Feng

[hczhai]

1. For standard fermionic SU2 operators C and D, the block2 code implemented the algorithm to adjust the order automatically to make it non-descending, so you can use adjust_order=True.
2. For a set of custom operators with various multiplicity and fermionic/non-fermionic properties, the above adjust_order=True algorithm is not supposed to work.
3. When adjust_order=False the non-descending is only required for groups of indices within each term, so for b.add_term("(C+D)0", [i, j, k, l], -(2 ** 0.5)) only i <= j and k <= l are required.

[cfengno1]

Dear Huanchen,

Thank you for your explanation!

I am still puzzled by the adjust_order parameters. In your example for SU(3) Heisenberg model https://block2.readthedocs.io/en/latest/tutorial/custom-hamiltonians.html#SU(3)-Heisenberg-Model, you also used adjust_order=true.

I am now defining spin operators for the tJ model:

       "z": np.array([[0, 0, 0], [0, 0.5, 0], [0, 0, -0.5]]),       # Sz
        "p": np.array([[0, 0, 0], [0, 0, 1], [0, 0, 0]]),       # S+
        "m": np.array([[0, 0, 0], [0, 0, 0], [0, 1, 0]]),       # S-

Do I need to adjust the order for these operators myself?

Best,
Feng

[hczhai]

In your example for SU(3) Heisenberg model, you also used adjust_order=true.

As indicated in the documentation, "We will only use Abelian symmetry groups for this problem". You can use adjust_order=True if you only have Abelian symmetry groups (like U1, U1Fermi), etc. In other words, you can use adjust_order=True if you do not have SU2 (which is non-Abelian) as a symmetry sub-group in driver.set_symmetry_groups.