I'm trying to measure expected values of Sx, Sy and Sz on each site for given Hamiltonian and graph. My code looks as follows:
#include "itensor/all.h"
using namespace itensor;
int
main(int argc, char* argv[])
{
int Nx = 8; // Nx can be understood as the size of an elementary cell
auto N = atoi(argv[1])*Nx;
auto sites = SpinOne(N,{"ConserveQNs=",false});
auto ampo = AutoMPO(sites);
// [Definitiion of my Hamiltonian]
auto H = toMPO(ampo);
auto state = InitState(sites);
for (int i = 1; i <= N; i++) {
state.set(i,"Up");
}
auto sweeps = Sweeps(10);
sweeps.maxm() = 50,100,200,300,400;
sweeps.cutoff() = 1E-10;
println(sweeps);
//
// Begin the DMRG calculation
// for the ground state
//
auto [en0,psi0] = dmrg(H,randomMPS(sites),sweeps,{"Quiet=",true});
//
// Print the final energies reported by DMRG
//
printfln("\n Ground State Energy = %.10f",en0);
//
// Measuring Sx, Sy & Sz
//
println("Ground state");
println("\nj | Sx | Sy | Sz | S_total |: ");
for( auto j : range1(N) ) {
// Re-gauge psi0 to get ready to measure at position j
psi0.position(j);
auto ket = psi0(j);
auto bra = dag(prime(ket,"Site"));
auto Sxjop = op(sites,"Sx",j);
auto Syjop = op(sites,"Sy",j);
auto Szjop = op(sites,"Sz",j);
//take an inner product
auto sxj = eltC(bra*Sxjop*ket);
auto syj = eltC(bra*Syjop*ket);
auto szj = eltC(bra*Szjop*ket);
println(j, " | ", sxj, " | ", syj, " | ", szj, " | ", sxj + syj + szj, " | ");
}
return 0;
}
When I run the above program I get:
Ground state
j | Sx | Sy | Sz |:
1 | (2.82726e-12,0) | (0,-2.64747e-28) | (-0.98518,0) |
[...]
However, if I only change the part:
[...]
auto Sxjop = op(sites,"Sx",j);
auto Syjop = op(sites,"Sy",j);
auto Szjop = op(sites,"Sz",j);
[...]
into:
[...]
auto Sxjop = op(sites,"Sx2",j);
auto Syjop = op(sites,"Sy2",j);
auto Szjop = op(sites,"Sz2",j);
[...]
I get:
Ground state
j | Sx^2 | Sy^2 | Sz^2 | S_total^2 |:
1 | (0.5,0) | (0.5,0) | (1,0) | (2,0) |
[...]
which makes sense, because I believe the square of total spin satisfies following rules:
$$
S{total}^2 = Sx^2 + Sy^2 + Sz^2
S_{total}^2 = S ( S + 1),
$$
so for the spin @@ S = 1 @@ I should get @@ S_{total}^2 = 2 @@, which holds.
But I have to measure the values Sx, Sy and Sz instead of their squares. Do you have any suggestions?