Hytd

source/module_esolver/esolver_ks_lcao.cpp

@@ -1106,7 +1106,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
         // ModuleRPA::DFT_RPA_interface rpa_interface(GlobalC::exx_info.info_global);
         // rpa_interface.rpa_exx_lcao().info.files_abfs = GlobalV::rpa_orbitals;
         RPA_LRI<TK, double> rpa_lri_double(GlobalC::exx_info.info_ri, GlobalC::exx_info.info_ewald);
-        rpa_lri_double.cal_postSCF_exx(*dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+        rpa_lri_double.cal_postSCF_exx(istep,
+                                       *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                        MPI_COMM_WORLD,
                                        this->kv,
                                        orb_);

source/module_hamilt_general/module_xc/exx_info.h

@@ -14,8 +14,7 @@ struct Exx_Info
 
         Conv_Coulomb_Pot_K::Ccp_Type ccp_type;
         double hybrid_alpha = 0.25;
-        double cam_alpha = 0.0;
-        double cam_beta = 0.0;
+        double hybrid_beta = 0.0;
         double hse_omega = 0.11;
         double mixing_beta_for_loop1 = 1.0;
 
@@ -43,8 +42,6 @@ struct Exx_Info
         Singular_Value::Fq_type fq_type;
         const bool& use_ewald;
 
-        double ewald_qdiv = 2;
-
         Exx_Info_Ewald(const Exx_Info::Exx_Info_Global& info_global) : use_ewald(info_global.use_ewald)
         {
         }
@@ -55,10 +52,11 @@ struct Exx_Info
     {
         const Conv_Coulomb_Pot_K::Ccp_Type& ccp_type;
         const double& hse_omega;
-        const double& cam_alpha;
-        const double& cam_beta;
+        const double& hybrid_alpha;
+        const double& hybrid_beta;
 
         bool real_number = false;
+        int Rcut_type = 1;
 
         double pca_threshold = 0;
         std::vector<std::string> files_abfs;
@@ -78,8 +76,8 @@ struct Exx_Info
         int abfs_Lmax = 0; // tmp
 
         Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
-            : ccp_type(info_global.ccp_type), hse_omega(info_global.hse_omega), cam_alpha(info_global.cam_alpha),
-              cam_beta(info_global.cam_beta)
+            : ccp_type(info_global.ccp_type), hse_omega(info_global.hse_omega), hybrid_alpha(info_global.hybrid_alpha),
+              hybrid_beta(info_global.hybrid_beta)
         {
         }
     };

source/module_hamilt_general/module_xc/xc_functional.cpp

@@ -18,14 +18,12 @@ 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::cam_alpha = 0.0;
-double XC_Functional::cam_beta = 0.0;
+double XC_Functional::hybrid_beta = 0.0;
 
-void XC_Functional::set_hybrid_alpha(const double alpha_in, const double cam_alpha_in, const double cam_beta_in)
+void XC_Functional::set_hybrid_alpha(const double alpha_in, const double beta_in)
 {
     hybrid_alpha = alpha_in;
-    cam_alpha = cam_alpha_in;
-    cam_beta = cam_beta_in;
+    hybrid_beta = beta_in;
 }
 
 double XC_Functional::get_hybrid_alpha()

source/module_hamilt_general/module_xc/xc_functional.h

@@ -66,7 +66,7 @@ class XC_Functional
     static void set_xc_type(const std::string xc_func_in);
 
     // For hybrid functional
-    static void set_hybrid_alpha(const double alpha_in, const double cam_alpha_in, const double cam_beta_in);
+    static void set_hybrid_alpha(const double alpha_in, const double beta_in);
     static double get_hybrid_alpha();
     /// Usually in exx caculation, the first SCF loop should be converged with PBE
     static void set_xc_first_loop(const UnitCell& ucell);
@@ -79,8 +79,7 @@ class XC_Functional
 
 	//exx_hybrid_alpha for mixing exx in hybrid functional:
 	static double hybrid_alpha;
-	static double cam_alpha;
-	static double cam_beta;
+	static double hybrid_beta;
 
 	public:
 	static std::vector<int> get_func_id() { return func_id; }

source/module_hamilt_general/module_xc/xc_functional_libxc.cpp

@@ -157,8 +157,8 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
             // and  0.0% short-range and 100.0% long-range exact exchange,
             // using the error function kernel.
 			add_func( XC_HYB_GGA_XC_LC_WPBE );	
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.cam_alpha,  //Fraction of Hartree-Fock exchange: 1.0
-				GlobalC::exx_info.info_global.cam_beta,  //Fraction of short-range exact exchange: -1.0
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,  //Fraction of Hartree-Fock exchange: 1.0
+				GlobalC::exx_info.info_global.hybrid_beta,  //Fraction of short-range exact exchange: -1.0
 				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.4
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}
@@ -168,8 +168,8 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
             // and  0.0% short-range and 100.0% long-range exact exchange,
             // using the error function kernel.
 			add_func( XC_HYB_GGA_XC_LRC_WPBE );	
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.cam_alpha,  //Fraction of Hartree-Fock exchange: 1.0
-				GlobalC::exx_info.info_global.cam_beta,  //Fraction of short-range exact exchange: -1.0
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,  //Fraction of Hartree-Fock exchange: 1.0
+				GlobalC::exx_info.info_global.hybrid_beta,  //Fraction of short-range exact exchange: -1.0
 				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.3
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}
@@ -179,8 +179,8 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
             // and 20.0% short-range and 100.0% long-range exact exchange,
             // using the error function kernel.
 			add_func( XC_HYB_GGA_XC_LRC_WPBEH );	
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.cam_alpha,  //Fraction of Hartree-Fock exchange: 1.0
-				GlobalC::exx_info.info_global.cam_beta,  //Fraction of short-range exact exchange: -0.8
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,  //Fraction of Hartree-Fock exchange: 1.0
+				GlobalC::exx_info.info_global.hybrid_beta,  //Fraction of short-range exact exchange: -0.8
 				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.2
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}
@@ -190,8 +190,8 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
             // and 100.0% short-range and 20.0% long-range exact exchange,
             // using the error function kernel.
 			add_func( XC_HYB_GGA_XC_CAM_PBEH);	
