anchor.diff:          ni = sys%nel
           do i = 1, sys%nel
               ims = (sys%basis%basis_fns(occ_list(i))%Ms+3)/2
  +            ! We want
  +            !   Gamma_totsym = Gamma_<D|o|D_i^a>
  +            !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
  +            !                = Gamma_i^-1 Gamma_op Gamma_a
  +            ! Thus we require
  +            !   Gamma_a = Gamma_i Gamma_op
               isym = sys%read_in%cross_product_sym_ptr(sys%read_in, sys%basis%basis_fns(occ_list(i))%sym, op_sym)
               if (symunocc(ims,isym) == 0) ni = ni - 1
           end do
  
positive.code: 
          ! In:
          !    sys: system object being studied.
          !    ij_sym: irreducible representation spanned by the (i,j) codensity.
          !        As symmetry is conserved in allowed excitations, this is also
          !        the irreducible representation spanned by the (a, b) codensity.
          !    a, b: unoccupied orbitals into which electrons are excited.
  
negative.code:                 ! Conserve symmetry (spatial and spin) in selecting a.
                  ims = (sys%basis%basis_fns(i)%Ms+3)/2
                  ! Assume op_sym is self-conjugate.
                  ! We want
                  !   Gamma_totsym = Gamma_<D|o|D_i^a>
                  !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
                  !                = Gamma_i^-1 Gamma_op Gamma_a
                  ! Thus we require
                  !   Gamma_a = Gamma_i Gamma_op
                  isym = sys%read_in%cross_product_sym_ptr(sys%read_in, sys%basis%basis_fns(i)%sym, op_sym)
                  if (symunocc(ims, isym) /= 0) then
                      ! Found i.  Now find a...
  