Adding an input keyword to set density thresholds in DFT calculations

source/module_elecstate/potentials/pot_xc.cpp

@@ -19,7 +19,6 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
     //----------------------------------------------------------
     //  calculate the exchange-correlation potential
     //----------------------------------------------------------
-
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
 #ifdef USE_LIBXC
@@ -44,4 +43,4 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
     ModuleBase::timer::tick("PotXC", "cal_v_eff");
 }
 
-} // namespace elecstate
+} // namespace elecstate

source/module_esolver/esolver_ks.cpp

@@ -189,6 +189,11 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     //! 4) it has been established that
     // xc_func is same for all elements, therefore
     // only the first one if used
+
+    // A little hacking of the code here. Adding the
+    // following line to set up the density threshold
+    // for DFT calculations.
+    XC_Functional::set_xc_ingred_thrs(PARAM.inp.xc_dens_thr);
     if (PARAM.inp.use_paw)
     {
         XC_Functional::set_xc_type(PARAM.inp.dft_functional);

source/module_hamilt_general/module_xc/xc_functional.cpp

@@ -14,21 +14,59 @@ XC_Functional::XC_Functional(){}
 
 XC_Functional::~XC_Functional(){}
 
+//double XC_Functional::dens_threshold;
+//double XC_Functional::zeta_threshold;
+//double XC_Functional::grho_threshold;
+//double XC_Functional::tau_threshold;
+
+double XC_Functional::dens_threshold = 1.0e-10;
+double XC_Functional::zeta_threshold = dens_threshold;
+double XC_Functional::grho_threshold = dens_threshold;
+double XC_Functional::tau_threshold = dens_threshold;
+
 std::vector<int> XC_Functional::func_id(1);
 int XC_Functional::func_type = 0;
 bool XC_Functional::use_libxc = true;
 double XC_Functional::hybrid_alpha = 0.25;
 
-void XC_Functional::set_hybrid_alpha(const double alpha_in)
+void XC_Functional::set_xc_ingred_thrs(double _thr_in)
 {
-    hybrid_alpha = alpha_in;
+    XC_Functional::dens_threshold = _thr_in;
+    XC_Functional::zeta_threshold = dens_threshold;
+    XC_Functional::grho_threshold = dens_threshold*1.0e-4;
+    XC_Functional::tau_threshold = dens_threshold*1.0e-4;
+}
+
+double XC_Functional::get_dens_threshold()
+{
+    return dens_threshold;
+}
+
+double XC_Functional::get_grho_threshold()
+{
+    return grho_threshold;
+}
+
+double XC_Functional::get_zeta_threshold()
+{
+    return zeta_threshold;
+}
+
+double XC_Functional::get_tau_threshold()
+{
+    return tau_threshold;
 }
 
 double XC_Functional::get_hybrid_alpha()
 {
     return hybrid_alpha;
 }
 
+void XC_Functional::set_hybrid_alpha(const double alpha_in)
+{
+    hybrid_alpha = alpha_in;
+}
+
 int XC_Functional::get_func_type()
 {
     return func_type;

source/module_hamilt_general/module_xc/xc_functional.h

@@ -65,6 +65,12 @@ class XC_Functional
 	static int get_func_type();
     static void set_xc_type(const std::string xc_func_in);
 
+    static void set_xc_ingred_thrs(double _thr_in);
+    static double get_dens_threshold();
+    static double get_grho_threshold();
+    static double get_zeta_threshold();
+    static double get_tau_threshold();
+
     // For hybrid functional
     static void set_hybrid_alpha(const double alpha_in);
     static double get_hybrid_alpha();
@@ -77,6 +83,11 @@ class XC_Functional
 	static int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga
 	static bool use_libxc;
 
+
+        static double dens_threshold;
+        static double grho_threshold;
+        static double zeta_threshold;
+        static double tau_threshold;
 	//exx_hybrid_alpha for mixing exx in hybrid functional:
 	static double hybrid_alpha;
 

source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp

@@ -191,8 +191,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 
 	}
 
-	const double epsr = 1.0e-6;
-	const double epsg = 1.0e-10;
+	const double epsr = dens_threshold;
+	const double epsg = grho_threshold;
 
 	double vtxcgc = 0.0;
 	double etxcgc = 0.0;
@@ -791,4 +791,4 @@ template void XC_Functional::grad_wfc<std::complex<double>, base_device::DEVICE_
     const ModulePW::PW_Basis_K* wfc_basis,
     const std::complex<double>* rhog,
     std::complex<double>* grad);
