Code: Select all
#To find neutral energy gap
import numpy as np
from tenpy.models.lattice import Chain
from tenpy.networks.site import BosonSite
from tenpy.models.model import CouplingModel, NearestNeighborModel, MPOModel
from tenpy.networks.mps import MPS
from tenpy.algorithms import dmrg
L=100
E_N_0=np.zeros(10) #Ground state energy
E_N_1=np.zeros(10) #1st Excited state
E_N_2=np.zeros(10) # second excited state
count=0 #index for above arrays
class BoseHubbardModelmuedge(CouplingModel,MPOModel): #fixing chemical potential at edges
def __init__(self,t,U,V,mu0,mu1,nmax,filling): #mu0 is chemical potential at leftmost site, mu1 is chemical potential at rightmost site
boson=BosonSite(Nmax=nmax,conserve='N',filling=filling)
lattice = Chain(L,boson,bc='open',bc_MPS='finite')
CouplingModel.__init__(self, lattice)
for u in range(len(self.lat.unit_cell)):
self.add_onsite(- U / 2., u, 'N')
self.add_onsite(U / 2., u, 'NN')
for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
self.add_coupling(-t, u1, 'Bd', u2, 'B', dx, plus_hc=True)
self.add_coupling(V, u1, 'N', u2, 'N', dx)
self.add_onsite_term(mu0,0,'N')
self.add_onsite_term(mu1,L-1,'N')
MPOModel.__init__(self, lattice, self.calc_H_MPO())
for VV in range(36,37,1):
V=VV/10
M = BoseHubbardModelmuedge(t=1,U=5,V=V,mu0=10,mu1=-10,nmax=2,filling=1) #"conserve":N
p=[]
for i in range(100):
p.append(1)
psi0=MPS.from_product_state(M.lat.mps_sites(),p,"finite")#creating an MPS state from a product state. The product state is all sites with single boson
dmrg_params = {"max_sweeps":20,"trunc_params": {"chi_max": 300, "svd_min": 1.e-10}}
H_MPO = M.calc_H_MPO(tol_zero = 1.0e-15)
states = dmrg_params['orthogonal_to']=[]
for i in range(3):
psii=psi0.copy()
results=dmrg.run(psii,M,dmrg_params)
states.append(psii)
psi_0=states[0] # ground state
psi_1=states[1] #first excited state
psi_2=states[2] #second excited state
E_N_0[count] = M.H_MPO.expectation_value(psi_0)
E_N_1[count] = M.H_MPO.expectation_value(psi_1)
E_N_2[count] = M.H_MPO.expectation_value(psi_2)
print("Ground state energy=",E_N_0[count])
print("First Excited state=",E_N_1[count])
print("Second Excited state=",E_N_2[count])
print("V=",V," : ","E_n=",E_N_1[count]-E_N_0[count])
print("-------------------------------------")
count=count+1
E_n=E_N_1-E_N_0 #neutral gap
print(E_n)
The output I am getting is
Code: Select all
Ground state energy= 193.12420756001185
First Excited state= 193.12420756001143
Second Excited state= 216.19939198650584
V= 3.6 : E_n= -4.263256414560601e-13