Expectation value of a local operator acting on a MPS

Hello everyone,
I am a newbie so I want to apologize first if this question is trivial. I am trying to compute the expectation value of a local operator on a MPS, and I would like to specify each component of the operator. I tried the following code but it is not working:

    auto s = Index(1,"s");
    auto j=1;
    psi.position(j);
    auto ket = psi(j);
    auto bra = dag(prime(ket,"Site"));
    auto psidag = dag(prime(psi));
    auto M0 = ITensor(s,prime(s)); 
    M0.set(s(1),prime(s)(1), +0.5); 
    auto P0 = elt(bra * M0 * ket); 
    println("P0 = ",P0); 

A follow-up question would be how to define an arbitrary local operator of arbitrary dimension and how to set all of his components.
Thanks to all the community,
Francesco

Hi Francesco,
Good question – so the Index s you used to make the operator shouldn’t be a new Index, but should be one of the indices of the MPS whose properties you are measuring. Say you have a site set called sites, then what you want to do is define

auto s = sites(j);

then I believe the rest of your code should work as you expect. Please let me know if it doesn’t!

For the other part of your question, when you say you want to define an operator of arbitary dimensions, could you clarify what you are asking? If the operator acts on a certain site of an MPS, its indices can only have the same dimension of that site, otherwise the action of that operator would not be well defined.

Thank you for your fast reply!
So I edited the code following your suggestion and the previous error disappeared but now I get the following error: “In .set, cannot set element with flux different from ITensor flux”.
Just for context, I have a spin 1/2 chain defined as:

    int N = 100;
    auto sites = SpinHalf(N); 
    auto state = InitState(sites);
    auto psi0 = randomMPS(state)    

Later I apply dmrg

auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",false});

Now what I am trying to do is to apply a 2x2 matrix ,with 1/2 as element (1,2) and 0 on all the other elements, to a single site of the chain.
Thank you again for your help :slight_smile:

Yes, so the error you are getting here is because by default, the site indices returned by the SpinHalf function are “quantum number conserving” and then our system will only let you construct operators with well-defined quantum numbers in this mode.

To disable this feature, construct your site set as:

auto sites = SpinHalf(N,{"ConserveQNs=",false});

then you should be able to construct any operator you want.

Thank you Miles!
I edited the code following your suggestion.
So is it the correct way to construct an operator O to initialize it with:

auto O = ITensor(s,prime(s));

And set the non zero elements O_ij with :

O.set(s(i),prime(s)(j),x);

Because now I am getting the following error:
error: ‘Itensor’ was not declared in this scope
auto M0 = Itensor(s,prime(s));
And I don’t really know why
Thank you,
Francesco

Yes that is the right way – glad you are figuring it out.

The reason for the error you are getting is that the ITensor type is spelled with a capital “T” so ITensor not Itensor.

Ups,my bad, I haven’t noticed that, thank you a lot !