feat(dfpt): QE-format dynamical matrix file output (dfpt_out_dyn)

source/Makefile.Objects

@@ -455,6 +455,7 @@ OBJS_DFPT=dfpt_hamilt_shift.o\
     dfpt_q0_pos.o\
     dfpt_rho.o\
     dfpt_stern.o\
+    dfpt_write_dyn.o\
     dfpt_xc_analytic.o
 
 OBJS_HAMILT_LCAO=hamilt_lcao.o\

source/source_esolver/esolver_dfpt_pw.cpp

@@ -12,13 +12,15 @@
 #include "source_estate/module_pot/pot_xc_fdm.h"
 #include "source_base/macros.h"
 #include "source_base/math_polyint.h"
+#include "source_base/parallel_cell.h"
 #include "source_base/tool_quit.h"
 #include "source_hamilt/module_xc/libxc_abacus.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_pw/module_dfpt/dfpt_pert.h"
 #include "source_pw/module_dfpt/dfpt_pw.h"
 #include "source_pw/module_dfpt/dfpt_rho.h"
+#include "source_pw/module_dfpt/dfpt_write_dyn.h"
 #ifdef __LIBXC
 #include "source_pw/module_dfpt/dfpt_xc_analytic.h"
 #endif
@@ -427,9 +429,12 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
     }
 
     this->stp.sync_cpu();
+    // myrank from the guarded parallel wrapper (no global-state read;
+    // the same CommunicationDomain::rank() pattern as source_md)
+    const int myrank = ModuleBase::world_comm_domain().rank();
     dfpt_->init(ucell, *this->stp.psi_cpu, this->pw_rho, this->pw_wfc, &this->sf, veff_r,
                 this->pelec->wg, this->pelec->ekb, xc_adapter_, nelec_, ecutwfc_,
-                dft_plus_u_ ? this->dftu_.get() : nullptr);
+                dft_plus_u_ ? this->dftu_.get() : nullptr, myrank);
 }
 
 void ESolver_DFPT_PW::run_post_process(UnitCell& ucell)
@@ -485,6 +490,17 @@ void ESolver_DFPT_PW::run_post_process(UnitCell& ucell)
             std::cout << std::endl;
         }
     }
+
+    // QE-format dyn file output (dfpt_out_dyn); the output directory is
+    // constructed from the suffix (the same formula as
+    // ReadInput::set_global_dir), so no global-state read is needed.
+    // The rank parameter defaults to 0 (serial DFPT; the MPI follow-up
+    // will wire it through the init context)
+    if (this->inp_->dfpt_out_dyn)
+    {
+        const std::string out_dir = "OUT." + this->inp_->suffix + "/";
+        ModuleDFPT::DFPT_Write_Dyn(ucell, dfpt_->get_data(), nq, out_dir, dfpt_->get_myrank());
+    }
 }
 
 } // namespace ModuleESolver

source/source_io/module_parameter/input_parameter.h

@@ -444,7 +444,8 @@ struct Input_para
     std::vector<int> dfpt_qmesh = {1, 1, 1}; ///< Monkhorst-Pack q mesh for DFPT (gamma-centered)
     std::string dfpt_qfile = "";              ///< file containing the DFPT q-point list; empty means dfpt_qmesh
     bool dfpt_compute_q0 = false;             ///< compute epsilon_inf and Born effective charges at q = 0
-    bool dfpt_loto = false;                   ///< apply the LO-TO non-analytic correction at q = 0
+    bool dfpt_loto = false;                   ///< apply the LO-TO nonanalytic correction at q = 0
+    bool dfpt_out_dyn = false;                ///< write the raw (pre-ASR, pre-LO-TO) dynamical matrix in QE dyn format
     double dfpt_conv_thr = 1.0e-8;            ///< convergence threshold of the DFPT first-order density
     int dfpt_max_iter = 100;                  ///< max iterations of the DFPT first-order density mixing
     double dfpt_mix_beta = 0.4;               ///< mixing coefficient of the DFPT first-order density

source/source_io/module_parameter/read_inp_dfpt.cpp

