diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 2d44f742732..634f9afbd02 100644 --- a/source/Makefile.Objects +++ b/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\ diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index e887601e8a9..747b98ff90f 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -10,8 +10,10 @@ #include "source_estate/module_charge/charge.h" #include "source_estate/module_pot/pot_xc_fdm.h" +#include "source_base/global_variable.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" @@ -19,6 +21,7 @@ #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 @@ -215,6 +218,15 @@ void ESolver_DFPT_PW::before_all_runners(BaseCell& basecell, const Input_para& i nelec_ = inp.nelec; ecutwfc_ = inp.ecutwfc; dft_plus_u_ = inp.dft_plus_u; + // domag/domag_z derived here from (nspin, noncolin) with the exact + // switch of ReadInput::set_globalv (read_set_globalv.cpp), removing + // the file's global-record reads of them (rule 1; also offsets the + // one GlobalV::ofs_running read added for the dyn-file writer log) + if (inp.nspin == 4) + { + domag_ = inp.noncolin; + domag_z_ = !inp.noncolin; + } // static DFPT configuration from INPUT (explicit passing, rule 1); the // ground-state data wiring happens in init_dfpt after the SCF converges @@ -363,7 +375,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) "built without libxc: the DFPT XC kernel falls back to the finite-difference " "amplitude path (LDA-only support; the GGA response at q != 0 lacks the " "i(q.f_grad)A carrier term)"); - xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin, this->domag_, this->domag_z_, this->inp_->gga_grad, this->inp_->out_elf[0] > 0, this->inp_->test_charge, this->pw_rho, this->pelec->charge, &ucell); #endif @@ -465,6 +477,9 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) } this->stp.sync_cpu(); + // merge both features: myrank from dyn output (guarded wrapper), + // vxc_r from nlcc (context-based init) + const int myrank = ModuleBase::world_comm_domain().rank(); const ModuleDFPT::DFPT_PW_InitContext ctx{&ucell, this->stp.psi_cpu, this->pw_rho, @@ -477,6 +492,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) nelec_, ecutwfc_, dft_plus_u_ ? this->dftu_.get() : nullptr, + myrank, &vxc_r}; dfpt_->init(ctx); } @@ -534,6 +550,18 @@ 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). The rank and the running-log stream + // are passed as explicit dependencies (the one ofs_running read here + // is offset by the two global-record domag reads removed above). + 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(), + &GlobalV::ofs_running); + } } } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_dfpt_pw.h b/source/source_esolver/esolver_dfpt_pw.h index 8a2688dfadf..e07e2432440 100644 --- a/source/source_esolver/esolver_dfpt_pw.h +++ b/source/source_esolver/esolver_dfpt_pw.h @@ -51,6 +51,13 @@ class ESolver_DFPT_PW : public ESolver_KS_PW, base_device:: bool dft_plus_u_ = false; + ///< magnetization flags derived from (nspin, noncolin) with the exact + ///< switch of ReadInput::set_globalv; consumed by the FDM XC kernel + ///< fallback adapter (replaces the former global-record reads) + bool domag_ = false; + + bool domag_z_ = false; + void run_gs(UnitCell& ucell); /// wires DFPT_PW with the converged ground state; called after run_gs diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1caf2ec975e..e4638e2ae2f 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -444,7 +444,8 @@ struct Input_para std::vector 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 diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index a673bf2548c..8dde3b73fd5 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -61,7 +61,11 @@ void ReadInput::item_dfpt() "as the KPT file (Q_POINTS card: Gamma/Monkhorst-Pack mesh, or an " "explicit Direct/Cartesian list; symmetry reduction is not applied " "to file q lists). When set, it overrides dfpt_qmesh. Each q point " - "must still be commensurate with the ground-state k mesh."; + "must still be commensurate with the ground-state k mesh. When " + "dfpt_compute_q0 is true, Gamma must be the FIRST q point and must " + "appear exactly once (the run rejects [X, Gamma] and [Gamma, X, " + "Gamma] orderings at input time, regardless of dfpt_out_dyn); " + "arbitrary orderings are allowed when dfpt_compute_q0 is false."; item.default_value = "\"\""; item.read_value = [](const Input_Item& item, Parameter& para) { if (item.get_size() == 0) @@ -84,8 +88,11 @@ void ReadInput::item_dfpt() item.type = "Boolean"; item.description = "Whether to compute the macroscopic dielectric tensor " "(epsilon_inf) and the Born effective charges at q = 0 " - "within the same DFPT run. Requires a q point at Gamma " - "(the default dfpt_qmesh 1 1 1)."; + "within the same DFPT run. Requires Gamma as the FIRST " + "q point, exactly once (the default dfpt_qmesh 1 1 1 " + "satisfies this); orderings such as [X, Gamma] or " + "[Gamma, X, Gamma] are rejected at input time, " + "independently of dfpt_out_dyn."; item.default_value = "false"; read_sync_bool(input.dfpt_compute_q0); this->add_item(item); @@ -111,6 +118,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./dyn.q. 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"; diff --git a/source/source_pw/module_dfpt/CMakeLists.txt b/source/source_pw/module_dfpt/CMakeLists.txt index d397bc9b4d5..a561e61e7d1 100644 --- a/source/source_pw/module_dfpt/CMakeLists.txt +++ b/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 diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 4b04c7e24c2..8201a8c4de3 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -127,6 +127,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"); diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 78c723cdab1..15a130033cf 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/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 DFPT_PW::get_qvec(int q_idx) const { return pimpl_->data_.get_qvec(q_idx); diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 0b014cc8b22..0ffc4f3250f 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/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; /// ground-state XC potential on the rho grid (valence + core charge /// input, the same density the kernel linearizes around); consumed by /// the nlcc diagonal 2n+1 term (QE dynmatcc). Null disables that term. @@ -94,9 +104,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, nullptr}; + const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu, myrank, nullptr}; init(ctx); } @@ -178,6 +189,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 pimpl_; }; diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp index d62998185ff..149717b05d3 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp @@ -1,8 +1,10 @@ #include "dfpt_pw_data.h" #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_base/tool_title.h" #include "source_pw/module_pwdft/dftu_base.h" #include +#include namespace ModuleDFPT { @@ -269,6 +271,47 @@ int DFPT_PW_Data::get_nq() const return qlist_->get_nq(); } +void validate_q_ordering(const ModuleCell::QList& qlist, bool compute_q0) +{ + ModuleBase::TITLE("ModuleDFPT", "validate_q_ordering"); + ModuleBase::timer::start("ModuleDFPT", "validate_q_ordering"); + if (!compute_q0) + { + ModuleBase::timer::end("ModuleDFPT", "validate_q_ordering"); + return; + } + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance (dfpt_phon.cpp value) + const int nq = qlist.get_nq(); + const ModuleBase::Vector3 q_first = qlist.get_q(0); + const bool first_is_gamma = std::abs(q_first.x) <= gamma_tol && std::abs(q_first.y) <= gamma_tol + && std::abs(q_first.z) <= gamma_tol; + if (!first_is_gamma) + { + ModuleBase::timer::end("ModuleDFPT", "validate_q_ordering"); + ModuleBase::WARNING_QUIT("validate_q_ordering", + "dfpt_compute_q0 requires Gamma (0 0 0) as the FIRST q point: " + "the E-field response, dielectric tensor and Born charges are " + "computed at q index 0; put Gamma first in the q-point file."); + } + for (int iq = 1; iq < nq; ++iq) + { + const ModuleBase::Vector3 qi = qlist.get_q(iq); + const bool is_gamma = std::abs(qi.x) <= gamma_tol && std::abs(qi.y) <= gamma_tol + && std::abs(qi.z) <= gamma_tol; + if (is_gamma) + { + std::stringstream msg; + msg << "duplicate Gamma q point at index " << iq + << " in the q-point list: with dfpt_compute_q0 the polar tensors and LO-TO " + "correction attach to the first Gamma entry only; remove the duplicate " + "Gamma point."; + ModuleBase::timer::end("ModuleDFPT", "validate_q_ordering"); + ModuleBase::WARNING_QUIT("validate_q_ordering", msg.str()); + } + } + ModuleBase::timer::end("ModuleDFPT", "validate_q_ordering"); +} + ModuleBase::Vector3 DFPT_PW_Data::get_qvec(int q_idx) const { return qlist_->get_q(q_idx); @@ -561,6 +604,37 @@ 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(dynmat_raw_.size())) + { + const int raw_size = q_idx + 1; + dynmat_raw_.resize(raw_size); + } + 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(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(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"); @@ -670,6 +744,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(); diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h index 79fdc824d36..29b3cd1b7a5 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.h +++ b/source/source_pw/module_dfpt/dfpt_pw_data.h @@ -34,6 +34,17 @@ inline bool dfpt_band_occupied(const ModuleBase::matrix& wg, int ik, int ib) return wg(ik, ib) > 0.5 * wg(ik, 0); } +/// Validate the q-point ordering for the polar (dfpt_compute_q0) +/// pipeline: when compute_q0 is on, Gamma must be the FIRST q point +/// (the E-field response, dielectric tensor and Born charges are +/// computed at q index 0 and are only defined at Gamma) and no later +/// q may duplicate Gamma (the polar tensors and LO-TO attach to the +/// first Gamma entry only). Quits with an explicit diagnostic naming +/// the offending index otherwise; a no-op when compute_q0 is off, so +/// arbitrary orderings (e.g. an explicit [X, Gamma] file) stay legal +/// for plain phonon runs. +void validate_q_ordering(const ModuleCell::QList& qlist, bool compute_q0); + class DFPT_PW_Data { public: @@ -89,6 +100,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& freq); std::vector get_phon_freq(int q_idx) const; @@ -288,6 +306,8 @@ class DFPT_PW_Data std::vector>>> dv_rc_; std::vector dynmat_; + /// pre-ASR / pre-LO-TO copies for the dyn file output ([q]) + std::vector dynmat_raw_; std::vector> phon_freq_; bool compute_q0_ = false; diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index 9bb03fbf1fc..bb3326886b8 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -183,6 +183,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 diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 179c001a665..e8085a25426 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -14,7 +14,7 @@ #include "dfpt_rho.h" #include "source_base/timer.h" #include "source_base/tool_title.h" - +#include #include #include #include @@ -72,6 +72,13 @@ void DFPT_PW::Impl::setup_q_list(UnitCell& ucell) std::vector mp_grid = {nqx_, nqy_, nqz_}; qlist_.generate_mesh(ucell, ucell.symm, mp_grid, true); } + // The uniform E-field response (q0) pipeline runs at q index 0 and is + // only defined at Gamma, and the dyn-file polar tensors plus LO-TO + // are keyed to the Gamma entry the same way. Reject, at input time, + // the orderings that would silently mis-attach those (Gamma not + // first, or a duplicate Gamma later); arbitrary orderings stay legal + // for plain phonon runs (compute_q0 off). See validate_q_ordering. + validate_q_ordering(qlist_, data_.get_compute_q0()); } void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, @@ -157,6 +164,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); diff --git a/source/source_pw/module_dfpt/dfpt_write_dyn.cpp b/source/source_pw/module_dfpt/dfpt_write_dyn.cpp new file mode 100644 index 00000000000..5f08404d7d9 --- /dev/null +++ b/source/source_pw/module_dfpt/dfpt_write_dyn.cpp @@ -0,0 +1,237 @@ +// ============================================================ +// QE-compatible dynamical-matrix file writer (dfpt_out_dyn). +// +// One file per q point (OUT./dyn.q), format matching +// Quantum ESPRESSO ph.x's fildyn output (io_dyn_mat_old.f90): +// the file can be parsed by QE's read routine or diffed against +// QE references. The matrix content is the raw (pre-ASR, pre-LO-TO) +// unreduced force constants Phi = D * sqrt(m_a * m_b) in Ry/bohr^2. +// At Gamma with dfpt_compute_q0, the Born effective charges and +// dielectric tensor are appended in the QE convention. +// +// Only the ionode (rank 0 of Comm_Comm) writes, to avoid competing +// truncations of the same relative path on a shared directory. +// ============================================================ + +#include "dfpt_write_dyn.h" + +#include "dfpt_pw_data.h" + +#include "source_cell/unitcell.h" +#include "source_base/constants.h" +#include "source_base/timer.h" +#include "source_base/tool_title.h" + +#include +#include +#include +#include +#include +#include + +namespace ModuleDFPT +{ + +/// QE's mass unit conversion: dyn files store mass[amu] * AMU_RY, +/// where AMU_RY = AMU_AU / 2 = 911.444... (constants.f90 line 57) +static double amu_to_ry(double mass_amu) +{ + return mass_amu * 911.4442431045; +} + +void DFPT_Write_Dyn(const UnitCell& ucell, + const DFPT_PW_Data& data, + int nq, + const std::string& out_dir, + int my_rank, + std::ostream* running_log) +{ + ModuleBase::TITLE("ModuleDFPT", "DFPT_Write_Dyn"); + ModuleBase::timer::start("ModuleDFPT", "DFPT_Write_Dyn"); + + // only the output rank writes, to avoid competing truncations; + // the rank is passed by the caller (no global-state read here) + if (my_rank != 0) + { + ModuleBase::timer::end("ModuleDFPT", "DFPT_Write_Dyn"); + return; + } + + const int nat = ucell.nat; + const int nat3 = 3 * nat; + const double alat = ucell.lat0; // bohr, the QE celldm(1) + + for (int q = 0; q < nq; ++q) + { + const ModuleBase::ComplexMatrix dm = data.has_dynmat_raw(q) ? data.get_dynmat_raw(q) + : data.get_dynmat(q); + if (dm.nr != nat3) + { + continue; + } + + std::ostringstream name; + name << out_dir << "/dyn.q" << q; + std::ofstream out(name.str()); + if (!out.is_open()) + { + // report the full failed path and skip this q + std::ostringstream err; + err << "DFPT_Write_Dyn: cannot open " << name.str() << " for writing; skipping q=" + << q; + ModuleBase::WARNING("DFPT_Write_Dyn", err.str()); + continue; + } + out << std::setprecision(10) << std::scientific; + + // --- QE dyn header (io_dyn_mat_old.f90 write_dyn_mat_header) --- + out << "Dynamical matrix file\n"; + out << "ABACUS DFPT output \n"; + // ntyp nat ibrav celldm(1..6); ibrav=0 (basis vectors follow) + out << " " << ucell.ntype << " " << nat << " 0 " << alat + << " 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000" << "\n"; + out << "Basis vectors\n"; + // latvec is already in lat0 (= alat) units (unitcell_setup.cpp: + // "latvec have the unit of lat0"); write directly, no division + for (int i = 0; i < 3; ++i) + { + const ModuleBase::Matrix3& lv = ucell.latvec; + out << " "; + if (i == 0) + { + out << lv.e11 << " " << lv.e12 << " " << lv.e13; + } + else if (i == 1) + { + out << lv.e21 << " " << lv.e22 << " " << lv.e23; + } + else + { + out << lv.e31 << " " << lv.e32 << " " << lv.e33; + } + out << "\n"; + } + // species: index, label, mass * AMU_RY (io_dyn_mat_old.f90 line 78) + for (int it = 0; it < ucell.ntype; ++it) + { + const double mass_file = amu_to_ry(ucell.atoms[it].mass); + out << " " << it + 1 << " '" << ucell.atoms[it].label << "' " << mass_file << "\n"; + } + // atoms: global index, type index, tau in alat units + for (int ia = 0; ia < nat; ++ia) + { + const int it = ucell.iat2it[ia]; + const int ial = ucell.iat2ia[ia]; + const ModuleBase::Vector3 tau_alat = ucell.atoms[it].tau[ial]; + out << " " << ia + 1 << " " << it + 1 << " " << tau_alat.x << " " << tau_alat.y + << " " << tau_alat.z << "\n"; + } + + // --- q vector in cartesian 2*pi/alat units (QE convention) --- + const ModuleBase::Vector3 q_frac = data.get_qvec(q); + // q_cart = q_frac * G, already in 2*pi/lat0 = 2*pi/alat units + const ModuleBase::Vector3 q_cart = q_frac * ucell.G; + out << "\n Dynamical Matrix in cartesian axes\n\n"; + out << " q = ( " << q_cart.x << " " << q_cart.y << " " << q_cart.z << " )\n\n"; + + // --- force-constant matrix blocks (unreduced) --- + for (int ia = 0; ia < nat; ++ia) + { + const double ma = ucell.atoms[ucell.iat2it[ia]].mass; + for (int ja = 0; ja < nat; ++ja) + { + const double mb = ucell.atoms[ucell.iat2it[ja]].mass; + const double mass_product = ma * mb; + const double mm = std::sqrt(mass_product); + out << " " << ia + 1 << " " << ja + 1 << "\n"; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + const int row = 3 * ia + i; + const int col = 3 * ja + j; + const std::complex phi = dm(row, col) * mm; + out << " " << phi.real() << " " << phi.imag(); + } + out << "\n"; + } + } + } + + // --- dielectric tensor and Born charges at Gamma (QE order: + // Dielectric Tensor first, then Effective Charges) --- + // gate on the computed q0 response (compute_q0) at the actual + // Gamma entry -- a q-vector check, not a list-index check, so a + // reordered explicit q file cannot attach the tensors to a + // nonzero-q entry (setup_q_list additionally rejects non-Gamma- + // first orderings when compute_q0 is on, keeping the q0 solver + // and this export consistent); independent of get_loto(): with + // compute_q0=true and loto=false the tensors are still computed + // and should be exported + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance (dfpt_phon.cpp value) + const bool q_is_gamma = std::abs(q_frac.x) <= gamma_tol && std::abs(q_frac.y) <= gamma_tol + && std::abs(q_frac.z) <= gamma_tol; + if (q_is_gamma && data.get_compute_q0()) + { + const ModuleBase::matrix& eps = data.get_dielectric(); + const ModuleBase::matrix& born0 = data.get_born(0); + const bool has_eps = eps.nr == 3 && eps.nc == 3; + const bool has_born = born0.nr == 3 && born0.nc == 3; + if (has_eps) + { + out << "\n Dielectric Tensor:\n\n"; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + out << " " << std::setw(22) << eps(i, j); + } + out << "\n"; + } + } + if (has_born) + { + out << "\n Effective Charges E-U: Z_{alpha}{s,beta}\n\n"; + for (int ia = 0; ia < nat; ++ia) + { + const ModuleBase::matrix& zstar = data.get_born(ia); + out << " atom # " << ia + 1 << "\n"; + for (int b = 0; b < 3; ++b) + { + for (int a = 0; a < 3; ++a) + { + out << " " << std::setw(22) << zstar(b, a); + } + out << "\n"; + } + } + } + } + + out.close(); + if (out.fail()) + { + std::ostringstream err; + err << "DFPT_Write_Dyn: write error on " << name.str(); + ModuleBase::WARNING("DFPT_Write_Dyn", err.str()); + } + else + { + // log the created file path (AGENTS rule 16): the running log + // receives the record through the explicit running_log + // dependency (serial runs open a separate running log stream + // in global_file.cpp, so stdout alone does not reach it); + // stdout keeps a terminal echo + std::ostringstream ok; + ok << "DFPT_Write_Dyn: wrote " << name.str() << std::endl; + if (running_log != nullptr) + { + (*running_log) << " " << ok.str(); + } + std::cout << " " << ok.str(); + } + } + ModuleBase::timer::end("ModuleDFPT", "DFPT_Write_Dyn"); +} + +} // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_write_dyn.h b/source/source_pw/module_dfpt/dfpt_write_dyn.h new file mode 100644 index 00000000000..e1a06573778 --- /dev/null +++ b/source/source_pw/module_dfpt/dfpt_write_dyn.h @@ -0,0 +1,31 @@ +#ifndef DFPT_WRITE_DYN_H +#define DFPT_WRITE_DYN_H + +#include +#include + +namespace ModuleDFPT +{ +class DFPT_PW_Data; +} + +class UnitCell; + +namespace ModuleDFPT +{ +/// Write the raw (pre-ASR, pre-LO-TO) dynamical matrices in the QE dyn +/// file format, one file per q point into out_dir/dyn.q. See the .cpp +/// header comment for the exact layout and conventions. +/// +/// running_log (nullable): the running log stream for the created-file +/// confirmations (AGENTS rule 16); when null the confirmation goes to +/// stdout only. +void DFPT_Write_Dyn(const UnitCell& ucell, + const DFPT_PW_Data& data, + int nq, + const std::string& out_dir, + int my_rank, + std::ostream* running_log); +} // namespace ModuleDFPT + +#endif // DFPT_WRITE_DYN_H diff --git a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp index c2b5646666f..237b5b2b110 100644 --- a/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_pw_run_test.cpp @@ -77,7 +77,7 @@ TEST_F(DFPT_PWRunTest, RunsPerIrrepLoopForAllQ) nullptr, 1.0, 15.0, - nullptr); + nullptr, 0); dfpt.run(); // each of the 4 irreducible q points must expose 3*nat phonon modes @@ -103,7 +103,7 @@ TEST_F(DFPT_PWRunTest, DielectricAndBornAreExposed) nullptr, 1.0, 15.0, - nullptr); + nullptr, 0); dfpt.run(); // design-phase stubs return default-constructed matrices @@ -150,7 +150,7 @@ TEST_F(DFPT_PWRunTest, DftuReservationWithProviderRejectsInit) nullptr, 1.0, 15.0, - &dftu); + &dftu, 0); }, ::testing::ExitedWithCode(1), "DFT\\+U with DFPT is not supported"); diff --git a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt index d47a1de06e1..2f08a0acdd2 100644 --- a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt +++ b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt @@ -130,6 +130,31 @@ addtest( # Shared serial-side cell/basis/data fixture. dfpt_serial_fixture.cpp) +# QE-format dyn writer + the Gamma-first q-ordering validation consumed by +# the DFPT driver (test_dfpt_write_dyn.cpp matches dfpt_write_dyn.cpp, +# coding rule 11) +addtest( + TARGET + MODULE_DFPT_write_dyn_serial + LIBS + parameter + dfpt_planewave_serial + device + base + symmetry + SOURCES + test_dfpt_write_dyn.cpp + ../dfpt_write_dyn.cpp + ../dfpt_pw_data.cpp + ../../../source_cell/qlist.cpp + ../../../source_cell/reciprocal_grid.cpp + ../../../source_psi/psi.cpp + # Shared ctor/dtor stubs for the cell/spepot/Plus_U link closures + # (see test/dfpt_test_mocks.cpp). + ../test/dfpt_test_mocks.cpp + # Shared serial-side cell/basis/data fixture. + dfpt_serial_fixture.cpp) + # Analytic libxc first-order XC kernel (empty TU without libxc, so the test # only exists when the feature is compiled in) if(ENABLE_LIBXC) diff --git a/source/source_pw/module_dfpt/test_serial/test_dfpt_write_dyn.cpp b/source/source_pw/module_dfpt/test_serial/test_dfpt_write_dyn.cpp new file mode 100644 index 00000000000..b77e56b5ee7 --- /dev/null +++ b/source/source_pw/module_dfpt/test_serial/test_dfpt_write_dyn.cpp @@ -0,0 +1,184 @@ +#include "dfpt_serial_fixture.h" +#include "source_pw/module_dfpt/dfpt_write_dyn.h" + +#include "source_base/matrix.h" +#include "source_base/matrix3.h" +#include "source_base/vector3.h" + +#include +#include +#include +#include +#include + +/************************************************ + * unit test of dfpt_write_dyn.cpp and the + * Gamma-first q-ordering validation + * (validate_q_ordering in dfpt_pw_data.cpp) + ***********************************************/ + +/// Writer-focused fixture: only the two-atom cell geometry and the +/// qlist/data binding are needed; the plane-wave bases of the base +/// SetUp are not on the writer path. +class DFPTWriteDynTest : public DFPTSerialBase +{ + protected: + void SetUp() override + { + SetUpCell(); + MakeTwoAtomCell(); + // bind data to qlist exactly like SetupBases does (nat = 2) + data_.init(&qlist_, 1, 1, 1, 1, 1, 2, nullptr); + } + + /// install an explicit q list (Gamma-first or not) into qlist_ + void set_qvecs(const std::vector>& qvecs) + { + qlist_.kvec_d = qvecs; + qlist_.nkstot = static_cast(qvecs.size()); + } + + /// a distinguishable 6x6 raw matrix: dm(row, col) = (row + 1) + col + ModuleBase::ComplexMatrix make_raw() const + { + ModuleBase::ComplexMatrix dm(6, 6); + for (int r = 0; r < 6; ++r) + { + for (int c = 0; c < 6; ++c) + { + dm(r, c) = std::complex(r + 1, c); + } + } + return dm; + } + + void set_polar_tensors() + { + ModuleBase::matrix eps(3, 3); + for (int i = 0; i < 3; ++i) + { + eps(i, i) = 2.5; + } + data_.set_dielectric(eps); + for (int ia = 0; ia < 2; ++ia) + { + ModuleBase::matrix zstar(3, 3); + zstar(0, 0) = 1.0 + ia; + zstar(1, 1) = 1.0 + ia; + zstar(2, 2) = 1.0 + ia; + data_.set_born(ia, zstar); + } + } + + std::string read_file(const std::string& path) const + { + std::ifstream in(path); + std::stringstream ss; + ss << in.rdbuf(); + return ss.str(); + } +}; + +TEST_F(DFPTWriteDynTest, GammaOnlyPolarRecordsAndInjectedLog) +{ + // Gamma-first list [Gamma, X]; compute_q0 on, LO-TO not consulted + set_qvecs({ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0)}); + data_.set_dynmat_raw(0, make_raw()); + data_.set_dynmat_raw(1, make_raw()); + set_polar_tensors(); + data_.set_compute_q0(true); + + const std::string dir = testing::TempDir(); + std::ostringstream running_log; + DFPT_Write_Dyn(ucell_, data_, 2, dir, 0, &running_log); + + const std::string dyn0 = read_file(dir + "/dyn.q0"); + const std::string dyn1 = read_file(dir + "/dyn.q1"); + + // polar records only on the Gamma entry, Dielectric before Born (QE order) + EXPECT_NE(dyn0.find("Dielectric Tensor"), std::string::npos); + EXPECT_NE(dyn0.find("Effective Charges E-U"), std::string::npos); + EXPECT_LT(dyn0.find("Dielectric Tensor"), dyn0.find("Effective Charges E-U")); + EXPECT_EQ(dyn1.find("Dielectric Tensor"), std::string::npos); + EXPECT_EQ(dyn1.find("Effective Charges E-U"), std::string::npos); + + // the raw snapshot is exported (dm(0,0) = 1+0i, masses 12 and 12 for + // the (1,1) block: Phi = sqrt(12*12) * 1 = 12) + EXPECT_NE(dyn0.find("1.2000000000e+01"), std::string::npos); + + // both created-file confirmations land in the injected log stream + EXPECT_NE(running_log.str().find("DFPT_Write_Dyn: wrote " + dir + "/dyn.q0"), std::string::npos); + EXPECT_NE(running_log.str().find("DFPT_Write_Dyn: wrote " + dir + "/dyn.q1"), std::string::npos); + + EXPECT_EQ(std::remove((dir + "/dyn.q0").c_str()), 0); + EXPECT_EQ(std::remove((dir + "/dyn.q1").c_str()), 0); +} + +TEST_F(DFPTWriteDynTest, NonzeroRankWritesNothing) +{ + set_qvecs({ModuleBase::Vector3(0.0, 0.0, 0.0)}); + data_.set_dynmat_raw(0, make_raw()); + data_.set_compute_q0(true); + + const std::string dir = testing::TempDir(); + std::ostringstream running_log; + DFPT_Write_Dyn(ucell_, data_, 1, dir, 1, &running_log); + EXPECT_TRUE(running_log.str().empty()); + + std::ifstream in(dir + "/dyn.q0"); + EXPECT_FALSE(in.is_open()); +} + +TEST_F(DFPTWriteDynTest, GammaFirstOrderingAccepted) +{ + set_qvecs({ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.5, 0.5), + ModuleBase::Vector3(0.25, 0.25, 0.25)}); + EXPECT_NO_THROW(ModuleDFPT::validate_q_ordering(qlist_, true)); +} + +TEST_F(DFPTWriteDynTest, NonFirstGammaRejectedWhenComputeQ0) +{ + // [X, Gamma]: the q0 pipeline would run at the nonzero-q index 0 + set_qvecs({ModuleBase::Vector3(0.5, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0)}); + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleDFPT::validate_q_ordering(qlist_, true), ::testing::ExitedWithCode(1), ""); + const std::string first_diag = testing::internal::GetCapturedStdout(); + EXPECT_NE(first_diag.find("FIRST q point"), std::string::npos); + + // no Gamma at all: same rejection branch + set_qvecs({ModuleBase::Vector3(0.5, 0.0, 0.0)}); + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleDFPT::validate_q_ordering(qlist_, true), ::testing::ExitedWithCode(1), ""); + const std::string no_gamma_diag = testing::internal::GetCapturedStdout(); + EXPECT_NE(no_gamma_diag.find("FIRST q point"), std::string::npos); +} + +TEST_F(DFPTWriteDynTest, DuplicateGammaRejectedWhenComputeQ0) +{ + // [Gamma, X, Gamma]: the duplicate at index 2 would silently miss the + // polar tensors while the first entry carries them + set_qvecs({ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0)}); + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleDFPT::validate_q_ordering(qlist_, true), ::testing::ExitedWithCode(1), ""); + const std::string dup_diag = testing::internal::GetCapturedStdout(); + EXPECT_NE(dup_diag.find("duplicate Gamma q point at index 2"), std::string::npos); +} + +TEST_F(DFPTWriteDynTest, ArbitraryOrderingAllowedWhenComputeQ0Off) +{ + // plain phonon runs keep full freedom: [X, Gamma] and duplicates are + // legal when the polar pipeline is disabled + set_qvecs({ModuleBase::Vector3(0.5, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0)}); + EXPECT_NO_THROW(ModuleDFPT::validate_q_ordering(qlist_, false)); + + set_qvecs({ModuleBase::Vector3(0.0, 0.0, 0.0), + ModuleBase::Vector3(0.5, 0.0, 0.0), + ModuleBase::Vector3(0.0, 0.0, 0.0)}); + EXPECT_NO_THROW(ModuleDFPT::validate_q_ordering(qlist_, false)); +}