fix bug and update code format in exx

PeizeLin · PeizeLin · commit 8f02084bb12d · 2025-05-29T11:43:47.000+08:00
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -183,12 +183,12 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
             // initialize 2-center radial tables for EXX-LRI
             if (GlobalC::exx_info.info_ri.real_number)
             {
-                this->exx_lri_double->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
+                this->exd->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
                 this->exd->exx_before_all_runners(this->kv, ucell, this->pv);
             }
             else
             {
-                this->exx_lri_complex->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
+                this->exc->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
                 this->exc->exx_before_all_runners(this->kv, ucell, this->pv);
             }
         }
diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h
@@ -21,21 +21,21 @@
 #include "module_exx_symmetry/symmetry_rotation.h"
 
 	class Parallel_Orbitals;
-	
+
 	template<typename T, typename Tdata>
 	class RPA_LRI;
 
 	template<typename T, typename Tdata>
 	class Exx_LRI_Interface;
 
-    namespace LR
-    {
-        template<typename T, typename TR>
-        class ESolver_LR;
+	namespace LR
+	{
+		template<typename T, typename TR>
+		class ESolver_LR;
 
-        template<typename T>
-        class OperatorLREXX;
-    }
+		template<typename T>
+		class OperatorLREXX;
+	}
 
 template<typename Tdata>
 class Exx_LRI
@@ -49,37 +49,37 @@ class Exx_LRI
 	using TatomR = std::array<double,Ndim>;		// tmp
 
 public:
-    Exx_LRI(const Exx_Info::Exx_Info_RI& info_in) :info(info_in) {}
-    Exx_LRI operator=(const Exx_LRI&) = delete;
-    Exx_LRI operator=(Exx_LRI&&);
+	Exx_LRI(const Exx_Info::Exx_Info_RI& info_in) :info(info_in) {}
+	Exx_LRI operator=(const Exx_LRI&) = delete;
+	Exx_LRI operator=(Exx_LRI&&);
 
-    void reset_Cs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs_in) { this->exx_lri.set_Cs(Cs_in, this->info.C_threshold); }
-    void reset_Vs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_in) { this->exx_lri.set_Vs(Vs_in, this->info.V_threshold); }
+	void reset_Cs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs_in) { this->exx_lri.set_Cs(Cs_in, this->info.C_threshold); }
+	void reset_Vs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_in) { this->exx_lri.set_Vs(Vs_in, this->info.V_threshold); }
 
-	void init(const MPI_Comm &mpi_comm_in, 
+	void init(const MPI_Comm &mpi_comm_in,
 			  const UnitCell &ucell,
-			  const K_Vectors &kv_in, 
+			  const K_Vectors &kv_in,
 			  const LCAO_Orbitals& orb);
 	void cal_exx_force(const int& nat);
-    void cal_exx_stress(const double& omega, const double& lat0);
+	void cal_exx_stress(const double& omega, const double& lat0);
 	void cal_exx_ions(const UnitCell& ucell, const bool write_cv = false);
-    void cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
+	void cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
 		const UnitCell& ucell,
-        const Parallel_Orbitals& pv,
-        const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr);
-    std::vector<std::vector<int>> get_abfs_nchis() const;
+		const Parallel_Orbitals& pv,
+		const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr);
+	std::vector<std::vector<int>> get_abfs_nchis() const;
 
 	std::vector< std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Hexxs;
-    double Eexx;
+	double Eexx;
 	ModuleBase::matrix force_exx;
 	ModuleBase::matrix stress_exx;
-    
+
 
 private:
 	const Exx_Info::Exx_Info_RI &info;
 	MPI_Comm mpi_comm;
 	const K_Vectors *p_kv = nullptr;
-    std::vector<double> orb_cutoff_;
+	std::vector<double> orb_cutoff_;
 
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
@@ -89,16 +89,16 @@ class Exx_LRI
 	RI::Exx<TA,Tcell,Ndim,Tdata> exx_lri;
 
 	void post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> &Hexxs_io ) const;
-    double post_process_Eexx(const double& Eexx_in) const;
+	double post_process_Eexx(const double& Eexx_in) const;
 
 	friend class RPA_LRI<double, Tdata>;
 	friend class RPA_LRI<std::complex<double>, Tdata>;
 	friend class Exx_LRI_Interface<double, Tdata>;
 	friend class Exx_LRI_Interface<std::complex<double>, Tdata>;
-    friend class LR::ESolver_LR<double, double>;
-    friend class LR::ESolver_LR<std::complex<double>, double>;
-    friend class LR::OperatorLREXX<double>;
-    friend class LR::OperatorLREXX<std::complex<double>>;
+	friend class LR::ESolver_LR<double, double>;
+	friend class LR::ESolver_LR<std::complex<double>, double>;
+	friend class LR::OperatorLREXX<double>;
+	friend class LR::OperatorLREXX<std::complex<double>>;
 };
 
 #include "Exx_LRI.hpp"
diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp
@@ -26,9 +26,9 @@
 #include <string>
 
 template<typename Tdata>
-void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in, 
+void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
 						  const UnitCell &ucell,
-						  const K_Vectors &kv_in, 
+						  const K_Vectors &kv_in,
 						  const LCAO_Orbitals& orb)
 {
 	ModuleBase::TITLE("Exx_LRI","init");
@@ -130,7 +130,7 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
 
 	// std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour.
-	const std::array<Tcell,Ndim> period_Vs = LRI_CV_Tools::cal_latvec_range<Tcell>(1+this->info.ccp_rmesh_times, ucell, orb_cutoff_);	
+	const std::array<Tcell,Ndim> period_Vs = LRI_CV_Tools::cal_latvec_range<Tcell>(1+this->info.ccp_rmesh_times, ucell, orb_cutoff_);
 	const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
 		list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
 
@@ -237,7 +237,7 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, R
 	}
 	this->Eexx = post_process_Eexx(this->Eexx);
 	this->exx_lri.set_symmetry(false, {});
