How to do non-interacting fermions?

ITensor Julia Questions
,
1

I am trying to run a calculation for non-interacting (spinful) fermions to check something, so I installed iTensor to check it out. This is the first time I’ve used the Julia version (and I don’t really know any Julia). I am getting an error that I can’t work out how to fix. My code is:


using ITensors

let

  N = 100
  sites = siteinds("Electron", N; conserve_qns=true)
  op = OpSum()
  t=1.0
  for j=1:N-1
      op += -t,"Cdagup",j,"Cup",j+1
      op += -t,"Cdagup",j+1,"Cup",j
      op += -t,"Cdagdn",j,"Cdn",j+1
      op += -t,"Cdagdn",j+1,"Cdn",j
  end

  H = MPO(op,sites)

  nsweeps = 30
  maxdim = [2,2,3,4,6,8,12,16,23,32,46,64,91,128,181,256,362,511,600,600,600,600,600,600,600,600,600,600,600,600] # gradually increase states kept
  cutoff = [1E-14] # desired truncation error
  noise = [0.0]

  states = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  print(states)
  psi0 = randomMPS(sites,states)

  print(psi0)

  energy,psi = dmrg(H,psi0; nsweeps, maxdim, cutoff, noise)

  return
end

But this gives an error:

ERROR: LoadError: ArgumentError: initial vector should not have norm zero
Stacktrace:
  [1] initialize(iter::KrylovKit.LanczosIterator{ProjMPO, ITensor, KrylovKit.ModifiedGramSchmidt2}; verbosity::Int64)
    @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/factorizations/lanczos.jl:171
  [2] initialize
    @ ~/.julia/packages/KrylovKit/r8GLV/src/factorizations/lanczos.jl:167 [inlined]
  [3] eigsolve(A::ProjMPO, x₀::ITensor, howmany::Int64, which::Symbol, alg::KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2, Float64})
    @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/lanczos.jl:11
  [4] #eigsolve#38
    @ ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/eigsolve.jl:200 [inlined]
  [5] eigsolve
    @ ~/.julia/packages/KrylovKit/r8GLV/src/eigsolve/eigsolve.jl:182 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:237 [inlined]
  [7] macro expansion
    @ ~/.julia/packages/TimerOutputs/RsWnF/src/TimerOutput.jl:253 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:236 [inlined]
  [9] macro expansion
    @ ./timing.jl:395 [inlined]
 [10] 
    @ ITensors ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:204
 [11] dmrg
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:156 [inlined]
 [12] #dmrg#1059
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:27 [inlined]
 [13] dmrg
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:20 [inlined]
 [14] #dmrg#1065
    @ ~/.julia/packages/ITensors/WMeVS/src/mps/dmrg.jl:391 [inlined]
 [15] top-level scope
    @ ~/dmrg.jl:29
 [16] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [17] top-level scope

I don’t know what this means, since the psi0 is initialized in the same way as the example codes. In desparation, looking at the differences between my code and the 2D hubbard model example on Github, I added the U term:

using ITensors

let

  N = 100
  sites = siteinds("Electron", N; conserve_qns=true)
  op = OpSum()
  t=1.0
  for j=1:N-1
      op += -t,"Cdagup",j,"Cup",j+1
      op += -t,"Cdagup",j+1,"Cup",j
      op += -t,"Cdagdn",j,"Cdn",j+1
      op += -t,"Cdagdn",j+1,"Cdn",j
  end

  U=0.01
  for n in 1:N
    op += U, "Nupdn", n
  end

  H = MPO(op,sites)

  nsweeps = 30
  maxdim = [2,2,3,4,6,8,12,16,23,32,46,64,91,128,181,256,362,511,600,600,600,600,600,600,600,600,600,600,600,600] # gradually increase states kept
  cutoff = [1E-14] # desired truncation error
  noise = [0.0]

  states = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  print(states)
  psi0 = randomMPS(sites,states)

  print(psi0)

  energy,psi = dmrg(H,psi0; nsweeps, maxdim, cutoff, noise)

  return
