# How do we implement custom global symmetries (generalized parity)?

+1 vote

I'm trying to write a generalized parity for a double-wire model of Fermions that conserves global number and the parity of the number of fermions in a given wire. Namely the model proposed in https://arxiv.org/pdf/1705.01786.pdf (equations 1-3) with the g term set to 0. The model has the usual nearest-neighbor tunneling, a nearest-neighbor potential energy term, and a cooper pair exchance term between two coupled quanutm wires. To implement this, I've mapped the two 1-d chains into a single 1-d chain with odd numbered sites corresponding to chain a and even numbered sites corresponding to chain b. The parity that is conserved is the parity of the number of fermions on one of the chains, say the odd number of sites. So I need to keep track of that parity as well as enforce global number conservation. The number conservation works fine, but the global parity does not.

I've attempted to adapt the Fermion.h file to include this generalized parity. The adapting the fermion.h file is based on some of the comments here, http://www.itensor.org/support/2030/3-component-hubbard-model?show=2030#q2030, for a 3-component Fermi-hubbard model.

Here is my attempt to generalize the fermion.h file names have been changed to avoid an error that was being thrown about redefining things like Fermion and FermionSite.
#pragma once

 #include "itensor/mps/siteset.h"
#include "itensor/util/str.h"

namespace itensor {
class FermionSiteAlt;
using FermionOddSitePar = BasicSiteSet<FermionSiteAlt>;
class FermionSiteAlt
{
Index s;
public:

FermionSiteAlt(Index I) : s(I) { }
FermionSiteAlt(Args const& args = Args::global())
{
auto ts = TagSet("Site,Fermion");
auto n = 1;
if(args.defined("SiteNumber"))
{
n = args.getInt("SiteNumber");
}
auto conserveQNs = args.getBool("ConserveQNs",true);
auto conserve_Nf = args.getBool("ConserveNf",conserveQNs);
auto conserve_Pfodd = args.getBool("ConservePfodd",false);
auto oddevenupdown = args.getBool("OddEvenUpDown",false);
if(not conserveQNs)
{
s = Index(2,ts);
}
else if(not oddevenupdown) //usual case
{
if(conserve_Nf) //usual case
{
if(not conserve_Pfodd) //usual case
{
s = Index(QN({"Nf",0,-1}),1,
QN({"Nf",1,-1}),1,Out,ts);
}
else // conserve_Pfodd not usual case
{
if( n%2==1 ) {
s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
QN({"Nf",1,-1}),1,QN({"Pfodd",1,-2}),1,
Out,ts);
}
if( n%2==0 ) {
s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
QN({"Nf",1,-1}),1,QN({"Pfodd",0,-2}),1,
Out,ts);
}
}
}
else
{
s = Index(QN({"Pf",0,-2}),1,
QN({"Pf",1,-2}),1,Out,ts);
}
}
else
{
auto q_emp = QN({"Sz",0},{"Nf",0,-1});
QN q_occ;
if(n%2==1) q_occ = QN({"Sz",+1},{"Nf",1,-1});
else       q_occ = QN({"Sz",-1},{"Nf",1,-1});

s = Index(q_emp,1,
q_occ,1,Out,ts);
}
}

Index
index() const { return s; }

IndexVal
state(std::string const& state)
{
if(state == "Emp" || state == "0")
{
return s(1);
}
else
if(state == "Occ" || state == "1")
{
return s(2);
}
else
{
throw ITError("State " + state + " not recognized");
}
return IndexVal{};
}

ITensor
op(std::string const& opname,
Args const& args) const
{
auto sP = prime(s);

auto Emp  = s(1);
auto EmpP = sP(1);
auto Occ  = s(2);
auto OccP = sP(2);

auto Op = ITensor(dag(s),sP);

if(opname == "N" || opname == "n")
{
Op.set(Occ,OccP,1);
}
else
if(opname == "C")
{
Op.set(Occ,EmpP,1);
}
else
if(opname == "Cdag")
{
Op.set(Emp,OccP,1);
}
else
if(opname == "A")
{
Op.set(Occ,EmpP,1);
}
else
{
Op.set(Emp,OccP,1);
}
else
if(opname == "F" || opname == "FermiPhase")
{
Op.set(Emp,EmpP,1);
Op.set(Occ,OccP,-1);
}
else
if(opname == "projEmp")
{
Op.set(Emp,EmpP,1);
}
else
if(opname == "projOcc")
{
Op.set(Occ,OccP,1);
}
else
{
throw ITError("Operator \"" + opname + "\" name not recognized");
}

return Op;
}

//
// Deprecated, for backwards compatibility
//

