anchor.diff:          end do
           
           if (allowed_excitation) then
  -            allowed_excitation = .false.
  -            do while (.not.allowed_excitation)
  -                call find_ia_mol(f, occ_list, i, a, allowed_excitation)
  +            ! We could wrap around find_ia_mol, but it's more efficient to have
  +            ! a custom generator instead.  The cost of an extra few lines is worth
  +            ! the speed...
  +
  +            do
  +                ! Select an occupied orbital at random.
  +                i = occ_list(int(genrand_real2()*nel)+1)
  +                ! Conserve symmetry (spatial and spin) in selecting a.
  +                ims = (basis_fns(i)%Ms+3)/2 
  +                isym = basis_fns(i)%sym
  +                if (symunocc(ims, isym) /= 0) then
  +                    ! Found i.  Now find a...
  +                        ! It's cheaper to draw additional random numbers than
  +                        ! decode the full list of unoccupied orbitals,
  +                        ! especially as the number of basis functions is usually
  +                        ! much larger than the number of electrons.
  +                        do
  +                            ind = int(nbasis_sym_spin(ims,isym)*genrand_real2())+1
  +                            a = sym_spin_basis_fns(ind,ims,isym)
  +                            if (.not.btest(f(bit_lookup(2,a)), bit_lookup(1,a))) exit
  +                        end do
  +                    exit
  +                end if
               end do
  +
           end if
   
       end subroutine choose_ia_mol
  
positive.code: 
      subroutine choose_ab_mol(f, sym, spin, symunocc, a, b, allowed_excitation)
  
          use basis, only: basis_length
          use system, only: nel, sym0, nsym
          use point_group_symmetry, only: cross_product_pg_sym
  
          integer(i0), intent(in) :: f(basis_length)
          integer, intent(in) :: sym, spin, symunocc(:,sym0:)
          integer, intent(out) :: a, b
          logical, intent(out) :: allowed_excitation
  
          integer :: isyma, isymb, imsa, imsb, ims_min, ims_max
  
          ! Is there a possible (a,b) pair which conserves symmetry once the (i,j)
          ! pair has been selected?
          ! We don't renormalise the probability generators for such instances, so
          ! need to return a null excitation.
          select case(spin)
          case(-2)
              ims_min = 1
              ims_max = 1
          case(0)
              ims_min = 1
              ims_max = 2
          case(2)
              ims_min = 2
              ims_max = 2
          end select
  
          allowed_excitation = .false.
          allowed: do isyma = sym0, nsym+(sym0-1)
              isymb = cross_product_pg_sym(isyma, sym)
              do imsa = ims_min, ims_max
                  if (spin == 0) then
                      imsb = mod(imsa,2)+1
                  else
                      imsb = imsa
                  end if
                  if  (spin /= 0 .and. isyma == isymb) then
                      if (symunocc(imsa, isyma) > 1) then
                            allowed_excitation = .true.
                            exit allowed
                      end if
                  else if  (symunocc(imsa,isyma) > 0 .and. symunocc(imsb,isymb) > 0) then
                        allowed_excitation = .true.
                        exit allowed
                    end if
                end do
          end do allowed
  
          if (allowed_excitation) then
              allowed_excitation = .false.
              do while (.not.allowed_excitation)
                  call find_ab_mol(f, sym, spin, a, b, allowed_excitation)
              end do
          end if
  
      end subroutine choose_ab_mol
  
negative.code:         !    occ_list: integer list of occupied spin-orbitals in the determinant.
          ! Out:
  
          use basis, only: basis_length, basis_fns, bit_lookup
          use point_group_symmetry, only: nbasis_sym_spin, sym_spin_basis_fns
          use system, only: nel, sym0
          
          use dSFMT_interface, only: genrand_real2
  
          integer(i0), intent(in) :: f(basis_length)
          integer, intent(in) :: occ_list(nel), symunocc(:,sym0:)
          integer, intent(out) :: i, a
          logical, intent(out) :: allowed_excitation
  
          integer :: ims, isym, ind
  
          ! Does this determinant have any possible single excitations?
          allowed_excitation = .false.
  