anchor.diff:          !    occ_list: integer list of occupied spin-orbitals in the determinant.
           ! Out:
   
  -        use basis, only: basis_length, basis_fns
  +        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
  +        integer :: ims, isym, ind
   
           ! Does this determinant have any possible single excitations?
           allowed_excitation = .false.
  
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:             shift = 0
              na = nbasis/2
          end select
  
          do
              ! We assume that the user is not so crazy that he/she is
              ! running a calculation where there exists no virtual 
  