Optimizing dmrg results for expected values of observables

DMRG and Numerical Methods
,
1

Hi,

consider I want dmrg to find a state with some specified value of a given observable. I think this should be possible in a very similar way to how dmrg looks for excited states, but can’t quite work it out.

Namely, I have a problem where my initial state has S_z = 0 and I know that the ground state is a singlet, S=0. There is a parameter regime where the triplet state S=1, S_z=0 is very close in energy, and thus dmrg converges into a state which is a mixture of the two, with \langle S^2 \rangle not exactly 0 (as expected for a singlet), but somewhat larger.
I would like to use my a priori physical knowledge and tell dmrg to look for the singlet states. Basically, I want to optimize \psi for

\langle \psi \vert H \vert \psi \rangle + \lambda ( \langle \psi \vert \hat{S}^2 \vert \psi \rangle - S(S+1) ),

where \lambda is some weight I would set empirically (similar to the excited state overlap weights), and S is the value of spin I want the final state to have.

One possibility would be to sum the three MPOs into

\hat{H} + \lambda \big( \hat{S}^2 - \hat{I} (S(S+1)) \big)

and run dmrg with this operator. How would one approach this?

Another option would require hacking the excited state calculation code, to make it compute the overlap in the second bracket at every step.

If possible, this approach would allow for optimisation of \psi for any a priori known observable, which might be generally useful, if it does not turn out to be computationally costly.
What would be the best approach?

Thanks, Luka

2

Thanks for a really interesting question. Yes this can be done and is not hard to set up with ITensor. There are two ingredients:

  1. our dmrg function has an interface where it can take multiple MPOs. These MPOs are understood by our code to be summed together. However, instead of actually summing the MPOs together, the code treats them as if they were summed and only performs the summation much later inside of the DMRG optimization step which is highly efficient.
  2. you will need to have an MPO form of the \mathbf{S}^2 operator which is not too hard to make

One option to make the MPO for \mathbf{S}^2 is to just expand the square and sum all of the N^2 or so different terms similar to how you can write out a sum using our OpSum/AutoMPO system for making Hamiltonians. However this will become inefficient for larger systems.

To go to larger systems, one can efficiently make the \mathbf{S}^2 MPO “by hand” using some common patterns for making MPOs, thinking of the MPO tensor as operator-valued matrices. Please let me know if you want to pursue this more sophisticated option and want some help with that.

3

Thank you for the extremely quick reply.
The approach indeed works, and performs very well! I think it even slightly improves convergence.

I already had the S^2 MPO implemented, thanks for offering help. For future reference, my implementation is inspired by this question from the old forum: To measure total spin angular momentum S^2 on MPS.

I added the option to subtract a constant from the MPO and to multiply the whole thing by a constant. The result is the operator \lambda ( \hat{S^2} - \mu \hat{I}) in MPO form, where \lambda and \mu are user inputs.

One more question though; how would I do the same for the excited states? It seems that I can only pass a vector of MPSs (of already computed states) OR a vector of MPOs, but not both.
One could create additional MPOs from the already computed states with the outer product \vert \psi \rangle \langle \psi \vert , but that does not sound terribly efficient. I found something that looks promising in the LocalMPO class (itensor/mps/localmpo.h):

//Use an MPS instead of an MPO. Equivalent to using an MPO
//of the outer product |Psi><Psi| but much more efficient.
//
LocalMPO(MPS const& Psi, 
         Args const& args = Args::global());

but could not quite figure out how to use it. Could one make an MPO equivalent to \vert \psi \rangle \langle \psi \vert using this?

4

Glad you got it working!

To extend this to be usable with DMRG for excited states, you would need to create a new class which merges the functionality of the LocalMPO_MPS class which is defined in mps/localmpo_mps.h together with the LocalMPOSet class defined in mps/localmposet.h.

LocalMPOSet handles the case of multiple MPOs which are summed during the eigensolver core step of DMRG. LocalMPO_MPS does something similar, but holds a collection of MPS that get applied to the MPS being optimized times the Weight parameter.

There’s no reason there couldn’t be a new type that applies both a collection of MPOs and then a collection of MPS to define the “projected Hamiltonian” or “local Hamiltonian” used in DMRG, we just didn’t think of that case. Let me know if you run into any issues while trying it.

