Add the framework of DFPT in plane wave basis

CMakeLists.txt

@@ -785,6 +785,7 @@ target_link_libraries(
   module_pwdft
   module_ofdft
   module_stodft
+  module_dfpt
   psi
   psi_initializer
   psi_overall_init

examples/18_md/02_lcao_output/INPUT

@@ -0,0 +1,38 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix                 Si_nhc_nvt
+calculation            md
+nbands                 20
+symmetry               0
+pseudo_dir             ../../../tests/PP_ORB
+orbital_dir            ../../../tests/PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc                10
+scf_thr                1e-5
+scf_nmax               100
+
+#Parameters (3.Basis)
+basis_type             lcao
+ks_solver              genelpa
+gamma_only             1
+
+#Parameters (4.Smearing)
+smearing_method        gaussian
+smearing_sigma         0.001
+
+#Parameters (5.Mixing)
+mixing_type            broyden
+mixing_beta            0.3
+chg_extrap             second-order
+
+#Parameters (6.MD)
+md_type                nvt
+md_nstep               10
+md_dt                  1
+md_tfirst              300
+md_tfreq               1
+md_tchain              1
+
+dft_functional       hse
+restart_save           1

examples/18_md/02_lcao_output/KPT

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

examples/18_md/02_lcao_output/STRU

@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+Si  28.085  Si_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Si_gga_8au_60Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+1.8897270    # 1 Angstrom = 1.8897270 bohr
+
+LATTICE_VECTORS
+5.43090      0.00000      0.00000
+0.00000      5.43090      0.00000
+0.00000      0.00000      5.43090
+
+ATOMIC_POSITIONS
+Direct 
+
+Si
+0.0
+2
+0.000	0.000	0.000	1	1	1
+0.5 0.5 0.5 1 1 1

source/source_cell/CMakeLists.txt

@@ -30,6 +30,8 @@ add_library(
     k_vector_utils.cpp
     sep.cpp
     sep_cell.cpp
+    qlist.cpp
+    qlist.h
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/qlist.cpp

@@ -0,0 +1,56 @@
+// ============================================================
+// This code is added by Mohan Chen on 2026-05-18.
+// This code is currently in the design phase and has not been
+// put into production yet. It may change in the future.
+// Please use this code with caution. Only developers who know
+// what they are doing should use this code.
+// ============================================================
+
+#include "qlist.h"
+
+namespace ModuleCell {
+
+QList::QList() {}
+
+QList::~QList() {}
+
+void QList::generate_mesh(UnitCell& ucell, ModuleSymmetry::Symmetry& symm,
+                          const std::vector<int>& mp_grid, bool use_irreps) {
+    (void)ucell;
+    (void)symm;
+    (void)mp_grid;
+    (void)use_irreps;
+
+    nq_ = 1;
+    qvec_.resize(nq_);
+    qvec_[0] = ModuleBase::Vector3<double>(0.0, 0.0, 0.0);
+
+    nirr_.resize(nq_);
+    nirr_[0] = 1;
+
+    irrep_modes_.resize(nq_);
+    irrep_modes_[0].resize(1);
+}
+
+void QList::read_from_file(const std::string& filename, UnitCell& ucell) {
+    (void)filename;
+    (void)ucell;
+}
+
+std::vector<int> QList::get_irrep_modes(int q_idx, int irrep_idx) const {
+    (void)q_idx;
+    (void)irrep_idx;
+    return std::vector<int>();
+}
+
+void QList::reduce(UnitCell& ucell, ModuleSymmetry::Symmetry& symm) {
+    (void)ucell;
+    (void)symm;
+}
+
+void QList::get_irreps(UnitCell& ucell, ModuleSymmetry::Symmetry& symm) {
+    (void)ucell;
+    (void)symm;
+}
+
+} // namespace ModuleCell

source/source_cell/qlist.h

@@ -0,0 +1,50 @@
+// ============================================================
+// This code is added by Mohan Chen on 2026-05-18.
+// This code is currently in the design phase and has not been
+// put into production yet. It may change in the future.
+// Please use this code with caution. Only developers who know
+// what they are doing should use this code.
+// ============================================================
+
+#ifndef QLIST_H
+#define QLIST_H
+
+#include "source_base/vector3.h"
+#include "module_symmetry/symmetry.h"
+#include "unitcell.h"
+#include <vector>
+
+namespace ModuleCell {
+
+class QList {
+public:
+    QList();
+    ~QList();
+
+    void generate_mesh(UnitCell& ucell, ModuleSymmetry::Symmetry& symm,
+                       const std::vector<int>& mp_grid, bool use_irreps);
+
+    void read_from_file(const std::string& filename, UnitCell& ucell);
+
+    int get_nq() const { return nq_; }
+
+    ModuleBase::Vector3<double> get_q(int idx) const { return qvec_[idx]; }
+
+    int get_nirr(int idx) const { return nirr_[idx]; }
+
+    std::vector<int> get_irrep_modes(int q_idx, int irrep_idx) const;
+
+private:
+    int nq_ = 0;
+    std::vector<ModuleBase::Vector3<double>> qvec_;
+    std::vector<int> nirr_;
+    std::vector<std::vector<std::vector<int>>> irrep_modes_;
+
+    void reduce(UnitCell& ucell, ModuleSymmetry::Symmetry& symm);
+
+    void get_irreps(UnitCell& ucell, ModuleSymmetry::Symmetry& symm);
+};
+
+} // namespace ModuleCell
+
+#endif // QLIST_H

source/source_esolver/CMakeLists.txt

@@ -13,6 +13,7 @@ list(APPEND objects
     esolver_of_interface.cpp
     esolver_of_tool.cpp
     pw_others.cpp
+    esolver_dfpt_pw.cpp
 )
 if(ENABLE_LCAO)
   list(APPEND objects

source/source_esolver/esolver_dfpt_pw.cpp

@@ -0,0 +1,82 @@
+// ============================================================
+// This code is added by Mohan Chen on 2026-05-18.
+// This code is currently in the design phase and has not been
+// put into production yet. It may change in the future.
+// Please use this code with caution. Only developers who know
+// what they are doing should use this code.
+// ============================================================
+
+#include "esolver_dfpt_pw.h"
+#include "source_base/tool_quit.h"
+
+namespace ModuleESolver {
+
+ESolver_DFPT_PW::ESolver_DFPT_PW() {
+    this->classname = "ESolver_DFPT_PW";
+    this->basisname = "PW";
+    gs_done_ = false;
+    dfpt_ = nullptr;
+}
+
+ESolver_DFPT_PW::~ESolver_DFPT_PW() {
+    if (dfpt_ != nullptr) {
+        delete dfpt_;
+        dfpt_ = nullptr;
+    }
+}
+
+void ESolver_DFPT_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "before_all_runners");
+
+    ESolver_KS_PW<std::complex<double>, base_device::DEVICE_CPU>::before_all_runners(ucell, inp);
+
+    init_dfpt(ucell);
+}
+
+void ESolver_DFPT_PW::runner(UnitCell& ucell, const int istep) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "runner");
+
+    if (!gs_done_) {
+        run_gs(ucell);
+        gs_done_ = true;
+    }
+
+    if (dfpt_ != nullptr) {
+        dfpt_->run();
+    }
+
+    run_post_process(ucell);
+}
+
+void ESolver_DFPT_PW::after_all_runners(UnitCell& ucell) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "after_all_runners");
+
+    ESolver_KS_PW<std::complex<double>, base_device::DEVICE_CPU>::after_all_runners(ucell);
+}
+
+void ESolver_DFPT_PW::run_gs(UnitCell& ucell) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "run_gs");
+
+    ESolver_KS_PW<std::complex<double>, base_device::DEVICE_CPU>::runner(ucell, 0);
+}
+
+void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "init_dfpt");
+
+    dfpt_ = new ModuleDFPT::DFPT_PW();
+
+//    dfpt_->init(ucell, *this->stp.psi, this->pelec->nelec, PARAM.inp.ecutwfc);
+
+    dfpt_->set_parameters("dfpt.in");
+
+    dfpt_->set_qmesh(1, 1, 1);
+
+    dfpt_->set_conv_thr(1e-8);
+    dfpt_->set_max_iter(100);
+}
+
+void ESolver_DFPT_PW::run_post_process(UnitCell& ucell) {
+    ModuleBase::TITLE("ESolver_DFPT_PW", "run_post_process");
+}
+
+} // namespace ModuleESolver