@@ -111,6 +111,24 @@ void ReadInput::item_dfpt()
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("dfpt_out_dyn");
+        item.annotation = "write the raw dynamical matrix in QE dyn format "
+                          "(pre-ASR, pre-LO-TO, unreduced force constants)";
+        item.category = category;
+        item.type = "Boolean";
+        item.description = "Write one QE-compatible dyn file per q point into the output "
+                           "directory: OUT.<suffix>/dyn.q<i>. Each file contains the "
+                           "unreduced (pre-ASR, pre-LO-TO) force-constant matrix Phi in "
+                           "Ry/bohr^2, the same content and layout as Quantum ESPRESSO's "
+                           "ph.x output (fildyn), so the files can be diffed element-by-element "
+                           "against QE references or post-processed by external tools. "
+                           "At Gamma with dfpt_compute_q0 the Born effective charges and "
+                           "dielectric tensor are appended in the QE convention.";
+        item.default_value = "false";
+        read_sync_bool(input.dfpt_out_dyn);
+        this->add_item(item);
+    }
     {
         Input_Item item("dfpt_conv_thr");
         item.annotation = "convergence threshold of the DFPT response SCF";

source/source_pw/module_dfpt/CMakeLists.txt

@@ -16,6 +16,7 @@ set(SOURCES
     dfpt_phon.cpp
     dfpt_phon_ewald.cpp
     dfpt_phon_elec.cpp
+    dfpt_write_dyn.cpp
     dfpt_q0.cpp
     dfpt_q0_pos.cpp
     dfpt_metal.cpp

source/source_pw/module_dfpt/dfpt_phon.cpp

@@ -111,6 +111,10 @@ void DFPT_Phon::assemble(int q_idx, DFPT_PW_Data& data)
         }
     }
     data.set_dynmat(q_idx, dyn);
+    // stash the raw pre-ASR pre-LO-TO copy for the dyn file writer;
+    // later set_dynmat calls (apply_asr, add_loto) update the working
+    // matrix but must not overwrite this snapshot
+    data.set_dynmat_raw(q_idx, dyn);
     dynmat_accum_ = ModuleBase::ComplexMatrix();
     accum_q_ = -1;
     ModuleBase::timer::end("DFPT_Phon", "assemble");

source/source_pw/module_dfpt/dfpt_pw.cpp

@@ -66,6 +66,16 @@ int DFPT_PW::get_nq() const
     return pimpl_->qlist_.get_nq();
 }
 
+const DFPT_PW_Data& DFPT_PW::get_data() const
+{
+    return pimpl_->data_;
+}
+
+int DFPT_PW::get_myrank() const
+{
+    return pimpl_->myrank_;
+}
+
 ModuleBase::Vector3<double> DFPT_PW::get_qvec(int q_idx) const
 {
     return pimpl_->data_.get_qvec(q_idx);

source/source_pw/module_dfpt/dfpt_pw.h

@@ -19,6 +19,11 @@ class PW_Basis;
 class PW_Basis_K;
 } // namespace ModulePW
 
