Computing correlation functions with mixed site types: unsure how to overload operator definitions

Hi,

I am currently working with a mixed site set consisting of spin-1/2 electrons and a bosonic spin-1/2. The implementation is as follows (copied from here, see also here), where the number of sites in the MPS is doubled and the electrons and spins are on neighboring sites.

using ITensors
using ITensors, ITensorMPS, ITensorInfiniteMPS

Nx=4
Ny=4
sites = Vector{Index}(undef,2*Nx*Ny)
for n=1:2*Nx*Ny
  nx=div(n - 1, Ny) + 1
  if iseven(nx) #This just places the spins on sites with even x
    sites[n] = Index([QN() => 2],"S=1/2,Site,n=$n")
  else
    sites[n] = siteind("Electron"; addtags="n=$n",conserve_qns=true,conserve_sz=false)
  end
end

Now, I wish to compute the charge-charge correlator, which in an ordinary spin-1/2 system is efficiently done with charge_corr=correlation_matrix(psi, "Ntot", "Ntot"). However, in the present case, “Ntot” is not defined for the site type “S=1/2”. My initial guess was to redefine the “S=1/2” site type by making a new file include("spinhalfkondo.jl"), identical to the one from ITensors save for this extra line at the very end:

function op(::OpName"Ntot", ::SiteType"S=1/2")
    return [
      0.0 0.0;
      0.0 0.0
    ]
  end

However, this solution does not work, throwing the error

ERROR: ArgumentError: Overload of "op" or "op!" functions not found for operator name "Ntot" and Index tags: ("S=1/2,Site,n=4",).

I’m not sure how to perform the requisite overload, given that I have constructed the site indices in an unusual manner. Any help would be appreciated.

For convenience, this code fragment generates a wavefunction compatible with the site set defined above.

  # Create starting state, for testing purposes
  state = map(CartesianIndices((Ny, 2*Nx))) do I
    y=I[1]
    x=I[2]
    if iseven(x) #i.e. if this is a spin index
      return iseven(div(x,2)) ⊻ iseven(y) ? "↓" : "↑"
    else #i.e. if this is a fermion index
      x_f = div(x,2)+1
      return (mod(x_f+y,filling)==0) ? "Up" : "0"
    end
  end
  display(state)

  psi = MPS(sites, state)

You need to have function ITensors.op( and it will work

Some other small comments:

  • state can be a bad variable name to use as it is also a function

  • Recall you can pass what sites you want in the correlation_matrix, so you could avoid this with annoying post processing with

    charge_corr=correlation_matrix(
        psi,
        "Ntot",
        "Ntot";
        sites = collect(1:length(psi))[hastags.(sites, "Electron")],
    )
    
  • Your state is a matrix, and when your state is converted inside MPS the code is doing

    [state(sites[j],state_[j]) for j=1:length(state_)]
    

    Just be sure that is what you intend

Thanks very much, the solution worked perfectly!

Just one quick question regarding the last point: I adapted that particular segment (with defining a state and then constructing the MPS) from the 2-D Hubbard DMRG example. Was there any concern with the implementation here? (I have benchmarked this code against exactly solvable things like a band insulator and things look OK, but if there was any potential problem that you had anticipated I would want to address it.)

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.