anchor.diff:  
                   ! 2b. Select orbitals to excite from
                   
  -                call choose_ij_ind(rng, sys, cs%occ_list, i_ind, j_ind, ij_spin)
  +                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 = cs%occ_list(i_ind)
  -                j = cs%occ_list(j_ind)
  +                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)
  
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:         use alias, only: select_weighted_value_prec
          use dSFMT_interface, only: dSFMT_t, get_rand_close_open
          use hamiltonian_data, only: hmatel_t
          use read_in_symmetry, only: cross_product_basis_read_in
          use search, only: binary_search
  
          type(sys_t), intent(in) :: sys
          type(excit_gen_data_t), intent(in) :: excit_gen_data
  