ITensor · tyler-bryson · Jun 8, 2018 · Jun 8, 2018 · Jun 21, 2018
diff --git a/itensor/fermion.h b/itensor/fermion.h
@@ -0,0 +1,137 @@
+//
+// Distributed under the ITensor Library License, Version 1.2
+//    (See accompanying LICENSE file.)
+//
+#ifndef __ITENSOR_FERMION_H
+#define __ITENSOR_FERMION_H
+#include "itensor/iqindex.h"
+#include "itensor/tensor/permutation.h"
+
+namespace itensor {
+
+// from two indexsets (possibly) sharing some indices, produce the permutations
+// of their labels (0-indexed) which would produce 'standard' orientation e.g.
+// if:      Ais   = J,K,I               Bis   = P,I,J,Q
+// then:    permA = (0,1),(1,0),(2,2)   permB = (0,2),(1,1),(2,0),(3,3)
+template <class index_type>
+void inline
+unriffle(IndexSetT<index_type> const& Ais,
+         IndexSetT<index_type> const& Bis,
+         Permutation & permA,
+         Permutation & permB)
+    {
+    auto rA = Ais.r();
+    auto rB = Bis.r();
+    int   n = 0;
+    int   k = 0;
+    permA   = Permutation(rA,-1);   //reinitialize the permutations
+    permB   = Permutation(rB,-1);   //with default values -1
+
+    for(auto i:range(rA))           //ascending i= 0,...,rA-1
+        {
+        for(auto j: range(rB))      //ascending j=0,...,rB-1
+            {
+            if(Ais[i]==Bis[j])
+                {
+                permA[i] = n;
+                permB[j] = n;
+                n++;
+                continue;
+                }
+            }
+        }
+
+    for(auto i: range(rA))
+        {
+        if(permA.dest(i)==-1)
+            {
+            permA[i]=k;
+            k++;
+            }
+        else
+            {
+            permA[i]=permA[i]+rA-n;
+            }
+        }
+
+    for(auto j:range(rB))
+        {
+        if(permB.dest(j)==-1)
+            {
+            permB[j]=n;
+            n++;
+            }
+        }
+    }
+template void unriffle(IQIndexSet const& Ais,
+                       IQIndexSet const& Bis,
+                       Permutation & permA,
+                       Permutation & permB);
+
+
+// Given a permutation P and an ordered subset of P's domain R,
+// returns the number of inversions in a sort of the image
+// of P precomposed with R:
+//      r0 r1 ... rn ---> P[r0] P[r1] ... P[rn]
+// When R is ordered, returns number of crossings in P restricted
+// to the domain R.
+// mod 2 this is the signature of P|R.
+int inline
+count_filtered_swaps(Permutation const& P, IntArray R)	//ask Miles, why fail for &R ?
+    {
+    auto rR = R.size();
+    for(auto i:range(rR))
+        {
+        R[i]=P[R[i]];
+        }
+
+	int swaps = 0;
+    auto Y = IntArray(P.size(),0);
+    for(int i=rR-1; i>=0; --i)
+        {
+        auto x = R[i];
+        for(auto j: range(x))
+            {
+            if(Y[j] != 0){swaps++;}
+            }
+        Y[x]=1;
+        }
+
+    return swaps;
+    }
+
+int inline
+total_swaps(Permutation const& P, IQIndexSet const I, IntArray const& block)
+    {
+    //probably add an assert I.size() == block.size() == P.size()
+    auto rI = I.r();
+    int num_swaps = 0;
+
+    //assumes ALL IQIndex have same total space of QN's
+    //use the first index to determine to determine where different
+    //'colors' of fermions can be.
+    auto i1 = I[0].qn(1);  //first index/sector of 0th IQIndex of B
+    for(auto c=1; c<=QNSize(); ++c)
+        {
+        if(not isActive(i1,c)) break;        //returns 1 until past dim of qn
+        if(not isFermionic(i1,c)) continue;  //skip non-fermionic
+        auto cib = IntArray();
+        for(auto i:range(rI))
+            {
+            //value of the quantum number of the c colored fermion
+            //within this block of the tensor
+            auto q = I[i].qn(block[i]+1)[c-1];
+            if(std::abs(q)%2==1)
+                {
+                cib.push_back(i);
+                }
+            }
+        num_swaps += count_filtered_swaps(P,cib);
+        }
+    return num_swaps;
+    }
+
+
+} //namespace itensor
+
+#endif
diff --git a/itensor/itdata/qdense.cc b/itensor/itdata/qdense.cc
@@ -11,6 +11,7 @@
 #include "itensor/itdata/qdense.h"
 #include "itensor/itdata/qutil.h"
 #include "itensor/util/print_macro.h"
