Fix minor bug of EXX operator and default value (#5252)

* allow calling contributeHR while type::k in opexx

* fix the default value of exx parameters

* move the test code to the right place

* fix test case

Author: maki49

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp

@@ -201,12 +201,14 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
 {
     ModuleBase::TITLE("OperatorEXX", "contributeHR");
     // 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 
-}
+    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 
+    if (this->add_hexx_type == Add_Hexx_Type::k) { return; }
+
     if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5)
     {
         // add H(R) normally
-        if (GlobalC::exx_info.info_ri.real_number) {
+        if (GlobalC::exx_info.info_ri.real_number)
+        {
             RI_2D_Comm::add_HexxR(
                 this->current_spin,
                 GlobalC::exx_info.info_global.hybrid_alpha,
@@ -215,7 +217,9 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
                 PARAM.globalv.npol,
                 *this->hR,
                 this->use_cell_nearest ? &this->cell_nearest : nullptr);
-        } else {
+        }
+        else
+        {
             RI_2D_Comm::add_HexxR(
                 this->current_spin,
                 GlobalC::exx_info.info_global.hybrid_alpha,
@@ -224,19 +228,20 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
                 PARAM.globalv.npol,
                 *this->hR,
                 this->use_cell_nearest ? &this->cell_nearest : nullptr);
-}
+        }
     }
-    if (PARAM.inp.nspin == 2) { this->current_spin = 1 - this->current_spin;
-}
+    if (PARAM.inp.nspin == 2) { this->current_spin = 1 - this->current_spin; }
 }
 
 template<typename TK, typename TR>
 void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
 {
     ModuleBase::TITLE("OperatorEXX", "constributeHR");
     // 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 
-}
+    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 
+
+    if (this->add_hexx_type == Add_Hexx_Type::R) { throw std::invalid_argument("Set Add_Hexx_Type::k sto call OperatorEXX::contributeHk()."); }
+
     if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5)
     {
         if (this->restart && this->two_level_step != nullptr)

source/module_io/read_input_item_exx_dftu.cpp

@@ -88,7 +88,7 @@ void ReadInput::item_exx()
         read_sync_string(input.exx_real_number);
         item.reset_value = [](const Input_Item& item, Parameter& para) {
             if (para.input.exx_real_number == "default")
-            {
+            {  // to run through here, the default value of para.input.exx_real_number should be "default"
                 if (para.input.gamma_only)
                 {
                     para.input.exx_real_number = "1";
@@ -181,7 +181,7 @@ void ReadInput::item_exx()
         read_sync_string(input.exx_ccp_rmesh_times);
         item.reset_value = [](const Input_Item& item, Parameter& para) {
             if (para.input.exx_ccp_rmesh_times == "default")
-            {
+            {   // to run through here, the default value of para.input.exx_ccp_rmesh_times should be "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);

source/module_io/write_vxc.hpp

@@ -284,10 +284,10 @@ void write_Vxc(const int nspin,
             vexxonly_op_ao.contributeHk(ik);
             std::vector<TK> vexx_k_mo = cVc(vexxonly_k_ao.get_hk(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d);
             e_orb_exx.emplace_back(orbital_energy(ik, nbands, vexx_k_mo, p2d));
+            // ======test=======
+            // exx_energy += all_band_energy(ik, vexx_k_mo, p2d, wg);
+            // ======test=======
         }
-        // ======test=======
-        // exx_energy += all_band_energy(ik, vexx_k_mo, p2d, wg);
-        // ======test=======
 #endif
         if (PARAM.inp.dft_plus_u)
         {

source/module_parameter/input_parameter.h

@@ -469,7 +469,7 @@ struct Input_para
     double exx_mixing_beta = 1.0;               ///< mixing_beta for outer-loop when exx_separate_loop=1
     double exx_lambda = 0.3;                    ///< used to compensate for divergence points at G=0 in the
                                                 ///< evaluation of Fock exchange using lcao_in_pw method
-    std::string exx_real_number = "0";          ///< exx calculated in real or complex
+    std::string exx_real_number = "default";          ///< exx calculated in real or complex
     double exx_pca_threshold = 0.0001;          ///< threshold to screen on-site ABFs in exx
     double exx_c_threshold = 0.0001;            ///< threshold to screen C matrix in exx
     double exx_v_threshold = 0.1;               ///< threshold to screen C matrix in exx
@@ -484,7 +484,7 @@ struct Input_para
                                                 ///< inequality
     double exx_cauchy_stress_threshold = 1e-07; ///< threshold to screen exx stress using Cauchy-Schwartz
                                                 ///< inequality
-    std::string exx_ccp_rmesh_times = "1";      ///< how many times larger the radial mesh required for
+    std::string exx_ccp_rmesh_times = "default";      ///< how many times larger the radial mesh required for
                                                 ///< calculating Columb potential is to that of atomic orbitals
     std::string exx_distribute_type = "htime";  ///< distribute type (assuming default as no specific value
                                                 ///< provided)

tests/integrate/392_NO_GO_LR_HF/INPUT

@@ -28,6 +28,7 @@ mixing_gg0 0
 
 dft_functional hf
 exx_real_number 1
+exx_ccp_rmesh_times 1
 
 lr_nstates 2
 xc_kernel hf