anchor.diff:          logical, intent(out) :: allowed_excitation
   
   
  -        integer ::  ij_spin
  +        integer ::  ij_spin, ij_sym, imsb, isyma, isymb
         
  -        integer :: a, b, i, j, a_ind, b_ind, i_ind, j_ind
  +        integer :: i_ind_ref, j_ind_ref, a_ind_ref, b_ind_cdet
  +        integer :: a_ind_rev_cdet, b_ind_rev_ref
  +        integer :: i_ref, j_ref, a_ref, b_ref, i_cdet, j_cdet, a_cdet, b_cdet
   
           integer :: nex
           integer :: cdet_store(sys%nel)
           integer :: ref_store(sys%nel)
           integer :: ii, jj, t
  +
  +        logical :: found
   
           ! 1. Select single or double.
           if (get_rand_close_open(rng) < excit_gen_data%pattempt_single) then
  
positive.code:         ! Out:
          !    i, j: orbitals in determinant from which two electrons are excited.
          !        Note that i,j are ordered such that i<j.
          ! [review] JSS: ij_sym has vanished.
          !    ij_sym: symmetry conjugate of the irreducible representation spanned by the codensity
          !        \phi_i*\phi_j. (We assume that ij is going to be in the bra of the excitation.)
          !    ij_spin: spin label of the combined ij codensity.
          !        ij_spin = -2   i,j both down
          !                =  0   i up and j down or vice versa
  
negative.code: 
                  ! 2b. Select orbitals to excite from
                  
                  call choose_ij_ind(rng, sys, cs%occ_list, i_ind_ref, j_ind_ref, ij_spin)
  
                  ! At this point we pretend we're the reference, and fix up mapping ref's orbitals to cdet's orbitals later.
                  i_ref = cs%occ_list(i_ind_ref)
                  j_ref = cs%occ_list(j_ind_ref)
  
                  ! We now need to select the orbitals to excite into which we do with weighting:
                  ! p(ab|ij) = p(a|i) p(b|j) + p(a|j) p(b|i)
  