Merge branch 'develop' into move_fft

Author: Critsium-xy

docs/advanced/input_files/input-main.md

@@ -2977,7 +2977,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
 - **Type**: Real
 - **Availability**: *[calculation](#calculation)==gen_opt_abfs*
 - **Description**: The threshold when solving for the zeros of spherical Bessel functions. A reasonable choice is 1e-12.
-- **Default**: 0
+- **Default**: 1E-12
 
 ### exx_real_number
 

source/Makefile.Objects

@@ -334,6 +334,8 @@ OBJS_HAMILT=hamilt_pw.o\
     op_exx_pw.o\
     ekinetic_pw.o\
     ekinetic_op.o\
+    exx_pw_ace.o\
+    exx_pw_pot.o\
     hpsi_norm_op.o\
     veff_pw.o\
     veff_op.o\

source/source_cell/pseudo.cpp

@@ -1,5 +1,6 @@
 #include "pseudo.h"
 #include "source_base/tool_title.h"
+#include <cstdint>
 
 pseudo::pseudo()
 {

source/source_esolver/esolver_ks_pw.cpp

@@ -589,10 +589,20 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
                 auto start = std::chrono::high_resolution_clock::now();
                 exx_helper.set_firstiter(false);
                 exx_helper.op_exx->first_iter = false;
+                double dexx = 0.0;
+                if (PARAM.inp.exx_thr_type == "energy")
+                {
+                    dexx = exx_helper.cal_exx_energy(this->kspw_psi);
+                }
                 exx_helper.set_psi(this->kspw_psi);
+                if (PARAM.inp.exx_thr_type == "energy")
+                {
+                    dexx -= exx_helper.cal_exx_energy(this->kspw_psi);
+                    // std::cout << "dexx = " << dexx << std::endl;
+                }
+                bool conv_ene = std::abs(dexx) < PARAM.inp.exx_ene_thr;
 
-                conv_esolver = exx_helper.exx_after_converge(iter);
-
+                conv_esolver = exx_helper.exx_after_converge(iter, conv_ene);
                 if (!conv_esolver)
                 {
                     auto duration = std::chrono::high_resolution_clock::now() - start;
@@ -602,6 +612,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
                     exx_helper.op_exx->first_iter = false;
                     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
                     update_pot(ucell, istep, iter, conv_esolver);
+                    exx_helper.iter_inc();
                 }
             }
         }
@@ -614,8 +625,9 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     //----------------------------------------------------------
     // 3) Print out electronic wavefunctions in pw basis
     //----------------------------------------------------------
-    if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
+    if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax)
     {
+        // conv_esolver == true has already been dealt with in after_scf
         ModuleIO::write_wfc_pw(GlobalV::KPAR,
                                GlobalV::MY_POOL,
                                GlobalV::MY_RANK,

source/source_hamilt/module_xc/exx_info.h

@@ -81,7 +81,7 @@ struct Exx_Info
     {
         int abfs_Lmax = 0; // tmp
         double ecut_exx = 60;
-        double tolerence = 1E-2;
+        double tolerence = 1E-12;
         double kmesh_times = 4;
     };
     Exx_Info_Opt_ABFs info_opt_abfs;

source/source_io/module_parameter/input_parameter.h

@@ -546,7 +546,7 @@ struct Input_para
                                                  ///< calculating Columb potential is to that of atomic orbitals
     int exx_opt_orb_lmax = 0;                    ///< the maximum l of the spherical Bessel functions for opt ABFs
     double exx_opt_orb_ecut = 0.0;               ///< the cut-off of plane wave expansion for opt ABFs
-    double exx_opt_orb_tolerence = 0.0;          ///< the threshold when solving for the zeros of spherical Bessel
+    double exx_opt_orb_tolerence = 1E-12;        ///< the threshold when solving for the zeros of spherical Bessel
                                                  ///< functions for opt ABFs
     bool exx_symmetry_realspace
         = true; ///< whether to reduce the real-space sector in when using symmetry=1 in EXX calculation
@@ -654,6 +654,8 @@ struct Input_para
     // EXX for planewave basis, rhx0820 2025-03-10
     bool exxace = true; // exxace, exact exchange for planewave basis, https://doi.org/10.1021/acs.jctc.6b00092
     bool exx_gamma_extrapolation = true; // gamma point extrapolation for exx, https://doi.org/10.1103/PhysRevB.79.205114
+    std::string exx_thr_type = "density"; // threshold type for exx outer loop, energy or density
+    double exx_ene_thr = 1e-5; // threshold for exx outer loop when exx_thr_type = energy
 
     // ====   #Parameters (23.XC external parameterization) ========
     /*

source/source_io/read_input_item_other.cpp

@@ -545,6 +545,32 @@ void ReadInput::item_others()
         read_sync_bool(input.exx_gamma_extrapolation);
         this->add_item(item);
     }
+    {
+        Input_Item item("exx_thr_type");
+        item.annotation = "threshold type for exx outer loop, energy or density";
+        read_sync_string(input.exx_thr_type);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            std::string thr_type = para.input.exx_thr_type;
+            std::transform(thr_type.begin(), thr_type.end(), thr_type.begin(), ::tolower);
+            if (thr_type != "energy" && thr_type != "density")
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "exx_thr_type should be energy or density");
+            }
+        };
+        this->add_item(item);
+    }
+    {
+        Input_Item item("exx_ene_thr");
+        item.annotation = "threshold for exx outer loop when exx_thr_type = energy";
+        read_sync_double(input.exx_ene_thr);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.exx_ene_thr <= 0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "exx_ene_thr must > 0");
+            }
+        };
+        this->add_item(item);
+    }
 
 }
 } // namespace ModuleIO

source/source_io/test/read_input_ptest.cpp

@@ -296,7 +296,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_DOUBLE_EQ(param.inp.rpa_ccp_rmesh_times, 10.0);
     EXPECT_EQ(param.inp.exx_opt_orb_lmax, 0);
     EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_ecut, 0.0);
-    EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_tolerence, 0.0);
+    EXPECT_DOUBLE_EQ(param.inp.exx_opt_orb_tolerence, 1E-12);
     EXPECT_FALSE(param.inp.noncolin);
     EXPECT_FALSE(param.inp.lspinorb);
     EXPECT_DOUBLE_EQ(param.inp.soc_lambda, 1.0);

source/source_io/test/support/INPUT

@@ -291,7 +291,7 @@ exx_v_grad_r_threshold         0 #
 exx_ccp_rmesh_times            default #
 exx_opt_orb_lmax               0 #
 exx_opt_orb_ecut               0 #
-exx_opt_orb_tolerence          0 #
+exx_opt_orb_tolerence          1E-12 #
 
 #Parameters (16.tddft)
 td_vext                        0 #add extern potential or not

source/source_lcao/module_deepks/deepks_basic.cpp

@@ -61,6 +61,15 @@ void DeePKS_domain::load_model(const std::string& model_file, torch::jit::script
     ModuleBase::TITLE("DeePKS_domain", "load_model");
     ModuleBase::timer::tick("DeePKS_domain", "load_model");
 
+    // check whether file exists
+    std::ifstream ifs(model_file.c_str());
+    if (!ifs)
+    {
+        ModuleBase::timer::tick("DeePKS_domain", "load_model");
+        ModuleBase::WARNING_QUIT("DeePKS_domain::load_model", "No model file named " + model_file + ", please check!");
+        return;
+    }
+    ifs.close();
     try
     {
         model = torch::jit::load(model_file);