Flying-dragon-boxing · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026
diff --git a/docs/advanced/elec_properties/band.md b/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***:
 

diff --git a/docs/advanced/elec_properties/hybrid_nscf.md b/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_PW_pw_hybrid_source
+OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_PW_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'
+```
diff --git a/docs/advanced/elec_properties/index.rst b/docs/advanced/elec_properties/index.rst
@@ -4,6 +4,7 @@ Electronic Properties and Outputs
 
 .. toctree::
    band
+   hybrid_nscf
    dos
    Mulliken
    potential

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -373,6 +373,8 @@ OBJS_GINT=batch_biggrid.o\
       unitcell_info.o\
 
 OBJS_HAMILT=hamilt_pw.o\
+    pw_hybrid_nscf.o\
+    pw_hybrid_source.o\
     td_history_pw.o\
     td_pw.o\
     hs_matrix_k.o\

diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp
@@ -1,4 +1,5 @@
 #include "esolver_ks_pw.h"
+#include "source_pw/module_pwdft/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);
 
+    pw::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,14 @@ 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;
+        auto* helper = static_cast<Exx_Helper<T, Device>*>(this->exx_helper);
+        pw::prepare_exx_nscf<T, Device>(inp, ucell, this->pw_wfc, *helper->op_exx, readin_dir);
+    }
+
     ModuleBase::timer::end("ESolver_KS_PW", "before_scf");
 }
 
@@ -315,7 +326,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);
     }

diff --git a/source/source_hamilt/module_xc/exx_info.cpp b/source/source_hamilt/module_xc/exx_info.cpp
@@ -114,7 +114,8 @@ bool init_general_exx_info(General_Exx_Info& info, const Input_para& inp)
                 {
                     info.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
                         {"alpha", ModuleBase::GlobalFunc::TO_STRING(erfc_alpha[i])},
-                        {"omega", ModuleBase::GlobalFunc::TO_STRING(inp.exx_erfc_omega[i])} }};
+                        {"omega", ModuleBase::GlobalFunc::TO_STRING(inp.exx_erfc_omega[i])},
+                        {"singularity_correction", inp.exx_singularity_correction} }};
                 }
             }
         }

diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp
@@ -466,14 +466,19 @@ void ReadInput::item_exx()
     {
         Input_Item item("exx_singularity_correction");
         item.annotation = "set the scheme of Coulomb singularity correction";
-        item.category = "Exact Exchange (LCAO)";
+        item.category = "Exact Exchange (Common)";
         item.type = "String";
-        item.description = R"(* spencer: see Phys. Rev. B 77, 193110 (2008).
-* revised_spencer: see Phys. Rev. Mater. 5, 013807 (2021). Set the scheme of Coulomb singularity correction.)";
+        item.description = R"(Scheme for Coulomb singularity / zero-transfer treatment in exact exchange.
+* limits: use the finite screened Coulomb-kernel limit, without an auxiliary correction. Default for HSE and cwp22 in both PW and LCAO. PW requires screened exchange and exx_gamma_extra=false; the latter defaults to false for this scheme.
+* gygi: retain the historical PW Gygi-Baldereschi auxiliary-function correction. PW only; default for unscreened PW hybrids. For independent-k screened hybrid NSCF this can cause discontinuities at source k points.
+* spencer: LCAO spherical truncation, Phys. Rev. B 77, 193110 (2008).
+* revised_spencer: LCAO, Phys. Rev. Mater. 5, 013807 (2021).
+* massidda / carrier: LCAO auxiliary singularity corrections.
+Use the same scheme for consistent SCF and NSCF bands. A source configuration mismatch produces a warning and does not prevent reading the existing SCF files.)";
         item.default_value = "default";
         item.unit = "";
         read_sync_string(input.exx_singularity_correction);
-        item.reset_value = [](const Input_Item& item, Parameter& para) {
+        item.reset_value = [this](const Input_Item& item, Parameter& para) {
             if (para.input.exx_singularity_correction == "default")
             {  
                 std::string& dft_functional = para.input.dft_functional;
@@ -490,13 +495,58 @@ void ReadInput::item_exx()
                     || dft_functional_lower == "lrc_wpbeh"
                     || dft_functional_lower == "cam_pbeh")
                 {
-                    para.input.exx_singularity_correction = "spencer";
+                    if (para.input.basis_type == "pw")
+                    {
+                        para.input.exx_singularity_correction = "gygi";
+                    }
+                    else
+                    {
+                        para.input.exx_singularity_correction = "spencer";
+                    }
                 }
                 else if (dft_functional_lower == "hse" || dft_functional_lower == "cwp22")
                 {
                     para.input.exx_singularity_correction = "limits";
                 }
             }
+            if (para.input.basis_type == "pw" && para.input.exx_singularity_correction == "limits")
+            {
+                const auto gamma_item = std::find_if(this->input_lists.begin(), this->input_lists.end(),
+                    [](const std::pair<std::string, Input_Item>& entry) {
+                        return entry.first == "exx_gamma_extra";
+                    });
+                if (gamma_item != this->input_lists.end() && !gamma_item->second.is_read())
+                {
+                    para.input.exx_gamma_extra = false;
+                }
+            }
+        };
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            const std::string& scheme = para.input.exx_singularity_correction;
+            if (para.input.basis_type == "pw" && scheme != "default")
+            {
+                if (scheme != "limits" && scheme != "gygi")
+                {
+                    ModuleBase::WARNING_QUIT("ReadInput", "PW exx_singularity_correction must be limits or gygi");
+                }
+                if (scheme == "limits")
+                {
+                    bool has_fock = false;
+                    for (const auto& alpha : para.input.exx_fock_alpha)
+                    {
+                        has_fock = has_fock || std::stod(alpha) != 0.0;
+                    }
+                    if (has_fock || para.input.exx_gamma_extra)
+                    {
+                        ModuleBase::WARNING_QUIT("ReadInput",
+                            "PW limits requires screened exchange and exx_gamma_extra=false");
+                    }
+                }
+            }
+            else if (para.input.basis_type != "pw" && scheme == "gygi")
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "exx_singularity_correction gygi is only supported for PW");
+            }
         };
         this->add_item(item);
     }

diff --git a/source/source_io/module_parameter/read_inp_other.cpp b/source/source_io/module_parameter/read_inp_other.cpp
@@ -886,8 +886,8 @@ When false (default), both the direction and magnitude of the magnetic moment ar
         item.annotation = "whether to perform gamma extrapolation in exxpw";
         item.category = "Exact Exchange (PW)";
         item.type = "Boolean";
-        item.description = "Whether to use the gamma point extrapolation method to calculate the Fock exchange operator. See https://doi.org/10.1103/PhysRevB.79.205114 for details. Should be set to true most of the time.";
-        item.default_value = "True";
+        item.description = "Whether to use the gamma point extrapolation method to calculate the Fock exchange operator. See https://doi.org/10.1103/PhysRevB.79.205114 for details. Defaults to false for PW exx_singularity_correction=limits (finite screened kernel), and true otherwise. Explicit true is incompatible with the PW limits scheme.";
+        item.default_value = "False for PW limits; True otherwise";
         item.unit = "";
         read_sync_bool(input.exx_gamma_extra);
         this->add_item(item);

diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp
@@ -151,7 +151,7 @@ void ReadInput::item_system()
         item.description = R"(Specify the type of calculation.
 
 * scf: perform self-consistent electronic structure calculations
-* nscf: perform non-self-consistent electronic structure calculations. A charge density file is required
+* nscf: perform non-self-consistent electronic structure calculations using a saved charge density. PW screened hybrids also require SCF binary wavefunctions and eig_occ.txt; see [hybrid NSCF](../elec_properties/hybrid_nscf.md).
 * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations
 * cell-relax: perform cell relaxation calculations
 * md: perform molecular dynamics simulations

diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp
@@ -1776,6 +1776,7 @@ TEST_F(InputTest, Item_test2)
     }
     { // exx_singluarity_correction
         auto it = find_label("exx_singularity_correction", readinput.input_lists);
+        TestParameters::input(param).basis_type = "lcao";
         TestParameters::input(param).exx_singularity_correction = "default";
         TestParameters::input(param).dft_functional = "HF";
         it->second.reset_value(it->second, param);
@@ -1790,6 +1791,27 @@ TEST_F(InputTest, Item_test2)
         TestParameters::input(param).dft_functional = "HSE";
         it->second.reset_value(it->second, param);
         EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "limits");
+
+        TestParameters::input(param).basis_type = "pw";
+        TestParameters::input(param).dft_functional = "PBE0";
+        TestParameters::input(param).exx_singularity_correction = "default";
+        it->second.reset_value(it->second, param);
+        EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "gygi");
+
+        TestParameters::input(param).dft_functional = "HSE";
+        TestParameters::input(param).exx_singularity_correction = "default";
+        TestParameters::input(param).exx_gamma_extra = true;
+        it->second.reset_value(it->second, param);
+        EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "limits");
+        EXPECT_FALSE(TestParameters::input(param).exx_gamma_extra);
+
+        auto gamma = find_label("exx_gamma_extra", readinput.input_lists);
+        gamma->second.str_values = {"true"};
+        TestParameters::input(param).exx_gamma_extra = true;
+        it->second.reset_value(it->second, param);
+        EXPECT_TRUE(TestParameters::input(param).exx_gamma_extra);
+        gamma->second.str_values.clear();
+        TestParameters::input(param).exx_gamma_extra = false;
     }
     { // exx_ccp_rmesh_times
         auto it = find_label("exx_ccp_rmesh_times", readinput.input_lists);

diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt
@@ -2,6 +2,8 @@
 
 list(APPEND objects
     hamilt_pw.cpp
+    pw_hybrid_nscf.cpp
+    pw_hybrid_source.cpp
     td_history_pw.cpp
     td_pw.cpp
     op_pw_ekin.cpp

diff --git a/source/source_pw/module_pwdft/exx_helper.cpp b/source/source_pw/module_pwdft/exx_helper.cpp
@@ -12,7 +12,7 @@ template <typename T, typename Device>
 void Exx_Helper<T, Device>::init(const UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& wg, const General_Exx_Info& exx_info)
 {
     if (inp.calculation != "scf" && inp.calculation != "relax"
-        && inp.calculation != "cell-relax" && inp.calculation != "md")
+        && inp.calculation != "cell-relax" && inp.calculation != "md" && inp.calculation != "nscf")
     {
         return;
     }
@@ -22,7 +22,7 @@ void Exx_Helper<T, Device>::init(const UnitCell& ucell, const Input_para& inp, c
         return;
     }
 
-    if (exx_info.separate_loop)
+    if (exx_info.separate_loop && inp.calculation != "nscf")
     {
         XC_Functional::set_xc_first_loop(ucell);
         this->set_firstiter();
@@ -36,7 +36,7 @@ void Exx_Helper<T, Device>::before_scf(void* p_hamilt, void* psi, const Input_pa
 {
     /// Return if not a valid calculation type
     if (inp.calculation != "scf" && inp.calculation != "relax"
-        && inp.calculation != "cell-relax" && inp.calculation != "md")
+        && inp.calculation != "cell-relax" && inp.calculation != "md" && inp.calculation != "nscf")
     {
         return;
     }
@@ -60,6 +60,11 @@ bool Exx_Helper<T, Device>::iter_finish(void* p_elec, Charge* p_charge, void* ps
                                         UnitCell& ucell, const Input_para& inp,
                                         bool& conv_esolver, int iter)
 {
+    if (inp.calculation == "nscf")
+    {
+        return false;
+    }
+
     /// Return if EXX is not enabled
     if (op_exx == nullptr)
     {