-			double parameter_hse[3] = { GlobalC::exx_info.info_global.cam_alpha,  //Fraction of Hartree-Fock exchange: 0.2
-				GlobalC::exx_info.info_global.cam_beta,  //Fraction of short-range exact exchange: 0.8
+			double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha,  //Fraction of Hartree-Fock exchange: 0.2
+				GlobalC::exx_info.info_global.hybrid_beta,  //Fraction of short-range exact exchange: 0.8
 				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.7
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp

@@ -353,7 +353,7 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
                                                                              &GlobalC::GridD);
             this->getOperator()->add(td_nonlocal);
         }
-        if ((elecstate::H_TDDFT_pw::stype == 2))
+        if (elecstate::H_TDDFT_pw::stype == 2)
         {
             Operator<TK>* td_pot_mixing = new TD_mixing_pot<OperatorLCAO<TK, TR>>(this->hsk,
                                                                            this->kv->kvec_d,

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp

@@ -242,7 +242,7 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
 template<typename TK, typename TR>
 void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
 {
-    ModuleBase::TITLE("OperatorEXX", "constributeHR");
+    ModuleBase::TITLE("OperatorEXX", "constributeHk");
     // Peize Lin add 2016-12-03
     if (PARAM.inp.calculation != "nscf" && this->two_level_step != nullptr && *this->two_level_step == 0 && !this->restart) { return; }  //in the non-exx loop, do nothing 
 

source/module_io/input_conv.cpp

@@ -367,8 +367,6 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_global.cal_exx = true;
         if (PARAM.inp.exx_use_ewald)
         {
-            GlobalC::exx_info.info_global.cam_alpha = std::stod(PARAM.inp.exx_cam_alpha);
-            GlobalC::exx_info.info_global.cam_beta = std::stod(PARAM.inp.exx_cam_beta);
             GlobalC::exx_info.info_global.use_ewald = true;
             GlobalC::exx_info.info_ewald.fq_type = Singular_Value::Fq_type(PARAM.inp.exx_fq_type);
             GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Ccp;
@@ -392,8 +390,6 @@ void Input_Conv::Convert()
             || dft_functional_lower == "cam_pbeh")
     {
         GlobalC::exx_info.info_global.cal_exx = true;
-        GlobalC::exx_info.info_global.cam_alpha = std::stod(PARAM.inp.exx_cam_alpha);
-        GlobalC::exx_info.info_global.cam_beta = std::stod(PARAM.inp.exx_cam_beta);
         if (PARAM.inp.exx_use_ewald)
         {
             GlobalC::exx_info.info_global.use_ewald = true;
@@ -413,15 +409,16 @@ void Input_Conv::Convert()
         // EXX case, convert all EXX related variables
         // GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.hybrid_alpha = std::stod(PARAM.inp.exx_hybrid_alpha);
+        GlobalC::exx_info.info_global.hybrid_beta = std::stod(PARAM.inp.exx_hybrid_beta);
         XC_Functional::set_hybrid_alpha(std::stod(PARAM.inp.exx_hybrid_alpha), 
-                                            std::stod(PARAM.inp.exx_cam_alpha), 
-                                            std::stod(PARAM.inp.exx_cam_beta));
+                                            std::stod(PARAM.inp.exx_hybrid_beta));
         GlobalC::exx_info.info_global.hse_omega = std::stod(PARAM.inp.exx_hse_omega);
         GlobalC::exx_info.info_global.separate_loop = PARAM.inp.exx_separate_loop;
         GlobalC::exx_info.info_global.hybrid_step = PARAM.inp.exx_hybrid_step;
         GlobalC::exx_info.info_global.mixing_beta_for_loop1 = PARAM.inp.exx_mixing_beta;
         GlobalC::exx_info.info_lip.lambda = PARAM.inp.exx_lambda;
 
+        GlobalC::exx_info.info_ri.Rcut_type = PARAM.inp.exx_spencer_type;
         GlobalC::exx_info.info_ri.real_number = std::stoi(PARAM.inp.exx_real_number);
         GlobalC::exx_info.info_ri.pca_threshold = PARAM.inp.exx_pca_threshold;
         GlobalC::exx_info.info_ri.C_threshold = PARAM.inp.exx_c_threshold;
@@ -435,7 +432,6 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_ri.cauchy_force_threshold = PARAM.inp.exx_cauchy_force_threshold;
         GlobalC::exx_info.info_ri.cauchy_stress_threshold = PARAM.inp.exx_cauchy_stress_threshold;
         GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times);
-        GlobalC::exx_info.info_ewald.ewald_qdiv = PARAM.inp.exx_ewald_qdiv;
 
         Exx_Abfs::Jle::Lmax = PARAM.inp.exx_opt_orb_lmax;
         Exx_Abfs::Jle::Ecut_exx = PARAM.inp.exx_opt_orb_ecut;

source/module_io/input_conv.h

@@ -47,7 +47,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
 {
     ModuleBase::TITLE("Input_Conv", "parse_expression");
     int count = 0;
-    std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)");
+    std::string pattern("([-+]?[0-9]+\\*[-+]?[0-9.]+|[-+]?[0-9,.]+)");
     std::vector<std::string> str;
     std::stringstream ss(fn);
     std::string section;
@@ -103,6 +103,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
         {
             int pos = sub_str.find("*");
             int num = stoi(sub_str.substr(0, pos));
+            assert(num>=0);
             T occ = stof(sub_str.substr(pos + 1, sub_str.size()));
             // std::vector<double> ocp_temp(num, occ);
             // const std::vector<double>::iterator dest = vec.begin() + count;

source/module_io/read_input.cpp

@@ -398,7 +398,7 @@ void ReadInput::write_txt_input(const Parameter& param, const std::string& filen
         {
             ofs << "\n#Parameters (14.vdW Correction)" << std::endl;
         }
-        else if (p_item->label == "exx_hybrid_alpha" || p_item->label == "exx_use_ewald")
+        else if (p_item->label == "exx_hybrid_alpha")
         {
             ofs << "\n#Parameters (15.exx)" << std::endl;
         }

source/module_io/read_input_item_exx_dftu.cpp

@@ -17,7 +17,9 @@ void ReadInput::item_exx()
                 std::string& dft_functional = para.input.dft_functional;
                 std::string dft_functional_lower = dft_functional;
                 std::transform(dft_functional.begin(), dft_functional.end(), dft_functional_lower.begin(), tolower);
-                if (dft_functional_lower == "hf")
+                if (dft_functional_lower == "hf" ||
+                    dft_functional_lower == "lc_pbe" || dft_functional_lower == "lc_wpbe" ||
+                    dft_functional_lower == "lrc_wpbe" || dft_functional_lower == "lrc_wpbeh")
                 {
                     para.input.exx_hybrid_alpha = "1";
                 }
@@ -26,6 +28,10 @@ void ReadInput::item_exx()
                 {
                     para.input.exx_hybrid_alpha = "0.25";
                 }
+                else if (dft_functional_lower == "cam_pbeh")
+                {
+                    para.input.exx_hybrid_alpha = "0.2";
+                }
                 else
                 { // no exx in scf, but will change to non-zero in
                     // postprocess like rpa
@@ -48,81 +54,48 @@ void ReadInput::item_exx()
         read_sync_int(input.exx_fq_type);
         this->add_item(item);
     }
-    {
-        Input_Item item("exx_ewald_qdiv");
-        item.annotation = "the order of q-divergence in auxiliary function";
-        read_sync_double(input.exx_ewald_qdiv);
-        this->add_item(item);
-    }
     {
         Input_Item item("exx_use_ewald");
         item.annotation = "if 1, use Ewald method to construct V matrix";
         read_sync_bool(input.exx_use_ewald);
         this->add_item(item);
     }
     {
-        Input_Item item("exx_cam_alpha");
-        item.annotation = "fraction of the full-range parts of Fock exchange in range-separated hybrid funtionals";
-        read_sync_string(input.exx_cam_alpha);
-        item.reset_value = [](const Input_Item& item, Parameter& para) {
-            if (para.input.exx_cam_alpha == "default")
-            {
-                std::string& dft_functional = para.input.dft_functional;
-                std::string dft_functional_lower = dft_functional;
-                std::transform(dft_functional.begin(), dft_functional.end(), dft_functional_lower.begin(), tolower);
-                if (dft_functional_lower == "hf" ||
-                    dft_functional_lower == "lc_pbe" || dft_functional_lower == "lc_wpbe" ||
-                    dft_functional_lower == "lrc_wpbe" || dft_functional_lower == "lrc_wpbeh")
-                {
-                    para.input.exx_cam_alpha = "1";
-                }
-                else if (dft_functional_lower == "cam_pbeh")
-                {
-                    para.input.exx_cam_alpha = "0.2";
-                }
-                else if (dft_functional_lower == "pbe0"
-                         || dft_functional_lower == "scan0")
-                {
-                    para.input.exx_cam_alpha = "0.25";
-                }
-                else
-                {
-                    para.input.exx_cam_alpha = "0";
-                }
-            }
-        };
+        Input_Item item("exx_spencer_type");
+        item.annotation = "set Spencer scheme type";
+        read_sync_int(input.exx_spencer_type);
         this->add_item(item);
     }
-        {
-        Input_Item item("exx_cam_beta");
-        item.annotation = "fraction of the short-range parts of Fock exchange in range-separated hybrid funtionals";
-        read_sync_string(input.exx_cam_beta);
+    {
+        Input_Item item("exx_hybrid_beta");
+        item.annotation = "another fraction of Fock exchange in range-separated hybrid funtionals";
+        read_sync_string(input.exx_hybrid_beta);
         item.reset_value = [](const Input_Item& item, Parameter& para) {
-            if (para.input.exx_cam_beta == "default")
+            if (para.input.exx_hybrid_beta == "default")
             {
                 std::string& dft_functional = para.input.dft_functional;
                 std::string dft_functional_lower = dft_functional;
                 std::transform(dft_functional.begin(), dft_functional.end(), dft_functional_lower.begin(), tolower);
                 if (dft_functional_lower == "lc_pbe" || dft_functional_lower == "lc_wpbe" ||
                     dft_functional_lower == "lrc_wpbe")
                 {
-                    para.input.exx_cam_beta = "-1";
+                    para.input.exx_hybrid_beta = "-1";
                 }
                 else if (dft_functional_lower == "lrc_wpbeh")
                 {
-                    para.input.exx_cam_beta = "-0.8";
+                    para.input.exx_hybrid_beta = "-0.8";
                 }
                 else if (dft_functional_lower == "cam_pbeh")
                 {
-                    para.input.exx_cam_beta = "0.8";
+                    para.input.exx_hybrid_beta = "0.8";
                 }
                 else if (dft_functional_lower == "hse")
                 {
-                    para.input.exx_cam_beta = "0.25";
+                    para.input.exx_hybrid_beta = "0.25";
                 }
                 else
                 {
-                    para.input.exx_cam_beta = "0";
+                    para.input.exx_hybrid_beta = "0";
                 }
             }
         };

source/module_io/test/read_input_ptest.cpp

@@ -266,9 +266,9 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.vdw_cutoff_period[1], 3);
     EXPECT_EQ(param.inp.vdw_cutoff_period[2], 3);
     EXPECT_EQ(std::stod(param.inp.exx_hybrid_alpha), 0.25);
-    EXPECT_EQ(std::stod(param.inp.exx_cam_alpha), 0.0);
-    EXPECT_EQ(std::stod(param.inp.exx_cam_beta), 0.0);
+    EXPECT_EQ(std::stod(param.inp.exx_hybrid_beta), 0.0);
     EXPECT_EQ(param.inp.exx_use_ewald, 0);
+    EXPECT_EQ(param.inp.exx_spencer_type, 0);
     EXPECT_EQ(param.inp.exx_fq_type, 1);
     EXPECT_EQ(param.inp.exx_real_number, "1");
     EXPECT_DOUBLE_EQ(param.inp.exx_hse_omega, 0.11);

source/module_io/test/support/INPUT

@@ -273,11 +273,11 @@ vdw_cutoff_period   3 3 3 #periods of periodic structure
 
 #Parameters (14.exx)
 exx_hybrid_alpha               default #
-exx_cam_alpha                  default #
-exx_cam_beta                   default #
+exx_hybrid_beta                   default #
 exx_hse_omega                  default #
 exx_separate_loop              1 #0 or 1
 exx_use_ewald                  0 
+exx_spencer_type                  0 
 exx_fq_type                    1 
 exx_hybrid_step                100 #
 exx_mixing_beta                1.0 #