Infinite DMRG with conserved charges

[pasa]

Good morning everyone,

I am again in need of your knowledge. I am trying to use infinite DMRG implementing a model which has conservation laws not yet implemented in Tenpy. So far I used either finite DMRG with conserved charges, or infinite DMRG with no charges. Those codes have worked well giving meaningful results.

First, to make it easier I started implementing the model \(H=\sum_j b_j b_{j+1} + H.c.\) which has a conserved quantity of the form \(Q=\sum_j (-1)^j n_j\) with \([b_i,b^\dagger_j]=\delta_{ij}\) and \(n_j=b^\dagger_j b_j\). This apparently runs fine with infinite DMRG and unit cell L=2. Let me share the full code. Here I am not entirely sure about one of the lines.

Code: Select all

from tenpy.algorithms import dmrg
from tenpy.models.hubbard import BoseHubbardModel
from tenpy.models.model import NearestNeighborModel as NNM
import pickle
import numpy as np
import time
from tenpy.linalg import np_conserved as npc
from tenpy.networks.site import BosonSite
import tenpy
from tenpy.networks.mps import MPS

class BoseModel(BoseHubbardModel):

    def init_sites(self, model_params):
        Nmax = model_params.get('Nmax', 2)
        
        L = model_params.get('L', 2)
        alpha_k = [1,-1]
        
        chinfo = npc.ChargeInfo([1], ['N*alpha'])
        
        sites = []
        for i in range(L):  # what is L here?
            site = BosonSite(Nmax=Nmax, conserve='N')
            leg = site.leg
            alpha = alpha_k[i]
            new_charge = np.hstack([leg.charges*alpha])
            new_leg = npc.LegCharge(chinfo, leg.slices, new_charge)
            site.change_charge(new_leg)
            sites.append(site)
        return sites

    def init_lattice(self, model_params):
        sites = self.init_sites(model_params)
        bc_MPS = model_params.get('bc_MPS', 'finite')
        bc = 'open' if bc_MPS == 'finite' else 'periodic'
        lat = tenpy.models.lattice.TrivialLattice(sites, bc_MPS=bc_MPS, bc=[bc])
        return lat

    def init_terms(self, model_params):
        
        J = model_params.get('J', 0.)
        mu = model_params.get('mu', 0.)
        L = self.lat.N_sites
        assert self.lat.bc_MPS == 'infinite', "indexing only works for finite!" #Not sure about this line!
        for u in range(L - 1):
            self.add_multi_coupling_term(J, [u, u+1], ['B', 'B'], ['Id']) #,plus_hc=True)
            self.add_multi_coupling_term(J, [u, u+1], ['Bd', 'Bd'], ['Id'])
        
        for u in range(L):
            self.add_onsite_term(-mu, u, 'N')
            



def iDMRG(J3,U,mu, n_max,chi,L,bcond,verbose=2):
    #Additional non-particle conserving term
    model_params = dict(L=L, J=J3, U=U, mu=mu, Nmax= n_max,   bc_MPS=bcond, conserve=['N*alpha'], verbose=verbose) 
    print("parameters: ", model_params)
    M = BoseModel(model_params)
    # print(M.all_coupling_terms().to_TermList())
    vec_j = np.zeros(n_max+1)
    vec_j[1] = 1
    product_state = [vec_j]*M.lat.N_sites #init state
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    
    chi_list = {}
    for n, nchi in enumerate(np.arange(5,60,10)):
        chi_list[10*n] = nchi
    chi_list[10*(n+1)] = chi

    dmrg_params = {
        'mixer': True,
        'mixer_params': { 'decay':2, 'amplitude':1.e-1, 'disable_after':50},
        'max_E_err': 1.e-10,
        'max_S_err': 1.e-8,
        'min_sweeps': 100,
        'max_sweeps': 2000,
        'lanczos_params': { 'N_max': 50},
        'chi_list': chi_list,
        'trunc_params': {
#            #'chi_max': D_bond,
            'svd_min': 1.e-10,
        },

        'start_env':30,
        'verbose':2,
    }
    eng = dmrg.TwoSiteDMRGEngine(psi,M, dmrg_params)
    E, psi = eng.run()
    
   
    return  psi, E
    