-#endif // __CUDA || __ROCM
+#endif // __CUDA || __ROCM

source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp

@@ -268,7 +268,7 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
 	const ModuleBase::matrix &v)
 {
 	assert(PARAM.inp.nspin==4);
-	constexpr double vanishing_charge = 1.0e-10;
+	const double vanishing_charge = XC_Functional::get_dens_threshold();
 	ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx);
 	for( int ir=0; ir<nrxx; ++ir )
 		v_nspin4(0,ir) = 0.5 * (v(0,ir)+v(1,ir));
@@ -287,4 +287,4 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
 	return v_nspin4;
 }
 
-#endif
+#endif

source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp

@@ -76,11 +76,13 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
     double vtxc = 0.0;
     ModuleBase::matrix v(nspin,nrxx);
 
+//    std::cout<< "hello from v_xc_libxc, threshold = " << XC_Functional::get_dens_threshold() << " " << XC_Functional::get_grho_threshold()<< std::endl; 
+
     for( xc_func_type &func : funcs )
     {
         // jiyy add for threshold
-        constexpr double rho_threshold = 1E-6;
-        constexpr double grho_threshold = 1E-10;
+        const double rho_threshold = XC_Functional::get_dens_threshold();
+        const double grho_threshold = XC_Functional::get_grho_threshold();
 
         xc_func_set_dens_threshold(&func, rho_threshold);
 
@@ -207,9 +209,9 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     std::vector<double> vtau   ( nrxx * nspin            );
     std::vector<double> vlapl  ( nrxx * nspin            );
 
-    constexpr double rho_th  = 1e-8;
-    constexpr double grho_th = 1e-12;
-    constexpr double tau_th  = 1e-8;
+    const double rho_th  = XC_Functional::get_dens_threshold();
+    const double grho_th = XC_Functional::get_grho_threshold();
+    const double tau_th  = XC_Functional::get_tau_threshold();
     // sgn for threshold mask
     std::vector<double> sgn( nrxx * nspin);
 #ifdef _OPENMP
@@ -392,4 +394,4 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) );
 }
 
-#endif
+#endif

source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp

@@ -1,6 +1,7 @@
 #ifdef USE_LIBXC
 
 #include "xc_functional_libxc.h"
+#include "xc_functional.h"
 
 #include <xc.h>
 #include <array>
