anchor.diff:                      im =  get_two_body_int_mol_nonzero(coulomb_ints_im, i, j, i, j, sys%basis%basis_fns)
                       hmatel = hmatel + cmplx(re, im, p)
                       if (sys%basis%basis_fns(i)%Ms == sys%basis%basis_fns(j)%Ms) then
  -                        re =  get_two_body_int_mol_nonzero(coulomb_ints, i, j, j, i, sys%basis%basis_fns)
  -                        im =  get_two_body_int_mol_nonzero(coulomb_ints_im, i, j, j, i, sys%basis%basis_fns)
  +                        if (sys%read_in%extra_exchange_integrals) then
  +                            re = get_two_body_exchange_pbc_int_nonzero(sys%read_in%additional_exchange_ints, &
  +                                        i, j, j, i, sys%basis%basis_fns)
  +                            im = get_two_body_exchange_pbc_int_nonzero(sys%read_in%additional_exchange_ints_imag, &
  +                                        i, j, j, i, sys%basis%basis_fns)
  +                        else
  +                            re =  get_two_body_int_mol_nonzero(coulomb_ints, i, j, j, i, sys%basis%basis_fns)
  +                            im =  get_two_body_int_mol_nonzero(coulomb_ints_im, i, j, j, i, sys%basis%basis_fns)
  +                        end if
                           hmatel = hmatel - cmplx(re, im, p)
                       end if
                   end do
  
positive.code:                     im = get_two_body_int_mol_nonzero(coulomb_ints_im, i, occ_list(iel), a, occ_list(iel), basis_fns)
                      hmatel%c = hmatel%c + cmplx(re, im, p)
                      if (basis_fns(occ_list(iel))%Ms == basis_fns(i)%Ms) then
                          re = get_two_body_int_mol_nonzero(coulomb_ints, i, occ_list(iel), occ_list(iel), a, basis_fns)
                          im = get_two_body_int_mol_nonzero(coulomb_ints_im, i, occ_list(iel), occ_list(iel), a, basis_fns)
                          hmatel%c = hmatel%c - cmplx(re, im, p)
                      end if
                  end if
  
negative.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, &
                                          get_two_body_exchange_pbc_int_nonzero
          use system, only: sys_t
          use hamiltonian_data, only: hmatel_t
  
  