Measure the properties of an MPS wavefunction

Here we will learn how to measure some properties of the Heisenberg model wavefunction by adding to the code from the basic DMRG example. Essentially we need three pieces to take an expectation value of some property: the operator corresponding to that property, the wavefunction, and a way to do the inner product of the operator with the wavefunction.

For example, consider trying to measure the z-component of the spin on some site j in the Heisenberg chain. This operator can be retrieved from the sites object that we defined last time, and can be saved as an ITensor:

auto Szjop = op(sites,"Sz",j);

Likewise we can obtain operators for the site j spin raising (S+) and lowering operators (S-), by calling op(sites,"Sp",j) and op(sites,"Sm",j) respectively. There is an efficient way to take the inner product of these local operators with the MPS wavefunction. To do this, we must first call the function

psi.position(j);

This step is very important. It "gauges" the MPS tensors to be orthogonal to the left and to the right of site j. Without this step, we would have to include all of the MPS tensors in our calculation to get the right answer. But with this gauging step we can just use the tensor at site j.

Getting the needed MPS site tensor is very simple:

auto ket = psi(j);

Here we have copied this tensor to a variable called ket to suggest its role as a Dirac ket in the expectation value we will compute.

To get the bra part of the Dirac bra-ket, we turn ket into a row vector and conjugate its entries. The way we do this is simple:

auto bra = dag(prime(ket,"Site"));

The prime function is what turns the ket into a row vector (because it will contract with the row index of our operator, which will be an index with a prime level of 1), and dag does the hermitian conjugation. The argument "Site" passed to prime tells it to prime only indices with the "Site" tag (physical indices).

Now we are ready to measure the expectation value of Sz by contracting the bra, operator, and ket. The call to elt(...) below converts the resulting scalar tensor tensor into a plain real number:

auto szj = elt(bra*Szjop*ket);

Let's do a slightly more complicated measurement that will involve an operator acting on on two sites j and j+1. The operator @@\vec{S}_j \cdot \vec{S}_{j+1}@@ will do just fine.
A useful way of rewriting this operator is

For example, one of the operators we'll need is

auto spm = 0.5*op(sites,"S+",j)*op(sites,"S-",j+1);

To represent the wavefunction for two sites, we simply contract together two site tensors:

auto bondket = psi(j)*psi(j+1);

The bondbra is made similarly to the bra from earlier:

auto bondbra = dag(prime(bondket,"Site"));

And expectation values are computed in the same way:

Real spm = elt(bondbra*spm*bondket);

Below is a complete code for measuring properties of the MPS wavefunction.
The code prints out the measured values of @@\langle S^z_j \rangle@@ and @@\langle {S}_j \cdot \vec{S}_{j+1} \rangle@@ .

 Download the sample code below

#include "itensor/all.h"

using namespace itensor;

int 
main()
    {
    int N = 100;

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

    auto ampo = AutoMPO(sites);
    for(int j = 1; j < N; ++j)
        {
        ampo += 0.5,"S+",j,"S-",j+1;
        ampo += 0.5,"S-",j,"S+",j+1;
        ampo +=     "Sz",j,"Sz",j+1;
        }
    auto H = toMPO(ampo);

    auto sweeps = Sweeps(5); //number of sweeps is 5
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;

    auto psi0 = randomMPS(sites);

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

    //
    // Measuring Sz
    //

    println("\nj Sz = ");
    for( auto j : range1(N) ) 
        {
        //re-gauge psi to get ready to measure at position j
        psi.position(j);

        auto ket = psi(j);
        auto bra = dag(prime(ket,"Site"));

        auto Szjop = op(sites,"Sz",j);

        //take an inner product 
        auto szj = elt(bra*Szjop*ket);
        printfln("%d %.12f",j,szj);
        }

    //
    // Measuring S.S
    //

    //Sum total S.S to check that it's 
    //equal to ground state energy
    Real totalSdS = 0.;

    println("\nj S.S = ");
    for( auto b : range1(N-1) ) 
        { 
        psi.position(b);

        auto bondket = psi(b)*psi(b+1);
        auto bondbra = dag(prime(bondket,"Site")); 

        auto zzop = op(sites,"Sz",b)*op(sites,"Sz",b+1); 
        auto pmop = 0.5*op(sites,"S+",b)*op(sites,"S-",b+1); 
        auto mpop = 0.5*op(sites,"S-",b)*op(sites,"S+",b+1); 

        auto zz = elt(bondbra*zzop*bondket);
        auto pm = elt(bondbra*pmop*bondket);
        auto mp = elt(bondbra*mpop*bondket);

        printfln("%d %.12f",b,zz+pm+mp);
        totalSdS += zz+pm+mp;
        }

    printfln("\nSum of S.S = %.12f",totalSdS);
    printfln("Ground state energy from DMRG = %.12f",energy);

    return 0;
    }