if __name__ == "__main__":
    L = 2
    U = 1.
    chi = 30
    n_max = 5
    bcond = 'infinite'
    J3 = 1.
    mu = 0.5
    
     
    psi, E= iDMRG(J3,U,mu, n_max,chi,L,bcond)
        

Nonetheless, the model I am actually interested in is this one \(H=\sum_j b_{j-1} b_j^\dagger b_{j+1}+ H.c.\) which has these two conserved quantities \( Q_1=\sum_j \alpha_j n_j\) with \( \alpha_j^1=[1,0,-1,-1,0,1] \) mod(6) and \( \alpha_j^2=[0,1,1,0,-1,-1] \) mod(6). Here the unit cell is L=6. The definition of model and charges is given in the following:

Code: Select all

from tenpy.algorithms import dmrg
from tenpy.models.hubbard import BoseHubbardModel
from tenpy.models.model import NearestNeighborModel as NNM
import pickle
import numpy as np
import time
from tenpy.linalg import np_conserved as npc
from tenpy.networks.site import BosonSite
import tenpy
from tenpy.networks.mps import MPS

class BoseModel(BoseHubbardModel):

    def init_sites(self, model_params):
        Nmax = model_params.get('Nmax', 2)
        
        L = model_params.get('L', 6)
        alpha1_k = [1,0,-1,-1,0,1]
        alpha2_k = [0,1,1,0,-1,-1]
        
        chinfo = npc.ChargeInfo([1,1], ['N*alpha1','N*alpha2'])
        
        sites = []
        for i in range(L):  # what is L here?
            site = BosonSite(Nmax=Nmax, conserve='N')
            leg = site.leg
            alpha1 = alpha1_k[i]
            alpha2 = alpha2_k[i]
            new_charge = np.hstack([leg.charges*alpha1,leg.charges*alpha2])
            new_leg = npc.LegCharge(chinfo, leg.slices, new_charge)
            site.change_charge(new_leg)
            sites.append(site)
        return sites

    def init_lattice(self, model_params):
        sites = self.init_sites(model_params)
        bc_MPS = model_params.get('bc_MPS', 'finite')
        bc = 'open' if bc_MPS == 'finite' else 'periodic'
        lat = tenpy.models.lattice.TrivialLattice(sites, bc_MPS=bc_MPS, bc=[bc])
        return lat

    def init_terms(self, model_params):
        
        J = model_params.get('J', 0.)
        mu = model_params.get('mu', 0.)
        muA = model_params.get('muA', 0.)
        muB = model_params.get('muB', 0.)
        L = self.lat.N_sites
        assert self.lat.bc_MPS == 'infinite', "indexing only works for finite!"
        for u in range(L - 2):
            print(u)
            self.add_multi_coupling_term(J, [u, u+1, u+2], ['B', 'Bd', 'B'], ['Id', 'Id']) #,plus_hc=True)
            self.add_multi_coupling_term(J, [u, u+1, u+2], ['Bd', 'B', 'Bd'], ['Id', 'Id'])
        
        alpha1_k = np.array([1,0,-1,-1,0,1])
        alpha2_k = np.array([0,1,1,0,-1,-1])
        
        for u in range(L):
            self.add_onsite_term(-mu, u, 'N')
            self.add_onsite_term(muA*alpha1_k, u, 'N')
            self.add_onsite_term(muB*alpha2_k, u, 'N')
            



