Electron number unfixed simulations

[cfengno1]

Hi Huanchen,

I would like to simulate a system without electron number conservation, i.e. grand canonical Hamiltonian, and I tested your example code for Hubbard model with Z2xZ2 symmetry https://block2.readthedocs.io/en/latest/tutorial/custom-hamiltonians.html#Custom-Symmetry-Groups.

However, I find that I always got a state with a fixed integer number of electrons throughout the DMRG calculations. I am wondering how to make non particle number eigenstate accessible?

Thank you so much!

Best,
Feng

[hczhai]

Hi, the typical procedure is adding a chemical potential term in your Hamiltonian.

[cfengno1]

Hi Huanchen,

Thank you for your reply!

I have added the chemical potential term. Actually, without the chemical term is the same as setting mu=0.

The thing is that I always get an integer number of electrons. I think my code does not access the superposition of states with different electron numbers. Physically speaking, I would like to target symmetry broken state that hosts superconductivity like in https://www.pnas.org/doi/abs/10.1073/pnas.2109978118

The following is my code:

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

L = 20
U = 8
N_ELEC = 8
TWO_SZ = 0
mu=-0.5

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

# 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
        "n": np.array([[0, 0, 0, 0], [0, 2, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), # N
    }
    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)
b.add_term("n", np.array([[i] for i in range(L)]).ravel(), -mu)

# [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=500, nroots=1)
energies = driver.dmrg(mpo, mps, n_sweeps=20, bond_dims=[500] * 10 + [1000] * 10,
    noises=[1e-1] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=20, iprint=1)

row = driver.get_npdm(mps, npdm_expr='n', fermionic_ops='', mask=[0])[0]
np.savetxt('Ne.dat', row, delimiter=',')

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)
n_expt = driver.expectation(mps,partile_n_mpo,mps)

Thank you so much for your help!

Best,
Feng

[hczhai]

To get fractional number of electrons you still need to adjust the Hamiltonian parameters, sweep schedules, pinning field, and MPS bond dimension.

The paper that you linked actually studied 2D systems. As the exact solution does not need to break symmetry, and DMRG is almost exact for 1D systems, broken-symmetry DMRG solution is unlikely to appear for 1D system unless you use very small bond dimension and carefully adjust the settings.

For your case the following setup may produce broken-symmetry solutions:

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

L = 20
U = 8
N_ELEC = 8
TWO_SZ = 0
mu=-0.15

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

# 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)
b.add_term("cd", np.array([[i, i] for i in range(L)]).ravel(), -mu)
b.add_term("CD", np.array([[i, i] for i in range(L)]).ravel(), -mu)

nroots = 10

# [Part C] Perform state-averaged DMRG
mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops="cdCD"), algo_type=MPOAlgorithmTypes.FastBipartite)
b = driver.expr_builder()

mps = driver.get_random_mps(tag="KET", bond_dim=500, nroots=nroots)

energies = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[14] * 10,
    noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-5] * 8, dav_max_iter=1000, iprint=1)

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)
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(nroots):
    n_expt = driver.expectation(kets[ir], n_mpo, kets[ir])
    print("Root = %d <E> = %20.15f <N> = %10.3f" % (ir, energies[ir], n_expt))

[cfengno1]

Hi Huanchen,

Thank you for the explanation!

I am thinking about applying pinning field first and then turn off the pinning field to get spontaneous symmetry breaking. In that case, I need to load the state from pinning field DMRG calculations as initial state for non-pinning DMRG calculations. You have a tutorial on restarting dmrg https://block2.readthedocs.io/en/latest/tutorial/restarting-dmrg.html, but will it use new Hamiltonian MPO and construct new environment tensors?

Best,
Feng

[hczhai]

When you load MPS it will only load the MPS and you can change MPO after reloading. The environment will be regenerated.

[cfengno1]

OK. Thanks!

A follow-up question the other way round: in very large scale simulations, it is very helpful to load the environment from the previous run. Can we do that when restart DMRG in block2?

Best,
Feng

[hczhai]

No, we do not support loading the environments when restarting.