Merge branch 'develop' into dev-1

Author: zgn-26714

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_base/module_device/device.h

@@ -65,8 +65,8 @@ std::string get_device_flag(const std::string& device,
 /**
  * @brief Get the rank of current node
  *        Note that GPU can only be binded with CPU in the same node
- * 
- * @return int 
+ *
+ * @return int
  */
 int get_node_rank();
 int get_node_rank_with_mpi_shared(const MPI_Comm mpi_comm = MPI_COMM_WORLD);
@@ -91,6 +91,14 @@ void record_device_memory(const Device* dev, std::ofstream& ofs_device, std::str
     return;
 }
 
+#if defined(__CUDA) || defined(__ROCM)
+template <>
+void print_device_info<base_device::DEVICE_GPU>(const base_device::DEVICE_GPU *ctx, std::ofstream &ofs_device);
+
+template <>
+void record_device_memory<base_device::DEVICE_GPU>(const base_device::DEVICE_GPU* dev, std::ofstream& ofs_device, std::string str, size_t size);
+#endif
+
 } // end of namespace information
 } // end of namespace base_device
 

source/source_base/module_device/output_device.cpp

@@ -190,7 +190,7 @@ void print_device_info<base_device::DEVICE_GPU>(
     ofs_device << "Detected " << deviceCount << " CUDA Capable device(s)\n";
   }
   int dev = 0, driverVersion = 0, runtimeVersion = 0;
-  cudaErrcheck(cudaSetDevice(dev));
+  cudaErrcheck(cudaGetDevice(&dev));
   cudaDeviceProp deviceProp;
   cudaErrcheck(cudaGetDeviceProperties(&deviceProp, dev));
   ofs_device << "\nDevice " << dev << ":\t " << deviceProp.name << std::endl;
@@ -429,7 +429,7 @@ void print_device_info<base_device::DEVICE_GPU>(
     ofs_device << "Detected " << deviceCount << " CUDA Capable device(s)\n";
   }
   int dev = 0, driverVersion = 0, runtimeVersion = 0;
-  hipErrcheck(hipSetDevice(dev));
+  hipErrcheck(hipGetDevice(&dev));
   hipDeviceProp_t deviceProp;
   hipErrcheck(hipGetDeviceProperties(&deviceProp, dev));
   ofs_device << "\nDevice " << dev << ":\t " << deviceProp.name << std::endl;

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);