Total spin operator for Electron sites system

DMRG and Numerical Methods
1

Hi there,

first of all thanks a lot for developing this great package for Julia, it’s really powerful and simple to use (from my limited experience)!

I’m considering a Hubbard model for a 1D ladder system (i.e. 2xN geometry), and want to constrain the solution to the singlet S^2=0 state, similar to what has been discussed here:
Optimizing dmrg results for expected values of observables
Essentially I want to implement the same operator as suggested in post 6
https://itensor.discourse.group/t/optimizing-dmrg-results-for-expected-values-of-observables/110/6
but I have a hard time understanding what is actually happening in that code.
I guess it is summing up all the local and cross-terms from expanding S^2 into the corresponding local operators?
If some kind person could give me some pointers how to adapt that code for Electron sites, I’d highly appreciate it.

As a follow up question, would you assume that adding such term to the Hamiltonian for a 2x8 system (max.bond dim. ~2000-4000) would still be a feasible calculation? (I guess I’ll just have to try).

Thanks so much,

2

Hi @steffen,

S^2 for S=1/2 particles

In the link you have posted, it is constructing the local tensors of the mpo directly. To understand why that works, consider that the a tensor in the bulk site i with link indices \alpha and \beta is:

W[i] = \begin{pmatrix} I & 0 & 0 & 0 & 0\\ S^2 & I & 2S_z & S^+ & S^-\\ S_z & 0 & I & 0 & 0\\ S^- & 0 & 0 & I & 0\\ S^+ & 0 & 0 & 0 & I \end{pmatrix}

Where the internal elements are matrices with indices i and i'. If you multiply W[i] * W[i+1] you obtain a new tensor where the second row would be:

S_1^2 + S_2^2 + 2 S_1^z S_2^z + S_1^+ S_2^- + S_1^-S_2^+

Which is exactly (\bm S_1 + \bm S_2)^2.
A recipe for constructing mpo’s this way is Finite State Automata (I’m sorry I haven’t found references for it). But more often one simply gets the hand of it practicing with simpler one, you may take a look here.

How to adapt it

Spin operators are already defined for electronic systems (SiteTypes Included with ITensor · ITensorMPS.jl) so I don’t think there is a need to adapt something.

Feasible calculation

It should make things easier, but haven’t never done it myself I’m not sure how it will turn out

1 Like
3

To continue Vince’s answer, check out the references here:

and Miles has nice notes below:

As a follow up question, would you assume that adding such term to the Hamiltonian for a 2x8 system (max.bond dim. ~2000-4000) would still be a feasible calculation? (I guess I’ll just have to try).

You can always pass MPOs individually (dmrg([H1,H2], psi0;...), and since this is a low “bond dimension” MPO it should not create too much more work compared your “energy” Hamiltonian

3 Likes
4

Adding to the answers above, a convenient way to make this operator is using our OpSum system. In some older posts, I had shown how to make it “manually” (by manipulating tensors) or I had said it might be too expensive to make using OpSum, but for system sizes that are reasonable (e.g. a few hundred sites) OpSum is fast enough that it’s not a problem.

The main thing needed is to define the “onsite” version of S^2_i = \vec{S}_i \cdot \vec{S}_i for the "Electron" site type as follows:

function ITensors.op(::OpName"S2",::SiteType"Electron")
  return 3/4*[1. 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
end

Once you put the function above somewhere in your code, you can then do the following to make the S^2 operator, meaning the following operator:

S^2 = (\sum_{j_1} \vec{S}_{j_2})\cdot(\sum_{j_2} \vec{S}_{j_2}) = \sum_j S^2_j + 2 \sum_{j_1 < j_2} \vec{S}_{j_1}\cdot\vec{S}_{j_2}

Here is the code:

  N = 50
  sites = siteinds("Electron",N)


  op = OpSum()
  for j1=1:N, j2=(j1+1):N
    op += 2,"Sz",j1,"Sz",j2
    op += "S+",j1,"S-",j2
    op += "S-",j1,"S+",j2
  end
  for j=1:N
    op += "S2",j
  end
  SdS = MPO(op,sites)
  @show linkdims(SdS)

I put the @show linkdims(SdS) at the end so you can verify that the bond dimension of the resulting MPO is always 5, independent of the system size.

Then as Ryan explained, you can “lazily” add the SdS MPO, say times some weighting factor, into your Hamiltonian by passing an array of MPOs into our DMRG code.

One precaution is that while this won’t change the cost of a single DMRG sweep very much (since the S^2 MPO has such a small bond dimension), you may need to do quite a bit more sweeps for the effect of this penalty term to really be “felt” by the DMRG code. The reason is something like that it can be a better ‘strategy’ for DMRG to ignore the total spin sector initially when lowering the energy, and only near the end when the state is quite precise does DMRG ‘notice’ that the total spin contribution is important to lower the energy all the way.

2 Likes
5

Thanks a lot for these great answers! Especially now with Miles’ simple code example it’s immediately clear (I got mentally blocked by the direct MPO construction).