I am trying to construct the following state of fermions:
which has a correlation matrix
using ITensors
s = siteinds("Fermion", 4);
ψ0 = productMPS(s, "Emp");
ampo24 = AutoMPO();
ampo24 += (1/sqrt(2), "Id", 1);
ampo24 += (1/sqrt(2), "c†", 2, "c†", 4);
ampo31 = AutoMPO();
ampo31 += (1/sqrt(2), "Id", 1);
ampo31 += (1/sqrt(2), "c†", 3, "c†", 1);
o24 = MPO(ampo24, s);
o31 = MPO(ampo31, s);
ψ = product(o31, product(o24, ψ0));
correlation_matrix(ψ, "c†", "c†")
This gives the correct result.
At the same time, I had expected (based on Matt’s response in [julia] TEBD gates involving fermions - ITensor Support Q&A) that it would be possible to do something like
o24_naive = op("Id", s[2])*op("Id", s[4]) + op("c†", s[2])*op("c†", s[4]);
o31_naive = op("Id", s[1])*op("Id", s[3]) + op("c†", s[3])*op("c†", s[1]);
ψ_naive = product(o31_naive, product(o24_naive, ψ0));
correlation_matrix(ψ_naive, "c†", "c†")
However, this returns
4×4 Matrix{Float64}:
0.0 0.0 -2.77556e-16 0.0
0.0 0.0 0.0 5.55112e-17
2.77556e-16 0.0 0.0 0.0
0.0 -5.55112e-17 0.0 0.0
Am I missing something crucial with this naive approach? eg, is it necessary to use enable_auto_fermion()
? Or is the best practice to always use OpSum/AutoMPO?