apply MPO to MPS that changes quantum number

Hi Miles and Matt,
I would like to construct an MPO that has a definite change of quantum number (e.g. a sum of C^dag_i C^dag_j terms) and apply it to an MPS that has this quantum number. One way I can do it is change MPS to dense which removes all quantum numbers and makes subsequent calculations very slow. Is there a better way to do this currently?

Best,
Shengtao

Hi Shengtao, thanks for the question. If the MPO causes a definite change to the quantum number, then it should be possible to represent it with block-sparse or QN ITensors and not need to pass to dense tensors. Did you try making it using our QN system and did you run into an error?

Hi Miles, thanks for the quick reply! I was doing it in the standard way using MPO(), here is a minimal example of it:

using ITensors

N = 2
sites=siteinds("Electron",N,conserve_qns=true)
psi=productMPS(sites,["Emp" for i=1:N])
ampo=OpSum()
ampo+="Cdagup",1,"Cdagdn",2
cdagmpo=MPO(ampo,sites)

On the last line it reports: ERROR: No block found with QN equal to QN((“Nf”,-2,-1),(“Sz”,0)). Is there a way to make this MPO and apply it to psi without removing the QNs?

Hi Shengtao,

Unfortunately it is a known issue that the OpSum to MPO constructor doesn’t work when the MPO has nonzero flux (in principle it could, but we haven’t gotten around to implementing it).

An alternative it to make the MPO manually, which isn’t too hard for this case. For example:

using ITensors

function single_site_op(sites, which_op, j)
  left_ops = "Id"
  right_ops = "Id"
  if has_fermion_string(which_op, sites[j])
    left_ops = "F"
  end
  ops = [n < j ? left_ops : (n > j ? right_ops : which_op) for n in 1:length(sites)]
  return [op(ops[n], sites[n]) for n in 1:length(sites)]
end

s = siteinds("Electron", 4; conserve_qns=true)

Cdagup1 = single_site_op(s, "Cdagup", 1)
Cdagdn2 = single_site_op(s, "Cdagdn", 2)
Cdag = MPO([apply(o1, o2) for (o1, o2) in zip(Cdagup1, Cdagdn2)])
@show flux(Cdag)

However, please double check that this gives the correct MPO, I haven’t thoroughly tested it.

Hi Matt,
Thanks and this works! Now I try to build on this and make an MPO for a sum of two CdagupCdagdn terms:

using ITensors

function single_site_op(sites, which_op, j)
  left_ops = "Id"
  right_ops = "Id"
  if has_fermion_string(which_op, sites[j])
    left_ops = "F"
  end
  ops = [n < j ? left_ops : (n > j ? right_ops : which_op) for n in 1:length(sites)]
  return [op(ops[n], sites[n]) for n in 1:length(sites)]
end

s = siteinds("Electron", 4; conserve_qns=true)

Cdagup1 = single_site_op(s, "Cdagup", 1)
Cdagdn2 = single_site_op(s, "Cdagdn", 2)
Cdagup3 = single_site_op(s, "Cdagup", 3)
Cdagdn4 = single_site_op(s, "Cdagdn", 4)
Cdagpt1 = MPO([apply(o1, o2) for (o1, o2) in zip(Cdagup1, Cdagdn2)])
Cdagpt2 = MPO([apply(o1, o2) for (o1, o2) in zip(Cdagup3, Cdagdn4)])

Cdag=Cdagpt1+Cdagpt2

The two individual MPOs were created successfully. However on the last line it reports:

