anchor.diff:                  do jel = iel+1, sys%nel
                       j = occ_list(jel)
                       hmatel = hmatel + get_two_body_int_mol_nonzero(coulomb_ints, i, j, i, j, sys%basis%basis_fns)
  -                    if (sys%basis%basis_fns(i)%Ms == sys%basis%basis_fns(j)%Ms) &
  -                                  hmatel = hmatel - get_two_body_int_mol_nonzero(coulomb_ints, i, j, j, i, sys%basis%basis_fns)
  +                    if (sys%basis%basis_fns(i)%Ms == sys%basis%basis_fns(j)%Ms) then
  +                        if (sys%read_in%extra_exchange_integrals) then
  +                            hmatel = hmatel - get_two_body_exchange_pbc_int_nonzero(sys%read_in%additional_exchange_ints, &
  +                                                    i, j, j, i, sys%basis%basis_fns)
  +                        else
  +                            hmatel = hmatel - get_two_body_int_mol_nonzero(coulomb_ints, i, j, j, i, sys%basis%basis_fns)
  +                        end if
  +                    end if
                   end do
               end do
           end associate
  
positive.code:         ! symmetry).  This is less safe that slater_condon1_mol but much faster
          ! as it allows symmetry checking to be skipped in the integral lookups.
  
          use molecular_integrals, only: get_one_body_int_mol_nonzero, get_two_body_int_mol_nonzero
          use system, only: sys_t
          use hamiltonian_data, only: hmatel_t
  
  
negative.code:                     if (occ_list(iel) /= i) &
                          hmatel = hmatel &
                              + get_two_body_int_mol_real(coulomb_ints, i, occ_list(iel), a, occ_list(iel), &
                                                      sys)
                          if (sys%read_in%extra_exchange_integrals) then
                              hmatel = hmatel - get_two_body_exchange_pbc_int_real(sys%read_in%additional_exchange_ints, &
                                                      i, occ_list(iel), occ_list(iel), a, sys)
                          else
                              hmatel = hmatel - get_two_body_int_mol_real(coulomb_ints, i, occ_list(iel), &
                                                  occ_list(iel), a, sys)
                          end if
                  end do
              end associate
  
  