Feature: support EXX PW (screened hybrid) NSCF

docs/advanced/elec_properties/band.md

@@ -1,5 +1,7 @@
 # Extracting Band Structure
 
+For plane-wave screened hybrid calculations, follow the [hybrid NSCF workflow](hybrid_nscf.md), which additionally requires the frozen SCF orbitals and source q mesh.
+
 In ABACUS, in order to obtain the eigenvalues of Hamiltonian, or generally called band structure, examples can be found in [examples/band](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/12_band).
 Similar to the [DOS case](https://abacus-rtd.readthedocs.io/en/latest/advanced/elec_properties/dos.html), one first needs to perform a ground-state energy calculation ***with one additional keyword "[out_chg](https://abacus-rtd.readthedocs.io/en/latest/advanced/input_files/input-main.html#out-chg)" in the INPUT file***:
 

docs/advanced/elec_properties/hybrid_nscf.md

@@ -0,0 +1,69 @@
+# PW screened hybrid NSCF
+
+Hybrid NSCF reads the charge density, binary `wfk*_pw.dat` wavefunctions and
+weighted occupations in `eig_occ.txt` from a converged SCF calculation. It keeps
+that source fixed while solving a new K-point list or additional target bands.
+No new file format is required.
+
+## SCF
+
+Use a complete, uniformly weighted source mesh and:
+
+```text
+calculation scf
+suffix hybrid
+basis_type pw
+dft_functional HSE
+smearing_method fixed
+symmetry -1
+kpar 1
+bndpar 1
+exxace false
+exx_gamma_extra false
+out_chg 1
+out_wfc_pw 2
+out_freq_ion 0
+cal_force 0
+cal_stress 0
+```
+
+## NSCF
+
+Keep the same structure, pseudopotentials and cutoffs. In a separate output
+directory, retain the shared settings above, replace `KPT` with the target
+points, and change:
+
+```text
+calculation nscf
+read_file_dir ../scf/OUT.hybrid
+init_wfc random
+nbands 20
+out_band 1
+out_wfc_pw 0
+out_chg 0
+```
+
+Supports CPU/CUDA and `nspin 1/2`; use one MPI rank for CUDA. ACE, SOC,
+unscreened exchange, forces and stress are unsupported in this workflow.
+Set `symmetry -1`: input `0` is accepted but PW EXX resets it to `-1` at runtime.
+
+HSE defaults to `exx_singularity_correction limits`, as in LCAO. Set `gygi` to
+retain the historical PW correction, which can produce jumps at source K points.
+Use matching SCF/NSCF settings for consistent bands. Do not use Methfessel-Paxton
+smearing (`mp`/`mp2`/`mp3`) for the source: it can write occupations outside
+`[0, k-point weight]`, which hybrid NSCF rejects; the non-negative schemes
+(`fixed`, `gauss`, `mv`, `fd`) are accepted. Configuration differences
+or missing `INPUT.info` produce warnings; missing or incompatible numerical
+source data produces an error.
+
+## Verification
+
+With MPI and LibXC enabled, substitute your configured executable/build paths:
+
+```bash
+cmake --build build --target MODULE_ESOLVER_pw_hybrid_source
+OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_ESOLVER_pw_hybrid_source$'
+# Frozen-source band regressions: tests/01_PW/nscf_hse_same and nscf_hse_path
+cd tests/01_PW
+OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a <abacus> -r '^nscf_hse'
+```

docs/advanced/elec_properties/index.rst

@@ -4,6 +4,7 @@ Electronic Properties and Outputs
 
 .. toctree::
    band
+   hybrid_nscf
    dos
    Mulliken
    potential

source/Makefile.Objects

@@ -334,6 +334,8 @@ OBJS_ESOLVER=esolver_factory.o\
     esolver_of_tool.o\
     esolver_dfpt_pw.o\
     pw_others.o\
+    pw_hybrid_nscf.o\
+    pw_hybrid_source.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\

source/source_esolver/CMakeLists.txt

@@ -3,6 +3,8 @@ list(APPEND objects
     esolver_ks.cpp
     esolver_fp.cpp
     esolver_ks_pw.cpp
+    pw_hybrid_nscf.cpp
+    pw_hybrid_source.cpp
     esolver_ks_pw_tddft.cpp
     esolver_ks_lcaopw.cpp
     esolver_sdft_pw.cpp

source/source_esolver/esolver_ks_pw.cpp

@@ -1,4 +1,5 @@
 #include "esolver_ks_pw.h"
+#include "pw_hybrid_nscf.h"
 
 #include "source_estate/elecstate_pw.h"
 #include "source_estate/module_charge/chg_symm.h"
@@ -82,6 +83,8 @@ void ESolver_KS_PW<T, Device>::before_all_runners(BaseCell& basecell, const Inpu
 
     ESolver_KS::before_all_runners(ucell, inp);
 
+    validate_hybrid_nscf(inp, this->general_exx_info_);
+
     //! setup and allocation for pelec, potentials, etc.
     elecstate::setup_estate_pw(ucell,
                                this->kv,
@@ -193,6 +196,13 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     //! Setup EXX helper for Hamiltonian and psi
     exx_helper->before_scf(this->p_hamilt, this->stp.template get_psi_t<T, Device>(), *this->inp_, this->general_exx_info_);
 
+    if (this->inp_->calculation == "nscf" && this->general_exx_info_.cal_exx)
+    {
+        const Input_para& inp = *this->inp_;
+        const std::string& readin_dir = PARAM.globalv.global_readin_dir;
+        ModuleESolver::prepare_exx_nscf<T, Device>(inp, ucell, this->pw_wfc, this->exx_helper, readin_dir);
+    }
+
     ModuleBase::timer::end("ESolver_KS_PW", "before_scf");
 }
 
@@ -315,7 +325,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     // Related to EXX
     bool cal_exx = general_exx_info_.cal_exx;
     double hybrid_alpha = general_exx_info_.hybrid_alpha;
-    if (cal_exx && !exx_helper->get_op_first_iter())
+    if (cal_exx && this->inp_->calculation != "nscf" && !exx_helper->get_op_first_iter())
     {
         this->pelec->set_exx(exx_helper->cal_exx_energy(this->stp.template get_psi_t<T, Device>()), cal_exx, hybrid_alpha);
     }