How to measure an operator that may locally change the quantum number?

sysy · March 21, 2024, 12:40pm

Dear itensor,

Thank you and your team for developing itensor.

Recently, I have encountered some issues while using itensor. I would like to measure the following operator,

\hat{\mathcal{O}}=\prod_{m=1}^{N_s}{\hat{o}_m}

It is a product operator acting on all lattice sites, where the matrix form of the local operator is as follows (the lattice sites I am using are electronic lattice sites):

\begin{array}{c} \left< 00 \right|\\ \left< \uparrow 0 \right|\\ \left< 0\downarrow \right|\\ \left< \uparrow \downarrow \right|\\ \end{array}\left( \begin{matrix} 0& 0& 0& -1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 1& 0& 0& 0\\ \end{matrix} \right)

This operator as a whole does not violate particle number conservation, but each local operator, when acting on a matrix product state, converts empty states and doubly occupied states into each other.

Therefore, if we measure it directly, the following error will be reported:

Trying to set element:
Index: (dim=4|id=689|“n=1,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2}), Val: 2
Index: (dim=4|id=689|“n=1,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2}), Val: 2
Element flux is: QN({“Nf”,0,-1},{“Lz”,0})
ITensor flux is: QN({“Nf”,-2,-1},{“Lz”,-2})
From line 327, file /home/customer/ys/itensor/itensor/itensor_impl.h

In .set, cannot set element with flux different from ITensor flux

The measurement code I am using is as follows:

auto psi2=psi;
for (int j1 = 1; j1 <= N; j1++)
{  
    auto PHop = op(sites_new, "O", j1);
    auto psitem = (PHop * psi2(j1));
    psitem.noPrime();
    psi2.set(j1, psitem);
}
auto PH = inner(psi2, psi);

the ‘O’ operator in the above code is defined as:

if(opname == "O")
            {
            Op.set(Em,UDP,-1); 
            Op.set(Up,UpP,+1);
            Op.set(Dn,DnP,+1);
            Op.set(UD,EmP,+1);
            }

The electronic lattice code I am using is as follows:

auto conserveQNs = args.getBool("ConserveQNs");
        auto conserveNf = args.getBool("ConserveNf");
        auto conserveLz = args.getBool("ConserveLz");
        
        if(conserveQNs || conserveNf || conserveLz)
            {
            if(conserveNf && conserveLz)
                {
                 s = Index(QN({"Nf",0,-1}, {"Lz",0}),1,
                          QN({"Nf",1,-1}, {"Lz",n}),1,
                          QN({"Nf",1,-1}, {"Lz",n}),1,
                          QN({"Nf",2,-1}, {"Lz",2*n}),1,Out,ts);
                          cout<<"s="<<endl;
                          print(s);
                }
            else if(conserveNf) // don't conserve Sz
                {
                s = Index(QN({"Nf",0,-1}),1,
                          QN({"Nf",1,-1}),1,
                          QN({"Nf",1,-1}),1,
                          QN({"Nf",2,-1}),1,Out,ts);
                }
            
            }
        else
            {
            s = Index(4,ts);
            }
        }

Is there any other way for me to measure this physical quantity?
Thank you for your help.

miles · March 21, 2024, 1:18pm

This kind of operator, where the whole operator conserves a quantity but each local or “factor” operator does not, would be a good candidate for using an MPO.

Since your operator is a product, then in your case the MPO would just have a bond dimension or virtual dimension of 1, so you can actually leave out the bond indices most likely. (I say “most likely” because many parts of ITensor are not too picky about whether the bond indices are there or not.)

So the first thing I’d suggest you try is to make an MPO and then set each tensor one by one to the operators you want. Then you can use the function named inner to take the expected value of this MPO.

Please let us know if you try that and run into any errors or technical obstacles.

sysy · March 21, 2024, 2:38pm

I am not sure if I fully understand what you said. But I tried it. Here is my code, it still gives an error, just like before.

auto ampoPH=MPO(N);
        for (int j1 = 1; j1 <= N; j1++)
        {
            auto ampoop=AutoMPO(sites);
            ampoop += 1,"O",j1 ;
            auto PHmpo=toMPO(ampoop);
            ampoPH=nmultMPO(ampoPH,PHmpo,{"MaxDim",500,"Cutoff",1E-8});
      
        }

        auto PH = inner(psi,ampoPH, psi);

sysy · March 22, 2024, 12:29am

I still haven’t succeeded. can you talk a little bit more about how to construct an mpo consisting of all the product of operators, and what does it mean to set each tensor one by one to the operator I want? Can you give me some sample codse to show me how it works?

miles · March 25, 2024, 7:19am

You can set each tensor of an MPO M by just accessing each tensor using M(j) (read-only) or M.ref(j) when you want to write to or modify that tensor.

Here is a working example of setting each tensor of an MPO to be the operator “Sz”:

#include "itensor/all.h"
#include "itensor/util/print_macro.h"
using namespace itensor;

int main()
    {
    auto N = 4;
    auto s = SpinHalf(N);
    auto M = MPO(N);

    for(int j : range1(N))
        {
        M.ref(j) = op(s,"Sz",j);
        }

    for(int j : range1(N))
        {
        PrintData(M(j));
        }

    return 0;
    }

