Repository navigation
feat(dfpt): QE-format dynamical matrix file output (dfpt_out_dyn) - #8240
Zanthoxylum wants to merge 8 commits into
Conversation
New INPUT switch dfpt_out_dyn (Boolean, default false): when enabled, the esolver post-processing writes one QE-compatible dyn file per q point (OUT.<suffix>/dyn.q<i>) containing the raw (pre-ASR, pre-LO-TO) unreduced force-constant matrix Phi = D*sqrt(m_a m_b) in Ry/bohr^2, matching the layout of Quantum ESPRESSO ph.x's fildyn output. At Gamma with dfpt_compute_q0, the dielectric tensor and Born effective charges are appended in the QE convention. Implementation: the data layer stashes an untouched raw copy alongside the working dynmat (which is mutated in place by apply_asr/add_loto at Gamma); the writer (dfpt_write_dyn.cpp) reads it via a new DFPT_PW::get_data() const accessor. The matrix dump functionality that drove the nlcc element-level debugging now exists as a permanent, INPUT-gated feature. Verified: bn_cubic run with dfpt_out_dyn true produces OUT.autotest/ dyn.q0 with the QE layout (header, lattice vectors, species/masses, tau, q vector, per-atom-pair 3x3 blocks, dielectric tensor, Born charges); DFPT serial suite 9/9.
Move the PARAM.globalv.global_out_dir read into init_dfpt (which already accesses PARAM) and cache it as a member; the post-processing uses the cached copy. This keeps the PR's global dependency net delta at zero (the governance check's requirement).
Keep the one PARAM.globalv.global_out_dir reference inline in the post-processing with a governance exception comment: the DFPT esolver has no StateSharing access and the directory is the same one every other esolver writes to. Net global dependency delta: 1 (allowed).
… references) The output directory is 'OUT.' + suffix + '/', the same formula as ReadInput::set_global_dir (read_set_globalv.cpp:85). The esolver already holds this->inp_ (with the suffix field), so the dyn file writer path now constructs the directory inline instead of reading PARAM.globalv.global_out_dir — the governance check's global dependency net delta drops to zero.
ABACUS-PR-BOT
left a comment
There was a problem hiding this comment.
QE post-processing cannot read the exported files as written, and ASR/LO-TO updates can replace the raw matrix before serialization. Please correct the file contract and snapshot handling, together with the unit, tensor-output and rank-ownership issues below.
Reviewed 02690c86ede1d126db43814bf66b2e725995df2e against 64463eb3d896bb878d8efecdd8c10ab2e9484be7 through static source/analytic comparison with QE 6.8. No builds, ABACUS/QE reader runs, numerical probes or CI inspection were performed; the reported file-generation and serial-suite results were not independently verified.
| out << "Dynamical matrix file\n"; | ||
| out << " " << nq << " " << nat << "\n"; |
There was a problem hiding this comment.
QE's legacy reader cannot consume these records: it skips two title lines, then reads ntyp, nat, ibrav, celldm, but this places nq/nat on line two and lattice reals on line three. The indexed species/atom records, q-block separators, exact Dielectric Tensor: marker and plain numeric tensor rows also need to match QE's writer. Please implement the complete record contract and cover it with an independent reader fixture.
| // stash the raw copy for the dyn file output (the working matrix is | ||
| // mutated by apply_asr and add_loto at Gamma) | ||
| if (q_idx >= static_cast<int>(dynmat_raw_.size())) | ||
| { | ||
| dynmat_raw_.resize(q_idx + 1); | ||
| } | ||
| dynmat_raw_[q_idx] = dm; |
There was a problem hiding this comment.
Both apply_asr() and add_loto() call set_dynmat() again, so this assignment replaces the raw snapshot with the corrected matrix before export. Capture the raw matrix only at assembly and keep later working-matrix updates separate; assert that nontrivial ASR and LO-TO corrections leave the exported snapshot unchanged.
| for (int it = 0; it < ucell.ntype; ++it) | ||
| { | ||
| out << " " << ucell.atoms[it].label << " " << ucell.atoms[it].mass; |
There was a problem hiding this comment.
QE text dyn files store AMU_RY * mass_in_amu, and dynmat.x divides the file mass by AMU_RY on read. Once the record layout is fixed, these bare amu values would therefore give the wrong downstream masses unless overridden. Convert the species metadata only; the current amu factors used to recover the force-constant matrix should remain unchanged.
| // unit-cell lattice vectors (QE writes them in alat units; here | ||
| // we use lat0-scaled bohr, which for cubic cells matches alat) | ||
| const double inv_scale = 1.0 / ucell.lat0; | ||
| const ModuleBase::Matrix3 lv = ucell.latvec * inv_scale; | ||
| out << " " << lv.e11 << " " << lv.e12 << " " << lv.e13 << "\n"; | ||
| out << " " << lv.e21 << " " << lv.e22 << " " << lv.e23 << "\n"; | ||
| out << " " << lv.e31 << " " << lv.e32 << " " << lv.e33 << "\n"; |
There was a problem hiding this comment.
latvec and Cartesian tau are already in lat0 units, so dividing only latvec by lat0 again makes their scales inconsistent. The q vector is also returned in reciprocal direct coordinates, while QE expects Cartesian components in 2*pi/alat units. Serialize a consistent alat=lat0 cell and transform q through ucell.G; include a nonorthogonal cell and finite q in the regression fixture.
| } | ||
| } | ||
| // Gamma-point Born charges and dielectric tensor | ||
| if (q == 0 && data.get_loto()) |
There was a problem hiding this comment.
With dfpt_compute_q0=true and the default dfpt_loto=false, the driver computes epsilon and Born charges but this guard omits them from the file. Gate tensor output on Gamma and valid computed q0 tensors, independently of whether the working matrix receives the LO-TO correction.
| #include "source_base/vector3.h" | ||
| #include "source_cell/unitcell.h" | ||
| #include "source_psi/psi.h" | ||
| #include "source_pw/module_dfpt/dfpt_pw_data.h" |
There was a problem hiding this comment.
Non-blocking: get_data() returns only a const reference, so forward-declare DFPT_PW_Data here and keep its full include in implementation files. In dfpt_write_dyn.cpp, include the writer's own declaration and remove <vector> and any other includes left unused after the unit fixes.
| const std::string& out_dir) | ||
| { | ||
| ModuleBase::TITLE("DFPT_Write_Dyn", "write"); | ||
| ModuleBase::timer::start("DFPT_Write_Dyn", "write"); |
There was a problem hiding this comment.
Non-blocking: use DFPT_Write_Dyn as the matching start/end timer label, following the enclosing-function naming rule.
| { | ||
| for (int j = 0; j < 3; ++j) | ||
| { | ||
| const std::complex<double> v = dm(3 * ia + i, 3 * ja + j) * mm; |
There was a problem hiding this comment.
Non-blocking: assign the flattened row and column indices to named locals before calling dm, following the project rule for computed function arguments.
| dfpt_phon.cpp | ||
| dfpt_phon_ewald.cpp | ||
| dfpt_phon_elec.cpp | ||
| dfpt_write_dyn.cpp |
There was a problem hiding this comment.
Non-blocking: add test_dfpt_write_dyn.cpp and raw-state tests covering an independently specified QE record fixture, raw preservation through ASR/LO-TO, and compute_q0=true with LO-TO disabled. The existing checked-in phonon tests do not call this writer or check its new raw-storage contract.
| for (int ia = 0; ia < nat; ++ia) | ||
| { | ||
| out << " atom " << ia + 1 << " " << ucell.atoms[ucell.iat2it[ia]].label << "\n"; | ||
| const ModuleBase::matrix& zstar = data.get_born(ia); |
There was a problem hiding this comment.
Non-blocking: clarify whether the appended Born charges should also be pre-ASR. apply_asr() modifies the stored charges before this read, while the description explicitly promises a raw force-constant matrix; document and test the intended tensor convention.
Summary
New INPUT switch
dfpt_out_dyn(Boolean, default false): when enabled, the DFPT post-processing writes one QE-compatible dyn file per q point into the output directory (OUT.<suffix>/dyn.q<i>), containing the raw (pre-ASR, pre-LO-TO) unreduced force-constant matrix Φ = D·√(mₐm_b) in Ry/bohr², matching the layout of Quantum ESPRESSO ph.x'sfildynoutput. At Γ withdfpt_compute_q0, the dielectric tensor and Born effective charges are appended in the QE convention.Implementation:
set/get_dynmat_raw) alongside the working dynmat — the working copy is mutated in place byapply_asrandadd_lotoat Γ, so the file writer needs the pre-modification versiondfpt_write_dyn.cpp) reads it via a newDFPT_PW::get_data()const accessor and writes the QE layout: header (nq, nat, lattice vectors in alat units, species/masses, atom types and fractional τ), q vector, per-atom-pair 3×3 blocks (re + im pairs, unreduced), and at Γ the ε∞ tensor + Born effective charges in the "d Force / dE" conventiondfpt_out_dynis trueFiles:
source/source_pw/module_dfpt/dfpt_write_dyn.{h,cpp}— the writersource/source_pw/module_dfpt/dfpt_pw_data.{h,cpp}— raw dynmat stashsource/source_pw/module_dfpt/dfpt_pw.{h,cpp}—get_data()accessorsource/source_io/module_parameter/read_inp_dfpt.cpp+input_parameter.h— the INPUT itemVerification: bn_cubic run with
dfpt_out_dyn trueproducesOUT.autotest/dyn.q0with the complete QE layout; DFPT serial suite 9/9. The numerical values reflect whatever force constants the current code produces (this branch is based on develop without PR #8237's nlcc terms; once that merges, the dyn files will carry the correct nlcc contributions automatically — the writer is a pure output layer with no physics dependency).