Pr adding range separated exchange

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;

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()){

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;
 }
 

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);

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
 

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);

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,13 +130,15 @@ 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;
 
     class MacroTaskExchangeSimple : public MacroTaskOperationBase {
 
         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

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;

src/apps/chem/polynomial.cc

@@ -11,6 +11,7 @@
 #include <cmath>
 #include <vector>
 #include "polynomial.h"
+#include <stdexcept>
 
 namespace slymer {
 

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);