anchor.diff:          integer, intent(out) :: i, a
           logical, intent(out) :: allowed_excitation
   
  -        integer :: ims, isym, ind
  +        integer :: imsa, isyma, ind
   
           ! Select an occupied orbital at random.
           i = occ_list(int(get_rand_close_open(rng)*sys%nel)+1)
   
           ! Conserve symmetry (spatial and spin) in selecting a.
  -        ims = (sys%basis%basis_fns(i)%Ms+3)/2
  +        imsa = (sys%basis%basis_fns(i)%Ms+3)/2
           ! We want
           !   Gamma_totsym = Gamma_<D|o|D_i^a>
           !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
           !                = Gamma_i^-1 Gamma_op Gamma_a
           ! Thus we require
           !   Gamma_a = Gamma_i Gamma_op
  -        isym = sys%read_in%cross_product_sym_ptr(sys%read_in, sys%basis%basis_fns(i)%sym,op_sym)
  -        ind = int(sys%read_in%pg_sym%nbasis_sym_spin(ims,isym)*get_rand_close_open(rng))+1
  -        if (sys%read_in%pg_sym%nbasis_sym_spin(ims,isym) == 0) then
  +        isyma = sys%read_in%cross_product_sym_ptr(sys%read_in, sys%basis%basis_fns(i)%sym,op_sym)
  +        ind = int(sys%read_in%pg_sym%nbasis_sym_spin(imsa,isyma)*get_rand_close_open(rng))+1
  +        if (sys%read_in%pg_sym%nbasis_sym_spin(imsa,isyma) == 0) then
               ! No orbitals with the correct symmetry.
               allowed_excitation = .false.
           else
  -            a = sys%read_in%pg_sym%sym_spin_basis_fns(ind,ims,isym)
  +            a = sys%read_in%pg_sym%sym_spin_basis_fns(ind,imsa,isyma)
               ! Is a already occupied in the determinant f?  If so, the excitation is
               ! not permitted.
               allowed_excitation = .not.btest(f(sys%basis%bit_lookup(2,a)), sys%basis%bit_lookup(1,a))
  
positive.code: 
          ni = sys%nel
          do i = 1, sys%nel
              ims = (sys%basis%basis_fns(occ_list(i))%Ms+3)/2
              ! We want
              !   Gamma_totsym = Gamma_<D|o|D_i^a>
              !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
              !                = Gamma_i^-1 Gamma_op Gamma_a
              ! Thus we require
              !   Gamma_a = Gamma_i Gamma_op
              isym = sys%read_in%cross_product_sym_ptr(sys%read_in, sys%basis%basis_fns(occ_list(i))%sym, op_sym)
              if (symunocc(ims,isym) == 0) ni = ni - 1
          end do
  
          ims = (sys%basis%basis_fns(a)%Ms+3)/2
          isym = sys%basis%basis_fns(a)%sym
          pgen = 1.0_p/(ni*symunocc(ims,isym))
  
      end function calc_pgen_single_mol
  
  
negative.code:                 ! Select an occupied orbital at random.
                  i = occ_list(int(get_rand_close_open(rng)*sys%nel)+1)
                  ! Conserve symmetry (spatial and spin) in selecting a.
                  imsa = (sys%basis%basis_fns(i)%Ms+3)/2
                  ! Assume op_sym is self-conjugate.
                  ! We want
                  !   Gamma_totsym = Gamma_<D|o|D_i^a>
  