source/source_esolver/esolver_dfpt_pw.h

@@ -0,0 +1,40 @@
+// ============================================================
+// This code is added by Mohan Chen on 2026-05-18.
+// This code is currently in the design phase and has not been
+// put into production yet. It may change in the future.
+// Please use this code with caution. Only developers who know
+// what they are doing should use this code.
+// ============================================================
+
+#ifndef ESOLVER_DFPT_PW_H
+#define ESOLVER_DFPT_PW_H
+
+#include "esolver_ks_pw.h"
+#include "source_pw/module_dfpt/dfpt_pw.h"
+
+namespace ModuleESolver {
+
+class ESolver_DFPT_PW : public ESolver_KS_PW<std::complex<double>, base_device::DEVICE_CPU> {
+public:
+    ESolver_DFPT_PW();
+    ~ESolver_DFPT_PW();
+
+    void before_all_runners(UnitCell& ucell, const Input_para& inp) override;
+    void runner(UnitCell& ucell, const int istep) override;
+    void after_all_runners(UnitCell& ucell) override;
+
+protected:
+    ModuleDFPT::DFPT_PW* dfpt_ = nullptr;
+
+    bool gs_done_ = false;
+
+    void run_gs(UnitCell& ucell);
+
+    void init_dfpt(UnitCell& ucell);
+
+    void run_post_process(UnitCell& ucell);
+};
+
+} // namespace ModuleESolver
+
+#endif // ESOLVER_DFPT_PW_H

source/source_pw/CMakeLists.txt

@@ -1,3 +1,4 @@
 add_subdirectory(module_pwdft)
 add_subdirectory(module_ofdft)
 add_subdirectory(module_stodft)
+add_subdirectory(module_dfpt)

source/source_pw/module_dfpt/CMakeLists.txt

@@ -0,0 +1,39 @@
+set(MODULE_NAME module_dfpt)
+
+set(SOURCES
+    dfpt_pw.cpp
+    dfpt_pw_data.cpp
+    dfpt_pert.cpp
+    dfpt_stern.cpp
+    dfpt_rho.cpp
+    dfpt_phon.cpp
+    dfpt_q0.cpp
+    dfpt_metal.cpp
+)
+
+set(HEADERS
+    dfpt_pw.h
+    dfpt_pw_data.h
+    dfpt_pert.h
+    dfpt_stern.h
+    dfpt_rho.h
+    dfpt_phon.h
+    dfpt_q0.h
+    dfpt_metal.h
+)
+
+add_library(${MODULE_NAME} ${SOURCES} ${HEADERS})
+
+target_include_directories(${MODULE_NAME}
+    PUBLIC
+        ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+target_link_libraries(${MODULE_NAME}
+    PUBLIC
+        base
+        cell
+        psi
+        elecstate
+        module_pwdft
+)