diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3d384199c5..250d8477a1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -785,6 +785,7 @@ target_link_libraries(
   module_pwdft
   module_ofdft
   module_stodft
+  module_dfpt
   psi
   psi_initializer
   psi_overall_init
diff --git a/examples/18_md/02_lcao_output/INPUT b/examples/18_md/02_lcao_output/INPUT
new file mode 100644
index 0000000000..7c29d6e135
--- /dev/null
+++ b/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
diff --git a/examples/18_md/02_lcao_output/KPT b/examples/18_md/02_lcao_output/KPT
new file mode 100644
index 0000000000..c289c0158a
--- /dev/null
+++ b/examples/18_md/02_lcao_output/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/examples/18_md/02_lcao_output/STRU b/examples/18_md/02_lcao_output/STRU
new file mode 100644
index 0000000000..4dcb19b678
--- /dev/null
+++ b/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
diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt
index 64b761e00c..3c1741ca2e 100644
--- a/source/source_cell/CMakeLists.txt
+++ b/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)
diff --git a/source/source_cell/qlist.cpp b/source/source_cell/qlist.cpp
new file mode 100644
index 0000000000..762f52bdad
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/source/source_cell/qlist.h b/source/source_cell/qlist.h
new file mode 100644
index 0000000000..6c48594c8c
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/source/source_esolver/CMakeLists.txt b/source/source_esolver/CMakeLists.txt
index 4a83a9e9a1..f0f91cc395 100644
--- a/source/source_esolver/CMakeLists.txt
+++ b/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
diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp
new file mode 100644
index 0000000000..9855fcac61
--- /dev/null
+++ b/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
diff --git a/source/source_esolver/esolver_dfpt_pw.h b/source/source_esolver/esolver_dfpt_pw.h
new file mode 100644
index 0000000000..59f88661c3
--- /dev/null
+++ b/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
diff --git a/source/source_pw/CMakeLists.txt b/source/source_pw/CMakeLists.txt
index d4dbd91af2..ee58a96433 100644
--- a/source/source_pw/CMakeLists.txt
+++ b/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)
diff --git a/source/source_pw/module_dfpt/CMakeLists.txt b/source/source_pw/module_dfpt/CMakeLists.txt
new file mode 100644
index 0000000000..a4b60dfb7a
--- /dev/null
+++ b/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
+)
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/README.md b/source/source_pw/module_dfpt/README.md
new file mode 100644
index 0000000000..c3133c1195
--- /dev/null
+++ b/source/source_pw/module_dfpt/README.md
@@ -0,0 +1,108 @@
+# DFPT-PW Module
+
+## Overview
+
+This module implements Density Functional Perturbation Theory (DFPT) for 
+plane-wave basis set in ABACUS. It allows calculation of phonon frequencies,
+dielectric tensor, Born effective charges, and related properties.
+
+**Note:** This code is currently in the design phase and has not been 
+put into production yet. It may change in the future. Please use with caution.
+
+## Directory Structure
+
+```
+module_dfpt/
+├── README.md           # This file
+├── dfpt_pw.h           # DFPT-PW main interface class
+├── dfpt_pw.cpp         # DFPT-PW implementation
+├── dfpt_pw_data.h      # PW-specific DFPT data container
+├── dfpt_pw_data.cpp
+├── dfpt_pert.h         # Perturbation construction
+├── dfpt_pert.cpp
+├── dfpt_stern.h        # Sternheimer equation solver
+├── dfpt_stern.cpp
+├── dfpt_rho.h          # First-order density handling
+├── dfpt_rho.cpp
+├── dfpt_phon.h         # Phonon/dynamical matrix
+├── dfpt_phon.cpp
+├── dfpt_q0.h           # q=0 special handling
+├── dfpt_q0.cpp
+├── dfpt_metal.h        # Metal system handling
+├── dfpt_metal.cpp
+└── CMakeLists.txt      # Build configuration
+```
+
+## Design Philosophy
+
+### 1. Separation of Concerns
+- **Data Layer**: `DFPT_PW_Data` stores all DFPT-related data
+- **Algorithm Layer**: Individual classes handle specific algorithms
+- **Interface Layer**: `DFPT_PW` provides a clean API to ESolver
+
+### 2. Encapsulation
+- All data members in `DFPT_PW_Data` are private
+- Access is through getter/setter methods
+- Pimpl idiom used to hide implementation details from ESolver
+
+### 3. Reusability
+- Uses existing ABACUS components: Psi, Charge_Mixing, Monkhorst-Pack
+- q-point management via `ModuleCell::QList`
+- Conjugate gradient via existing HSolver
+
+### 4. KISS Principle
+- Short, descriptive function and variable names
+- Minimal dependencies between components
+- Clear separation of PW-specific code
+
+## Module Dependencies
+
+```
+DFPT_PW
+    ├── DFPT_PW_Data     # Data container
+    ├── DFPT_Pert        # Perturbation construction
+    ├── DFPT_Stern       # Sternheimer solver
+    ├── DFPT_Rho         # Density handling
+    ├── DFPT_Phon        # Phonon calculation
+    ├── DFPT_Q0          # q=0 special handling
+    ├── DFPT_Metal       # Metal system handling
+    └── ModuleCell::QList # q-point management
+```
+
+## Key Features
+
+1. **Monochromatic Perturbation**: Handles q≠0 perturbations
+2. **q=0 Specialization**: Computes dielectric tensor, Born charges, LO-TO splitting
+3. **Metal System Support**: Handles smearing, Fermi level correction
+4. **Phonon Calculation**: Assembles dynamical matrix, computes frequencies
+5. **Symmetry Support**: Uses irreducible representations for efficiency
+
+## Usage
+
+```cpp
+// In ESolver
+ModuleDFPT::DFPT_PW dfpt;
+dfpt.init(ucell, psi, nelec, ecutwfc);
+dfpt.set_qmesh(4, 4, 4);
+dfpt.set_conv_thr(1e-8);
+dfpt.run();
+
+// Get results
+std::vector<double> freq = dfpt.get_phonon_freq(q_idx);
+ModuleBase::matrix eps = dfpt.get_dielectric_tensor();
+```
+
+## Development Status
+
+- **Phase**: Design phase
+- **Author**: Mohan Chen
+- **Date**: 2026-05-18
+- **Status**: Not yet production-ready
+
+## Future Work
+
+1. Implement core Sternheimer solver
+2. Add proper error handling
+3. Complete unit tests
+4. Optimize parallelization
+5. Add LCAO support (separate module)
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_metal.cpp b/source/source_pw/module_dfpt/dfpt_metal.cpp
new file mode 100644
index 0000000000..3131a5186d
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_metal.cpp
@@ -0,0 +1,62 @@
+// ============================================================
+// 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 "dfpt_metal.h"
+
+namespace ModuleDFPT {
+
+DFPT_Metal::DFPT_Metal() {}
+
+DFPT_Metal::~DFPT_Metal() {}
+
+void DFPT_Metal::init(double sigma, const std::string& smearing_type) {
+    sigma_ = sigma;
+    smearing_type_ = smearing_type;
+}
+
+void DFPT_Metal::dfdeps(const ModuleBase::matrix& eig, double efermi, 
+                         ModuleBase::matrix& dfdeps) {
+    (void)eig;
+    (void)efermi;
+    (void)dfdeps;
+}
+
+void DFPT_Metal::compute_dmu(int q_idx, const psi::Psi<std::complex<double>>& psi,
+                             const ModuleBase::matrix& wg, const ModuleBase::matrix& dfdeps,
+                             DFPT_PW_Data& data) {
+    (void)q_idx;
+    (void)psi;
+    (void)wg;
+    (void)dfdeps;
+    (void)data;
+}
+
+void DFPT_Metal::compute_drho_metal(int q_idx, const psi::Psi<std::complex<double>>& psi,
+                                    const ModuleBase::matrix& wg, const ModuleBase::matrix& dfdeps,
+                                    double dmu, DFPT_PW_Data& data) {
+    (void)q_idx;
+    (void)psi;
+    (void)wg;
+    (void)dfdeps;
+    (void)dmu;
+    (void)data;
+}
+
+double DFPT_Metal::fd_dfdeps(double e, double efermi) {
+    (void)e;
+    (void)efermi;
+    return 0.0;
+}
+
+double DFPT_Metal::gauss_dfdeps(double e, double efermi) {
+    (void)e;
+    (void)efermi;
+    return 0.0;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_metal.h b/source/source_pw/module_dfpt/dfpt_metal.h
new file mode 100644
index 0000000000..fa5469d49f
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_metal.h
@@ -0,0 +1,46 @@
+// ============================================================
+// 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 DFPT_METAL_H
+#define DFPT_METAL_H
+
+#include "dfpt_pw_data.h"
+#include "source_psi/psi.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Metal {
+public:
+    DFPT_Metal();
+    ~DFPT_Metal();
+    
+    void init(double sigma, const std::string& smearing_type);
+    
+    void dfdeps(const ModuleBase::matrix& eig, double efermi, 
+                ModuleBase::matrix& dfdeps);
+    
+    void compute_dmu(int q_idx, const psi::Psi<std::complex<double>>& psi,
+                     const ModuleBase::matrix& wg, const ModuleBase::matrix& dfdeps,
+                     DFPT_PW_Data& data);
+    
+    void compute_drho_metal(int q_idx, const psi::Psi<std::complex<double>>& psi,
+                            const ModuleBase::matrix& wg, const ModuleBase::matrix& dfdeps,
+                            double dmu, DFPT_PW_Data& data);
+
+private:
+    double sigma_ = 0.0;
+    std::string smearing_type_ = "gaussian";
+    
+    double fd_dfdeps(double e, double efermi);
+    
+    double gauss_dfdeps(double e, double efermi);
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_METAL_H
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pert.cpp b/source/source_pw/module_dfpt/dfpt_pert.cpp
new file mode 100644
index 0000000000..4c086cbf5c
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pert.cpp
@@ -0,0 +1,64 @@
+// ============================================================
+// 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 "dfpt_pert.h"
+
+namespace ModuleDFPT {
+
+DFPT_Pert::DFPT_Pert() {}
+
+DFPT_Pert::~DFPT_Pert() {}
+
+void DFPT_Pert::init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, 
+                     ModulePW::PW_Basis_K* pw_wfc, Structure_Factor& sf) {
+    ucell_ = &ucell;
+    pw_rho_ = pw_rho;
+    pw_wfc_ = pw_wfc;
+    sf_ = &sf;
+}
+
+void DFPT_Pert::build_dv(int q_idx, int atom_idx, int dir, DFPT_PW_Data& data) {
+    (void)q_idx;
+    (void)atom_idx;
+    (void)dir;
+    (void)data;
+}
+
+void DFPT_Pert::apply_dv(int q_idx, int k_idx, const psi::Psi<std::complex<double>>& psi, 
+                         DFPT_PW_Data& data) {
+    (void)q_idx;
+    (void)k_idx;
+    (void)psi;
+    (void)data;
+}
+
+void DFPT_Pert::build_efield(const ModuleBase::Vector3<double>& field, DFPT_PW_Data& data) {
+    (void)field;
+    (void)data;
+}
+
+void DFPT_Pert::dVloc_dtau(int atom_idx, int dir, const ModuleBase::Vector3<double>& q, 
+                           std::vector<std::complex<double>>& dv) {
+    (void)atom_idx;
+    (void)dir;
+    (void)q;
+    (void)dv;
+}
+
+void DFPT_Pert::dVnl_dtau(int atom_idx, int dir, const ModuleBase::Vector3<double>& q,
+                          const psi::Psi<std::complex<double>>& psi, int k_idx,
+                          std::vector<std::vector<std::complex<double>>>& dv_psi) {
+    (void)atom_idx;
+    (void)dir;
+    (void)q;
+    (void)psi;
+    (void)k_idx;
+    (void)dv_psi;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pert.h b/source/source_pw/module_dfpt/dfpt_pert.h
new file mode 100644
index 0000000000..e5e3661fd7
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pert.h
@@ -0,0 +1,52 @@
+// ============================================================
+// 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 DFPT_PERT_H
+#define DFPT_PERT_H
+
+#include "dfpt_pw_data.h"
+#include "source_cell/unitcell.h"
+#include "source_psi/psi.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+#include "source_pw/module_pwdft/structure_factor.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Pert {
+public:
+    DFPT_Pert();
+    ~DFPT_Pert();
+    
+    void init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, 
+              ModulePW::PW_Basis_K* pw_wfc, Structure_Factor& sf);
+    
+    void build_dv(int q_idx, int atom_idx, int dir, DFPT_PW_Data& data);
+    
+    void apply_dv(int q_idx, int k_idx, const psi::Psi<std::complex<double>>& psi, 
+                  DFPT_PW_Data& data);
+    
+    void build_efield(const ModuleBase::Vector3<double>& field, DFPT_PW_Data& data);
+
+private:
+    UnitCell* ucell_ = nullptr;
+    ModulePW::PW_Basis* pw_rho_ = nullptr;
+    ModulePW::PW_Basis_K* pw_wfc_ = nullptr;
+    Structure_Factor* sf_ = nullptr;
+    
+    void dVloc_dtau(int atom_idx, int dir, const ModuleBase::Vector3<double>& q, 
+                    std::vector<std::complex<double>>& dv);
+    
+    void dVnl_dtau(int atom_idx, int dir, const ModuleBase::Vector3<double>& q,
+                   const psi::Psi<std::complex<double>>& psi, int k_idx,
+                   std::vector<std::vector<std::complex<double>>>& dv_psi);
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_PERT_H
diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp
new file mode 100644
index 0000000000..8659237d67
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_phon.cpp
@@ -0,0 +1,71 @@
+// ============================================================
+// 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 "dfpt_phon.h"
+
+namespace ModuleDFPT {
+
+DFPT_Phon::DFPT_Phon() {}
+
+DFPT_Phon::~DFPT_Phon() {}
+
+void DFPT_Phon::init(UnitCell& ucell) {
+    ucell_ = &ucell;
+}
+
+void DFPT_Phon::assemble(int q_idx, DFPT_PW_Data& data) {
+    int nat = ucell_->nat;
+    ModuleBase::matrix dynmat(3 * nat, 3 * nat);
+    dynmat.zero_out();
+    
+    ModuleBase::Vector3<double> q = data.get_qvec(q_idx);
+    ion_ion(q, dynmat);
+    electron(q_idx, data, dynmat);
+    
+    data.set_dynmat(q_idx, dynmat);
+}
+
+void DFPT_Phon::diagonalize(int q_idx, DFPT_PW_Data& data) {
+    ModuleBase::matrix dynmat = data.get_dynmat(q_idx);
+    int nat = ucell_->nat;
+    
+    std::vector<double> freq(3 * nat, 0.0);
+    for (int i = 0; i < 3 * nat; ++i) {
+        freq[i] = static_cast<double>(i);
+    }
+    
+    data.set_phon_freq(q_idx, freq);
+}
+
+void DFPT_Phon::add_loto(DFPT_PW_Data& data) {
+    (void)data;
+}
+
+bool DFPT_Phon::check_sum_rule(int q_idx, DFPT_PW_Data& data) const {
+    (void)q_idx;
+    (void)data;
+    return true;
+}
+
+void DFPT_Phon::ion_ion(const ModuleBase::Vector3<double>& q, ModuleBase::matrix& dyn) {
+    (void)q;
+    (void)dyn;
+}
+
+void DFPT_Phon::electron(int q_idx, DFPT_PW_Data& data, ModuleBase::matrix& dyn) {
+    (void)q_idx;
+    (void)data;
+    (void)dyn;
+}
+
+void DFPT_Phon::ewald_sum(const ModuleBase::Vector3<double>& q, ModuleBase::matrix& dyn) {
+    (void)q;
+    (void)dyn;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h
new file mode 100644
index 0000000000..6e343b8cd0
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_phon.h
@@ -0,0 +1,47 @@
+// ============================================================
+// 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 DFPT_PHON_H
+#define DFPT_PHON_H
+
+#include "dfpt_pw_data.h"
+#include "source_cell/unitcell.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Phon {
+public:
+    DFPT_Phon();
+    ~DFPT_Phon();
+    
+    void init(UnitCell& ucell);
+    
+    void assemble(int q_idx, DFPT_PW_Data& data);
+    
+    void diagonalize(int q_idx, DFPT_PW_Data& data);
+    
+    void add_loto(DFPT_PW_Data& data);
+    
+    bool check_sum_rule(int q_idx, DFPT_PW_Data& data) const;
+
+private:
+    UnitCell* ucell_ = nullptr;
+    
+    double ewald_alpha_ = 0.0;
+    double ewald_rcut_ = 0.0;
+    
+    void ion_ion(const ModuleBase::Vector3<double>& q, ModuleBase::matrix& dyn);
+    
+    void electron(int q_idx, DFPT_PW_Data& data, ModuleBase::matrix& dyn);
+    
+    void ewald_sum(const ModuleBase::Vector3<double>& q, ModuleBase::matrix& dyn);
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_PHON_H
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp
new file mode 100644
index 0000000000..3937e8516d
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pw.cpp
@@ -0,0 +1,112 @@
+// ============================================================
+// 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 "dfpt_pw.h"
+#include "dfpt_pw_data.h"
+#include "dfpt_pert.h"
+#include "dfpt_stern.h"
+#include "dfpt_rho.h"
+#include "dfpt_phon.h"
+#include "dfpt_q0.h"
+#include "dfpt_metal.h"
+#include "source_cell/qlist.h"
+
+namespace ModuleDFPT {
+
+class DFPT_PW::Impl {
+public:
+    Impl() {}
+    ~Impl() {}
+    
+    DFPT_PW_Data data_;
+    DFPT_Pert pert_;
+    DFPT_Stern stern_;
+    DFPT_Rho rho_;
+    DFPT_Phon phon_;
+    DFPT_Q0 q0_;
+    DFPT_Metal metal_;
+    ModuleCell::QList qlist_;
+    
+    psi::Psi<std::complex<double>> gs_psi_;
+    UnitCell* ucell_ = nullptr;
+    double nelec_ = 0.0;
+    double ecutwfc_ = 0.0;
+    
+    int nqx_ = 1, nqy_ = 1, nqz_ = 1;
+    double conv_thr_ = 1e-8;
+    int max_iter_ = 100;
+};
+
+DFPT_PW::DFPT_PW() : pimpl_(new Impl()) {}
+
+DFPT_PW::~DFPT_PW() {
+    delete pimpl_;
+}
+
+void DFPT_PW::init(UnitCell& ucell, const psi::Psi<std::complex<double>>& psi,
+                   double nelec, double ecutwfc) {
+    pimpl_->ucell_ = &ucell;
+    pimpl_->gs_psi_ = psi;
+    pimpl_->nelec_ = nelec;
+    pimpl_->ecutwfc_ = ecutwfc;
+    
+    std::vector<int> mp_grid = {pimpl_->nqx_, pimpl_->nqy_, pimpl_->nqz_};
+    pimpl_->qlist_.generate_mesh(ucell, ucell.symm, mp_grid, true);
+    
+    int nq = pimpl_->qlist_.get_nq();
+    int nk = psi.get_nk();
+    int nbands = psi.get_nbands();
+    int npw_max = psi.get_current_ngk();
+    int nrxx = 0;
+    int nspin = 1;
+    int nat = ucell.nat;
+    
+    pimpl_->data_.init(&pimpl_->qlist_, nk, nbands, npw_max, nrxx, nspin, nat);
+}
+
+void DFPT_PW::run() {
+    int nq = pimpl_->qlist_.get_nq();
+    for (int q_idx = 0; q_idx < nq; ++q_idx) {
+        pimpl_->phon_.assemble(q_idx, pimpl_->data_);
+        pimpl_->phon_.diagonalize(q_idx, pimpl_->data_);
+    }
+}
+
+std::vector<double> DFPT_PW::get_phonon_freq(int q_idx) const {
+    return pimpl_->data_.get_phon_freq(q_idx);
+}
+
+ModuleBase::matrix DFPT_PW::get_dielectric_tensor() const {
+    return pimpl_->data_.get_dielectric();
+}
+
+ModuleBase::matrix DFPT_PW::get_born_charges(int atom_idx) const {
+    return pimpl_->data_.get_born(atom_idx);
+}
+
+void DFPT_PW::set_parameters(const std::string& param_file) {
+    (void)param_file;
+}
+
+void DFPT_PW::set_qmesh(int nqx, int nqy, int nqz) {
+    pimpl_->nqx_ = nqx;
+    pimpl_->nqy_ = nqy;
+    pimpl_->nqz_ = nqz;
+}
+
+void DFPT_PW::set_conv_thr(double thr) {
+    pimpl_->conv_thr_ = thr;
+    pimpl_->data_.set_conv_thr(thr);
+}
+
+void DFPT_PW::set_max_iter(int max_iter) {
+    pimpl_->max_iter_ = max_iter;
+    pimpl_->data_.set_max_iter(max_iter);
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h
new file mode 100644
index 0000000000..b6ddee13eb
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pw.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 DFPT_PW_H
+#define DFPT_PW_H
+
+#include <string>
+#include <vector>
+#include "source_cell/unitcell.h"
+#include "source_psi/psi.h"
+
+namespace ModuleDFPT {
+
+class DFPT_PW {
+public:
+    DFPT_PW();
+    ~DFPT_PW();
+    
+    void init(UnitCell& ucell, const psi::Psi<std::complex<double>>& psi, 
+              double nelec, double ecutwfc);
+    
+    void run();
+    
+    std::vector<double> get_phonon_freq(int q_idx) const;
+    
+    ModuleBase::matrix get_dielectric_tensor() const;
+    
+    ModuleBase::matrix get_born_charges(int atom_idx) const;
+    
+    void set_parameters(const std::string& param_file);
+    
+    void set_qmesh(int nqx, int nqy, int nqz);
+    
+    void set_conv_thr(double thr);
+    
+    void set_max_iter(int max_iter);
+    
+private:
+    class Impl;
+    Impl* pimpl_;
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_PW_H
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp
new file mode 100644
index 0000000000..b0cb27c498
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp
@@ -0,0 +1,174 @@
+// ============================================================
+// 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 "dfpt_pw_data.h"
+
+namespace ModuleDFPT {
+
+DFPT_PW_Data::DFPT_PW_Data() {}
+
+DFPT_PW_Data::~DFPT_PW_Data() {
+    clean();
+}
+
+void DFPT_PW_Data::init(ModuleCell::QList* qlist, int nk, int nbands, int npw_max, 
+                        int nrxx, int nspin, int nat) {
+    qlist_ = qlist;
+    nk_ = nk;
+    nbands_ = nbands;
+    npw_max_ = npw_max;
+    nrxx_ = nrxx;
+    nspin_ = nspin;
+    nat_ = nat;
+    
+    allocate_memory();
+    is_initialized_ = true;
+}
+
+void DFPT_PW_Data::clean() {
+    deallocate_memory();
+    is_initialized_ = false;
+}
+
+int DFPT_PW_Data::get_nq() const {
+    return qlist_->get_nq();
+}
+
+ModuleBase::Vector3<double> DFPT_PW_Data::get_qvec(int q_idx) const {
+    return qlist_->get_q(q_idx);
+}
+
+int DFPT_PW_Data::get_nirr(int q_idx) const {
+    return qlist_->get_nirr(q_idx);
+}
+
+std::vector<int> DFPT_PW_Data::get_irrep_modes(int q_idx, int irrep) const {
+    return qlist_->get_irrep_modes(q_idx, irrep);
+}
+
+void DFPT_PW_Data::set_dpsi(int q_idx, int k_idx, int band_idx, 
+                             const std::vector<std::complex<double>>& psi) {
+    (void)q_idx;
+    (void)k_idx;
+    (void)band_idx;
+    (void)psi;
+}
+
+std::vector<std::complex<double>> DFPT_PW_Data::get_dpsi(int q_idx, int k_idx, int band_idx) const {
+    (void)q_idx;
+    (void)k_idx;
+    (void)band_idx;
+    return std::vector<std::complex<double>>();
+}
+
+psi::Psi<std::complex<double>>& DFPT_PW_Data::get_dpsi_obj(int q_idx) {
+    static psi::Psi<std::complex<double>> dummy;
+    (void)q_idx;
+    return dummy;
+}
+
+void DFPT_PW_Data::set_drho_r(int q_idx, int spin, const std::vector<double>& rho) {
+    (void)q_idx;
+    (void)spin;
+    (void)rho;
+}
+
+std::vector<double> DFPT_PW_Data::get_drho_r(int q_idx, int spin) const {
+    (void)q_idx;
+    (void)spin;
+    return std::vector<double>();
+}
+
+void DFPT_PW_Data::set_drho_g(int q_idx, int spin, const std::vector<std::complex<double>>& rho) {
+    (void)q_idx;
+    (void)spin;
+    (void)rho;
+}
+
+std::vector<std::complex<double>> DFPT_PW_Data::get_drho_g(int q_idx, int spin) const {
+    (void)q_idx;
+    (void)spin;
+    return std::vector<std::complex<double>>();
+}
+
+void DFPT_PW_Data::set_dv_r(int q_idx, int spin, const std::vector<double>& v) {
+    (void)q_idx;
+    (void)spin;
+    (void)v;
+}
+
+std::vector<double> DFPT_PW_Data::get_dv_r(int q_idx, int spin) const {
+    (void)q_idx;
+    (void)spin;
+    return std::vector<double>();
+}
+
+void DFPT_PW_Data::set_dynmat(int q_idx, const ModuleBase::matrix& dm) {
+    if (q_idx >= static_cast<int>(dynmat_.size())) {
+        dynmat_.resize(q_idx + 1);
+    }
+    dynmat_[q_idx] = dm;
+}
+
+ModuleBase::matrix DFPT_PW_Data::get_dynmat(int q_idx) const {
+    if (q_idx < static_cast<int>(dynmat_.size())) {
+        return dynmat_[q_idx];
+    }
+    return ModuleBase::matrix();
+}
+
+void DFPT_PW_Data::set_phon_freq(int q_idx, const std::vector<double>& freq) {
+    if (q_idx >= static_cast<int>(phon_freq_.size())) {
+        phon_freq_.resize(q_idx + 1);
+    }
+    phon_freq_[q_idx] = freq;
+}
+
+std::vector<double> DFPT_PW_Data::get_phon_freq(int q_idx) const {
+    if (q_idx < static_cast<int>(phon_freq_.size())) {
+        return phon_freq_[q_idx];
+    }
+    return std::vector<double>();
+}
+
+void DFPT_PW_Data::set_dielectric(const ModuleBase::matrix& eps) {
+    dielectric_ = eps;
+}
+
+ModuleBase::matrix DFPT_PW_Data::get_dielectric() const {
+    return dielectric_;
+}
+
+void DFPT_PW_Data::set_born(int atom_idx, const ModuleBase::matrix& z) {
+    if (atom_idx >= static_cast<int>(born_.size())) {
+        born_.resize(atom_idx + 1);
+    }
+    born_[atom_idx] = z;
+}
+
+ModuleBase::matrix DFPT_PW_Data::get_born(int atom_idx) const {
+    if (atom_idx < static_cast<int>(born_.size())) {
+        return born_[atom_idx];
+    }
+    return ModuleBase::matrix();
+}
+
+void DFPT_PW_Data::allocate_memory() {
+    dynmat_.resize(get_nq());
+    phon_freq_.resize(get_nq());
+    born_.resize(nat_);
+}
+
+void DFPT_PW_Data::deallocate_memory() {
+    dynmat_.clear();
+    phon_freq_.clear();
+    born_.clear();
+    residuals_.clear();
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h
new file mode 100644
index 0000000000..64aaa0bd0d
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_pw_data.h
@@ -0,0 +1,123 @@
+// ============================================================
+// 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 DFPT_PW_DATA_H
+#define DFPT_PW_DATA_H
+
+#include "source_base/matrix.h"
+#include "source_base/vector3.h"
+#include "source_psi/psi.h"
+#include "source_cell/qlist.h"
+#include <vector>
+#include <complex>
+
+namespace ModuleDFPT {
+
+class DFPT_PW_Data {
+public:
+    DFPT_PW_Data();
+    ~DFPT_PW_Data();
+    
+    void init(ModuleCell::QList* qlist, int nk, int nbands, int npw_max, 
+              int nrxx, int nspin, int nat);
+    
+    void clean();
+    
+    int get_nq() const;
+    ModuleBase::Vector3<double> get_qvec(int q_idx) const;
+    int get_nirr(int q_idx) const;
+    std::vector<int> get_irrep_modes(int q_idx, int irrep) const;
+    
+    void set_dpsi(int q_idx, int k_idx, int band_idx, 
+                  const std::vector<std::complex<double>>& psi);
+    std::vector<std::complex<double>> get_dpsi(int q_idx, int k_idx, int band_idx) const;
+    psi::Psi<std::complex<double>>& get_dpsi_obj(int q_idx);
+    
+    void set_drho_r(int q_idx, int spin, const std::vector<double>& rho);
+    std::vector<double> get_drho_r(int q_idx, int spin) const;
+    void set_drho_g(int q_idx, int spin, const std::vector<std::complex<double>>& rho);
+    std::vector<std::complex<double>> get_drho_g(int q_idx, int spin) const;
+    
+    void set_dv_r(int q_idx, int spin, const std::vector<double>& v);
+    std::vector<double> get_dv_r(int q_idx, int spin) const;
+    
+    void set_dynmat(int q_idx, const ModuleBase::matrix& dm);
+    ModuleBase::matrix get_dynmat(int q_idx) const;
+    void set_phon_freq(int q_idx, const std::vector<double>& freq);
+    std::vector<double> get_phon_freq(int q_idx) const;
+    
+    void set_dielectric(const ModuleBase::matrix& eps);
+    ModuleBase::matrix get_dielectric() const;
+    void set_born(int atom_idx, const ModuleBase::matrix& z);
+    ModuleBase::matrix get_born(int atom_idx) const;
+    
+    void set_compute_q0(bool flag) { compute_q0_ = flag; }
+    bool get_compute_q0() const { return compute_q0_; }
+    void set_loto(bool flag) { loto_ = flag; }
+    bool get_loto() const { return loto_; }
+    
+    void set_is_metal(bool flag) { is_metal_ = flag; }
+    bool get_is_metal() const { return is_metal_; }
+    void set_dmu(double dmu) { dmu_ = dmu; }
+    double get_dmu() const { return dmu_; }
+    
+    void set_max_iter(int iter) { max_iter_ = iter; }
+    int get_max_iter() const { return max_iter_; }
+    void set_conv_thr(double thr) { conv_thr_ = thr; }
+    double get_conv_thr() const { return conv_thr_; }
+    void set_current_iter(int iter) { current_iter_ = iter; }
+    int get_current_iter() const { return current_iter_; }
+    void set_converged(bool flag) { converged_ = flag; }
+    bool get_converged() const { return converged_; }
+    
+    void add_residual(double r) { residuals_.push_back(r); }
+    std::vector<double> get_residuals() const { return residuals_; }
+    
+private:
+    ModuleCell::QList* qlist_ = nullptr;
+    
+    int nk_ = 0;
+    int nbands_ = 0;
+    int npw_max_ = 0;
+    int nrxx_ = 0;
+    int nspin_ = 1;
+    int nat_ = 0;
+    
+    std::vector<psi::Psi<std::complex<double>>> dpsi_;
+    
+    std::vector<std::vector<std::vector<double>>> drho_r_;
+    std::vector<std::vector<std::vector<std::complex<double>>>> drho_g_;
+    
+    std::vector<std::vector<std::vector<double>>> dv_r_;
+    
+    std::vector<ModuleBase::matrix> dynmat_;
+    std::vector<std::vector<double>> phon_freq_;
+    
+    bool compute_q0_ = false;
+    bool loto_ = false;
+    ModuleBase::matrix dielectric_;
+    std::vector<ModuleBase::matrix> born_;
+    
+    bool is_metal_ = false;
+    double dmu_ = 0.0;
+    
+    int max_iter_ = 100;
+    double conv_thr_ = 1e-8;
+    int current_iter_ = 0;
+    bool converged_ = false;
+    std::vector<double> residuals_;
+    
+    bool is_initialized_ = false;
+    
+    void allocate_memory();
+    void deallocate_memory();
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_PW_DATA_H
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_q0.cpp b/source/source_pw/module_dfpt/dfpt_q0.cpp
new file mode 100644
index 0000000000..b49f9c54b2
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_q0.cpp
@@ -0,0 +1,42 @@
+// ============================================================
+// 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 "dfpt_q0.h"
+
+namespace ModuleDFPT {
+
+DFPT_Q0::DFPT_Q0() {}
+
+DFPT_Q0::~DFPT_Q0() {}
+
+void DFPT_Q0::init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, 
+                   ModulePW::PW_Basis_K* pw_wfc) {
+    ucell_ = &ucell;
+    pw_rho_ = pw_rho;
+    pw_wfc_ = pw_wfc;
+}
+
+void DFPT_Q0::compute_eps(const psi::Psi<std::complex<double>>& psi, 
+                          const ModuleBase::matrix& wg, DFPT_PW_Data& data) {
+    (void)psi;
+    (void)wg;
+    (void)data;
+}
+
+void DFPT_Q0::compute_born(const psi::Psi<std::complex<double>>& psi, DFPT_PW_Data& data) {
+    (void)psi;
+    (void)data;
+}
+
+void DFPT_Q0::pos_matrix(const psi::Psi<std::complex<double>>& psi,
+                         std::vector<std::vector<ModuleBase::Vector3<std::complex<double>>>>& r_mat) {
+    (void)psi;
+    (void)r_mat;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_q0.h b/source/source_pw/module_dfpt/dfpt_q0.h
new file mode 100644
index 0000000000..f96d4d67af
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_q0.h
@@ -0,0 +1,44 @@
+// ============================================================
+// 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 DFPT_Q0_H
+#define DFPT_Q0_H
+
+#include "dfpt_pw_data.h"
+#include "source_cell/unitcell.h"
+#include "source_psi/psi.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Q0 {
+public:
+    DFPT_Q0();
+    ~DFPT_Q0();
+    
+    void init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, 
+              ModulePW::PW_Basis_K* pw_wfc);
+    
+    void compute_eps(const psi::Psi<std::complex<double>>& psi, 
+                     const ModuleBase::matrix& wg, DFPT_PW_Data& data);
+    
+    void compute_born(const psi::Psi<std::complex<double>>& psi, DFPT_PW_Data& data);
+    
+private:
+    UnitCell* ucell_ = nullptr;
+    ModulePW::PW_Basis* pw_rho_ = nullptr;
+    ModulePW::PW_Basis_K* pw_wfc_ = nullptr;
+    
+    void pos_matrix(const psi::Psi<std::complex<double>>& psi,
+                    std::vector<std::vector<ModuleBase::Vector3<std::complex<double>>>>& r_mat);
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_Q0_H
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp
new file mode 100644
index 0000000000..4cf239adb5
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_rho.cpp
@@ -0,0 +1,53 @@
+// ============================================================
+// 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 "dfpt_rho.h"
+
+namespace ModuleDFPT {
+
+DFPT_Rho::DFPT_Rho() {}
+
+DFPT_Rho::~DFPT_Rho() {
+    if (mixer_ != nullptr) {
+        delete mixer_;
+        mixer_ = nullptr;
+    }
+}
+
+void DFPT_Rho::init(int nspin, int nrxx, ModulePW::PW_Basis* pw_rho, 
+                    ModulePW::PW_Basis_K* pw_wfc, const std::string& mix_type, 
+                    double mix_beta) {
+    nspin_ = nspin;
+    nrxx_ = nrxx;
+    pw_rho_ = pw_rho;
+    pw_wfc_ = pw_wfc;
+    (void)mix_type;
+    (void)mix_beta;
+}
+
+void DFPT_Rho::compute_drho(const psi::Psi<std::complex<double>>& psi, 
+                            const ModuleBase::matrix& wg, int q_idx, 
+                            DFPT_PW_Data& data) {
+    (void)psi;
+    (void)wg;
+    (void)q_idx;
+    (void)data;
+}
+
+void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) {
+    (void)q_idx;
+    (void)data;
+}
+
+double DFPT_Rho::get_residual(int q_idx, DFPT_PW_Data& data) const {
+    (void)q_idx;
+    (void)data;
+    return 0.0;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h
new file mode 100644
index 0000000000..9790e58382
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_rho.h
@@ -0,0 +1,51 @@
+// ============================================================
+// 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 DFPT_RHO_H
+#define DFPT_RHO_H
+
+#include "dfpt_pw_data.h"
+#include "source_psi/psi.h"
+#include "source_estate/module_charge/charge_mixing.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_basis/module_pw/pw_basis_k.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Rho {
+public:
+    DFPT_Rho();
+    ~DFPT_Rho();
+    
+    void init(int nspin, int nrxx, ModulePW::PW_Basis* pw_rho, 
+              ModulePW::PW_Basis_K* pw_wfc, const std::string& mix_type, 
+              double mix_beta);
+    
+    void compute_drho(const psi::Psi<std::complex<double>>& psi, 
+                      const ModuleBase::matrix& wg, int q_idx, 
+                      DFPT_PW_Data& data);
+    
+    void mix_drho(int q_idx, DFPT_PW_Data& data);
+    
+    double get_residual(int q_idx, DFPT_PW_Data& data) const;
+
+private:
+    int nspin_ = 1;
+    int nrxx_ = 0;
+    ModulePW::PW_Basis* pw_rho_ = nullptr;
+    ModulePW::PW_Basis_K* pw_wfc_ = nullptr;
+    
+    Charge_Mixing* mixer_ = nullptr;
+    
+    std::vector<std::vector<std::vector<double>>> drho_in_;
+    std::vector<std::vector<std::vector<double>>> drho_out_;
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_RHO_H
diff --git a/source/source_pw/module_dfpt/dfpt_stern.cpp b/source/source_pw/module_dfpt/dfpt_stern.cpp
new file mode 100644
index 0000000000..14fb3e47f6
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_stern.cpp
@@ -0,0 +1,65 @@
+// ============================================================
+// 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 "dfpt_stern.h"
+
+namespace ModuleDFPT {
+
+DFPT_Stern::DFPT_Stern() {}
+
+DFPT_Stern::~DFPT_Stern() {}
+
+void DFPT_Stern::init(int nk, int nbands, int npw_max, const ModuleBase::matrix& eig, 
+                      const ModuleBase::matrix& wg, double alpha) {
+    nk_ = nk;
+    nbands_ = nbands;
+    npw_max_ = npw_max;
+    eig_ = eig;
+    wg_ = wg;
+    alpha_ = alpha;
+}
+
+void DFPT_Stern::solve(const psi::Psi<std::complex<double>>& psi, 
+                       const std::vector<std::complex<double>>& dv_psi,
+                       int q_idx, int k_idx, int band_idx, double omega, 
+                       DFPT_PW_Data& data) {
+    (void)psi;
+    (void)dv_psi;
+    (void)q_idx;
+    (void)k_idx;
+    (void)band_idx;
+    (void)omega;
+    (void)data;
+}
+
+void DFPT_Stern::apply_pv(const std::vector<std::complex<double>>& x, 
+                          std::vector<std::complex<double>>& px) {
+    (void)x;
+    (void)px;
+}
+
+void DFPT_Stern::apply_op(const std::vector<std::complex<double>>& x,
+                          std::vector<std::complex<double>>& y,
+                          int k_idx, int band_idx) {
+    (void)x;
+    (void)y;
+    (void)k_idx;
+    (void)band_idx;
+}
+
+void DFPT_Stern::cg_solve(const std::vector<std::complex<double>>& b,
+                          std::vector<std::complex<double>>& x,
+                          int k_idx, int band_idx, double& residual) {
+    (void)b;
+    (void)x;
+    (void)k_idx;
+    (void)band_idx;
+    (void)residual;
+}
+
+} // namespace ModuleDFPT
\ No newline at end of file
diff --git a/source/source_pw/module_dfpt/dfpt_stern.h b/source/source_pw/module_dfpt/dfpt_stern.h
new file mode 100644
index 0000000000..38305e2526
--- /dev/null
+++ b/source/source_pw/module_dfpt/dfpt_stern.h
@@ -0,0 +1,53 @@
+// ============================================================
+// 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 DFPT_STERN_H
+#define DFPT_STERN_H
+
+#include "dfpt_pw_data.h"
+#include "source_psi/psi.h"
+
+namespace ModuleDFPT {
+
+class DFPT_Stern {
+public:
+    DFPT_Stern();
+    ~DFPT_Stern();
+    
+    void init(int nk, int nbands, int npw_max, const ModuleBase::matrix& eig, 
+              const ModuleBase::matrix& wg, double alpha);
+    
+    void solve(const psi::Psi<std::complex<double>>& psi, 
+               const std::vector<std::complex<double>>& dv_psi,
+               int q_idx, int k_idx, int band_idx, double omega, 
+               DFPT_PW_Data& data);
+    
+private:
+    int nk_ = 0;
+    int nbands_ = 0;
+    int npw_max_ = 0;
+    double alpha_ = 1.0;
+    
+    ModuleBase::matrix eig_;
+    ModuleBase::matrix wg_;
+    
+    void apply_pv(const std::vector<std::complex<double>>& x, 
+                  std::vector<std::complex<double>>& px);
+    
+    void apply_op(const std::vector<std::complex<double>>& x,
+                  std::vector<std::complex<double>>& y,
+                  int k_idx, int band_idx);
+    
+    void cg_solve(const std::vector<std::complex<double>>& b,
+                  std::vector<std::complex<double>>& x,
+                  int k_idx, int band_idx, double& residual);
+};
+
+} // namespace ModuleDFPT
+
+#endif // DFPT_STERN_H
\ No newline at end of file