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

Question

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");
      ts.addTags("n="+str(n));
      }
    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
    if(opname == "Adag")
        {
        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.

Comments

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!

---

Excellent. Thanks for keeping things updated!

---

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 -

---

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.

---

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.

Answer 1

(Marking this question as answered - see discussion above.)