anchor.diff:                  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>
  +                !   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:         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* 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
  