Merge branch 'develop' into toolchain-202501

Author: QuantumMisaka

examples/hse/pw_Si2/INPUT

@@ -0,0 +1,32 @@
+INPUT_PARAMETERS
+pseudo_dir              ../../../tests/PP_ORB
+orbital_dir             ../../../tests/PP_ORB
+nbands                  4
+nspin                   1
+calculation             scf
+basis_type              pw
+ks_solver               dav
+ecutwfc                 50
+scf_thr                 1e-9
+scf_nmax                100
+gamma_only              0
+symmetry                -1
+smearing_method         fixed
+mixing_type             broyden
+mixing_beta             0.7
+
+dft_functional          hse
+
+# init_wfc file
+# init_chg file
+
+pseudo_mesh            1
+pseudo_rcut            10
+
+# out_wfc_pw 1
+# out_chg 1
+
+exx_hybrid_alpha        0.25
+
+cal_stress 1
+# cal_force 1

examples/hse/pw_Si2/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+3 3 3 0 0 0

examples/hse/pw_Si2/STRU

@@ -0,0 +1,24 @@
+ATOMIC_SPECIES
+Si	1	Si_ONCV_PBE_FR-1.1.upf					#Pseudopotentials are downloaded from http://www.quantum-simulation.org/potentials/sg15_oncv/upf/
+
+NUMERICAL_ORBITAL
+orb_Si.dat
+
+LATTICE_CONSTANT
+1.889766
+
+LATTICE_VECTORS
+0.0	2.708337	2.708337
+2.708337	0.0	2.708337
+2.708337	2.708337	0.0
+
+
+
+ATOMIC_POSITIONS
+Direct
+
+Si
+0
+2
+0.125	0.125	0.125	0 0 0
+0.875	0.875	0.875	0 0 0

source/Makefile.Objects

@@ -50,6 +50,7 @@ VPATH=./src_global:\
 ./module_hamilt_pw/hamilt_stodft:\
 ./module_hamilt_pw/hamilt_pwdft/operator_pw:\
 ./module_hamilt_pw/hamilt_pwdft/kernels:\
+./module_hamilt_pw/hamilt_pwdft/module_exx_helper:\
 ./module_hamilt_pw/hamilt_stodft/kernels:\
 ./module_hamilt_lcao/module_hcontainer:\
 ./module_hamilt_lcao/hamilt_lcaodft:\
@@ -222,6 +223,7 @@ OBJS_ELECSTAT=elecstate.o\
     elecstate_pw_sdft.o\
     elecstate_pw_cal_tau.o\
     elecstate_op.o\
+    elecstate_tools.o\
     efield.o\
     gatefield.o\
     potential_new.o\
@@ -294,6 +296,7 @@ OBJS_HAMILT=hamilt_pw.o\
     hamilt_sdft_pw.o\
     operator.o\
     operator_pw.o\
+    op_exx_pw.o\
     ekinetic_pw.o\
     ekinetic_op.o\
     hpsi_norm_op.o\
@@ -305,6 +308,7 @@ OBJS_HAMILT=hamilt_pw.o\
     meta_op.o\
     velocity_pw.o\
     radial_proj.o\
+    exx_helper.o\
 
 OBJS_HAMILT_OF=kedf_tf.o\
     kedf_vw.o\
@@ -670,6 +674,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     sto_stress_pw.o\
     stress_func_cc.o\
     stress_func_ewa.o\
+    stress_func_exx.o\
     stress_func_gga.o\
     stress_func_mgga.o\
     stress_func_har.o\

source/module_elecstate/CMakeLists.txt

@@ -4,6 +4,7 @@ list(APPEND objects
     elecstate_energy.cpp
     elecstate_exx.cpp
     elecstate_print.cpp
+    elecstate_tools.cpp
     elecstate_pw.cpp
     elecstate_pw_sdft.cpp
     elecstate_pw_cal_tau.cpp

source/module_elecstate/elecstate.cpp

@@ -15,43 +15,6 @@ const double* ElecState::getRho(int spin) const
     return &(this->charge->rho[spin][0]);
 }
 
-void ElecState::fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec)
-{
-    assert(nbands > 0);
-    assert(nelec > 0.0);
-
-    const double ne_thr = 1.0e-5;
-
-    const int num = this->klist->get_nks() * nbands;
-    if (num != ocp_kb.size())
-    {
-        ModuleBase::WARNING_QUIT("ElecState::fixed_weights",
-                                 "size of occupation array is wrong , please check ocp_set");
-    }
-
-    double num_elec = 0.0;
-    for (int i = 0; i < ocp_kb.size(); ++i)
-    {
-        num_elec += ocp_kb[i];
-    }
-
-    if (std::abs(num_elec - nelec) > ne_thr)
-    {
-        ModuleBase::WARNING_QUIT("ElecState::fixed_weights",
-                                 "total number of occupations is wrong , please check ocp_set");
-    }
-
-    for (int ik = 0; ik < this->wg.nr; ++ik)
-    {
-        for (int ib = 0; ib < this->wg.nc; ++ib)
-        {
-            this->wg(ik, ib) = ocp_kb[ik * this->wg.nc + ib];
-        }
-    }
-    this->skip_weights = true;
-
-    return;
-}
 
 
 void ElecState::init_nelec_spin()
