diff --git a/src/apps/chem/PNO.cpp b/src/apps/chem/PNO.cpp
index b34333c9922..812aa58ffc9 100644
--- a/src/apps/chem/PNO.cpp
+++ b/src/apps/chem/PNO.cpp
@@ -2620,14 +2620,21 @@ PNOPairs PNO::load_pnos(PNOPairs& pairs) const {
 	for (ElectronPairIterator it = pit(); it; ++it) {
 		vector_real_function_3d tmp;
 		std::string name = pairs.name(it);
-		try {
-			load_function(world, tmp, name);
-			msg << "Found PNOs for " << name << " with rank " << tmp.size()
-					<< "\n";
-		} catch (...) {
-			msg << "failed to find PNO " << name
-					<< " initialize as zero function\nn";
-		}
+        bool exists = archive::ParallelInputArchive<archive::BinaryFstreamInputArchive>::exists(world, name.c_str());
+        if (exists) {
+            try {
+                load_function(world, tmp, name);
+                msg << "Found PNOs for " << name << " with rank " << tmp.size()
+                    << "\n";
+            } catch (...) {
+                msg << "failed to find PNO " << name
+                    << " initialize as zero function\nn";
+            }
+        }
+        else {
+            msg << "failed to find PNO " << name
+                << " initialize as zero function\nn";
+        }
 		pno_ij[it.ij()] = tmp;
 	}
 	return pairs;
diff --git a/src/apps/chem/SCF.cc b/src/apps/chem/SCF.cc
index 70732fc4225..8c94cc290af 100644
--- a/src/apps/chem/SCF.cc
+++ b/src/apps/chem/SCF.cc
@@ -1718,7 +1718,7 @@ vecfuncT SCF::apply_potential(World & world, const tensorT & occ,
 	enl = 0.0;
 
 	// compute the local DFT potential for the MOs
-	if (xc.is_dft() && !(xc.hf_exchange_coefficient() == 1.0)) {
+	if (xc.is_dft()) {
 		START_TIMER(world);
 
 		XCOperator<double,3> xcoperator(world,this,ispin,param.dft_deriv());
@@ -1739,7 +1739,7 @@ vecfuncT SCF::apply_potential(World & world, const tensorT & occ,
 
 	END_TIMER(world, "V*psi");
 	print_meminfo(world.rank(), "V*psi");
-	if (xc.hf_exchange_coefficient()) {
+	if (xc.hf_exchange_coefficient() && xc.get_hf_yukawa_separation() > 1e-6) {
 		START_TIMER(world);
 		//            vecfuncT Kamo = apply_hf_exchange(world, occ, amo, amo);
 		Exchange<double,3> K=Exchange<double,3>(world,this,ispin).set_symmetric(true);
@@ -1753,8 +1753,27 @@ vecfuncT SCF::apply_potential(World & world, const tensorT & occ,
 			exchf *= 2.0;
 		gaxpy(world, 1.0, Vpsi, -xc.hf_exchange_coefficient(), Kamo);
 		Kamo.clear();
-		END_TIMER(world, "HF exchange");
+
 		exc = exchf * xc.hf_exchange_coefficient() + exc;
+
+        const double yukawa_separation_parameter = xc.get_hf_yukawa_separation();
+        if (yukawa_separation_parameter < 1e4) {
+            Exchange<double,3> K_short_range=
+                    Exchange<double,3>(world,this,ispin, yukawa_separation_parameter).set_symmetric(true);
+            vecfuncT K_short_range_amo=K_short_range(amo);
+            tensorT excv_short_range = inner(world, K_short_range_amo, amo);
+            double exchf_short_range = 0.0;
+            for (long i = 0; i < amo.size(); ++i) {
+                exchf_short_range += 0.5 * excv_short_range[i] * occ[i];
+            }
+            if (!xc.is_spin_polarized())
+                exchf_short_range *= 2.0;
+            gaxpy(world, 1.0, Vpsi, xc.hf_exchange_coefficient(), K_short_range_amo);
+            K_short_range_amo.clear();
+
+            exc = exchf_short_range * xc.hf_exchange_coefficient() + exc;
+        }
+        END_TIMER(world, "HF exchange");
 	}
 	// need to come back to this for psp - when is this used?
 	if (param.pure_ae()){
diff --git a/src/apps/chem/SCFOperators.cc b/src/apps/chem/SCFOperators.cc
index ba9d74cd093..a6d632f1c86 100644
--- a/src/apps/chem/SCFOperators.cc
+++ b/src/apps/chem/SCFOperators.cc
@@ -664,11 +664,13 @@ void XCOperator<T, NDIM>::prep_xc_args_response(const real_function_3d &dens_pt,
 
 /// ctor with a conventional calculation
 template<typename T, std::size_t NDIM>
-Exchange<T,NDIM>::Exchange(World& world, const SCF *calc, const int ispin) : impl(new Exchange<T,NDIM>::ExchangeImpl(world,calc,ispin)) {};
+Exchange<T,NDIM>::Exchange(World& world, const SCF *calc, const int ispin, const double mu)
+        : impl(new Exchange<T,NDIM>::ExchangeImpl(world,calc,ispin,mu)) {}
 
 /// ctor with a nemo calculation
 template<typename T, std::size_t NDIM>
-Exchange<T,NDIM>::Exchange(World& world, const Nemo *nemo, const int ispin) : impl(new Exchange<T,NDIM>::ExchangeImpl(world,nemo,ispin)) {};
+Exchange<T,NDIM>::Exchange(World& world, const Nemo *nemo, const int ispin, const double mu)
+        : impl(new Exchange<T,NDIM>::ExchangeImpl(world,nemo,ispin, mu)) {}
 
 /// apply the exchange operator on a vector of functions
 
@@ -682,12 +684,13 @@ std::vector<Function<T,NDIM>> Exchange<T,NDIM>::operator()(const std::vector<Fun
 };
 
 template<typename T, std::size_t NDIM>
-Exchange<T,NDIM>& Exchange<T,NDIM>::set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1) {
+Exchange<T,NDIM>& Exchange<T,NDIM>::set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1,
+                                                   const double mu1) {
     if (not impl) {
         World& world=bra.front().world();
         impl.reset(new Exchange<T,NDIM>::ExchangeImpl(world));
     }
-    impl->set_parameters(bra,ket,lo1);
+    impl->set_parameters(bra,ket,lo1,mu1);
     return *this;
 }
 
diff --git a/src/apps/chem/SCFOperators.h b/src/apps/chem/SCFOperators.h
index 8c975f25cae..c0f2d6fe047 100644
--- a/src/apps/chem/SCFOperators.h
+++ b/src/apps/chem/SCFOperators.h
@@ -122,10 +122,10 @@ class Exchange : public SCFOperatorBase<T,NDIM> {
     Exchange(std::shared_ptr<MacroTaskQ> taskq) : SCFOperatorBase<T, NDIM>(taskq) {}
 
     /// ctor with a conventional calculation
-    Exchange(World& world, const SCF *calc, const int ispin);
+    Exchange(World& world, const SCF *calc, int ispin, double mu = 0.0);
 
     /// ctor with a nemo calculation
-    Exchange(World& world, const Nemo *nemo, const int ispin);
+    Exchange(World& world, const Nemo *nemo, int ispin, double mu = 0.0);
 
     std::string info() const {return "K";}
 
@@ -142,7 +142,7 @@ class Exchange : public SCFOperatorBase<T,NDIM> {
         return *this;
     }
 
-    Exchange& set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1);
+    Exchange& set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1, const double mu1 = 0.0);
 
     Function<T, NDIM> operator()(const Function<T, NDIM>& ket) const {
         vecfuncT vket(1, ket);
diff --git a/src/apps/chem/commandlineparser.h b/src/apps/chem/commandlineparser.h
index db9b395ad68..23893184925 100644
--- a/src/apps/chem/commandlineparser.h
+++ b/src/apps/chem/commandlineparser.h
@@ -5,7 +5,10 @@
 #ifndef MADNESS_COMMANDLINEPARSER_H
 #define MADNESS_COMMANDLINEPARSER_H
 
-#include<map>
+#include <algorithm>
+#include <map>
+#include <sstream>
+#include <vector>
 namespace madness {
 /// very simple command line parser
 
diff --git a/src/apps/chem/exchangeoperator.cc b/src/apps/chem/exchangeoperator.cc
index 5dc4393991f..f6c80e66ce0 100644
--- a/src/apps/chem/exchangeoperator.cc
+++ b/src/apps/chem/exchangeoperator.cc
@@ -11,8 +11,8 @@ namespace madness {
 
 
 template<typename T, std::size_t NDIM>
-Exchange<T, NDIM>::ExchangeImpl::ExchangeImpl(World& world, const SCF *calc, const int ispin)
-        : world(world), symmetric_(false), lo(calc->param.lo()) {
+Exchange<T, NDIM>::ExchangeImpl::ExchangeImpl(World& world, const SCF *calc, const int ispin, const double mu)
+        : world(world), symmetric_(false), lo(calc->param.lo()), mu(mu) {
     if (ispin == 0) { // alpha spin
         mo_ket = convert<double, T, NDIM>(world, calc->amo);        // deep copy necessary if T==double_complex
     } else if (ispin == 1) {  // beta spin
@@ -23,8 +23,8 @@ Exchange<T, NDIM>::ExchangeImpl::ExchangeImpl(World& world, const SCF *calc, con
 
 template<typename T, std::size_t NDIM>
 Exchange<T, NDIM>::ExchangeImpl::ExchangeImpl(World& world, const Nemo *nemo,
-                            const int ispin) // @suppress("Class members should be properly initialized")
-        : ExchangeImpl(world, nemo->get_calc().get(), ispin) {
+                            const int ispin, const double mu) // @suppress("Class members should be properly initialized")
+        : ExchangeImpl(world, nemo->get_calc().get(), ispin, mu) {
 
     if (ispin == 0) { // alpha spin
         mo_ket = convert<double, T, NDIM>(world,
@@ -92,7 +92,7 @@ Exchange<T, NDIM>::ExchangeImpl::K_macrotask_efficient(const vecfuncT& vf, const
 
     // the result is a vector of functions living in the universe
     const long nresult = vf.size();
-    MacroTaskExchangeSimple xtask(nresult, lo, mul_tol, is_symmetric());
+    MacroTaskExchangeSimple xtask(nresult, mu, lo, mul_tol, is_symmetric());
     vecfuncT Kf;
 
     // deferred execution if a taskq is provided by the user
@@ -119,12 +119,12 @@ std::vector<Function<T, NDIM> > Exchange<T, NDIM>::ExchangeImpl::K_small_memory(
     const long nocc = mo_ket.size();
     const long nf = vket.size();
     vecfuncT Kf = zero_functions_compressed<T, NDIM>(world, nf);
-    auto poisson = set_poisson(world, lo);
+    auto yukawa = set_yukawa(world, mu, lo);
 
     for (int i = 0; i < nocc; ++i) {
         vecfuncT psif = mul_sparse(world, mo_bra[i], vket, mul_tol); /// was vtol
         truncate(world, psif);
-        psif = apply(world, *poisson.get(), psif);
+        psif = apply(world, *yukawa.get(), psif);
         truncate(world, psif);
         psif = mul_sparse(world, mo_ket[i], psif, mul_tol); /// was vtol
         gaxpy(world, 1.0, Kf, 1.0, psif);
@@ -137,8 +137,8 @@ template<typename T, std::size_t NDIM>
 std::vector<Function<T, NDIM> > Exchange<T, NDIM>::ExchangeImpl::K_large_memory(const vecfuncT& vket,
                                                                   const double mul_tol) const {    // Larger memory algorithm ... use i-j sym if psi==f
 
-    auto poisson = set_poisson(world, lo);
-    vecfuncT result = compute_K_tile(world, mo_bra, mo_ket, vket, poisson, is_symmetric(), mul_tol);
+    auto yukawa = set_yukawa(world, mu, lo);
+    vecfuncT result = compute_K_tile(world, mo_bra, mo_ket, vket, yukawa, is_symmetric(), mul_tol);
     truncate(world, result);
     return result;
 }
@@ -234,8 +234,8 @@ Exchange<T, NDIM>::ExchangeImpl::MacroTaskExchangeSimple::compute_offdiagonal_ba
     mul1_timer += long((cpu1 - cpu0) * 1000l);
 
     cpu0 = cpu_time();
-    auto poisson = set_poisson(subworld, lo);
-    vecfuncT Nij = apply(subworld, *poisson.get(), orbital_product_flat);
+    auto yukawa = set_yukawa(subworld, mu, lo);
+    vecfuncT Nij = apply(subworld, *yukawa.get(), orbital_product_flat);
     truncate(subworld, Nij);
     cpu1 = cpu_time();
     apply_timer += long((cpu1 - cpu0) * 1000l);
diff --git a/src/apps/chem/exchangeoperator.h b/src/apps/chem/exchangeoperator.h
index 0ee87c57c62..44da0204867 100644
--- a/src/apps/chem/exchangeoperator.h
+++ b/src/apps/chem/exchangeoperator.h
@@ -53,25 +53,27 @@ class Exchange<T,NDIM>::ExchangeImpl {
     ExchangeImpl(World& world) : world(world) {}
 
     /// ctor with a conventional calculation
-    ExchangeImpl(World& world, const SCF *calc, const int ispin) ;
+    ExchangeImpl(World& world, const SCF *calc, int ispin, double mu = 0.0) ;
 
     /// ctor with a nemo calculation
-    ExchangeImpl(World& world, const Nemo *nemo, const int ispin);
+    ExchangeImpl(World& world, const Nemo *nemo, int ispin, double mu = 0.0);
 
     /// set the bra and ket orbital spaces, and the occupation
 
     /// @param[in]	bra		bra space, must be provided as complex conjugate
     /// @param[in]	ket		ket space
-    void set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1) {
+    void set_parameters(const vecfuncT& bra, const vecfuncT& ket, const double lo1, const double mu1 = 0.0) {
         mo_bra = copy(world, bra);
         mo_ket = copy(world, ket);
         lo = lo1;
+        mu = mu1;
     }
 
     std::string info() const {return "K";}
 
-    static auto set_poisson(World& world, const double lo, const double econv = FunctionDefaults<3>::get_thresh()) {
-        return std::shared_ptr<real_convolution_3d>(CoulombOperatorPtr(world, lo, econv));
+    static auto set_yukawa(World& world, const double mu, const double lo,
+                           const double econv = FunctionDefaults<3>::get_thresh()) {
+        return std::shared_ptr<real_convolution_3d>(YukawaOperatorPtr(world, mu, lo, econv));
     }
 
     /// apply the exchange operator on a vector of functions
@@ -128,6 +130,7 @@ class Exchange<T,NDIM>::ExchangeImpl {
     bool symmetric_ = false;      /// is the exchange matrix symmetric? K phi_i = \sum_k \phi_k \int \phi_k \phi_i
     vecfuncT mo_bra, mo_ket;    ///< MOs for bra and ket
     double lo = 1.e-4;
+    double mu = 0.0;
     long printlevel = 0;
     double mul_tol = 0.0;
 
@@ -135,6 +138,7 @@ class Exchange<T,NDIM>::ExchangeImpl {
 
         long nresult;
         double lo = 1.e-4;
+        double mu = 0.0;
         double mul_tol = 1.e-7;
         bool symmetric = false;
 
@@ -186,8 +190,9 @@ class Exchange<T,NDIM>::ExchangeImpl {
         };
 
     public:
-        MacroTaskExchangeSimple(const long nresult, const double lo, const double mul_tol, const bool symmetric)
-                : nresult(nresult), lo(lo), mul_tol(mul_tol), symmetric(symmetric) {
+        MacroTaskExchangeSimple(const long nresult, const double mu,
+                                const double lo, const double mul_tol, const bool symmetric)
+                : nresult(nresult), mu(mu), lo(lo), mul_tol(mul_tol), symmetric(symmetric) {
             partitioner.reset(new MacroTaskPartitionerExchange(symmetric));
         }
 
@@ -264,9 +269,9 @@ class Exchange<T,NDIM>::ExchangeImpl {
         ) const {
             double mul_tol = 0.0;
             double symmetric = true;
-            auto poisson = Exchange<double, 3>::ExchangeImpl::set_poisson(subworld, lo);
-            return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, poisson, symmetric,
-                                                     mul_tol);
+            auto yukawa = Exchange<double, 3>::ExchangeImpl::set_yukawa(subworld, mu, lo);
+            return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, yukawa, symmetric,
+                                                                   mul_tol);
         }
 
         /// compute a batch of the exchange matrix, with non-identical ranges
@@ -282,9 +287,9 @@ class Exchange<T,NDIM>::ExchangeImpl {
                                                     const vecfuncT& vf_batch) const {
             double mul_tol = 0.0;
             double symmetric = false;
-            auto poisson = Exchange<double, 3>::ExchangeImpl::set_poisson(subworld, lo);
-            return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, poisson, symmetric,
-                                                     mul_tol);
+            auto yukawa = Exchange<double, 3>::ExchangeImpl::set_yukawa(subworld, mu, lo);
+            return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, yukawa, symmetric,
+                                                                   mul_tol);
         }
 
         /// compute a batch of the exchange matrix, with non-identical ranges
diff --git a/src/apps/chem/nemo.cc b/src/apps/chem/nemo.cc
index 983b5f20d17..cd6fe3db993 100644
--- a/src/apps/chem/nemo.cc
+++ b/src/apps/chem/nemo.cc
@@ -252,9 +252,16 @@ std::shared_ptr<Fock<double,3>> Nemo::make_fock_operator() const {
     fock->add_operator("J",std::make_shared<Coulomb<double,3> >(J));
     fock->add_operator("V",std::make_shared<Nuclear<double,3> >(world,this));
     fock->add_operator("T",std::make_shared<Kinetic<double,3> >(world));
-    if (calc->xc.hf_exchange_coefficient()>0.0) {
+    const double yukawa_separation = calc->xc.get_hf_yukawa_separation();
+    if (calc->xc.hf_exchange_coefficient()>0.0 && yukawa_separation > 1e-6) {
         Exchange<double,3> K=Exchange<double,3>(world,this,ispin).set_symmetric(false);
         fock->add_operator("K",{-1.0,std::make_shared<Exchange<double,3>>(K)});
+
+        if (yukawa_separation < 1e4) {
+            Exchange<double,3> K_short_range =
+                    Exchange<double,3>(world,this,ispin, yukawa_separation).set_symmetric(false);
+            fock->add_operator("K_short_range",{1.0,std::make_shared<Exchange<double,3>>(K_short_range)});
+        }
     }
     if (calc->xc.is_dft()) {
         XCOperator<double,3> xcoperator(world,this,ispin);
@@ -573,6 +580,10 @@ void Nemo::compute_nemo_potentials(const vecfuncT& nemo, vecfuncT& psi,
 		vecfuncT& Unemo) const {
 
     {
+        if (calc->xc.get_hf_yukawa_separation() < 1e4) {
+            throw std::invalid_argument("solve_cphf cannot yet handle range separation parameters less than 1e4!");
+        }
+
         timer t(world);
         real_function_3d vcoul;
         int ispin = 0;
@@ -1203,6 +1214,10 @@ vecfuncT Nemo::solve_cphf(const size_t iatom, const int iaxis, const Tensor<doub
 
     print("\nsolving nemo cphf equations for atom, axis",iatom,iaxis);
 
+    if (calc->xc.get_hf_yukawa_separation() < 1e4) {
+        throw std::invalid_argument("solve_cphf cannot yet handle range separation parameters less than 1e4!");
+    }
+
     vecfuncT xi=guess;
     // guess for the perturbed MOs
     const vecfuncT nemo=calc->amo;
diff --git a/src/apps/chem/polynomial.cc b/src/apps/chem/polynomial.cc
index 07c9427ab0d..9982a4ed25c 100644
--- a/src/apps/chem/polynomial.cc
+++ b/src/apps/chem/polynomial.cc
@@ -11,6 +11,7 @@
 #include <cmath>
 #include <vector>
 #include "polynomial.h"
+#include <stdexcept>
 
 namespace slymer {
 
diff --git a/src/apps/chem/test_dft.cc b/src/apps/chem/test_dft.cc
index 8e983131ecb..873afa8d0d9 100644
--- a/src/apps/chem/test_dft.cc
+++ b/src/apps/chem/test_dft.cc
@@ -32,6 +32,7 @@
 //#define WORLD_INSTANTIATE_STATIC_TEMPLATES
 #include <madness.h>
 #include <chem/SCFOperators.h>
+#include <chem/xcfunctional.h>
 
 using namespace madness;
 
@@ -108,6 +109,37 @@ int test_slater_exchange(World& world) {
 
 }
 
+int test_get_hf_yukawa_separation(World& world) {
+
+    const bool spin_polarized=false;
+    const std::string default_xc_data =
+            "GGA_X_ITYH_PBE .75 EXTERNAL_PARAMETERS {OMEGA: 0.2, kappa: 1.2} GGA_C_PBE 1. HF_X .25";
+    XCfunctional default_xc;
+    default_xc.initialize(default_xc_data, spin_polarized, world);
+
+    const double default_yukawa_separation = default_xc.get_hf_yukawa_separation();
+    print("default_yukawa_separation:", default_yukawa_separation);
+
+    const double expected_default_yukawa_separation = 1e9;
+    const double default_error = default_yukawa_separation-expected_default_yukawa_separation;
+    const double thresh = 1e-6;
+    int result=0;
+    if (check_err(default_error,thresh,"Default Yukawa separation error")) result+=1;
+
+    const std::string xc_data =
+            "GGA_X_ITYH_PBE .75 EXTERNAL_PARAMETERS {kappa: 1.2, _omega: 200} GGA_C_PBE 1. HF_X .25 YUKAWA_SEPARATION 0.7";
+    XCfunctional xc;
+    xc.initialize(xc_data, spin_polarized, world);
+
+    const double yukawa_separation = xc.get_hf_yukawa_separation();
+    print("yukawa_separation:", yukawa_separation);
+
+    const double expected_yukawa_separation = 0.7;
+    const double error = yukawa_separation-expected_yukawa_separation;
+    if (check_err(error,thresh,"Yukawa separation error")) result+=1;
+    return result;
+}
+
 int test_external_parameters(World& world) {
 
     real_function_3d dens=real_factory_3d(world).f(slater2).truncate_on_project();
@@ -215,6 +247,7 @@ int main(int argc, char** argv) {
     int result=0;
 
     result += test_slater_exchange(world);
+    result += test_get_hf_yukawa_separation(world);
     result += test_external_parameters(world);
     result += test_throw_if_external_parameters_are_specified_before_functionals(world);
     result += test_throw_if_external_parameters_missing_opening_brace(world);
diff --git a/src/apps/chem/xcfunctional.h b/src/apps/chem/xcfunctional.h
index 91517735e67..bff34afe0e6 100644
--- a/src/apps/chem/xcfunctional.h
+++ b/src/apps/chem/xcfunctional.h
@@ -97,6 +97,7 @@ class XCfunctional {
 
     bool spin_polarized;        ///< True if the functional is spin polarized
     double hf_coeff;            ///< Factor multiplying HF exchange (+1.0 gives HF)
+    double hf_yukawa_separation; ///< Range-separation parameter for the
     double rhomin, rhotol;      ///< See initialize and munge*
     double ggatol;              ///< See initialize and munge*
 
@@ -283,6 +284,12 @@ class XCfunctional {
         return hf_coeff;
     }
 
+    /// Returns the value of the Yukawa range-separation parameter for HF exchange
+    double get_hf_yukawa_separation() const
+    {
+        return hf_yukawa_separation;
+    }
+
     /// Computes the energy functional at given points
 
     /// This uses the convention that the total energy is
diff --git a/src/apps/chem/xcfunctional_ldaonly.cc b/src/apps/chem/xcfunctional_ldaonly.cc
index 948c2553441..b9b4e99ed02 100644
--- a/src/apps/chem/xcfunctional_ldaonly.cc
+++ b/src/apps/chem/xcfunctional_ldaonly.cc
@@ -17,12 +17,12 @@ int x_uks_s__(double *ra, double *rb, double *f, double *dfdra, double *dfdrb);
 int c_uks_vwn5__(double *ra, double *rb, double *f, double *dfdra, double *dfdrb);
 
 XCfunctional::XCfunctional() : hf_coeff(0.0) {
-    rhotol=1e-7; rhomin=1e-12; // default values
+    rhotol=1e-7; rhomin=1e-12; hf_yukawa_separation = 1e9; // default values
 }
 
 void XCfunctional::initialize(const std::string& input_line, bool polarized,
         World& world, bool verbose) {
-    rhotol=1e-7; rhomin=1e-12; // default values
+    rhotol=1e-7; rhomin=1e-12; hf_yukawa_separation = 1e9; // default values
 
     spin_polarized = polarized;
 
@@ -41,6 +41,9 @@ void XCfunctional::initialize(const std::string& input_line, bool polarized,
         else if (token == "RHOTOL") {
             s >> rhotol;
         }
+        else if (name == "YUKAWA_SEPARATION") {
+            line >> hf_yukawa_separation;
+        }
         else if (token == "HF") {
             hf_coeff = 1.0;
             found_valid_token = true;
diff --git a/src/apps/chem/xcfunctional_libxc.cc b/src/apps/chem/xcfunctional_libxc.cc
index 35f12474930..1c78cfd2685 100644
--- a/src/apps/chem/xcfunctional_libxc.cc
+++ b/src/apps/chem/xcfunctional_libxc.cc
@@ -60,6 +60,7 @@ XCfunctional::XCfunctional() : hf_coeff(0.0) {
     ggatol=1.e-4;
     nderiv=0;
     spin_polarized=false;
+    hf_yukawa_separation = 1e9;
 }
 
 void XCfunctional::initialize(const std::string& input_line, bool polarized,
@@ -77,6 +78,7 @@ void XCfunctional::initialize(const std::string& input_line, bool polarized,
 
     nderiv = 0;
     hf_coeff = 0.0;
+    hf_yukawa_separation = 1e9;
     funcs.clear();
 
     if (printit) print("\nConstruct XC Functional from LIBXC Library");
@@ -107,6 +109,8 @@ void XCfunctional::initialize(const std::string& input_line, bool polarized,
             line >> rhotol;
         } else if (name == "GGATOL") {
             line >> ggatol;
+        } else if (name == "YUKAWA_SEPARATION") {
+            line >> hf_yukawa_separation;
         } else if (name == "EXTERNAL_PARAMETERS") {
             /* EXTERNAL_PARAMETERS are functional specific parameters which can be set in Libxc.
              * In each Libxc functional file (e.g. src/gga_x_pbe.c), there are lines that look as follows:
diff --git a/src/apps/moldft/input b/src/apps/moldft/input
index e34470b5cf6..a49006a38d5 100644
--- a/src/apps/moldft/input
+++ b/src/apps/moldft/input
@@ -3,6 +3,7 @@ dft
   maxsub 5
 end
   #xc "GGA_X_PBE .75 GGA_C_PBE 1. HF_X .25"
+  #xc "GGA_X_SFAT 1.0 EXTERNAL_PARAMETERS {omega: 4.0} GGA_C_LYPR 1. EXTERNAL_PARAMETERS {omega: 2.4} HF_X 1.0 YUKAWA_SEPARATION 4.0"
   #xc "GGA_X_PBE .75 EXTERNAL_PARAMETERS {kappa: 0.7, mu: 0.3} GGA_C_PBE 1. EXTERNAL_PARAMETERS {gamma: 0.05, beta: 0.08} HF_X .25"
   # The external parameters are defined in the source files of each Libxc functional. For the previous input line,
   # they are defined in src/gga_x_pbe.c and src/gga_c_pbe.c.
diff --git a/src/madness/mra/gfit.h b/src/madness/mra/gfit.h
index d785adbab48..c133d7a7770 100644
--- a/src/madness/mra/gfit.h
+++ b/src/madness/mra/gfit.h
@@ -66,6 +66,13 @@ class GFit {
 		return fit;
 	}
 
+    /// return a fit for the Yukawa potential, f(r) = exp(-\mu r)/r
+    static GFit YukawaFit(double mu, double lo, double hi, double eps, bool prnt=false) {
+        GFit fit=BSHFit(mu,lo,hi,eps/(4.0*constants::pi),prnt);
+        fit.coeffs_.scale(4.0*constants::pi);
+        return fit;
+    }
+
 	/// return a fit for the bound-state Helmholtz function
 
 	/// the BSH function is defined by
diff --git a/src/madness/mra/operator.h b/src/madness/mra/operator.h
index 80e177de394..799d014fb21 100644
--- a/src/madness/mra/operator.h
+++ b/src/madness/mra/operator.h
@@ -1607,6 +1607,56 @@ namespace madness {
         return new SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
     }
 
+    /// Factory function generating separated kernel for convolution with exp(-mu*r)/r in 3D.
+    static
+    inline
+    SeparatedConvolution<double,3> YukawaOperator(World& world,
+                                                   double mu,
+                                                   double lo,
+                                                   double eps,
+                                                   const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+
+                                                   int k=FunctionDefaults<3>::get_k())
+    {
+        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
+        double hi = cell_width.normf(); // Diagonal width of cell
+        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
+
+        GFit<double,3> fit=GFit<double,3>::YukawaFit(mu,lo,hi,eps,false);
+        Tensor<double> coeff=fit.coeffs();
+        Tensor<double> expnt=fit.exponents();
+
+
+        if (bc(0,0) == BC_PERIODIC) {
+            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
+        }
+        return SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
+    }
+
+
+    /// Factory function generating separated kernel for convolution with exp(-mu*r)/r in 3D.
+    static
+    inline
+    SeparatedConvolution<double,3>* YukawaOperatorPtr(World& world,
+                                                       double mu,
+                                                       double lo,
+                                                       double eps,
+                                                       const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
+                                                       int k=FunctionDefaults<3>::get_k())
+    {
+        const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
+        double hi = cell_width.normf(); // Diagonal width of cell
+        if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
+        GFit<double,3> fit=GFit<double,3>::YukawaFit(mu,lo,hi,eps,false);
+        Tensor<double> coeff=fit.coeffs();
+        Tensor<double> expnt=fit.exponents();
+
+        if (bc(0,0) == BC_PERIODIC) {
+            fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
+        }
+        return new SeparatedConvolution<double,3>(world, coeff, expnt, bc, k);
+    }
+
 
     /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
     template <std::size_t NDIM>