+#include "itensor/fermion.h"
 
 using std::vector;
 using std::move;
@@ -409,6 +410,7 @@ add(PlusEQ<IQIndex> const& P,
         }
     }
 
+
 template<typename TA, typename TB>
 void
 doTask(PlusEQ<IQIndex> const& P,
@@ -459,10 +461,15 @@ doTask(Contract<IQIndex>& Con,
     STOP_TIMER(33)
     auto& C = *nd;
 
+    //--------------------------------------------------------------------//
+    Permutation pL,pR;
+    unriffle(Con.Lis,Con.Ris,pL,pR);
+    //--------------------------------------------------------------------//
+
     //Function to execute for each pair of
     //contracted blocks of A and B
     auto do_contract = 
-        [&Con,&Lind,&Rind,&Cind]
+        [&Con,&Lind,&Rind,&Cind,&pL,&pR]    // [&Con,&Lind,&Rind,&Cind]
         (DataRange<const VA> ablock, Labels const& Ablockind,
          DataRange<const VB> bblock, Labels const& Bblockind,
          DataRange<VC>       cblock, Labels const& Cblockind)
@@ -484,8 +491,33 @@ doTask(Contract<IQIndex>& Con,
 
         //Compute cref += aref*bref
         START_TIMER(2)
-        contract(aref,Lind,bref,Rind,cref,Cind,1.,1.);
+
+        //----------------------------------------------------------------//
+        Print(Global::debug1());
+        if(Global::debug1())
+            {
+            //compute number of fermion over fermion crossings
+            auto Aswaps = total_swaps(pL,Con.Lis,Ablockind);
+            auto Bswaps = total_swaps(pR,Con.Ris,Bblockind);
+            auto swaps  = Aswaps+Bswaps;
+            auto fermion_sign = std::pow(-1,swaps);
+            if(fermion_sign < 0){println("-1");}
+
+            contract(aref,Lind,bref,Rind,cref,Cind,fermion_sign,1.);
+            }
+        else
+            {
+            contract(aref,Lind,bref,Rind,cref,Cind,1.,1.);
+            }
+
+
+
+        //----------------------------------------------------------------//
+
+
+
         STOP_TIMER(2)
+
         };
 
     START_TIMER(20)

diff --git a/itensor/mps/sites/spinless.h b/itensor/mps/sites/spinless.h
@@ -27,7 +27,7 @@ class SpinlessSite
         auto q_occ = QN("Nf=",1);
         if(not conserve_Nf) q_occ = QN("Pf=",1);
         s = IQIndex{nameint("Spinless ",n),
-            Index(nameint("Emp ",n),1,Site),QN(),
+            Index(nameint("Emp ",n),1,Site),QN({0,-1}),
             Index(nameint("Occ ",n),1,Site),q_occ};
         }
 

diff --git a/itensor/tensor/permutation.h b/itensor/tensor/permutation.h
@@ -27,6 +27,9 @@ struct Permutation
 
     Permutation(size_type size);
 
+    //TYLER
+    Permutation(size_type size, size_type val);
+
     size_type
     size() const { return store_.size(); }
 
@@ -67,6 +70,13 @@ Permutation(size_type size)
         store_[n] = n;
     }
 
+//TYLER: for constructing permutations.
+inline Permutation::
+Permutation(size_type size, size_type val)
+  : store_(size,val)
+     {
+     //intentionally left blank
+     }
 
 Permutation inline
 inverse(Permutation const& P)

diff --git a/itensor/tensor/range.h b/itensor/tensor/range.h
@@ -14,8 +14,7 @@ namespace itensor {
 //Parent type for all ranges
 class RangeType { };
 
-template<typename index_type,
-         size_t start = 0>
+template<typename index_type, size_t start = 0>
 class RangeT;
 
 template<typename range_type>