@@ -64,143 +27,6 @@ void ElecState::init_nelec_spin()
     }
 }
 
-
-void ElecState::calculate_weights()
-{
-    ModuleBase::TITLE("ElecState", "calculate_weights");
-    if (this->skip_weights)
-    {
-        return;
-    }
-
-    const int nbands = this->ekb.nc;
-    const int nks = this->ekb.nr;
-
-    if (!Occupy::use_gaussian_broadening && !Occupy::fixed_occupations)
-    {
-        if (PARAM.globalv.two_fermi)
-        {
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[0],
-                             this->ekb,
-                             this->eferm.ef_up,
-                             this->wg,
-                             0,
-                             this->klist->isk);
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[1],
-                             this->ekb,
-                             this->eferm.ef_dw,
-                             this->wg,
-                             1,
-                             this->klist->isk);
-            // ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
-        }
-        else
-        {
-            // -1 means don't need to consider spin.
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             PARAM.inp.nelec,
-                             this->ekb,
-                             this->eferm.ef,
-                             this->wg,
-                             -1,
-                             this->klist->isk);
-        }
-    }
-    else if (Occupy::use_gaussian_broadening)
-    {
-        if (PARAM.globalv.two_fermi)
-        {
-            double demet_up = 0.0;
-            double demet_dw = 0.0;
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[0],
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef_up,
-                             demet_up,
-                             this->wg,
-                             0,
-                             this->klist->isk);
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[1],
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef_dw,
-                             demet_dw,
-                             this->wg,
-                             1,
-                             this->klist->isk);
-            this->f_en.demet = demet_up + demet_dw;
-        }
-        else
-        {
-            // -1 means is no related to spin.
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             PARAM.inp.nelec,
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef,
-                             this->f_en.demet,
-                             this->wg,
-                             -1,
-                             this->klist->isk);
-        }
-#ifdef __MPI
-        const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-        Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, this->f_en.demet);
-#endif
-    }
-    else if (Occupy::fixed_occupations)
-    {
-        ModuleBase::WARNING_QUIT("calculate_weights", "other occupations, not implemented");
-    }
-
-    return;
-}
-
-
-void ElecState::calEBand()
-{
-    ModuleBase::TITLE("ElecState", "calEBand");
-    // calculate ebands using wg and ekb
-    double eband = 0.0;
-#ifdef _OPENMP
-#pragma omp parallel for collapse(2) reduction(+ : eband)
-#endif
-    for (int ik = 0; ik < this->ekb.nr; ++ik)
-    {
-        for (int ibnd = 0; ibnd < this->ekb.nc; ibnd++)
-        {
-            eband += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
-        }
-    }
-    this->f_en.eband = eband;
-
-#ifdef __MPI
-    const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, this->f_en.eband);
-#endif
-    return;
-}
-
-
 void ElecState::init_scf(const int istep, 
                          const UnitCell& ucell,
                          const Parallel_Grid& pgrid,

source/module_elecstate/elecstate.h

@@ -72,12 +72,11 @@ class ElecState
         return;
     }
 
-    // calculate wg from ekb
-    virtual void calculate_weights();
+
 
     // use occupied weights from INPUT and skip calculate_weights
     // mohan updated on 2024-06-08
-    void fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec);
+    
 
     // if nupdown is not 0(TWO_EFERMI case),
     // nelec_spin will be fixed and weights will be constrained
@@ -132,12 +131,8 @@ class ElecState
     bool vnew_exist = false;
     void cal_converged();
     void cal_energies(const int type);
-#ifdef __EXX
-#ifdef __LCAO
     void set_exx(const double& Eexx);
     void set_exx(const std::complex<double>& Eexx);
-#endif //__LCAO
-#endif //__EXX
 
     double get_hartree_energy();
     double get_etot_efield();
@@ -167,11 +162,11 @@ class ElecState
     ModuleBase::matrix wg;  ///< occupation weight for each k-point and band
 
   public:
-    // calculate ebands for all k points and all occupied bands
-    void calEBand();
 
     bool skip_weights = false;
 };
+    
+
 
 } // namespace elecstate
 #endif

source/module_elecstate/elecstate_exx.cpp

@@ -3,8 +3,6 @@
 namespace elecstate
 {
 
-#ifdef __EXX
-#ifdef __LCAO
 /// @brief calculation if converged
 /// @date Peize Lin add 2016-12-03
 void ElecState::set_exx(const double& Eexx)
@@ -17,7 +15,5 @@ void ElecState::set_exx(const double& Eexx)
     }
     return;
 }
-#endif //__LCAO
-#endif //__EXX
 
 }