about MixedSiteSet for electron and spin

ITensor Julia Questions
,
1

Dear there,
I want to use dmrg to solve a Hamiltonian which has couplied spin 1/2 operator and electron operator on same site.
For example: H= a \sum_n S_x(n)C^\dag_{\uparrow}(n)C_{\uparrow}(n) + b\sum_n S_z(n).

After some searching, I found I need use MixedSiteSet. But so far, I only find a previous example for a Holstein-type new site set , which is:
using Holstein = MixedSiteSet(ElectronSite,BosonSite);

And I saw in that answer, it mentions “the odd numbered sites will be electron sites and the even numbered sites will be boson sites”. But in my case, I need electron and spin on same site. I’m not sure if this kind of MixedSiteSet already exits in Itensor or how could I construct this?

Thank you so much!!

1 Like
2

Hi, thanks for the question. What version of ITensor are you planning to use, the Julia or C++ version? I ask because the bit of code you mentioned looks like C+ code. However, I would strongly recommend you use the Julia version for this kind of mixed site type Hamiltonian because it’s much easier to do there.

Secondly, when you say you “need electron and spin on the same site” it’s important to realize that the physical sites and MPS sites don’t have to be exactly the same as each other. Even though in your Hamiltonian there is a spin on the same site as an electron, in practice in the code you can always map each term of the Hamiltonian into a system where the spins are on different sites from the electrons. Then a term like \sum_n S^z_n remains a local, on-site term while a term like \sum_n S^x_n c^\dagger_{\uparrow n} c_{\uparrow n} becomes a nearest-neighbor term in your code.

4

Thanks for your reply! I use Julia! (Julia Version 1.7.2, ITensors v0.3.18). But I’m not very sure why we can map each term of H into a system where spin are on different sites from electron? like why it’s equivalent to the original H…? And then do we still need define a MixedSiteSet somehow?

Thank you so much!

5

As to why it’s possible to map one Hamiltonian to another, as long as the Hilbert or state space on which each version of the Hamiltonian acts is the same size and as long as both versions of the Hamiltonians map equivalent states into each other, then all properties (energy, expected values of appropriately mapped operators) will also be the same. Another way to see this is that the idea that the spin (in your model) is on the same site as the electron is just a sort of arbitrary mathematical concept: why couldn’t we just as well say the spin is next to the electron and they are grouped into a larger concept called a “site”? But the real mathematical reason is the first one I mentioned.

As an example of how to make your Hamiltonian simpler, you could use the following one:

H' = \sum_{j \in \text{odd}} a C^\dagger_{\uparrow j } C_{\uparrow j} S^x_{j+1} + \sum_{j \in text{even}} S^z_j

where I’ve put the electrons on sites 1,3,5,… and the spins on sites 2,4,6,…

This alternative Hamiltonian should have all the same relevant properties (appropriately defined) as the one you wrote but could be easier to code and more efficient to solve numerically since the sites are made of things we already have in ITensor (an “Electron” site followed by a “S=1/2” site and then another “Electron” site etc.)

Finally, regarding how to actually set up a calculation with a mixture of site types, we have a code example about this in the documentation here:
https://itensor.github.io/ITensors.jl/stable/examples/DMRG.html#DMRG-Calculation-with-Mixed-Local-Hilbert-Space-Types
Please take a close look at it and see if it gives you enough information (in your case if you use the H' I showed above, you would replace the “S=1/2” and “S=1” site types with “Electron” and “S=1/2”).

6

Sorry for late response! I just got a chance to carefully tried this way, and it does exactly what I want! Thank you so much for the detailed explanation!

1 Like
7

Really glad you got it to work!

8

I’ve one more quick question! Is there a way to set up some symmetry for electron and spin separately? I’m trying to use conserve_nf=true or conserve_qns=true for electron(even site) and conserve_sz=true for spin(odd site), and my sites definition is as in the code example:
sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N;)
I tried several ways this morning but always get error messages…
Thank you!

9

There should be a way to do that. What error messages did you get?

10

I get several different error messages. First when I try to add conserve for both electron and spin, I tried:
sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true,conserve_nf=true)
ERROR: LoadError: MethodError: no method matching space(::SiteType{S=1/2}; conserve_qns=true, conserve_nf=true)
which I feel is like a syntext error, saying I should mention two conserve condition in other places or format…

Then I try to only add conserve_qns=true for electron:
sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true)
ERROR: LoadError: In setindex!, the element (1, 2) of ITensor:
Dim 1: (dim=2|id=103|“S=1/2,Site,n=1”)’
1: QN(“Sz”,1) => 1
2: QN(“Sz”,-1) => 1
Dim 2: (dim=2|id=103|“S=1/2,Site,n=1”)
1: QN(“Sz”,1) => 1
2: QN(“Sz”,-1) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}
2×2
Block(2, 1)
[2:2, 1:1]
0.5
you are trying to set is in a block with flux QN(“Sz”,2), which is different from the flux QN(“Sz”,-2) of the other blocks of the ITensor. You may be trying to create an ITensor that does not have a well defined quantum number flux.

11

