anchor.diff:          integer :: ref_store(sys%nel)
           integer :: ii, jj, t
   
  -        logical :: found
  +        logical :: found, a_found
   
           ! 1. Select single or double.
           if (get_rand_close_open(rng) < excit_gen_data%pattempt_single) then
  
positive.code:                 ! [review] - JSS: is it really worth constructing this (O(N) time) rather than just going for a binary (or even
                  ! [review] - JSS: linear) search on the cumulative table directly?
                  ! [reply] - VAN: can you please clarify how you would do what you say?
                  allocate(ia_weights(1:sys%nvirt_beta), stat=ierr)
                  call check_allocate('ia_weights', sys%nvirt_beta, ierr)
                  call create_weighted_excitation_list_ptr(sys, i, 0, cdet%unocc_list_beta, sys%nvirt_beta, ia_weights, &
                                                           ia_weights_tot)
                  ! Use the alias method to select i with the appropriate probability
                  a_ind = select_weighted_value(rng, sys%nvirt_beta, ia_weights, ia_weights_tot)
                  a = cdet%unocc_list_beta(a_ind) 
              else
                  allocate(ia_weights(1:sys%nvirt_alpha), stat=ierr)
                  call check_allocate('ia_weights', sys%nvirt_alpha, ierr)
  
                  call create_weighted_excitation_list_ptr(sys, i, 0, cdet%unocc_list_alpha, sys%nvirt_alpha, ia_weights, &
                                                           ia_weights_tot)
                  ! Use the alias method to select i with the appropriate probability
                  a_ind = select_weighted_value(rng, sys%nvirt_alpha, ia_weights, ia_weights_tot)
                  a = cdet%unocc_list_alpha(a_ind) 
              end if
  
              ! Given i,j,a construct the weights of all possible b
              ! This requires that total symmetry and spin are conserved.
              ! The symmetry of b (isymb) is given by 
              ! (sym_i* x sym_j* x sym_a)* = sym_b
              ! (at least for Abelian point groups)
              isymb = sys%read_in%sym_conj_ptr(sys%read_in, &
                      sys%read_in%cross_product_sym_ptr(sys%read_in, ij_sym, sys%basis%basis_fns(a)%sym))
              ! Ms_i + Ms_j = Ms_a + Ms_b => Ms_b = Ms_i + Ms_j - Ms_a
              ! Ms_k is +1 if up and -1 if down but imsb is +2 if up and +1 if down, 
              ! therefore a conversion is necessary.
              imsb = (ij_spin-sys%basis%basis_fns(a)%Ms+3)/2
              
              if (sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb) > 0) then 
                  allocate(jb_weights(1:sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb)), stat=ierr)
                  call check_allocate('jb_weights', sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb), ierr)
                  call create_weighted_excitation_list_ptr(sys, j, a, sys%read_in%pg_sym%sym_spin_basis_fns(:,imsb,isymb), &
                                          sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb), jb_weights, jb_weights_tot)
              else
                  jb_weights_tot = 0.0_p
              end if
  
              ! Test whether at least one possible b given i,j,a exists.
              ! Note that we did not need a btest for orbital a because we only considered
              ! virtual orbitals there.
  
              if (jb_weights_tot > 0.0_p) then
                  ! Use the alias method to select b with the appropriate probability
                  b_ind = select_weighted_value(rng, sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb), jb_weights, jb_weights_tot)
                  b = sys%read_in%pg_sym%sym_spin_basis_fns(b_ind,imsb,isymb)
  
negative.code:             j = cs%occ_list(i)  ! The elec we're looking at
              if (sys%basis%basis_fns(j)%Ms == -1) then ! beta
                  nv = sys%nvirt_beta
                  if (nv > 0) then
                      call create_weighted_excitation_list_ptr(sys, j, 0, cs%virt_list_beta, nv, cs%ia_weights(:,i), &
                                                              cs%ia_weights_tot(i))
                      call generate_alias_tables(nv, cs%ia_weights(:,i), cs%ia_weights_tot(i), cs%ia_aliasU(:,i), &
                                                 cs%ia_aliasK(:,i))
                  end if
                  do bsym = sys%sym0_tot, sys%sym_max_tot
                      if (sys%read_in%pg_sym%nbasis_sym_spin(1,bsym) > 0) then
                          call create_weighted_excitation_list_ptr(sys, j, 0, sys%read_in%pg_sym%sym_spin_basis_fns(:,1,bsym), &
                              sys%read_in%pg_sym%nbasis_sym_spin(1,bsym), cs%jb_weights(:,bsym,i), cs%jb_weights_tot(bsym,i))
                          call generate_alias_tables(sys%read_in%pg_sym%nbasis_sym_spin(1,bsym), cs%jb_weights(:,bsym,i), &
                              cs%jb_weights_tot(bsym,i), cs%jb_aliasU(:,bsym,i), cs%jb_aliasK(:,bsym,i))
                      end if
                  end do
              else ! alpha
                  nv = sys%nvirt_alpha
                  if (nv > 0) then
                      call create_weighted_excitation_list_ptr(sys, j, 0, cs%virt_list_alpha, nv, cs%ia_weights(:,i), &
                                                               cs%ia_weights_tot(i))
                      call generate_alias_tables(nv, cs%ia_weights(:,i), cs%ia_weights_tot(i), cs%ia_aliasU(:,i), &
                                                 cs%ia_aliasK(:,i))
                  end if
                  do bsym = sys%sym0_tot, sys%sym_max_tot
                      if (sys%read_in%pg_sym%nbasis_sym_spin(2,bsym) > 0) then
                          call create_weighted_excitation_list_ptr(sys, j, 0, sys%read_in%pg_sym%sym_spin_basis_fns(:,2,bsym), &
                              sys%read_in%pg_sym%nbasis_sym_spin(2,bsym), cs%jb_weights(:,bsym,i), cs%jb_weights_tot(bsym,i))
                          call generate_alias_tables(sys%read_in%pg_sym%nbasis_sym_spin(2,bsym), cs%jb_weights(:,bsym,i), &
                              cs%jb_weights_tot(bsym,i), cs%jb_aliasU(:,bsym,i), cs%jb_aliasK(:,bsym,i))
                      end if
                  end do
              end if
          end do
  