Wavefunction projection into fixed particle number subsectors

Question

Hi Miles,

iTensor is amazing! I haven’t played with it since we were all in Sherbrooke and I’m blown away by the examples and how simple it is to get something useful up and running. I think my question arises from a lack of understanding of IQIndex and IQMPS.

I’m thinking about spinless fermions with nearest neighbor interactions, i.e.:

N = 8;
auto sites = Spinless(N,{"ConserveNf",true});
AutoMPO ampo(sites);
for(int j = 1; j < N;++j)
{
    ampo += -1.0,"Cdag",j,"C",j+1;
    ampo += -1.0,"Cdag",j+1,"C",j;
ampo += 0.765,"N",j+1,"N",j;
}
auto H = IQMPO(ampo);

After getting the ground state as per the https://github.com/ITensor/ITensor/blob/master/sample/iqdmrg.cc sample code:

auto state = InitState(sites);
 for(int i = 1; i <= N; ++i)
 {
     if(i%2 == 0) state.set(i,"Emp");
     else         state.set(i,"Occ");
}
auto psi = IQMPS(state);

auto sweeps = Sweeps(5);
sweeps.maxm() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 2;
sweeps.noise() = 1E-7,1E-8,0.0;

auto energy = dmrg(psi,H,sweeps,"Quiet");

I’m trying to figure out if it is possible to compute the entanglement entropy for “each local particle number block” i.e. I know that the reduced density matrix must commute with the particle number operator for the spatial subregion n, so ρA is block diagonal in n.

Is there a way to use IQIndex or 'Proj' to get the wavefunction projected into a fixed number of particles subsector in the subregion (corresponding to 1/2 the sites)?

Thanks in advance for any help!

Best,
Adrian.