anchor.diff:          !    <D|H|D>, the diagonal Hamiltonian matrix element involving D for
           !    systems defined by integrals read in from an FCIDUMP file.
   
  -        use molecular_integrals, only: get_one_body_int_mol_nonzero, get_two_body_int_mol_nonzero
  +        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
   
           complex(p) :: hmatel
  
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:                 hmatel = get_one_body_int_mol(one_e_ints, one_e_ints_im, i, a, sys)
  
                  do iel = 1, sys%nel
                      if (occ_list(iel) /= i) then
                          hmatel = hmatel &
                              + get_two_body_int_mol(coulomb_ints, coulomb_ints_im, 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_complex(sys%read_in%additional_exchange_ints, &
                                          sys%read_in%additional_exchange_ints_imag, i, occ_list(iel), occ_list(iel), a, sys)
  
                          else
                              hmatel = hmatel &
                                  - get_two_body_int_mol(coulomb_ints, coulomb_ints_im, i, occ_list(iel), occ_list(iel), a, &
                                                      sys)
                          end if
                      end if
                  end do
              end associate
  
  