anchor.diff:                  if (occ_list(iel) /= i) then
                       hmatel%r = hmatel%r &
                                   + get_two_body_int_mol_nonzero(coulomb_ints, i, occ_list(iel), a, occ_list(iel), basis_fns)
  -                    if (basis_fns(occ_list(iel))%Ms == basis_fns(i)%Ms) &
  -                        hmatel%r = hmatel%r &
  +                    if (basis_fns(occ_list(iel))%Ms == basis_fns(i)%Ms) then
  +                        if (sys%read_in%extra_exchange_integrals) then
  +                            hmatel%r = hmatel%r - get_two_body_exchange_pbc_int_nonzero(sys%read_in%additional_exchange_ints, &
  +                                                    i, occ_list(iel), occ_list(iel), a, basis_fns)
  +                        else
  +                            hmatel%r = hmatel%r &
                                       - get_two_body_int_mol_nonzero(coulomb_ints, i, occ_list(iel), occ_list(iel), a, basis_fns)
  +                        end if
  +                    end if
                   end if
               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
  
  