# How to measure the average value of the magnetic moment Mz in the z-axis direction using metts

Hello everyone, I have encountered a problem in learning the metts method, and I hope to get your help.
The original code is as follows:

``````using ITensors
using Printf

#=

This example code implements the minimally entangled typical thermal state (METTS).
- "Minimally entangled typical quantum states at finite temperature", Steven R. White,
Phys. Rev. Lett. 102, 190601 (2009)
and arxiv:0902.4475 (https://arxiv.org/abs/0902.4475)
- "Minimally entangled typical thermal state algorithms", E M Stoudenmire and Steven R White,
New Journal of Physics 12, 055026 (2010) https://doi.org/10.1088/1367-2630/12/5/055026

=#

function ITensors.op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
h =
1 / 2 * op("S+", s1) * op("S-", s2) +
1 / 2 * op("S-", s1) * op("S+", s2) +
op("Sz", s1) * op("Sz", s2)
return exp(τ * h)
end

"""
Given a Vector of numbers, returns
the average and the standard error
(= the width of distribution of the numbers)
"""
function avg_err(v::Vector)
N = length(v)
avg = v[1] / N
avg2 = v[1]^2 / N
for j in 2:N
avg += v[j] / N
avg2 += v[j]^2 / N
end
return avg, √((avg2 - avg^2) / N)
end

function main(; N=10, cutoff=1E-8, δτ=0.1, beta=2.0, NMETTS=3000, Nwarm=10)

# Make an array of 'site' indices
s = siteinds("S=1/2", N)

# Make gates (1,2),(2,3),(3,4),...
gates = ops([("expτSS", (n, n + 1), (τ=-δτ / 2,)) for n in 1:(N - 1)], s)
# Include gates in reverse order to complete Trotter formula
append!(gates, reverse(gates))

# Make y-rotation gates to use in METTS collapses
Ry_gates = ops([("Ry", n, (θ=π / 2,)) for n in 1:N], s)

# Arbitrary initial state
psi = randomMPS(s)

# Make H for measuring the energy
terms = OpSum()
for j in 1:(N - 1)
terms += 1 / 2, "S+", j, "S-", j + 1
terms += 1 / 2, "S-", j, "S+", j + 1
terms += "Sz", j, "Sz", j + 1
end
H = MPO(terms, s)

# Make τ_range and check δτ is commensurate
τ_range = δτ:δτ:(beta / 2)
if norm(length(τ_range) * δτ - beta / 2) > 1E-10
error("Time step δτ=\$δτ not commensurate with beta/2=\$(beta/2)")
end

energies = Float64[]

for step in 1:(Nwarm + NMETTS)
if step <= Nwarm
println("Making warmup METTS number \$step")
else
println("Making METTS number \$(step-Nwarm)")
end

# Do the time evolution by applying the gates
for τ in τ_range
psi = apply(gates, psi; cutoff)
normalize!(psi)
end

# Measure properties after >= Nwarm
if step > Nwarm
energy = inner(psi', H, psi)
push!(energies, energy)
@printf("  Energy of METTS %d = %.4f\n", step - Nwarm, energy)
a_E, err_E = avg_err(energies)
@printf(
"  Estimated Energy = %.4f +- %.4f  [%.4f,%.4f]\n",
a_E,
err_E,
a_E - err_E,
a_E + err_E
)
end

# Measure in X or Z basis on alternating steps
if step % 2 == 1
psi = apply(Ry_gates, psi)
samp = sample!(psi)
new_state = [samp[j] == 1 ? "X+" : "X-" for j in 1:N]
else
samp = sample!(psi)
new_state = [samp[j] == 1 ? "Z+" : "Z-" for j in 1:N]
end
psi = productMPS(s, new_state)
end

return nothing
end

``````

I found the average code for measuring energy:：

``````   energy = inner(psi', H, psi)
push!(energies, energy)
@printf("  Energy of METTS %d = %.4f\n", step - Nwarm, energy)
a_E, err_E = avg_err(energies)
``````

Now I want to measure the average value of magnetization in the z-axis direction Mz, how should I measure it, I hope to get your help.
Thank you.

You can measure the expected values of local operators such as S^z by using the `expect` function:
https://itensor.github.io/ITensors.jl/dev/examples/MPSandMPO.html#Expected-Value-of-Local-Operators