To make other tensors for each site besides the “Sz” operator, you can make individual ITensors yourself with the indices s(j) and prime(s(j)) for the operator acting on site j. Let me know if you have questions about how to do that, but mainly you can use the .set function on an ITensor to do this e.g.

    auto T = ITensor(s(1),prime(s(1)));
    T.set(1,1, -0.3);
    PrintData(T);

(By the way, this is all much easier to do in the Julia version of ITensor, so I would encourage you to use that version instead especially if you are a beginner or are new to the ITensor approach. If you are combining ITensor with another C++ code, then I would understand using the C++ version.)

sysy · March 25, 2024, 10:12am

Hello, I tried the following code,

    auto M = MPO(N);
        for (int j : range1(N))
        {
            auto T = ITensor(sites(j), prime(sites(j)));
            T.set(1, 4, -1);
            T.set(4, 1, 1);
            T.set(2, 2, 1);
            T.set(3, 3, 1);
            M.ref(j) = T;
        }
        for (int j : range1(N))
        {
            PrintData(M(j));
        }
        auto PH = inner(psi, M, psi);

and it still reported the error like this:
Mismatched QN Index from set 1 (dim=4|id=312|“n=1,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2})
Mismatched QN Index from set 2 (dim=4|id=312|“n=1,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2})
From line 870, file itensor.cc

Mismatched QN Index arrows

sysy · March 29, 2024, 2:46am

I’m not sure if I fully understand what you’re saying, but these codes still haven’t worked.

miles · March 29, 2024, 3:17am

I did not realize you were conserving quantum numbers, so I simplified the example I gave to be casual about the use of Index arrows. For QN-conserving tensors, it is important to set the arrow direction (covariance or contravariance) of indices in a particular way.

I believe if you change the construction of the tensor to be

auto T = ITensor(sites(j), dag(prime(sites(j))));

then it should work.

See this documentation page for more information about arrows.

sysy · March 29, 2024, 4:27am

sorry, it still reports another error:
Trying to set element:
Index: (dim=4|id=215|“n=1,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2}), Val: 4
Index: (dim=4|id=215|“n=1,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Lz”,0})
2: 1 QN({“Nf”,1,-1},{“Lz”,1})
3: 1 QN({“Nf”,1,-1},{“Lz”,1})
4: 1 QN({“Nf”,2,-1},{“Lz”,2}), Val: 1
Element flux is: QN({“Nf”,2,-1},{“Lz”,2})
ITensor flux is: QN({“Nf”,-2,-1},{“Lz”,-2})
From line 327, file /home/customer/ys/itensor/itensor/itensor_impl.h

In .set, cannot set element with flux different from ITensor flux

the code I am using is:

 auto M = MPO(N);
        for (int j : range1(N))
        {
            auto T = ITensor(sites(j), dag(prime(sites(j))));
            T.set(1, 4, -1);
             T.set(4, 1, 1);
             T.set(2, 2, 1);
             T.set(3, 3, 1); 
            M.ref(j) = T;
        }
        for (int j : range1(N))
        {
            PrintData(M(j));
        }
        auto PH = inner(psi, M, psi);

miles · March 29, 2024, 5:14am

When you get that error, it is saying that the operator you are trying to make is not consistent with the quantum numbers you are trying to keep track of. The only kinds of operators or tensors you are allowed to make in our system are those which change the quantum numbers by a “well-defined” amount. (More formally, transform under a single irreducible representation of the symmetry group.)

Does your operator have that property?

An example for a spin system would be that, if you were trying to conserve the total “Sz” quantum number, then you are allowed to define the operators S^z, S^+, and S^-. But you cannot make the operators S^x or S^y because these are a superposition of S^+ and S^-.

sysy · April 1, 2024, 2:27am

·The operator I want to define indeed has the property you mentioned. It should be a superposition of creating a spin-up and spin-down electron simultaneously and annihilating a spin-up and spin-down electron simultaneously. So in itensor, there is no way to measure this operator unless I turn off particle number conservation, right?

miles · April 1, 2024, 2:47am

Correct, you cannot measure it while keeping particle conservation on. In fact, in most cases where users want to measure operators like this, the expectation value would be zero when using a QN-conserving state.

Turning off conservation would definitely work and let you measure this.

Another option is to break the operator into “equivariant” terms (terms which have a well-defined quantum number changing behavior). So if you take the example of c + c^\dagger you can use linearity to measure it this way:

\bra{\psi}c + c^\dagger\ket{\psi} = \bra{\psi}c\ket{\psi} + \bra{\psi}c^\dagger\ket{\psi}

and just measure the two expectation values separately then add the results. (You can see why in the above example, if \ket{\psi} has a well-defined total particle number the answer will be zero though.)

Topic		Replies	Views
Measure an operator that is a product of many local operators? ITensor C++ Questions cpp , dmrg , mps , mpo , fermions	1	78	May 16, 2024
Expectation value of a local operator acting on a MPS ITensor C++ Questions	6	703	July 26, 2022
Applying an MPO to an MPS ITensor Julia Questions	1	262	November 13, 2023
How to give local measurement for MPS and sample the measurement bitstrings from the calculated probabilities. DMRG and Numerical Methods	15	1207	June 16, 2023
About operations on a state ITensor Julia Questions	4	86	October 13, 2024

How to measure an operator that may locally change the quantum number?

Related topics