FermionSiteAlt(int n, Args const& args = Args::global())
{
*this = FermionSiteAlt({args,"SiteNumber=",n});
}

};

//
// Deprecated, for backwards compatability
//

// Commented out
//// Rebels don't care about backwards compatability. Burn it down!
// using SpinlessSite = FermionSiteAlt;
// using Spinless = BasicSiteSet<SpinlessSite>;

} //namespace itensor


Continued in first comment below - character limit prevented posting both sets of code.

commented by (480 points)
Here is the code to actually compute the MPSs:

#include "itensor/all.h"
#include "itensor/FermionSiteOddSitesParity.h"
#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <map>
#include <random>
#include <chrono>
#include "itensor/util/print_macro.h"

using namespace itensor;

int main()
{

//
// Set up model parameters
//

int PhysicalLatticeLength = 5 ;
int numberOfParticles = 5;

int NumberFermionicSpecies = 2;
int L = PhysicalLatticeLength*NumberFermionicSpecies;

int Lp1by2m1 = (L+1)/2 - 1 ;

println("L is \n");
std::cout << L << " \n ";

auto minW = 1.0 ;
auto maxW = 1.0 ;
auto deltaW = 0.5 ;
int     NumberWIndices = floor((maxW - minW)/ deltaW + 1.01);

int NumberOfStates = 9; // 5

double EnergySpectrum[NumberWIndices][NumberOfStates];
double randphysicalonebodypot[PhysicalLatticeLength];
double randonebodypotential[L];
double overlapMatrix[NumberWIndices][NumberOfStates][NumberOfStates];

double globalParityArray[NumberWIndices][NumberOfStates];

int   globalNumberTemp;
std::vector<int> globalNumberVector;

// Build the space, conserving Fermionic particle number
auto psites = FermionOddSitePar(L,{"ConserveNf",true, "ConservePfodd", true});
// auto psites = Fermion(L,{"ConserveNf",true});

// Parameters for Iemini et al. Wire Model

double t = 1.0 ; // Tunneling sets energy & time scale
double U  = 0.5;
double G = 0.0;
// W will be swept over

//
// Build Hamiltonian and initial MPS state to stick into DMRG
//

for( auto WIndex: range(NumberWIndices) )
{    // start Delta Inex loop

auto W = minW + deltaW*WIndex;

println("About to begin the calculation for W =   \n" );
std::cout << W << " \n ";

// Create the operators for building the Hamiltonian
auto ampo = AutoMPO(psites);

// Nearest Neighbor Couplings
for( auto ii:range1(PhysicalLatticeLength-1))
{
// Tunneling for all for both wires species
ampo += -t, "Cdag", NumberFermionicSpecies*ii-1, "C", NumberFermionicSpecies*ii+1; // Wire a
ampo += -t, "Cdag", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii-1;      // Wire a
ampo += -t, "Cdag", NumberFermionicSpecies*ii-0, "C", NumberFermionicSpecies*ii+2; // Wire b
ampo += -t, "Cdag", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii-0; // Wire b

// Potential Term of adjascent sites
ampo += -2*U, "N", NumberFermionicSpecies*ii-1, "N", NumberFermionicSpecies*ii+1;
ampo += -2*U, "N", NumberFermionicSpecies*ii-0, "N", NumberFermionicSpecies*ii+2;

// Nearest-Neighbor Cooper Pair exchange between wires
ampo += W, "Cdag", NumberFermionicSpecies*ii-0, "Cdag", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii-1;
ampo += W, "Cdag", NumberFermionicSpecies*ii-1, "Cdag", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii-0;

}

// Single-site Terms
for( auto ii:range1(PhysicalLatticeLength))
{

// One-body term to set zero
ampo += 2*U, "N", NumberFermionicSpecies*ii-1; // na(j)
ampo += 2*U, "N", NumberFermionicSpecies*ii-0; // nb(j)

}

// Generate the Hamiltonian
auto H = toMPO(ampo);

// Generating an initial state over the space specified above

auto state = InitState(psites);
// Set the initial state to be an MPS over the set of sites call psites
//// With defined QNs from FermionSiteOddSitesParity

// Set QNs for total particle number and fixed-parity
for( int ll : range1(L))
{
state.set(ll,"0");
}
for( int ll : range1(numberOfParticles) )
{
state.set(2*ll-1,"1");
}

// Generate the initial MPS to actually use
auto psi0 = state;

//
// Set DMRG properties
//

for( auto ii : range1(L) ) if( ii%2==1 ) state.set(ii,"Occ");
// Set the number of DMRG sweeps
auto sweeps = Sweeps(8);
// Set the maximum bond dimension per sweep
// Setting it low initially allows the computation to be much faster
sweeps.maxdim() = 80, 80, 100, 150, 150,  250, 350, 400;
// Set the cutoff for the bond dimension
sweeps.cutoff() = 1E-6, 1E-6, 1E-8, 1E-8, 1E-10, 1E-10, 1E-10, 1E-11;
// Set the number of Davidson iterations per sweep
sweeps.niter() =  4, 4, 4, 3, 3, 3, 2, 2;
// Set the noise to add to the MPS itself each sweep
sweeps.noise() = 1E-5, 1E-5, 1E-6, 1E-7, 1E-8, 1E-10, 1E-10, 1E-11;

//
//  Run the DMRG
//
int nstates = NumberOfStates;
std::vector<MPS> states;
std::vector<double> energies;

//     int nn = 0 ;
auto [energy,psi] = dmrg(H,states,psi0,sweeps,"Silent");
//// Calculate the energy and state
//// Well, that was short
states.push_back(psi);
energies.push_back(energy);

// Can save ITensor objects to file, but not generic C data arrays - not great?
writeToFile("psig", psi);

EnergySpectrum[WIndex][0] = energy;

for(int nn = 1; nn < nstates; nn++){
auto [energye,psie] = dmrg(H,states,psi0,sweeps,{"Silent=", false, "Weight=", 5.0});

states.push_back(psie);
energies.push_back(energye);
EnergySpectrum[WIndex][nn] = energye;
}

//
//  Post-process the MPS for parity
//

for( int nn : range(NumberOfStates) )
{
psi = states[nn];
psi.noPrime();

// 'gauge' the MPS to site 1
//  // any 'position' between 1 and L, inclusive, would work here
psi.position(1);
psidag = dag(psi);
for( int kk =1; kk < L; ++kk){
auto kkSiteInd = siteIndex(psi, kk) ;
if(kk%2==1) psidag.prime(kkSiteInd) ;
}

psi.prime(lii_1);
auto Par = psi(1);
Par *= op(psites,"F",1) ; // op(psites,"Id",1) - 2*op(psites,"N", 1);
Par *= psidag(1);

for(int kk = 2; kk < L; ++kk)
{
if( kk%2==1 ){ // odd sites, calculate 1-2N
// Need to count the number of Fermions between sites ii and jj
//// including phase for Jordan-Wigner string
Par *= psi(kk);
Par *= op(psites,"F",kk); // this is the operator we need!
// op(psites,"Id",kk) - 2*op(psites,"N", kk);
// J-W string - need to add the "F" argument in between
//// to get proper phase with Fermions
Par *= psidag(kk);
std::cout << " Computed parity operator contribution for an odd site ";
} else { // even sites, use an identity operator
// Need to count the number of Fermions between sites ii and jj
//// including phase for Jordan-Wigner string
Par *= psi(kk);
// Par *= op(psites,"F",kk); // Dont need for commuting operators?
// J-W string - need to add the "F" argument in between
//// to get proper phase with Fermions
Par *= psidag(kk);
std::cout << " Computed parity operator contribution for an even site ";
} // else an even site
std::cout << " Computed parity operator contribution for site " << kk  << " \n " << " \n ";
} // over site kk

// Final site
// index linking jj to jj+1 and use an identity operator
Par *= prime(psi(L),lL);
// Like above for site ii, prime again at site jj (this time on the RHS)
// Par *= op(psites,"F",L);
Par *= psidag(L);
std::cout << " Computed parity operator contribution for site " << L  << " \n ";

std::cout << " About to print the operator Par" << " \n ";
Print(Par);
std::cout << " \n " << " Printed the operator Par" << " \n ";

double Parityresult = elt(Par); // use eltC if expecting complex
globalParityArray[WIndex][nn] = Parityresult;

} // nn state of interest loop

for( int nn : range(NumberOfStates) )
{
auto psinn=states[nn];
for(int mm : range(NumberOfStates))
{
auto psimm=states[mm] ;
overlapMatrix[WIndex][nn][mm]=inner(psinn,psimm);
} // over mm
} // over nn

println("Energies at this parameter set are \n" );
for( int nn : range(NumberOfStates) )
{
std::cout << EnergySpectrum[WIndex][nn]  << " \n ";
}

println("Global Parities of Eigenstates at this parameter set are \n" );
for( int nn : range(NumberOfStates) )
{
std::cout << globalParityArray[WIndex][nn]  << " \n ";
}

} // WIndex ending loop

string parityFileName = "ParityDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
std::ofstream paritySpectrumFile(parityFileName); // ("EnergySpectrumFile.txt");
paritySpectrumFile << NumberWIndices << NumberOfStates << "\n";
for (int i = 0; i < NumberWIndices; i++)
{
for (int n = 0; n < NumberOfStates; n++)
{
paritySpectrumFile << globalParityArray[i][n] << " ";
}
paritySpectrumFile << "\n";
}

string energyFileName = "EnergiesDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
std::ofstream EnergySpectrumFile(energyFileName); // ("EnergySpectrumFile.txt");
EnergySpectrumFile << NumberWIndices << NumberOfStates << "\n";
for (int i = 0; i < NumberWIndices; i++)
{
for (int n = 0; n < NumberOfStates; n++)
{
EnergySpectrumFile << EnergySpectrum[i][n] << " ";
}
EnergySpectrumFile << "\n";
}

string overlapFileName = "overlapMatrixDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
std::ofstream overlapFile(overlapFileName);
overlapFile << NumberWIndices*NumberOfStates << " " << L << "\n";
for (int i = 0; i < NumberWIndices; i++)
{
for (int nn = 0; nn < NumberOfStates; nn++)
{
for (int mm = 0; mm < NumberOfStates; mm++)
{
overlapFile << overlapMatrix[i][nn][mm] << " ";
} // j
overlapFile << "\n";
} // n
}    // i

return 0;
}

