Converting from SAnySU2 to SZ

[Shovan-Physics]

Hi, suppose I am working with the SU(2) t-J model as shown in the documentation, using the symmetry type SAnySU2 and target=Q(N_ELEC, TWO_S, TWO_S): https://block2.readthedocs.io/en/latest/tutorial/custom-hamiltonians.html#SU(2)-t-J-Model.

Now I want to convert the ground state to the symmetry type SZ where the basis states are labeled by N_ELEC and TWO_SZ. I'm not sure how to do this. mps_change_to_sz only works in the SU2 mode. There is also mps_change_symm(mps, tag, target), but I'm not sure what the proper syntax for target is and how to specify that I want TWO_SZ = TWO_S, for example.

[hczhai]

You can use the following script to do the translation from SAnySU2 to SAnySZ, which requires a few modifications compared with the one provided in the documentation:

1. Add AbelianPG symmetry subgroup in the driver.set_symmetry_groups.
2. Use an alternative argument list in driver.initialize_system to set the spin/projected spin explicitly via spin.
3. Add argument add_ident=False in driver.get_mpo.

After switching to the symmetry mode SAnySZ, use driver.initialize_system to set the desired TWO_SZ = TWO_S via the spin parameter. Then the MPO and MPS can be translated from SAnySU2 to SAnySZ using driver.mpo_change_symm and driver.mps_change_symm. The translated objects should only be used the SAnySZ symmetry mode.

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 = 2
TWO_SZ = 2

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

driver.set_symmetry_groups("U1Fermi", "SU2", "SU2", "AbelianPG")
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, 1, 1, 0), 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, n_elec=N_ELEC, spin=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, add_ident=False)
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))

# [Part D] Set states and matrix representation of operators in SAnySZ mode
driver.symm_type = SymmetryTypes.SAnySZ
driver.set_symmetry_groups("U1Fermi", "U1", "AbelianPG")
Q = driver.bw.SX

site_basis_u1, site_ops_u1 = [], []

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

driver.initialize_system(n_sites=L, n_elec=N_ELEC, spin=TWO_SZ, hamil_init=False)
driver.ghamil = driver.get_custom_hamiltonian(site_basis_u1, site_ops_u1)

# [Part E] Translate objects for using in SAnySZ mode
uket = driver.mps_change_symm(mps, 'U1-KET', driver.target)
umpo = driver.mpo_change_symm(mpo, 'U1-MPO')
impo = driver.get_identity_mpo()

print('E(U1)  = ', driver.expectation(uket, umpo, uket) / driver.expectation(uket, impo, uket))

[Shovan-Physics]

Thanks a lot for the script. However, if I run the same code with LX, LY = 10, 2, the two energies do not match. I get

DMRG energy =  -14.731348089307268 (per site =  -0.736567)
E(U1)  =  -9.991014746241166

Is it resulting from some truncation of the converted MPS to a particular bond dimension? How do we get the bond dimension of a given MPS, and is there a fix to the mismatch?

[hczhai]

Thanks for the feedback. The energy mismatch can be fixed by adding a line

driver.align_mps_center(mps, ref=0)

before

uket = driver.mps_change_symm(mps, 'U1-KET', driver.target)

The max bond dimension of an MPS can be obtained using

print('Max D(SU2) = ', mps.info.get_max_bond_dimension())
print('Max D(U1)  = ', uket.info.get_max_bond_dimension())

[Shovan-Physics]

Great, thanks!