Two impurity Anderson model not converging, even for U = 0 (help a beginner)

DMRG and Numerical Methods
1

Hello,

I have a convergence problem that I think must be simple to spot since I am by no means a DMRG specialist. I’ve been trying to get some simple expectation values for the following model:

H = \sum_{\sigma} V_{g} (\hat{n}_{1\sigma} + \hat{n}_{6\sigma} ) + V \sum_{\sigma} (c^{\dagger}_{1\sigma} c_{2\sigma} + c^{\dagger}_{5\sigma} c_{6\sigma}+ H.c.) + U \hat{n}_{1\uparrow}\hat{n}_{1\downarrow} + U \hat{n}_{6\uparrow}\hat{n}_{6\downarrow} - t \sum_{i = 2}^{4}\sum_{\sigma} (c^{\dagger}_{i\sigma} c_{i+1\sigma}+ H.c.)

It has only six sites, so I can compare these results with exact diagonalization. Every time I run the code, I get a different value for my expectation values. Since it is a one-dimensional chain with only two interacting sites, I think DMRG should be able to handle it perfectly. Which means I did something wrong. Would anyone have any idea how to correct my error?

I am writing the code here since it’s minimal.

Kind regards,
Ana Luiza

function chi(; N, U, t, V, Vg, h)
  i1 = 1
  i2 = 4

  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 10, 20, 40, 800)
  setcutoff!(sweeps, 1e-6)
  setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 1e-8,1e-8,1e-8,1e-8,1e-8,1e-8,1e-8,1e-8, 1e-11)

  sites = siteinds("Electron", N+2)
  hop = fill(-t, (N-1))
  ampo = OpSum()
  for i in collect(1:N-1)
    ampo += hop[i], "Cdagup", i, "Cup", i+1
    ampo += hop[i], "Cdagup", i+1, "Cup", i
    ampo += hop[i], "Cdagdn", i, "Cdn", i+1
    ampo += hop[i], "Cdagdn", i+1, "Cdn", i
  end
    
  ampo += Vg, "Nup", N+1
  ampo += Vg, "Ndn", N+1
    
  ampo += Vg + h, "Nup", N+2
  ampo += Vg - h, "Ndn", N+2
        
  ampo += 0, "Nupdn", N+1
  ampo += U, "Nupdn", N+2
    
  ampo += V, "Cdagup", i1, "Cup", N+1
  ampo += V, "Cdagup", N+1, "Cup", i1
  ampo += V, "Cdagdn", i1, "Cdn", N+1
  ampo += V, "Cdagdn", N+1, "Cdn", i1

  ampo += V, "Cdagup", i2, "Cup", N+2
  ampo += V, "Cdagup", N+2, "Cup", i2
  ampo += V, "Cdagdn", i2, "Cdn", N+2
  ampo += V, "Cdagdn", N+2, "Cdn", i2
    
  H = MPO(ampo, sites)
  psi0 = randomMPS(sites, 10)

  energy, psi = dmrg(H, psi0, sweeps)
    
  bondket = psi[N+1]*psi[N+2]
  bondbra = dag(prime(bondket,"Site"))
    
  zzop = op(sites,"Sz",N+1)*op(sites,"Sz",N+2)
  Szcorr = scalar(bondbra*zzop*bondket)
  print(Szcorr)

  
end
2

Hi Ana,
Looking at your code, I can’t see anything majorly wrong. But I did notice this: you have a term ampo += 0, "Nupdn", N+1. Did you mean to include a term with a coefficient of zero in your Hamiltonian?

Also I notice that some details of your code seem different from the mathematical expression for H at the top, such as V_g acting on sites 1 and 6 at the top but seemingly different sites from that in your code.

As to why you get different answers every time, I don’t know. Sometimes that can happen because the Hamiltonian as programmed is not actually Hermitian. Other reasons can be that it is a very difficult Hamiltonian to converge and you may need to do more sweeps, although normally 20 sweeps is enough for a small system size.