anchor.diff:  
   contains
   
  -    pure function calc_pgen_hub_k(ij_sym, ispin, f, unocc_alpha, unocc_beta) result(pgen)
  +    pure function calc_pgen_hub_k(ij_sym, f, unocc_alpha, unocc_beta) result(pgen)
   
           ! Calculate the generation probability of a given excitation for the
           ! Hubbard model in momentum space.  The Hubbard model is a special case
  
positive.code:         use symmetry, only: inv_sym, sym_table
  
          real(dp) :: pgen
          integer, intent(in) :: ij_sym, ispin
          integer(i0), intent(in) :: f(basis_length)
          integer, intent(in) :: unocc_alpha(nvirt_alpha), unocc_beta(nvirt_beta)
  
  
negative.code: 
          ! pgen then follows immediately.
  
          ! exciting from alpha, beta orbitals.
          do a = 1, nvirt_alpha
              b = 2*sym_table(ab_sym, unocc_beta(a)/2) ! b is a beta orbital.
              b_pos = bit_lookup(1,b)
              b_el = bit_lookup(2,b)
              if (btest(f(b_el), b_pos)) forbidden_excitations = forbidden_excitations + 1
          end do
          do a = 1, nvirt_beta
              b = 2*sym_table(ab_sym, (unocc_beta(a)+1)/2) + 1 ! b is an alpha orbital.
              b_pos = bit_lookup(1,b)
              b_el = bit_lookup(2,b)
              if (btest(f(b_el), b_pos)) forbidden_excitations = forbidden_excitations + 1
          end do
          mallow = nsites - nbeta
  
          ! N^2/4 is the number of ways of choosing i,j.
          ! An additional factor of two arises from the symmetry in the order of
          ! choosing a and b.
          pgen = 8.0_dp/(nel*nel*(mallow - forbidden_excitations))
  
      end function calc_pgen_hub_k
  
  