QN randomMPS not meeting requested linkdim

[emstoudenmire]

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

[sujay-kazi]

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

[emstoudenmire]

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.