In `svd`, the left or right indices are empty (the indices of `A` are (((dim=16|id=880|"Electron,Site,n=1") <Out>
 1: QN(("Nf",-2,-1),("Sz",0)) => 1
 2: QN(("Nf",-1,-1),("Sz",-1)) => 2
 3: QN(("Nf",-1,-1),("Sz",1)) => 2
 4: QN(("Nf",0,-1),("Sz",-2)) => 1
 5: QN(("Nf",0,-1),("Sz",0)) => 4
 6: QN(("Nf",0,-1),("Sz",2)) => 1
 7: QN(("Nf",1,-1),("Sz",-1)) => 2
 8: QN(("Nf",1,-1),("Sz",1)) => 2
 9: QN(("Nf",2,-1),("Sz",0)) => 1,)), but the input indices are (Index{Vector{Pair{QN, Int64}}}[(dim=16|id=880|"Electron,Site,n=1") <Out>
 1: QN(("Nf",-2,-1),("Sz",0)) => 1
 2: QN(("Nf",-1,-1),("Sz",-1)) => 2
 3: QN(("Nf",-1,-1),("Sz",1)) => 2
 4: QN(("Nf",0,-1),("Sz",-2)) => 1
 5: QN(("Nf",0,-1),("Sz",0)) => 4
 6: QN(("Nf",0,-1),("Sz",2)) => 1
 7: QN(("Nf",1,-1),("Sz",-1)) => 2
 8: QN(("Nf",1,-1),("Sz",1)) => 2
 9: QN(("Nf",2,-1),("Sz",0)) => 1])). For now, this is not supported. You may have accidentally input the wrong indices.

If I change conserve_qn to false then the addition of MPO works. Am I doing this correctly?

Best,
Shengtao

Hi Shengtao,

This is caused by a (now fixed) bug in the MPO code that MPO functions like contraction and addition weren’t working if the MPOs didn’t have explicit link indices. If you update to ITensors 0.3.12 it should now work.

Here is a complete and slightly simplified code that should work with ITensors 0.3.12:

using ITensors

function op_mpo(sites, which_op, j)
  left_ops = "Id"
  right_ops = "Id"
  if has_fermion_string(which_op, sites[j])
    left_ops = "F"
  end
  ops = [n < j ? left_ops : (n > j ? right_ops : which_op) for n in 1:length(sites)]
  return MPO([op(ops[n], sites[n]) for n in 1:length(sites)])
end

s = siteinds("Electron", 4; conserve_qns=true)

Cdagup1 = op_mpo(s, "Cdagup", 1)
Cdagdn2 = op_mpo(s, "Cdagdn", 2)
Cdagup3 = op_mpo(s, "Cdagup", 3)
Cdagdn4 = op_mpo(s, "Cdagdn", 4)
Cdag12 = apply(Cdagup1, Cdagdn2)
Cdag34 = apply(Cdagup3, Cdagdn4)
Cdag = Cdag12 + Cdag34
@show flux(Cdag)

Hi Matt,
Thanks and glad to know the bug has been fixed. But it seems that I encounter another bug when trying to apply it to my really large system:

using ITensors

function op_mpo(sites, which_op, j)
  left_ops = "Id"
  right_ops = "Id"
  if has_fermion_string(which_op, sites[j])
    left_ops = "F"
  end
  ops = [n < j ? left_ops : (n > j ? right_ops : which_op) for n in 1:length(sites)]
  return MPO([op(ops[n], sites[n]) for n in 1:length(sites)])
end

N=260  #When N>250 last line doesn't work
s = siteinds("Electron", N; conserve_qns=true)
Cdagup = op_mpo(s, "Cdagup", 5)
Cdagdn = op_mpo(s, "Cdagdn", 10)
Cdagmpo = apply(Cdagup, Cdagdn)

When N>250 , ITensor reports the following error for the last line:

ERROR: BoundsError: attempt to access 0-element Vector{Pair{QN, Int64}} at index [1]

EDIT: Actually it’s when N>256 it will throw that error.

Hi Shengtao,

Sorry for the slow response. The issue you are seeing there is that the norm of the MPO diverges as you make the system larger (try printing norm(Cdagmpo)). This causes problems for the truncation we are doing in the MPO contraction code (the norm of the MPO gets collected on a single tensor which we are truncating, so the norm of that tensor starts to diverge and causes numerical issues). I’m working on a fix/workaround for that.

However, Miles recently fixed the original issue you were seeing (with constructing the MPO, so code like this should work now:

using ITensors

N = 2
sites = siteinds("Electron",N,conserve_qns=true)
psi = productMPS(sites,["Emp" for i=1:N])
os = OpSum()
os += "Cdagup", 1, "Cdagdn", 2
cdagmpo = MPO(os, sites)

Please try updating to the latest version of ITensors.jl and see if that works now.

Cheers,
Matt

Hi Matt,
Actually I’ve noticed that change and this works perfectly now. Thank you and Miles for the effort to improve this!

Best,
Shengtao

1 Like