end

and suddenly it works! But I want non-interacting fermions. Setting U=0 brings back the original problem. In fact even U=1e-3 has problems, but U=1e-2 works. But I want U=0. How do I do non-interacting fermions?

2

Thanks for letting us know. I find this behavior a bit surprising too. But I think I figured out the reason, which is that for whatever reason, the product state \ket{\uparrow\downarrow\uparrow\downarrow \cdots} has exactly zero expectation value with H. So then it’s not a good initial state to use with DMRG, and that can explain the zero norm error you got (which came from inside KrylovKit, the library we use to solve the “mini diagonalization” or Krylov eigensolver step inside of each step of DMRG).

A fix for this which might be what you need is to pass the linkdims argument to randomMPS to have it make an MPS which has a bond dimension greater than 1.

psi0 = randomMPS(sites,states; linkdims=10)

When I make this change, I find the error goes away. Please let me know if the error persists for you.

3

Thanks Miles, yes that works for me.

But (somewhat coincidentally) that state is actually fairly similar to what I want to use as the initial state; not exactly the same, but similar properties. What I am trying to do is confirm the behavior of 2-site DMRG in some settings similar to the calculations in https://journals.aps.org/prl/supplemental/10.1103/PhysRevLett.130.246402, eg figs S-6 (Electron), S-9 (spinless fermion) etc, with specific initial states.

It seems strange that a state that has \langle \psi_0 \vert H \vert \psi_0 \rangle = 0 leads to an error - that shouldn’t normally be a problem with Lanczos, it just means that the first diagonal element, \alpha_0, of the sub-matrix is zero. H \vert \psi_0 \rangle is not a zero vector (and even if it was, it just means that the Lanczos has found an eigenvalue and can stop immediately - that should be the case in the calculation from fig S-9, for example).

As an alternative to the U term or changing the initial state, I tried adding

  for j=1:N
    op += 1.0,"I",j
  end

and this also works fine, with the original product initial state. EDIT: sorry, I messed up – I forgot to remove the ; linkdims=10 part. Shifting the Hamiltonian by a constant doesn’t work. But that makes it even more mysterious!

The location of the error in lanczos.jl is checking the initial norm \beta_0, but this shouldn’t be zero. I added

  psi0 = randomMPS(sites,states)
  print(norm(psi0))

and this prints 1.0.

4

Hi Miles,

Just a quick update. I don’t understand at all what psi0 = random(sites,states;linkdims=10) is doing, or how to interpret the output of print(psi0), but as far as I can tell the calculation is doing what I want anyway (fingers crossed).

With psi0 = random(sites,states;linkdims=10) I saw it crash once, with

After sweep 1 energy=-94.43063223563396  maxlinkdim=2 maxerr=4.81E-01 time=8.115
After sweep 2 energy=-93.29194763551637  maxlinkdim=2 maxerr=5.00E-01 time=0.118
ERROR: LoadError: ArgumentError: initial vector should not have norm zero
Stacktrace:
  [1] initialize(iter::KrylovKit.LanczosIterator{ProjMPO, ITensor, KrylovKit.ModifiedGramSchmidt2}; verbosity::Int64)
    @ KrylovKit ~/.julia/packages/KrylovKit/r8GLV/src/factorizations/lanczos.jl:171
...

but since then its been working OK, which is a low enough probability of failure for me :wink:

Just a check: this calculation is using standard 2-site DMRG? Real fp64 arithmetic? I’ve installed MKL using the instructions at GitHub - JuliaLinearAlgebra/MKL.jl: Intel MKL linear algebra backend for Julia, is there a way to verify that iTensor is actually using it? It didn’t change the CPU time before & after, but I’ve got MKL set up as the default BLAS anyway on my Ubuntu box, so it was probably already using it. Finally, is there an option for single-site DMRG? (I don’t need it, but just for an interesting comparison…)

Edit: is there a way to find out how many Hamiltonian matrix-vector multiplies it is doing per iteration?

5