The energies printed are noticably out of order, and only some have a defined odd-site (single wire) parity.

Energies at this parameter set are:
-2.2028
-2.22071
-1.46725
-2.21413
-2.20282
-1.95795
-1.95795
-1.62416
-1.69117

Global Parity Measurements of Eigenstates are
-1
1
1
-0.997684
0.997674
0.999984
-0.999984
0.859696
-0.678753

Of course, the energies can be tweaked by playing with the term that is proportional to the overlap between the states, but it doesn't seem to be conserving the parity. Any insight would be greatly appreciated.
commented by (480 points)
edited
After moving
for( auto ii : range1(L) ) if( ii%2==1 ) state.set(ii,"Occ");
to just before assigning
auto psi0 = state; ,
the issue seems to be that the states simply aren't keeping track of the parity.

Without the noise term, the energies are wildly off, but the variant of parity is conserved. With the noise term on, the energies seem correct but the parity is not conserved. Using the randomMPS then redoing the symmetrization of the initial MPS labeled "psi0" doesn't do anything to help. The issue is I think in the modified fermion.h file (called FermionSiteOddSitesParity.h here, a not great filename but it's there).

The main problem is constructing a parity operator that only acts on half of the lattice non-trivially:
// conserves total particle number and the parity of the number of particles on an odd-indexed sites
{
if( n%2==1 ) {
s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
QN({"Nf",1,-1}),1,QN({"Pfodd",1,-2}),1,
Out,ts);
}
if( n%2==0 ) {
s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
QN({"Nf",1,-1}),1,QN({"Pfodd",0,-2}),1,
Out,ts);
}
}
commented by (480 points)
I attempted to make an asantz from the states that are of the form (1,0,1,0,...,1,0,0,...,0), (0,0,1,0,...,1,0,0,0,0), ... (0,..., 0,1,0,1,0) with particles only on odd sites with desired total number of particles with starting positions that go until the final odd site is filled. When I print this MPS:
- Sites not on the edges have things like this:
ITensor ord=3:
1: 1 QN({"Nf",0,-1})
(dim=4|id=450|"n=13,Site,Fermion") <Out>
1: 1 QN({"Nf",0,-1})
2: 1 QN({"Pfodd",0,-2})
3: 1 QN({"Nf",1,-1})
4: 1 QN({"Pfodd",1,-2})
1: 2 QN({"Nf",0,-1})
{norm=1.41 (QDense Real)}