@@ -12,9 +13,12 @@ void XC_Functional_Libxc::gcxc_libxc(
 {
     sxc = v1xc = v2xc = 0.0;
 
-    constexpr double small = 1.e-6;
-    constexpr double smallg = 1.e-10;
-    if (rho <= small || grho < smallg)
+//    constexpr double small = 1.e-6;
+//    constexpr double smallg = 1.e-10;
+
+//    std::cout << "hello from XC_Functional_Libxc, threshold = " << XC_Functional::get_dens_threshold() << std::endl;
+
+    if (rho <= XC_Functional::get_dens_threshold() || grho < XC_Functional::get_grho_threshold())
     {
         return;
     }
@@ -49,8 +53,9 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
     {
         if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
         {
-            constexpr double rho_threshold = 1E-6;
-            constexpr double grho_threshold = 1E-10;
+//            std::cout << "hello libxc gcxc_spin_libxc, threshold: " << XC_Functional::get_dens_threshold() << std::endl;
+            const double rho_threshold = XC_Functional::get_dens_threshold();
+            const double grho_threshold = XC_Functional::get_grho_threshold();
             std::array<double,2> sgn = {1.0, 1.0};
             if(func.info->kind==XC_CORRELATION)
             {
@@ -76,4 +81,4 @@ void XC_Functional_Libxc::gcxc_spin_libxc(
     XC_Functional_Libxc::finish_func(funcs);
 }
 
-#endif
+#endif

source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp

@@ -5,6 +5,7 @@
 #ifdef USE_LIBXC
 
 #include "xc_functional_libxc.h"
+#include "xc_functional.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include <array>
 
@@ -67,8 +68,8 @@ void XC_Functional_Libxc::tau_xc_spin(
     {
         if( func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA)
         {
-            constexpr double rho_threshold = 1E-6;
-            constexpr double grho_threshold = 1E-10;
+            const double rho_threshold = XC_Functional::get_dens_threshold();
+            const double grho_threshold = XC_Functional::get_grho_threshold();
             std::array<double,2> sgn = {1.0, 1.0};
             if(func.info->kind==XC_CORRELATION)
             {
@@ -98,4 +99,4 @@ void XC_Functional_Libxc::tau_xc_spin(
     XC_Functional_Libxc::finish_func(funcs);
 }
 
-#endif
+#endif

source/module_hamilt_general/module_xc/xc_functional_vxc.cpp

@@ -40,7 +40,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc(
     // in Rydeberg unit, so * 2.0.
     double e2 = 2.0;
 
-    double vanishing_charge = 1.0e-10;
+    double vanishing_charge = dens_threshold;
 
     if (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z))
     {

source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp

@@ -40,12 +40,10 @@ void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc,
     // USE kinds
     // implicit none
     // real rho, grho, sx, sc, v1x, v2x, v1c, v2c;
-    const double small = 1.e-6;
-    const double smallg = 1.e-10;
     double s,v1,v2;
     sxc = v1xc = v2xc = 0.0;
-
-    if (rho <= small || grho < smallg)
+    std::cout << "threshold for densities " << dens_threshold << " " << rho << " " << grho<< std::endl;
+    if (rho <= dens_threshold || grho < grho_threshold)
     {
         return;
     }
@@ -127,7 +125,7 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double
     // derivatives of exchange wr. grho
 
     // parameter :
-    double small = 1.e-10;
+    double small = dens_threshold;
     double sxup, sxdw;
     int iflag;
 
@@ -139,6 +137,7 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double
     v1xdw = 0.00; v2xdw = 0.00;
     sxup  = 0.00; sxdw  = 0.00;
 
+
     if (rho <= small)
     {
         return;
@@ -245,14 +244,17 @@ void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc,
     // derivatives of correlation wr. grho
 
     // parameter :
-    double small = 1.0e-10;
-	double epsr = 1.0e-6;
+    double small = dens_threshold;
+    double epsr = zeta_threshold;
 
     double x;
 
     sc = 0.00;
     v1cup = 0.00; v1cdw = 0.00;
     v2c = 0.00;
+
+//    std::cout << "hello from gcc_spin, threshold = " << dens_threshold << " " << grho_threshold << std::endl;
+
     if (std::abs(zeta) - 1.0 > small || rho <= small || sqrt(std::abs(grho)) <= small)
     {
         return;

source/module_hamilt_general/module_xc/xc_functional_wrapper_xc.cpp

@@ -121,4 +121,4 @@ void XC_Functional::xc_spin(const double &rho, const double &zeta,
         vxcdw += vdw;
 	}
 	return;
-}
+}

source/module_io/read_input_item_elec_stru.cpp

@@ -541,6 +541,19 @@ void ReadInput::item_elec_stru()
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("xc_dens_thr");
+        item.annotation = "threshold for density and other dft functional ingredients";
+        read_sync_double(input.xc_dens_thr);
+        item.check_value = [] (const Input_Item& item, const Parameter& para){
+            if (para.input.xc_dens_thr <= 0.0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "xc_dens_thr should be positive");
+            }
+        std::cout<< "resetting xc_dens_thr " << para.input.xc_dens_thr<< std::endl;
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("scf_ene_thr");
         item.annotation = "total energy error threshold";

source/module_parameter/input_parameter.h

@@ -112,6 +112,7 @@ struct Input_para
     bool gamma_only = false;   ///< for plane wave.
     int scf_nmax = 100;        ///< number of max elec iter
     double scf_thr = -1.0;     ///< \sum |rhog_out - rhog_in |^2
+    double xc_dens_thr = 1.0e-10; ///< threshold for density and other DFT ingredients to toss away in DFT calculation.
     double scf_ene_thr = -1.0; ///< energy threshold for scf convergence, in eV
     int scf_thr_type = -1;     ///< type of the criterion of scf_thr, 1: reci drho, 2: real drho
     bool final_scf = false;    ///< whether to do final scf