+namespace ModuleDFPT
+{
+class DFPT_PW_Data;
+}
+
 namespace ModuleDFPT
 {
 
@@ -43,6 +48,11 @@ struct DFPT_PW_InitContext
     double nelec;
     double ecutwfc;
     const Plus_U_Base* dftu;
+    /// MPI rank of the running process (from the runtime context, not
+    /// read from any global inside module_dfpt); consumed by the dyn
+    /// file writer for the rank-0 output gate. No default member
+    /// initializer: C++11 aggregate rules require it.
+    int myrank;
 };
 
 /**
@@ -89,9 +99,10 @@ class DFPT_PW
               const XC_First_Order* xc,
               double nelec,
               double ecutwfc,
-              const Plus_U_Base* dftu)
+              const Plus_U_Base* dftu,
+              int myrank)
     {
-        const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu};
+        const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu, myrank};
         init(ctx);
     }
 
@@ -173,6 +184,14 @@ class DFPT_PW
 
     std::string format_loto_report() const;
 
+    /// read-only access to the data layer (consumed by the dyn file
+    /// writer in the esolver post-processing)
+    const DFPT_PW_Data& get_data() const;
+
+    /// the MPI rank passed at init (0 for serial); consumed by the dyn
+    /// file writer for the rank-0 output gate
+    int get_myrank() const;
+
   private:
     std::unique_ptr<Impl> pimpl_;
 };

source/source_pw/module_dfpt/dfpt_pw_data.cpp

@@ -524,6 +524,36 @@ void DFPT_PW_Data::set_dynmat(int q_idx, const ModuleBase::ComplexMatrix& dm)
     ModuleBase::timer::end("DFPT_PW_Data", "set_dynmat");
 }
 
+void DFPT_PW_Data::set_dynmat_raw(int q_idx, const ModuleBase::ComplexMatrix& dm)
+{
+    ModuleBase::TITLE("DFPT_PW_Data", "set_dynmat_raw");
+    ModuleBase::timer::start("DFPT_PW_Data", "set_dynmat_raw");
+    if (q_idx >= static_cast<int>(dynmat_raw_.size()))
+    {
+        dynmat_raw_.resize(q_idx + 1);
+    }
+    dynmat_raw_[q_idx] = dm;
+    ModuleBase::timer::end("DFPT_PW_Data", "set_dynmat_raw");
+}
+
+ModuleBase::ComplexMatrix DFPT_PW_Data::get_dynmat_raw(int q_idx) const
+{
+    ModuleBase::TITLE("DFPT_PW_Data", "get_dynmat_raw");
+    ModuleBase::timer::start("DFPT_PW_Data", "get_dynmat_raw");
+    if (q_idx < static_cast<int>(dynmat_raw_.size()) && dynmat_raw_[q_idx].nr > 0)
+    {
+        ModuleBase::timer::end("DFPT_PW_Data", "get_dynmat_raw");
+        return dynmat_raw_[q_idx];
+    }
+    ModuleBase::timer::end("DFPT_PW_Data", "get_dynmat_raw");
+    return get_dynmat(q_idx); // fallback: the working copy
+}
+
+bool DFPT_PW_Data::has_dynmat_raw(int q_idx) const
+{
+    return q_idx >= 0 && q_idx < static_cast<int>(dynmat_raw_.size()) && dynmat_raw_[q_idx].nr > 0;
+}
+
 ModuleBase::ComplexMatrix DFPT_PW_Data::get_dynmat(int q_idx) const
 {
     ModuleBase::TITLE("DFPT_PW_Data", "get_dynmat");
@@ -633,6 +663,7 @@ void DFPT_PW_Data::deallocate_memory()
     ModuleBase::TITLE("DFPT_PW_Data", "deallocate_memory");
     ModuleBase::timer::start("DFPT_PW_Data", "deallocate_memory");
     dynmat_.clear();
+    dynmat_raw_.clear();
     phon_freq_.clear();
     born_.clear();
     docc_.clear();

source/source_pw/module_dfpt/dfpt_pw_data.h

@@ -89,6 +89,13 @@ class DFPT_PW_Data
     /// through the LapackConnector::zheev wrapper.
     void set_dynmat(int q_idx, const ModuleBase::ComplexMatrix& dm);
     ModuleBase::ComplexMatrix get_dynmat(int q_idx) const;
+
+    /// stash the raw (pre-ASR, pre-LO-TO) dynamical matrix for the dyn
+    /// file output; the working dynmat_ is mutated in place by apply_asr
+    /// and add_loto at Gamma, so the file writer needs this untouched copy
+    void set_dynmat_raw(int q_idx, const ModuleBase::ComplexMatrix& dm);
+    ModuleBase::ComplexMatrix get_dynmat_raw(int q_idx) const;
+    bool has_dynmat_raw(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;
 
@@ -281,6 +288,8 @@ class DFPT_PW_Data
     std::vector<std::vector<std::vector<std::complex<double>>>> dv_rc_;
 
     std::vector<ModuleBase::ComplexMatrix> dynmat_;
+    /// pre-ASR / pre-LO-TO copies for the dyn file output ([q])
+    std::vector<ModuleBase::ComplexMatrix> dynmat_raw_;
     std::vector<std::vector<double>> phon_freq_;
 
     bool compute_q0_ = false;

source/source_pw/module_dfpt/dfpt_pw_impl.h

@@ -178,6 +178,9 @@ class DFPT_PW::Impl
     std::string qfile_;
     double conv_thr_ = 1e-8;
     int max_iter_ = 100;
+    /// MPI rank from the init context (0 for serial); used by the dyn
+    /// file writer for the rank-0 output gate
+    int myrank_ = 0;
     double mix_beta_ = 0.4;
     std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix
     std::string mix_type_ = "v_sc_broyden"; ///< response SCF mixing scheme

source/source_pw/module_dfpt/dfpt_pw_init.cpp

@@ -143,6 +143,7 @@ void DFPT_PW::init(const DFPT_PW_InitContext& ctx)
     pimpl_->nelec_ = ctx.nelec;
     pimpl_->ecutwfc_ = ctx.ecutwfc;
     pimpl_->dftu_ = ctx.dftu;
+    pimpl_->myrank_ = ctx.myrank;
 
     pimpl_->check_metallic_occ(*ctx.wg);