@ianmcc a few comments:

  1. using LinearAlgebra; BLAS.get_config() shows your BLAS/LAPACK backend. If it shows MKL is enabled, that’s what ITensors.jl will be using.
  2. By default it uses real, double precision floating point arithmetic. You can check with ITensors.scalartype(psi0) and ITensors.scalartype(H), or eltype(psi0[1]), eltype(H[1]), etc. to check the element types of each tensor of the MPS/MPO. You can change to single precision with psi0 = randomMPS(Float32, sites, states; linkdims=10) and H = MPO(Float32, op, sites).
  3. You can just write maxdim = [2,2,3,4,6,8,12,16,23,32,46,64,91,128,181,256,362,511,600], it will repeat the last value up to the number of sweeps.
  4. psi0 = randomMPS(sites, states; linkdims=10) essentially uses a random circuit applied to the product state defined by states (though a bit smarter than that, but basically it is based off of applying random unitaries to the specified product state). It is probably more helpful to look at the individual tensors, for example with @show psi0[1] to show the data of first MPS tensor.
  5. Regarding the number of effective Hamiltonian matrix-vector multiplications, that is controlled with eigsolve_krylovdim (see the docs for the dmrg function: DMRG · ITensors.jl, which you can also access at the Julia prompt by executing julia> ? dmrg). By default it is set to 3, which means a maximum of 2 matrix-vector multiplications (since the initial vector is included as part of the Krylov space). Unfortunately we don’t have an easy way to save or print how many matvecs were actually performed in practice. That information is available from KrylovKit.eigsolve in the numops field of the ConvergenceInfo object it outputs (see Eigenvalue problems · KrylovKit.jl), and it would be relatively easy to access that information through ITensors observer system (Observer System for DMRG · ITensors.jl), but we would need to modify the dmrg code to make that work.

By default dmrg from ITensors.jl uses 2-site DMRG, and it is a bit hard-coded to that right now though it would be pretty easy to allow users to use 1-site DMRG since the infrastructure is in place to do that internally.

ITensorTDVP.jl is a package that generalizes the ITensors.dmrg function to provide a common framework for doing MPS-based DMRG, TDVP, linear equation solving, DMRG-X, variational fitting to perform MPO*MPS contraction, etc. You can use 1-site DMRG in that package by using ITensorTDVP.dmrg with nsite=1, though we don’t have a good subspace expansion built into that package yet so you would likely need to use 2-site DMRG and then switch to 1-site.

There are open PRs for subspace expansion: Pull requests · ITensor/ITensorTDVP.jl · GitHub but those got stalled since development efforts have been focused on the ITensorNetworks.jl package, where we are developing a general tree tensor network solver framework that includes TTN-based DMRG, TDVP, linear equation solving, DMRG-X, variational fitting to perform TTN*TTN contractions etc., and we are actively working on including state-of-the-art subspace expansion to perform 1-site versions of those algorithms on any TTN structure.

Regarding issues around the Krylov solver for that starting state, that does seem pretty strange. We are using Jutho Haegeman’s wonderful KrylovKit.jl for the Krylov eigensolver in our DMRG code. A suggestion for investigating that issue would be to try to reproduce the linear map (i.e. the DMRG effective Hamiltonian) that is being passed to the Krylov solver, in isolation from DMRG, and investigate why it might be hitting that issue in KrylovKit.eigsolve, and maybe raise an issue over at KrylovKit.jl.

I’d also like to make a shameless plug for the package ITensorGaussianMPS.jl, which implements an algorithm I developed with Steve White for constructing MPS of free fermion states, which you might find useful.

6

Thanks Matthew, for your detailed reply!

I was using iTensor to confirm the behavior of DMRG in some specific circumstance, which I have done, so I won’t be able to follow up with the Krylov solver issue. But the original code is a reproducer, if anyone is interested. The final calculation didn’t make the cut to be included in a figure, but I’ll acknowledge iTensor in the paper. Incidentally, it is on subspace expansion, so look out for it early next week (I intend to submit to the arxiv on Friday).

1 Like