Open emstoudenmire opened 4 years ago
Just for added clarity, I have put the code here as well. This was code that I wrote up to attempt to answer this question on the Support Q&A (http://itensor.org/support/2140/implementing-local-gauge-symmetries):
ferm_full = QN("Sz",-1) => 1
ferm_vac = QN("Sz",0) => 1
antiferm_full = QN("Sz",+1) => 1
antiferm_vac = QN("Sz",0) => 1
even_link_p1 = QN("Sz",+2) => 1 # p1 = "positive 1"
even_link_0 = QN("Sz",0) => 1
even_link_n1 = QN("Sz",-2) => 1 # n1 = "negative 1"
odd_link_p1 = QN("Sz",-2) => 1
odd_link_0 = QN("Sz",0) => 1
odd_link_n1 = QN("Sz",+2) => 1
function Kogut_Susskind_sites(num_spatial_sites)
sites = Index[]
for i = 1:num_spatial_sites
push!(sites, Index(ferm_full,ferm_vac; tags="Site,S=1/2,n=$(4*i-3)"))
push!(sites, Index(odd_link_p1,odd_link_0,odd_link_n1; tags="Site,S=1,n=$(4*i-2)"))
push!(sites, Index(antiferm_full,antiferm_vac; tags="Site,S=1/2,n=$(4*i-1)"))
push!(sites, Index(even_link_p1,even_link_0,even_link_n1; tags="Site,S=1,n=$(4*i)"))
end
return sites
end
function Schwinger_Hamiltonian(sites, num_spatial_sites, x, μ)
N = num_spatial_sites * 2
ampo_H = AutoMPO()
for j=0:N-1
add!(ampo_H, 1, "Sz2", 2*j+2)
add!(ampo_H, (μ/2)*(-1)^j * 2, "Sz", 2*j+1)
if j != N-1
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*j+1, "S-", 2*j+2, "S-", 2*j+3)
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*j+3, "S+", 2*j+2, "S-", 2*j+1)
end
end
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*N-1, "S-", 2*N, "S-", 1)
add!(ampo_H, x * 2*2/sqrt(2), "S+", 1, "S+", 2*N, "S-", 2*N-1)
H = MPO(ampo_H, sites)
return H
end
function Schwinger_DMRG(num_spatial_sites)
sites = Kogut_Susskind_sites(num_spatial_sites)
H = Schwinger_Hamiltonian(sites, num_spatial_sites, x, μ)
sweeps = Sweeps(10)
maxdim!(sweeps, 10,20,100,100,200)
cutoff!(sweeps, 1E-10)
# construct initial state as completely empty
init_state = ["Up" for n = 1:4*num_spatial_sites] # just to initialize
for n = 1:4*num_spatial_sites
if n % 4 == 1
init_state[n] = "Dn" # empty fermion sites
elseif n % 4 == 3
init_state[n] = "Up" # empty antifermion states
else # n % 4 == 2 || n % 4 == 0
init_state[n] = "Up" # flux = 0 at all links
end
end
psi0 = randomMPS(sites, init_state, 100) # replace 100 with some value
energy, psi = dmrg(H, psi0, sweeps)
println("Ground state energy = $energy")
return energy
end
I have run "Schwinger_DMRG(10)" (i.e. 10 spatial sites) for every value of the bond dimension from 2 to 10,000 (replacing the 100 value in the randomMPS function), and it produced the same error every time, but if I use a value of 1, then no error is thrown.
I have also included a copy of the error message below:
ERROR: LoadError: MPS center bond dim less than requested
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] randomizeMPS!(::ITensors.MPS, ::Array{ITensors.Index,1}, ::Int64) at C:\Users\user\.julia\packages\ITensors\tgmxZ\src\mps\mps.jl:114
[3] randomMPS at C:\Users\user\.julia\packages\ITensors\tgmxZ\src\mps\mps.jl:208 [inlined]
[4] Schwinger_DMRG(::Int64) at d:\Julia Files\Schwinger_periodic.jl:166
[5] top-level scope at d:\Julia Files\Schwinger_periodic.jl:174
[6] include(::Module, ::String) at .\Base.jl:377
[7] include at c:\Users\user\.vscode\extensions\julialang.language-julia-0.15.40\scripts\debugger\debugger.jl:1 [inlined]
[8] run_request(::Base.PipeEndpoint, ::Main.VSCodeDebugger.DebuggerState, ::String, ::String) at c:\Users\user\.vscode\extensions\julialang.language-julia-0.15.40\scripts\debugger\debugger_requests.jl:9
[9] startdebug(::String) at c:\Users\user\.vscode\extensions\julialang.language-julia-0.15.40\scripts\debugger\debugger.jl:87
[10] top-level scope at c:\Users\user\.vscode\extensions\julialang.language-julia-0.15.40\scripts\debugger\run_debugger.jl:8
[11] include(::Module, ::String) at .\Base.jl:377
[12] exec_options(::Base.JLOptions) at .\client.jl:288
[13] _start() at .\client.jl:484
in expression starting at d:\Julia Files\Schwinger_periodic.jl:174
Noting a related issue shared by a user in this forum post.
Briefly, the issue is that for certain initial states using "Fermion" sites (representing "split" up and down electrons), randomMPS
always returns a product state even when a larger bond dimension is requested, and even though when keeping the sites together as "Electrons" a larger bond dimension is possible for the randomMPS
algorithm.
Apparently sometimes the
randomMPS
function which takes an initial state vector is not meeting the requested linkdim value. This was reported for an input state which was a S=1/2 Neel state and a linkdim of 4 in this forum post: http://itensor.org/support/2186/questions-random-initial-conservation-randommps-linkdim