anchor.diff:  
           ! Conserve symmetry (spatial and spin) in selecting a.
           ims = (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
  
positive.code:         do i = 1, sys%nel
              ims = (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
              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
  
negative.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
  