depending on the site index, the dimension of the in/out of the first and last on that list vary in dimension, but always of the form 1: i QN({"Nf",0,-1}) .

So it seems like the issue is in the definition of the quantum numbers of my altered fermion.h file (called FermionSiteOddSitesParity.h here). Will try messing with that some.
commented by (480 points)
edited
Editing this final comment. I commented out some code that I *thought* was just a redundancy in forming the initial ansatz (different from the version posted above), and it seems to have worked. My code is slightly different from above. I'll organize all this Kitaev Wire stuff into a document once I get it working smoothly (and tested some).
commented by (70.1k points)
Hi Jared,
Thanks for being patient with my slow reply on this. So from your last comment above, did you find a solution now?

Looking at your custom site type file, it seems your solution is the one I would suggest: defining the same "Nf" quantum number for both types of sites, but for the "a" sites also defining a "Pf" parity quantum number. (Technically you don't have to define or include it at all for the "b" sites, because QN objects automatically "merge" their QN values and any that are missing are assumed to just have the value zero. You can try that out by adding two individual QN objects yourself if you want.)

My only suggestion looking at your code is that you can probably remove or simplify a lot of the if-else cases implementing user-defined options because I'm guessing you only want one of those cases for your system. The reason the library code is some complex is that we wanted to offer a lot of different options for people.

Miles
commented by (480 points)
Hi Miles,

Thanks for taking the time to look at this. I think the issue with the alternative site definition was that I was attempting to define two independent quantum numbers, but since the phase and number operators are both tied to the same variable of whether or not that site is occupied. Once I combined them and made it so the onsite index was binary, it seemed to fix it.

Removing the extra if-thens is a good idea for this system, I don't need this to get bigger and uglier with each new problem.

It does seem that the code only works when I set up the ansatz using the state language then generate an MPS from the state, which from the documentation I thought was antiquated (the InitState documentation gives a warning about the comments being for old code: http://itensor.org/docs.cgi?page=classes/initstate&vers=cppv3).

There was one more thing that was somewhat difficult for me to understand. From the fermion.h file, it seems that setting either to "occupied" or to "1" should give the same effect, but when I used both (as a redundancy), the code didn't preserve the abstracted parity, but I think when I did this for the Kitaev wire model, it worked fine with the redundancy.
commented by (70.1k points)
I see, so that page is only saying that a certain interface to InitState is antiquated. (Perhaps I should remove that page now too, because very few people may be even using an older version.) But the InitState system is necessary to use for systems with QN conservation (using its more up-to-date interface) because MPS are required to be initialized with well defined quantum numbers, so we can't just construct a 'generic' MPS for you just from the site indices without more information. Here is an example of a code using the newest interface for InitState:

(in the Julia version we've been able to just replace an InitState with an array of strings or numbers, and might pursue that simpler design in the future for the C++ version).

Regarding your other question, when making an InitState for Fermion sites, yes the strings "Occ" and "1" should have exactly the same meaning and effect. The string "occupied" is not one of the recognized inputs, though. I just updated the docs to reflect this: http://itensor.org/docs.cgi?vers=cppv3&page=classes/fermion

Miles
commented by (480 points)
Oh, thank you for clarifying that about INitState. I don't know if that page needs to be removed (or if anyone uses a version of ITensor before version 2 these days), but if you added a note saying what you did in the first paragraph, it would remove the confusion.

I think I did use "Occ" or at least I did in the sample code above, but I don't think there's any need to write the same thing twice like I had before since it just invites errors.

Thanks again, Miles.

Jared
commented by (70.1k points)
Good suggestions about the docs. I think what's happening there is that even though those pages aren't linked to or reachable by following links on the ITensor site, Google is indexing them anyway so they are showing up in search results and confusing people. I just updated that InitState one and removed some others that were out of date (and again aren't actually linked to on the site but can still be accessed directly if you find them in a search).

I agree it might be time to take down the ITensor v2 docs!
commented by (480 points)
Excellent. Thanks for keeping things updated!
commented by (70.1k points)
Also, if you don't mind, I'm collecting feedback from users here and there about their choice of language (C++ versus Julia). What are your current reasons for using the C++ version of ITensor versus the Julia one? Just more familiarity with C++ and/or already using that from before or another reason? Thanks -
commented by (480 points)
I don't mind at all. I'm not very familiar at all with Julia, but more or less assumed C++ would simply compute faster (of course, even if that's true at a fundamental level, user-written code is likely the bigger factor). There two other minor reasons, like C(++) is a more general purpose language, and since I'm uncertain on the probability of a post-doc, having some background in C is a more profitable bet. A more unique-to-me reason is that I work part-time with an experimental rock physics group rather than TA in my department, and use OpenFOAM some for that - OpenFOAM uses C with C++ comments,  so using it for ITensor allows me to consolidate my syntax frustrations.
commented by (70.1k points)
Thanks for that feedback! Yes it definitely makes sense given all those reasons. I'm thinking of adding a small FAQ to the ITensor site to help people weigh this choice.

Just for your info, Julia is just as fast of a language as C++, even for very low-level code. This is because Julia is a compiled language too, and all the types of variables in Julia can be known to the compiler, so just like in C++. There are even BLAS libraries now being written in pure Julia that rival Fortran BLAS! That being said, because of its greater flexibility as a language, it is easier to accidentally write slow code in Julia (slower than the comparable C++ I mean), though one can catch these slowdown by profiling the code.