anchor.diff:          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_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:         ! 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_(<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
  