Finally, the code you are referencing in LocalMPO is something you shouldn’t need to change or copy. What that does is let a LocalMPO object represent the operator | \Psi \rangle \langle \Psi| for some state | \Psi \rangle without needing to actually square an MPS which would be inefficient. Inside of LocalMPO_MPS it holds the MPS by storing them in an LocalMPO using that feature of the LocalMPO class. In hindsight, a better design would have been to make a separate class for that feature.

5

Hi,

I’m new to the Julia version and couldn’t find an example where one creates a custom (QN-conserving) MPO by hand. I’m particularly interested in the total spin S^2 of a spin-1/2 system, which we know can be written as a compact MPO. I was looking to adapt your older answer for C++: To measure total spin angular momentum S^2 on MPS - ITensor Support Q&A. However, some built-in functions and index notations seem to have changed and I haven’t got it to work so far. I was wondering if you could help with adapting the code for the Julia version or point me to some documentation.

Thanks!
Shovan

6

Hi Shovan, Yes the modern julia version is quite different. Here is a julia version for S2 and S=1/2 sites:

using ITensors

function makeS2(sites)
    N = length(sites)
    S2 = MPO(sites)
    qns=[QN("Sz",0)=>3,QN("Sz",-2)=>1,QN("Sz",2)=>1]
    col= Index(qns,"Link,l=0")
    for n in 1:N
        row = dag(col)
        col = Index(qns,"Link,l=$n")
        s=sites[n]
     
        W = ITensor(row,col,dag(s),prime(s))
       
        W += op(s,"Id") * onehot(row=>1,col=>1)
        W += op(s,"Id") * onehot(row=>2,col=>2)

        W += op(s,"S2") * onehot(row=>2,col=>1) 

        W += 2*op(s,"Sz") *onehot(row=>2,col=>3) 
        W += op(s,"Id") * onehot(row=>3,col=>3) 
        W += op(s,"Sz") * onehot(row=>3,col=>1)

        W += op(s,"S+") * onehot(row=>2,col=>4) 
        W += op(s,"Id") * onehot(row=>4,col=>4) 
        W += op(s,"S-") * onehot(row=>4,col=>1) 

        W += op(s,"S-") * onehot(row=>2,col=>5) 
        W += op(s,"Id") * onehot(row=>5,col=>5) 
        W += op(s,"S+") * onehot(row=>5,col=>1)

        #
        #  Terminate the MPO by effectively capping with l/r vectors.
        #
        if n==1
            W*=onehot(dag(row)=>5)
        end
        if n==N
            W*=onehot(dag(col)=>1)
        end

        S2[n] = W
    end

    return S2
end

sites = siteinds("S=1/2", 5,conserve_qns=true)
S2=makeS2(sites)

Next time please start a new topic. The discourse system is telling me the topic is solved and I should not be replying :slight_smile:

Anyway I hope this helps.

2 Likes
7

Hey,
one more thing I learned since this question was opened and might be worth documenting; do you want to implement an algorithm for enforcing total spin (like what the original question is about) or simply measure it once?
If latter is the case, I find it is easier and sometimes computationally less expensive to simply calculate the spin-spin correlation matrix \langle S_i \cdot S_j \rangle. The sum over all of its elements gives \langle S^2 \rangle, while partial sums can be used to extract information about the total magnetic moments of certain subspaces.

1 Like
8

Thanks a lot! Couldn’t have asked for a better answer.

Will start a new topic next time.

1 Like
9

Thanks for the suggestion. Indeed, right now we only want to measure it, which can also be found from all the two-site correlations. However, as it turns out we were not computing these correlations.

If there is a way to enforce total spin quantum number and only minimize the energy within a given sector that would be even more useful. But as far as I’m aware ITensor doesn’t yet have such SU(2) sitesets?

10

Enforcing total spin is precisely what this post is about - see the first three entries. (At least in the cpp implementation) one simply provides the S^2 MPO as well as the Hamiltonian to the dmrg() call, and the code treats them as if they were summed together.
Take care though, in my experience some playing around with the relative weight and other parameters is necessary to get reliable results.

11

@JanReimers Just want to point out one typo. For n=1, the row should be 2 not 5.

@PavesicL Ah nice, that’s a good idea. Somehow I missed it on first reading. Is it also there in the Julia version?
In our case we actually want to find out what is the total spin in the ground state. I suppose one can still effectively minimize within the different total spin blocks by adding appropriate chemical potentials, but one is still searching in the full Hilbert space. So there’s probably no advantage over just finding the ground state without adding any S^2 term.