anchor.diff:          ! Conserve symmetry (spatial and spin) in selecting a.
           imsa = (sys%basis%basis_fns(i)%Ms+3)/2
           ! We want
  -        !   Gamma_totsym = Gamma_<D|o|D_i^a>
  +        !   Gamma_totsym = Gamma_(<D|o|D_i^a>)
  +        !                = Gamma_(<D|op|a^†_a a_i D>)
           !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
           !                = Gamma_i^-1 Gamma_op Gamma_a
           ! Thus we require
  
positive.code:             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>
              !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
              !                = Gamma_i^-1 Gamma_op Gamma_a
              ! Thus we require
  
negative.code:         do i = 1, sys%nel
              imsa = (sys%basis%basis_fns(occ_list(i))%Ms+3)/2
              ! We want
              !   Gamma_totsym = Gamma_(<D|o|D_i^a>)
              !                = Gamma_(<D|op|a^†_a a_i D>)
              !                = Gamma_D* Gamma_op Gamma_i^-1 Gamma_a Gamma_D
              !                = Gamma_i^-1 Gamma_op Gamma_a
              ! Thus we require
  