def iDMRG(J3,U,mu,muA, muB, n_max,chi,L,bcond,verbose=2):
    #Additional non-particle conserving term
    model_params = dict(L=L, J=J3, U=U, mu=mu, muA=muA, muB=muB, Nmax= n_max,   bc_MPS=bcond, conserve=['N*alpha1','N*alpha2'], verbose=verbose) 
    print("parameters: ", model_params)
    M = BoseModel(model_params)
    # print(M.all_coupling_terms().to_TermList())
    vec_j = np.zeros(n_max+1)
    vec_j[1] = 1
    product_state = [vec_j]*M.lat.N_sites #init state
    psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
    
    chi_list = {}
    for n, nchi in enumerate(np.arange(5,60,10)):
        chi_list[10*n] = nchi
    chi_list[10*(n+1)] = chi

    dmrg_params = {
        'mixer': True,
        'mixer_params': { 'decay':2, 'amplitude':1.e-1, 'disable_after':50},
        'max_E_err': 1.e-10,
        'max_S_err': 1.e-8,
        'min_sweeps': 100,
        'max_sweeps': 2000,
        'lanczos_params': { 'N_max': 50},
        'chi_list': chi_list,
        'trunc_params': {
#            #'chi_max': D_bond,
            'svd_min': 1.e-10,
        },

        'start_env':30,
        'verbose':2,
    }
    eng = dmrg.TwoSiteDMRGEngine(psi,M, dmrg_params)
    E, psi = eng.run()
  
   
    return  psi, E
    
if __name__ == "__main__":
    L = 6
    U = 1.
    chi = 30
    n_max = 5
    bcond = 'infinite'
    J3 = 1.
    sign = +1
    mu = 0.5
    muA = 1.
    muB = 0.
    
    
    psi, E= iDMRG(J3,U,mu,muA, muB, n_max,chi,L,bcond)

Nonetheless, I now get the following error

Code: Select all

H_MPO = self.calc_H_MPO()
  File "/home//TeNPy/tenpy/models/model.py", line 1640, in calc_H_MPO
    ot.remove_zeros(tol_zero)
  File "/home//TeNPy/tenpy/networks/terms.py", line 389, in remove_zeros
    if abs(term[op]) < tol_zero:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
/home/pablo/TeNPy/tenpy/tools/params.py:237: UserWarning: unused options for config BoseModel:
['U', 'conserve', 'verbose']
  warnings.warn(msg.format(keys=sorted(unused), name=self.name))

Could anyone tell me where the error is? Again, I think it has to be something related to the MPS and lattice boundary conditions but my several attempts did not work out. It might just be something stupid I might be missing.


Thank you very much.

Best,
Pablo

[Johannes]

Just to make sure I understand: when you say \([b_i,b^\dagger_j]=i\delta_{ij}\) , was the \(i\) supposed to be there??? I guess not, correct?
In the Hamiltonians you wrote, each b is acting on the same site j - I think you mean j, j+1 and j+2 (for 3 sites), correct?

That said, in your models, you don't couple the unit cells for infinite MPS, thus effectively simulating independent copies of the finite system of just L sites.
To add couplings between the unit cells, you can just change the for u in range(L - 1): and for u in range(L-2) to run until L for 'infinite' MPS (but not for finite).
In other words something like:

Code: Select all

N_terms = L-1 if bc == 'finite' else L
for u in range(N_terms):
    self.add_multi_coupling_term(J, [u, u+1, u+2], ['B', 'Bd', 'B'], ['Id', 'Id'] ,plus_hc=True)

Note that you had the wrong operators ['B', 'B', 'B'] in your example!

Finally, the error message you currently get is due to the onsite terms which you have: you forgot a [u] selecting the entries of alpha1_k and alpha2_k.

[pasa]

Dear Johannes,

thank you very much for your quick answer.

Indeed, all your first corrections are typos including the

Code: Select all

 ['B','B','B'] 

in the code. I have modified them in the original post to avoid confusions. This time I checked I did not introduce new ones .

I was exactly missing the couplings between different terms when using infinite boundary conditions. Thank you very much for helping to find that out.

Best,
Pablo