anchor.diff:          ! Conserve symmetry (spatial and spin) in selecting the fourth orbital.
           ! Ms_i + Ms_j = Ms_a + Ms_b (Ms_i = -1,+1)
           ! => Ms_b = Ms_i + Ms_j - Ms_a
  -        ims = (spin-sys%basis%basis_fns(a)%Ms+3)/2
  +        imsb = (spin-sys%basis%basis_fns(a)%Ms+3)/2
           ! (sym_i* x sym_j* x sym_a)* = sym_b
           ! (at least for Abelian point groups)
  -        isym = sys%read_in%sym_conj_ptr(sys%read_in, &
  +        isymb = sys%read_in%sym_conj_ptr(sys%read_in, &
                   sys%read_in%cross_product_sym_ptr(sys%read_in, sym, sys%basis%basis_fns(a)%sym))
   
  -        if (sys%read_in%pg_sym%nbasis_sym_spin(ims,isym) == 0) then
  +        if (sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb) == 0) then
               ! No orbitals with the correct symmetry.
               allowed_excitation = .false.
  -        else if (spin /= 0 .and. isym == sys%basis%basis_fns(a)%sym .and. &
  -                    sys%read_in%pg_sym%nbasis_sym_spin(ims,isym) == 1) then
  +        else if (spin /= 0 .and. isymb == sys%basis%basis_fns(a)%sym .and. &
  +                    sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb) == 1) then
               allowed_excitation = .false.
           else
               do
  -                ind = int(sys%read_in%pg_sym%nbasis_sym_spin(ims,isym)*get_rand_close_open(rng))+1
  -                b = sys%read_in%pg_sym%sym_spin_basis_fns(ind,ims,isym)
  +                ind = int(sys%read_in%pg_sym%nbasis_sym_spin(imsb,isymb)*get_rand_close_open(rng))+1
  +                b = sys%read_in%pg_sym%sym_spin_basis_fns(ind,imsb,isymb)
                   if (b /= a) exit
               end do
   
  
positive.code:             !                = 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) then
                  allowed_excitation = .true.
                  exit
              end if
  
negative.code:         integer, intent(out) :: i, a
          logical, intent(out) :: allowed_excitation
  
          integer :: imsa, isyma, ind
  
          ! Does this determinant have any possible single excitations?
          allowed_excitation = .false.
          do i = 1, sys%nel
              imsa = (sys%basis%basis_fns(occ_list(i))%Ms+3)/2
              ! In principle here we should have (Gamma_i* Gamma_op)*.  We'll assume Gamma_op*=Gamma_op
              ! We want
              !   Gamma_totsym = Gamma_<D|o|D_i^a>
  