I see, so here are the answers:

  • for the first code, when you pass conserve_nf the "S=1/2" type doesn’t accept that keyword argument so it’s leading to a function calling error. (We should probably fix that!) But as long as you are ok with all of the QN’s being conserved, then passing conserve_qns=true should be enough.
  • for the second code, you pass conserve_qns=true which is great and then the error is certainly not coming from that line of code but from a different line. What is the line of code producing the error (based on the backtrace you get)?

Could you please share a minimal working example of your code that reproduces the error?

13

I see! I attach my code that reproduces the second error. The line produces the error is H = MPO(ampo, sites). And I’ve thought my model more carefully and realized the symmetry I need is particle number and spin SU(2) symmtry (which I think is equivalent to conserve_qns=true?) for electron site, and No symmtry for spin1/2 site. I’m not very sure, if I write sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true) , does conserve_qns=true only applys to electron or both?

A minimal Code:

using ITensors
using HDF5

	let
	  N = 200 
      Npart = 100
	  t = 1.0
      hz=5
      hx=5
	  error=1E-8

	  sites = siteinds(n->isodd(n) ? "S=1/2" : "Electron", N; conserve_qns=true)

	  ampo = OpSum()
	  # hopping term
	  for i in 2:2:(N - 2)
        ampo += -t, "Cdagup", i, "Cup", i + 2
        ampo += -t, "Cdagup", i + 2, "Cup", i
        ampo += -t, "Cdagdn", i, "Cdn", i + 2
        ampo += -t, "Cdagdn", i + 2, "Cdn", i
      end
      # pseudospin -hz\sum_i Sz(i) at odd sites
	  for i in 1:2:N-2
	    ampo += -hz, "Sz", i
	  end
      # coupling term  -hx\sum_i Sx(i)N_up(i+1)  
      for i in 1:2:N-1
        ampo += -hx, "Sx", i, "Nup", i+1
	    ampo += -hx, "Sx", i, "Ndn", i+1
	  end
	  

      println(ampo)
	  H = MPO(ampo, sites)
	  

	  state = ["Emp" for n in 1:N]
	  p = Npart
	  for i in N:-2:2
	    if p > i/2
	      println("Doubly occupying site $i")
	      state[i] = "UpDn"
	      p -= 2
	    elseif p > 0
	      println("Singly occupying site $i")
	      state[i] = (isodd(i÷2) ? "Up" : "Dn")
	      p -= 1
	    end
	  end
	  for i in N-1:-2:1
	    state[i]=(isodd(i÷2) ? "Up" : "Dn")
	  end
	  # Initialize wavefunction to be bond dimension 10 random MPS with number of particles the same as `state`
      psi0 = productMPS(sites, state)


	  # Start DMRG calculation: 
	  sweeps = Sweeps(50)
	  setmaxdim!(sweeps,100)
	  setcutoff!(sweeps, error)
	  @show sweeps
	  
	  energy, psi = dmrg(H, psi0, sweeps)
    
      f = h5open("hx$(hx)hz$(hz)_1E-8.h5","w")
      write(f,"PGS",psi)
      close(f)
	  println("\nGround State Energy = $energy")
	end
14

I see. So your Hamiltonian has terms that include the “Sx” operator which would not conserve total Sz of the spins. At a glance I would conclude that it’s therefore not true that the symmetry your model conserves includes the spin SU(2) symmetry of the spin sites. (Unless I have missed something and the “Sx” operators combine in some way to preserve the spin symmetry?)

15

Thanks for your reply! Sorry I may not describe my model clear enough, just to make it clearer, the spin operator in my Hamiltonian is not electron spin, they are some effective pseudospin. So only C^\dag and C are electron operators, while S_x and S_z are pseudospin operators. The picture is like, I have two chains, one electron chain and one pseudospin chain, then I add some site-to-site coupling between them.

And as you said, my pseudospin doens’t have SU(2) symmetry, but my electrons have SU(2). So what I mean is, my electron sites has usual number conservation and SU(2) symmetry, but my pseudospin sites don’t have any symmetry. But I’m not very sure how to set this in Itensor?

Thanks!

16

I see, now I understand the symmetry you want to conserve.

There isn’t a simple way to set up the case you want using the siteinds function, but here is a workaround you can use:

sites = [isodd(n) ? Index([QN("Sz",0) => 2];tags="S=1/2,Site,n=$n") : siteind("Electron";addtags="n=$n",conserve_qns=true) for n=1:N]

What it does is to just manually define the Index objects for the “S=1/2” sites (which is what siteinds does internally anyway, it’s just that the keyword options that are allowed don’t cover the case you want). Also here I am using a bit of a hack of turning off QN conservation for those sites by always setting the QN to zero, but right now a limitation of ITensor is that if you are using QNs all indices must have them. In the future we plan to let you have a mix of indices, some with QNs and some without.

Please see if that works and gives correct results – I would suggest doing tests like:

  • computing the density of every electron site and adding it up to check its total value
  • computing the spin of every electron site and checking the total similarly
  • making plots of local properties to check that they look reasonable
  • comparing your energies to exact calculations for small systems
1 Like