-	ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec");	
+	ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec");
 }
 
 template<typename Tdata>
@@ -283,11 +283,6 @@ void Exx_LRI<Tdata>::cal_exx_force(const int& nat)
 	ModuleBase::TITLE("Exx_LRI","cal_exx_force");
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_force");
 
-    if (!this->exx_lri.flag_finish.D)
-    {
-        ModuleBase::WARNING_QUIT("Force_Stress_LCAO", "Cannot calculate EXX force when the first PBE loop is not converged.");
-    }
-
 	this->force_exx.create(nat, Ndim);
 	for(int is=0; is<PARAM.inp.nspin; ++is)
 	{
diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h
@@ -44,12 +44,19 @@ class Exx_LRI_Interface
     void read_Hexxs_cereal(const std::string& file_name);
 
     std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& get_Hexxs() const { return this->exx_ptr->Hexxs; }
-    
+
     double& get_Eexx() const { return this->exx_ptr->Eexx; }
 
+    // Processes in ESolver_KS_LCAO
+    /// @brief Exx_LRI::init()
+    void init(const MPI_Comm &mpi_comm,
+              const UnitCell &ucell,
+              const K_Vectors &kv,
+              const LCAO_Orbitals& orb);
+
     // Processes in ESolver_KS_LCAO
     /// @brief in before_all_runners: set symmetry according to irreducible k-points
-    /// since k-points are not reduced again after the variation of the cell and exx-symmetry must be consistent with k-points. 
+    /// since k-points are not reduced again after the variation of the cell and exx-symmetry must be consistent with k-points.
     /// In the future, we will reduce k-points again during cell-relax, then this setting can be moved to `exx_beforescf`.
     void exx_before_all_runners(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv);
 
@@ -58,23 +65,23 @@ class Exx_LRI_Interface
 
     /// @brief in eachiterinit:  do DM mixing and calculate Hexx when entering 2nd SCF
     void exx_eachiterinit(const int istep,
-                          const UnitCell& ucell, 
+                          const UnitCell& ucell,
                           const elecstate::DensityMatrix<T, double>& dm/**< double should be Tdata if complex-PBE-DM is supported*/,
-                          const K_Vectors& kv, 
+                          const K_Vectors& kv,
                           const int& iter);
 
     /// @brief in hamilt2density: calculate Hexx and Eexx
     void exx_hamilt2density(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter);
 
     /// @brief in iter_finish: write Hexx, do something according to whether SCF is converged
-    void exx_iter_finish(const K_Vectors& kv, 
+    void exx_iter_finish(const K_Vectors& kv,
                          const UnitCell& ucell,
-                         hamilt::Hamilt<T>& hamilt, 
-                         elecstate::ElecState& elec, 
+                         hamilt::Hamilt<T>& hamilt,
+                         elecstate::ElecState& elec,
                          Charge_Mixing& chgmix,
-                         const double& scf_ene_thr, 
-                         int& iter, 
-                         const int istep, 
+                         const double& scf_ene_thr,
+                         int& iter,
+                         const int istep,
                          bool& conv_esolver);
     /// @brief: in do_after_converge: add exx operators; do DM mixing if seperate loop
     bool exx_after_converge(const UnitCell& ucell,
@@ -86,6 +93,13 @@ class Exx_LRI_Interface
                             const int& istep,
                             const double& etot,
                             const double& scf_ene_thr);
+
+    /// @brief: in cal_exx_force: Exx_LRI::cal_exx_force()
+    void cal_exx_force(const double& omega, const double& lat0);
+
+    /// @brief: in cal_exx_stress: Exx_LRI::cal_exx_stress()
+    void cal_exx_stress(const double& omega, const double& lat0);
+
     int two_level_step = 0;
     double etot_last_outer_loop = 0.0;
     elecstate::DensityMatrix<T, double>* dm_last_step;
@@ -95,6 +109,16 @@ class Exx_LRI_Interface
 
     bool exx_spacegroup_symmetry = false;
     ModuleSymmetry::Symmetry_rotation symrot_;
+
+    struct Flag_Finish
+    {
+        bool init = false;
+        bool ions = false;
+        bool elec = false;
+        bool force = false;
+        bool stress = false;
+    };
+    Flag_Finish flag_finish;
 };
 
 #include "Exx_LRI_interface.hpp"
diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp
diff --git a/source/module_ri/Mix_Matrix.cpp b/source/module_ri/Mix_Matrix.cpp

Original file line number	Diff line number	Diff line change
`@@ -183,12 +183,12 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa`
`183`	`183`	`// initialize 2-center radial tables for EXX-LRI`
`184`	`184`	`if (GlobalC::exx_info.info_ri.real_number)`
`185`	`185`	`{`
`186`		`- this->exx_lri_double->init(MPI_COMM_WORLD, ucell, this->kv, orb_);`
	`186`	`+ this->exd->init(MPI_COMM_WORLD, ucell, this->kv, orb_);`
`187`	`187`	`this->exd->exx_before_all_runners(this->kv, ucell, this->pv);`
`188`	`188`	`}`
`189`	`189`	`else`
`190`	`190`	`{`
`191`		`- this->exx_lri_complex->init(MPI_COMM_WORLD, ucell, this->kv, orb_);`
	`191`	`+ this->exc->init(MPI_COMM_WORLD, ucell, this->kv, orb_);`
`192`	`192`	`this->exc->exx_before_all_runners(this->kv, ucell, this->pv);`
`193`	`193`	`}`
`194`	`194`	`}`