adding a parameter for change SCF thresholds

fakeone · fakeone · commit 1b0f451b311e · 2024-12-02T11:43:36.000+08:00
diff --git a/source/module_elecstate/potentials/pot_xc.cpp b/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
diff --git a/source/module_esolver/esolver_ks.cpp b/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);
diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp
@@ -14,21 +14,27 @@ 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;
-double XC_Functional::default_dft_threshold=1.0e-10;
-double XC_Functional::dens_threshold=1.0e-10;
-double XC_Functional::zeta_threshold=1.0e-10;
-double XC_Functional::grho_threshold=1.0e-10;
-double XC_Functional::tau_threshold=1.0e-10;
 
-void XC_Functional::set_dft_ingred_thrs(double d_thr, double z_thr, double g_thr, double t_thr)
-{  dens_threshold = d_thr;
-   zeta_threshold = z_thr;
-   grho_threshold = g_thr;
-   tau_threshold = t_thr;
+void XC_Functional::set_xc_ingred_thrs(double _thr_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()
diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h
@@ -65,7 +65,7 @@ class XC_Functional
 	static int get_func_type();
     static void set_xc_type(const std::string xc_func_in);
 
-    static void set_dft_ingred_thrs(double d_thr, double z_thr, double g_thr, double t_thr);
+    static void set_xc_ingred_thrs(double _thr_in);
     static double get_dens_threshold();
     static double get_grho_threshold();
     static double get_zeta_threshold();
@@ -88,7 +88,6 @@ class XC_Functional
         static double grho_threshold;
         static double zeta_threshold;
         static double tau_threshold;
-        static double default_dft_threshold;
 	//exx_hybrid_alpha for mixing exx in hybrid functional:
 	static double hybrid_alpha;
 
diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/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
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp b/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
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/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
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/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
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/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
diff --git a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp b/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))
     {
diff --git a/source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp b/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;
diff --git a/source/module_hamilt_general/module_xc/xc_functional_wrapper_xc.cpp b/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;
-}
+}
diff --git a/source/module_io/read_input_item_elec_stru.cpp b/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";
diff --git a/source/module_parameter/input_parameter.h b/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

Original file line number	Diff line number	Diff line change
`@@ -19,7 +19,6 @@ void PotXC::cal_v_eff(const Chargeconst chg, const UnitCellconst ucell, Module`
`19`	`19`	`//----------------------------------------------------------`
`20`	`20`	`// calculate the exchange-correlation potential`
`21`	`21`	`//----------------------------------------------------------`
`22`		`-`
`23`	`22`	`if (XC_Functional::get_func_type() == 3 \|\| XC_Functional::get_func_type() == 5)`
`24`	`23`	`{`
`25`	`24`	`#ifdef USE_LIBXC`
`@@ -44,4 +43,4 @@ void PotXC::cal_v_eff(const Chargeconst chg, const UnitCellconst ucell, Module`
`44`	`43`	`ModuleBase::timer::tick("PotXC", "cal_v_eff");`
`45`	`44`	`}`
`46`	`45`
`47`		`-} // namespace elecstate`
	`46`	`+} // namespace elecstate`
Original file line number	Diff line number	Diff line change
`@@ -268,7 +268,7 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(`
`268`	`268`	`const ModuleBase::matrix &v)`
`269`	`269`	`{`
`270`	`270`	`assert(PARAM.inp.nspin==4);`
`271`		`- constexpr double vanishing_charge = 1.0e-10;`
	`271`	`+ const double vanishing_charge = XC_Functional::get_dens_threshold();`
`272`	`272`	`ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx);`
`273`	`273`	`for( int ir=0; ir<nrxx; ++ir )`
`274`	`274`	`v_nspin4(0,ir) = 0.5 * (v(0,ir)+v(1,ir));`
`@@ -287,4 +287,4 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(`
`287`	`287`	`return v_nspin4;`
`288`	`288`	`}`
`289`	`289`
`290`		`-#endif`
	`290`	`+#endif`
Original file line number	Diff line number	Diff line change
`@@ -1,6 +1,7 @@`
`1`	`1`	`#ifdef USE_LIBXC`
`2`	`2`
`3`	`3`	`#include "xc_functional_libxc.h"`
	`4`	`+#include "xc_functional.h"`
