Hi,
I am recently to reproduce the charge density wave (CDW) state of a one-dimensional extended Bose-Hubbard model with open boundary condition. The model can be written as,
H=-t\sum_{ i}(b_i^\dag b_{i+1}+ h.c.)+\frac{U}{2}\sum_in_i(n_i-1)+V\sum_{i}n_in_{i+1},
According to Ref.[ Phys. Rev. B 61 , 12474(2000), url: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.61.12474], a CDW density profile can be seen at t/U=0.1, V/U=0.4 and particle filling 1/2 (See Fig.6 of the reference).
However, when I reproduce this result with total latice size L=100, I got a different density profile which have attached as following:
1 0.9124231273975183
2 0.14970860238647662
3 0.835638249019661
4 0.1815491194311886
5 0.8095297013017156
6 0.19752402219695384
7 0.7944770447783213
8 0.2084615825914796
9 0.7832714507898253
10 0.2175051554962943
11 0.7735159302715544
12 0.22592959282080677
13 0.76417425961646
14 0.23435322543735432
15 0.754727223902288
16 0.24311001953992215
17 0.7448909450185272
18 0.2523919529973155
19 0.7345033001932816
20 0.26231198650512755
21 0.7234718333265456
22 0.272934943886998
23 0.711747452635697
24 0.28429391711369256
25 0.6993101299072558
26 0.296399593544292
27 0.6861606086111528
28 0.30924584360623236
29 0.6723153101596289
30 0.3228132819701113
31 0.6578030265956872
32 0.33707166280292333
33 0.6426626403870428
34 0.35198159253882744
35 0.6269414507572181
36 0.3674958402693919
37 0.610693870838916
38 0.38356040170281824
39 0.5939803395228962
40 0.40011542338878675
41 0.5768663730939148
42 0.41709603867523465
43 0.5594216975939265
44 0.43443315599188825
45 0.5417194285504348
46 0.45205422269525203
47 0.5238352798074614
48 0.46988397756436084
49 0.5058467863170942
50 0.48784520097806666
51 0.487832538990757
52 0.505859463391281
53 0.4698714252347378
54 0.5238478772132479
55 0.45204187300950827
56 0.5417318525391266
57 0.43442109897398196
58 0.5594338569744906
59 0.4170843591190376
60 0.5768781816435744
61 0.4001041997533515
62 0.5939917167515864
63 0.3835497054720258
64 0.610704742979596
65 0.3674857352808265
66 0.6269517514966265
67 0.35197213376992875
68 0.6426723120061396
69 0.3370628942195659
70 0.6578120223101307
71 0.3228052372354159
72 0.6723235929203218
73 0.3092385450318608
74 0.6861681523639005
75 0.2963930526173842
76 0.6993169197860075
77 0.2842881337282399
78 0.7117534850206948
79 0.2729299062462747
80 0.7234771173360137
81 0.26230767511189923
82 0.734507855885176
83 0.2523883331345254
84 0.7448948006531395
85 0.24310705178269798
86 0.754730417606847
87 0.2343508666590108
88 0.7641768319826925
89 0.2259277883752413
90 0.7735179338028265
91 0.21750384420180302
92 0.7832729462962496
93 0.20846069076218043
94 0.7944780954181354
95 0.1975234371111525
96 0.8095304273224601
97 0.18154881186565536
98 0.8356386748208199
99 0.14970848815911933
100 0.9124233291074241
It shows a site-dependent oscillation ampltudide, which is different with my expection that the oscillation ampltuide should not changed in the bulk. So my question is: Am I missed something or did I make some mistakes? Any suggestions are much more appreciated.
Best,
Chen-Rong
A minimal code is attached.
using ITensors,ITensorMPS
mutable struct DemoObserver <: AbstractObserver
energy_tol::Float64
last_energy::Float64
DemoObserver(energy_tol=0.0) = new(energy_tol,1000.0)
end
function ITensors.checkdone!(o::DemoObserver;kwargs...)
sw = kwargs[:sweep]
energy = kwargs[:energy]
if abs(energy-o.last_energy) < o.energy_tol
println("Stopping DMRG after sweep $sw")
return true
end
# Otherwise, update last_energy and keep going
o.last_energy = energy
return false
end
d=6
function ITensors.op(::OpName"N^2", ::SiteType"Boson", d::Int;)
o = zeros(d, d)
o[2, 2] = 1
o[3, 3] = 4
o[4, 4] = 9
o[5, 5] = 16
o[6, 6] = 25
return o
end
include(get(ARGS, 1, "input.jl"))
sweeps = Sweeps(nsweep, sweeps_args)
sites = siteinds("Boson", N;conserve_qns=true,dim=6)
os = OpSum()
for b in 1:(N - 1)
os .+= -Jb, "Adag", b, "A", b + 1
os .+= -Jb, "A", b, "Adag", b + 1
end
for i in 1:N
os .+= Ub/2, "N^2", i
os .+= -Ub/2, "N", i
end
for i in 1:N-1
os .+= V, "N", i,"N", i+1
end
H = MPO(os, sites)
state=[isodd(n) ? "0" : "1" for n in 1:N ]
psi0 = MPS(sites, state)
# Check total number of particles:
@show flux(psi0)
obs = DemoObserver(etol)
# Start DMRG calculation:
energy, psi = dmrg(H, psi0,sweeps,observer=obs)
nmea = expect(psi, "N")
println("Density:")
for j in 1:N
println("$j $(nmea[j])")
end
println()
println("\nGround State Energy = $energy")
The input file is,
etol = 1E-11
N = 100
Jb = 0.1
Ub = 1.0
V = 0.4
nsweep = 800
sweeps_args = [
"maxdim" "mindim" "cutoff" "noise"
500 20 1e-12 1e-6
500 20 1e-12 1e-6
1000 20 1e-12 1e-6
1000 20 1e-12 1e-6
1000 20 1e-12 1e-6
1000 20 1e-12 1e-8
1800 20 1e-12 1e-8
1800 20 1e-12 1e-8
1800 20 1e-12 1e-8
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
1800 20 1e-12 0
]