`4`	`5`
`5`	`6`	`#include <xc.h>`
`6`	`7`	`#include <array>`
`@@ -12,9 +13,12 @@ void XC_Functional_Libxc::gcxc_libxc(`
`12`	`13`	`{`
`13`	`14`	`sxc = v1xc = v2xc = 0.0;`
`14`	`15`
`15`		`- constexpr double small = 1.e-6;`
`16`		`- constexpr double smallg = 1.e-10;`
`17`		`- if (rho <= small \|\| grho < smallg)`
	`16`	`+// constexpr double small = 1.e-6;`
	`17`	`+// constexpr double smallg = 1.e-10;`
	`18`	`+`
	`19`	`+// std::cout << "hello from XC_Functional_Libxc, threshold = " << XC_Functional::get_dens_threshold() << std::endl;`
	`20`	`+`
	`21`	`+ if (rho <= XC_Functional::get_dens_threshold() \|\| grho < XC_Functional::get_grho_threshold())`
`18`	`22`	`{`
`19`	`23`	`return;`
`20`	`24`	`}`
`@@ -49,8 +53,9 @@ void XC_Functional_Libxc::gcxc_spin_libxc(`
`49`	`53`	`{`
`50`	`54`	`if( func.info->family == XC_FAMILY_GGA \|\| func.info->family == XC_FAMILY_HYB_GGA)`
`51`	`55`	`{`
`52`		`- constexpr double rho_threshold = 1E-6;`
`53`		`- constexpr double grho_threshold = 1E-10;`
	`56`	`+// std::cout << "hello libxc gcxc_spin_libxc, threshold: " << XC_Functional::get_dens_threshold() << std::endl;`
	`57`	`+ const double rho_threshold = XC_Functional::get_dens_threshold();`
	`58`	`+ const double grho_threshold = XC_Functional::get_grho_threshold();`
`54`	`59`	`std::array<double,2> sgn = {1.0, 1.0};`
`55`	`60`	`if(func.info->kind==XC_CORRELATION)`
`56`	`61`	`{`
`@@ -76,4 +81,4 @@ void XC_Functional_Libxc::gcxc_spin_libxc(`
`76`	`81`	`XC_Functional_Libxc::finish_func(funcs);`
`77`	`82`	`}`
`78`	`83`
`79`		`-#endif`
	`84`	`+#endif`
Original file line number	Diff line number	Diff line change
`@@ -5,6 +5,7 @@`
`5`	`5`	`#ifdef USE_LIBXC`
`6`	`6`
`7`	`7`	`#include "xc_functional_libxc.h"`
	`8`	`+#include "xc_functional.h"`
`8`	`9`	`#include "module_hamilt_pw/hamilt_pwdft/global.h"`
`9`	`10`	`#include <array>`
`10`	`11`
`@@ -67,8 +68,8 @@ void XC_Functional_Libxc::tau_xc_spin(`
`67`	`68`	`{`
`68`	`69`	`if( func.info->family == XC_FAMILY_MGGA \|\| func.info->family == XC_FAMILY_HYB_MGGA)`
`69`	`70`	`{`
`70`		`- constexpr double rho_threshold = 1E-6;`
`71`		`- constexpr double grho_threshold = 1E-10;`
	`71`	`+ const double rho_threshold = XC_Functional::get_dens_threshold();`
	`72`	`+ const double grho_threshold = XC_Functional::get_grho_threshold();`
`72`	`73`	`std::array<double,2> sgn = {1.0, 1.0};`
`73`	`74`	`if(func.info->kind==XC_CORRELATION)`
`74`	`75`	`{`
`@@ -98,4 +99,4 @@ void XC_Functional_Libxc::tau_xc_spin(`
`98`	`99`	`XC_Functional_Libxc::finish_func(funcs);`
`99`	`100`	`}`
`100`	`101`
`101`		`-#endif`
	`102`	`+#endif`
Original file line number	Diff line number	Diff line change
`@@ -40,7 +40,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc(`
`40`	`40`	`// in Rydeberg unit, so * 2.0.`
`41`	`41`	`double e2 = 2.0;`
`42`	`42`
`43`		`- double vanishing_charge = 1.0e-10;`
	`43`	`+ double vanishing_charge = dens_threshold;`
`44`	`44`
`45`	`45`	`if (PARAM.inp.nspin == 1 \|\| ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z))`
`46`	`46`	`{`