diff --git a/source/Makefile.Objects b/source/Makefile.Objects index f393a5eb31..2de92aedfb 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -970,6 +970,7 @@ OBJS_SRCPW=h_ewald_pw.o\ dftu_base.o\ dftu_base_io.o\ dftu_base_io_read.o\ + dftu_base_io_write.o\ dftu_pw.o\ dftu_pw_tools.o\ yukawa_screening.o\ diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index bb07dc8894..2622e59307 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -678,7 +678,6 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // Read globalv values once into locals so the repeated DFT+U call sites // below do not each re-enter the global dependency surface. const std::string& global_out_dir = PARAM.globalv.global_out_dir; - const int npol = PARAM.globalv.npol; const bool gamma_only_local = PARAM.globalv.gamma_only_local; finish_dftu_lcao(conv_esolver, this->inp_->dft_plus_u, @@ -690,7 +689,6 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& hamilt_lcao, global_out_dir, this->inp_->nspin, - npol, gamma_only_local, istep, iter, @@ -720,15 +718,14 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // overwrite occ_mat.txt with the latest occupation matrix and the // actual charge-density residual; ESolver_KS::iter_finish() above has // computed drho for this electronic step. + const DFTU_BASE::OccmatStepInfo occmat_step{istep, + iter, + this->scf_thr, + this->drho}; DFTU_BASE::write_latest_occmat(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, - this->scf_thr, - this->drho, + occmat_step, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); @@ -738,14 +735,9 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& DFTU_BASE::append_ion_step_snapshot(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, + occmat_step, conv_esolver, true, - this->scf_thr, - this->drho, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index b77f95c9b0..1a10a37c31 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -249,8 +249,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const init_occmat_cfg, this->kv.isk.data(), this->inp_->nspin, - this->inp_->kpar, - PARAM.globalv.npol); + this->inp_->kpar); // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState if (this->inp_->dft_plus_u) @@ -380,22 +379,20 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int this->inp_->scf_nmax, this->inp_->out_occ_mat, this->inp_->dft_plus_u}; + const DFTU_BASE::OccmatStepInfo occmat_step{istep, + iter, + this->scf_thr, + this->drho}; // Read globalv values once into locals so the repeated DFT+U call sites // below do not each re-enter the global dependency surface. const std::string& global_out_dir = PARAM.globalv.global_out_dir; - const int npol = PARAM.globalv.npol; const bool latest_ready = (iter > 1 || istep > 0); if (latest_ready) { DFTU_BASE::write_latest_occmat(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, - this->scf_thr, - this->drho, + occmat_step, occmat_cfg, DFTU_BASE::SOC_LAYOUT_PAULI); } @@ -408,14 +405,9 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int DFTU_BASE::append_ion_step_snapshot(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, + occmat_step, conv_esolver, occmat_ready, - this->scf_thr, - this->drho, occmat_cfg, DFTU_BASE::SOC_LAYOUT_PAULI); diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index 562766505f..437043ef60 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -144,6 +144,7 @@ if(ENABLE_MPI) ../../../source_pw/module_pwdft/dftu_base.cpp ../../../source_pw/module_pwdft/dftu_base_io.cpp ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp @@ -213,6 +214,7 @@ if(ENABLE_MPI) ../../../source_pw/module_pwdft/dftu_base.cpp ../../../source_pw/module_pwdft/dftu_base_io.cpp ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp diff --git a/source/source_lcao/setup_dftu_lcao.cpp b/source/source_lcao/setup_dftu_lcao.cpp index 77fe5e0f5d..5e3f7f33c3 100644 --- a/source/source_lcao/setup_dftu_lcao.cpp +++ b/source/source_lcao/setup_dftu_lcao.cpp @@ -49,7 +49,6 @@ void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, @@ -82,8 +81,7 @@ void finish_dftu_lcao(const bool conv_esolver, DFTU_LCAO::cal_energy_correction(*dftu_ptr, ucell, PARAM.inp.nspin); } } - DFTU_BASE::output(*dftu_ptr, ucell, global_out_dir, nspin, npol, - istep, iter, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + DFTU_BASE::prepare_ion_step_output(global_out_dir, istep, iter, occmat_cfg); /// use the converged occupation matrix for next MD/Relax SCF calculation if (conv_esolver) @@ -103,7 +101,6 @@ template void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, @@ -119,7 +116,6 @@ template void finish_dftu_lcao>(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, diff --git a/source/source_lcao/setup_dftu_lcao.h b/source/source_lcao/setup_dftu_lcao.h index ddd1a79bcd..cd85efa00b 100644 --- a/source/source_lcao/setup_dftu_lcao.h +++ b/source/source_lcao/setup_dftu_lcao.h @@ -56,7 +56,6 @@ void init_dftu_lcao(int dft_plus_u, * @param hamilt_lcao Hamiltonian LCAO object * @param global_out_dir Output directory for occ_mat.txt and occ_matg{#}.txt * @param nspin Number of spin channels (1, 2, or 4) - * @param npol Number of polarizations * @param gamma_only_local Whether only the Gamma point is used for LCAO * @param istep Ionic-step index, starting from 0 * @param iter Electronic-iteration index, starting from 1 @@ -74,7 +73,6 @@ void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index dab7babe6d..2b0af1ec74 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -36,6 +36,7 @@ list( dftu_base.cpp dftu_base_io.cpp dftu_base_io_read.cpp + dftu_base_io_write.cpp dftu_pw.cpp dftu_pw_tools.cpp yukawa_screening.cpp diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 0bb70b999a..35e33c4581 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -184,7 +184,7 @@ void Plus_U_Base::init_base(UnitCell& cell, DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, readin_fn, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); #endif this->set_occmat_ready(); @@ -216,7 +216,7 @@ void Plus_U_Base::init_base(UnitCell& cell, DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, readin_fn, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); #endif this->set_occmat_ready(); } diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index a8f680b33e..fa8023dfaf 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -2,219 +2,26 @@ #include "source_cell/unitcell.h" #include "source_pw/module_pwdft/dftu_base.h" -#include "source_base/constants.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" -#include "source_base/module_external/lapack_connector.h" #include "source_base/parallel_common.h" #include "source_base/parallel_global.h" #include "source_base/timer.h" #include "source_main/version.h" -#include -#include -#include #include #include #include #include #include -// local helpers for eigenvalue calculation -// migrated from dftu_base.cpp, mohan 2025-11-08 -namespace -{ - -inline void JacobiRotate(std::vector>& A, int p, int q, int n) -{ - if (std::abs(A[p][q]) > 1e-10) - { - double r = (A[q][q] - A[p][p]) / (2.0 * A[p][q]); - double t = 0.0; - if (r >= 0) - { - t = 1.0 / (r + sqrt(1.0 + r * r)); - } - else - { - t = -1.0 / (-r + sqrt(1.0 + r * r)); - } - double c = 1.0 / sqrt(1.0 + t * t); - double s = t * c; - - A[p][p] -= t * A[p][q]; - A[q][q] += t * A[p][q]; - A[p][q] = A[q][p] = 0.0; - - for (int k = 0; k < n; k++) - { - if (k != p && k != q) - { - double Akp = c * A[k][p] - s * A[k][q]; - double Akq = s * A[k][p] + c * A[k][q]; - A[k][p] = A[p][k] = Akp; - A[k][q] = A[q][k] = Akq; - } - } - } -} - -inline std::vector CalculateEigenvalues(std::vector>& A, int n) -{ - std::vector eigenvalues(n); - while (true) - { - int p = 0, q = 1; - for (int i = 0; i < n; i++) - { - for (int j = i + 1; j < n; j++) - { - if (std::abs(A[i][j]) > std::abs(A[p][q])) - { - p = i; - q = j; - } - } - } - - if (std::abs(A[p][q]) < 1e-10) - { - for (int i = 0; i < n; i++) - { - eigenvalues[i] = A[i][i]; - } - break; - } - - JacobiRotate(A, p, q, n); - } - return eigenvalues; -} - -/// @brief Extract the four Pauli-component blocks of one SOC correlated shell. -/// -/// The blocks, each (2l+1)x(2l+1) stored row-major with index m0*m+m1, are -/// the charge channel and the sigma_x/y/z spin channels. The PW storage keeps -/// the four blocks contiguous; the LCAO storage is a real symmetric matrix in -/// spin basis, whose spin-off-diagonal imaginary channel is not stored. -/// -/// @param dftu DFT+U object holding the occupation matrix -/// @param iat atom index -/// @param l angular-momentum channel -/// @param layout storage layout of the nspin == 4 occupation matrix -/// @param blocks [out] four Pauli-component blocks -inline void extract_soc_pauli_blocks(const Plus_U_Base& dftu, - const int iat, - const int l, - const DFTU_BASE::OccmatSocLayout layout, - std::vector>& blocks) -{ - const int m = 2 * l + 1; - const int m2 = m * m; - blocks.assign(4, std::vector(m2, 0.0)); - - const OccupationMatrix& occmat = dftu.occmat(); - if (layout == DFTU_BASE::SOC_LAYOUT_PAULI) - { - const ModuleBase::matrix& occ = occmat.mat(iat, l, 0); - for (int is = 0; is < 4; ++is) - { - for (int k = 0; k < m2; ++k) - { - blocks[is][k] = occ.c[is * m2 + k]; - } - } - } - else // SOC_LAYOUT_SPIN_BASIS_REAL - { - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double suu = occmat.get(iat, l, 0, m0, m1); - const double sud = occmat.get(iat, l, 0, m0, m1 + m); - const double sdu = occmat.get(iat, l, 0, m0 + m, m1); - const double sdd = occmat.get(iat, l, 0, m0 + m, m1 + m); - blocks[0][k] = suu + sdd; - blocks[1][k] = sud + sdu; - // The spin-basis storage drops imaginary parts, so the - // sigma_y channel has no contribution. - blocks[2][k] = 0.0; - blocks[3][k] = suu - sdd; - } - } - } -} - -/// @brief Diagonalize the SOC occupation matrix reconstructed in spinor space. -/// -/// The Hermitian spinor matrix is built from the Pauli blocks as -/// n_upup = (b0 + b3)/2, n_dndn = (b0 - b3)/2, -/// n_updown = (b1 + i*b2)/2. -/// Its 2*(2l+1) real eigenvalues are computed with LAPACK zheev. -/// -/// @param blocks four Pauli-component blocks -/// @param m 2l+1 -/// @return eigenvalues sorted ascending by zheev -inline std::vector calculate_spinor_eigenvalues( - const std::vector>& blocks, - const int m) -{ - const int n = 2 * m; - ModuleBase::ComplexMatrix spinor_n(n, n); - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double b0 = blocks[0][k]; - const double b1 = blocks[1][k]; - const double b2 = blocks[2][k]; - const double b3 = blocks[3][k]; - - spinor_n(m0, m1) = std::complex(0.5 * (b0 + b3), 0.0); - spinor_n(m0, m1 + m) = std::complex(0.5 * b1, 0.5 * b2); - spinor_n(m0 + m, m1) = std::complex(0.5 * b1, -0.5 * b2); - spinor_n(m0 + m, m1 + m) = std::complex(0.5 * (b0 - b3), 0.0); - } - } - - std::vector eigenvalues(n, 0.0); - std::vector> work(1, std::complex(0.0, 0.0)); - std::vector rwork(std::max(1, 3 * n - 2), 0.0); - int info = 0; - int lwork = -1; - - LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), - work.data(), lwork, rwork.data(), &info); - if (info != 0) - { - ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", - "LAPACK zheev workspace query failed"); - } - lwork = std::max(1, static_cast(work[0].real())); - work.resize(lwork); - LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), - work.data(), lwork, rwork.data(), &info); - if (info != 0) - { - ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", - "LAPACK zheev diagonalization failed"); - } - return eigenvalues; -} - -} // namespace - - -namespace +namespace DFTU_BASE { // Write the STRU-style provenance header shared by occ_mat.txt and the // per-ionic-step snapshot files: ABACUS version, local timestamp and the // 1-based relaxation step this file belongs to. -void write_provenance_header(std::ostream& os, const int istep) +static void write_provenance_header(std::ostream& os, const int istep) { std::time_t now = std::time(nullptr); std::tm tm_local; @@ -223,14 +30,28 @@ void write_provenance_header(std::ostream& os, const int istep) std::strftime(time_buf, sizeof(time_buf), "%Y-%m-%d %H:%M:%S", &tm_local); os << "# ABACUS version: " << VERSION << "\n"; os << "# Written at " << time_buf << "\n"; - os << "# RELAX STEP " << istep + 1 << "\n"; + os << "# RELAX STEP " << istep + 1 << "\n\n"; } -} // namespace - - -namespace DFTU_BASE +static void validate_output_configuration(const char* context, + const int istep, + const int iter, + const OccmatOutputCfg& cfg) { + if (istep < 0) + { + ModuleBase::WARNING_QUIT(context, "istep must be >= 0"); + } + if (iter < 1) + { + ModuleBase::WARNING_QUIT(context, "iter must be >= 1"); + } + if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) + { + ModuleBase::WARNING_QUIT(context, + "invalid occupation-matrix output frequency configuration"); + } +} bool is_ion_step_output_step(const int istep, const OccmatOutputCfg& cfg) { @@ -280,11 +101,10 @@ std::string find_first_existing_file(const std::string& dir, /// element. void local_occup_bcast(const UnitCell& ucell, OccupationMatrix& occ, - const std::vector& l_channel, - int nspin, - int npol) + const std::vector& l_channel) { ModuleBase::TITLE("DFTU_BASE", "local_occup_bcast"); + const int nspin = occ.nspin(); for (int T = 0; T < ucell.ntype; T++) { @@ -292,43 +112,31 @@ void local_occup_bcast(const UnitCell& ucell, { continue; } - + const int l = l_channel[T]; for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); - const int L = l_channel[T]; - - for (int l = 0; l <= ucell.atoms[T].nwl; l++) + if (nspin == 1 || nspin == 2) { - if (l != l_channel[T]) - { - continue; - } - - if (nspin == 1 || nspin == 2) + for (int spin = 0; spin < 2; spin++) { - for (int spin = 0; spin < 2; spin++) - { - Parallel_Common::bcast_double(occ.mat(iat, l, spin).c, - occ.mat(iat, l, spin).nr * occ.mat(iat, l, spin).nc); - } - } - else if (nspin == 4) // SOC - { - Parallel_Common::bcast_double(occ.mat(iat, l, 0).c, - occ.mat(iat, l, 0).nr * occ.mat(iat, l, 0).nc); + Parallel_Common::bcast_double(occ.mat(iat, l, spin).c, + occ.mat(iat, l, spin).nr * occ.mat(iat, l, spin).nc); } } + else if (nspin == 4) // SOC + { + Parallel_Common::bcast_double(occ.mat(iat, l, 0).c, + occ.mat(iat, l, 0).nr * occ.mat(iat, l, 0).nc); + } } } - return; } #endif -void prepare_ion_step_file(const std::string& global_out_dir, - const int istep, - const OccmatOutputCfg& cfg) +static void prepare_ion_step_file(const std::string& global_out_dir, + const int istep) { // Only the root process creates the per-ionic-step file. if (GlobalV::MY_RANK != 0) @@ -348,59 +156,23 @@ void prepare_ion_step_file(const std::string& global_out_dir, } // Provenance header shared with occ_mat.txt. write_provenance_header(ofs_ion_step, istep); - ofs_ion_step.close(); - return; } - -void output(const Plus_U_Base& dftu, - const UnitCell& ucell, - const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) +void prepare_ion_step_output(const std::string& global_out_dir, + const int istep, + const int iter, + const OccmatOutputCfg& cfg) { - ModuleBase::TITLE("DFTU_BASE", "output"); - - if (istep < 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", "istep must be >= 0"); - } - if (iter < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", "iter must be >= 1"); - } - if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", - "invalid occupation-matrix output frequency configuration"); - } - - // The per-type U values and the per-atom occupation matrices are not - // dumped into the running log here: for large cells (thousands of - // atoms) this floods the log at every electronic step. The dedicated - // occupation-matrix files (occ_mat.txt, occ_matg{#}.txt) gated by - // out_occ_mat carry the same information instead. - - // occ_mat.txt is not written here: its section records drho, which - // is only available after the electronic solve. write_latest_occmat() - // overwrites it at the iter_finish stage. - - // At the first electronic step of an output ionic step, (re)create the - // per-ionic-step file and write its header. The LCAO path enters here - // directly; the PW path additionally prepares the file from - // iter_init_dftu_pw(), because it skips output() at istep 0 / iter 1. - // Electronic-step sections are appended later by append_ion_step_snapshot(). + ModuleBase::TITLE("DFTU_BASE", "prepare_ion_step_output"); + validate_output_configuration("DFTU_BASE::prepare_ion_step_output", + istep, + iter, + cfg); const bool ion_step_output = is_ion_step_output_step(istep, cfg); if (cfg.out_occ_mat && ion_step_output && iter == 1) { - prepare_ion_step_file(global_out_dir, istep, cfg); + prepare_ion_step_file(global_out_dir, istep); } - - return; } // Write one electronic-step section shared by occ_mat.txt and the @@ -408,20 +180,16 @@ void output(const Plus_U_Base& dftu, // charge-density convergence threshold, the actual charge-density // residual and the full occupation matrices. When occmat_ready is false, // the matrix body is an "N/A" placeholder. -void write_snapshot_section(std::ostream& os, - const Plus_U_Base& dftu, - const UnitCell& ucell, - const int iter, - const double scf_thr, - const double drho, - const bool occmat_ready, - const int nspin, - const int npol, - const OccmatSocLayout soc_layout) +static void write_snapshot_section(std::ostream& os, + const Plus_U_Base& dftu, + const UnitCell& ucell, + const OccmatStepInfo& step, + const bool occmat_ready, + const OccmatSocLayout soc_layout) { - os << "# Electronic step " << iter << "\n"; - os << "# scf_thr " << std::scientific << std::setprecision(8) << scf_thr << "\n"; - os << "# drho " << std::scientific << std::setprecision(8) << drho << "\n"; + os << "# Electronic step " << step.iter << "\n"; + os << "# scf_thr " << std::scientific << std::setprecision(8) << step.scf_thr << "\n"; + os << "# drho " << std::scientific << std::setprecision(8) << step.drho << "\n"; // The PW path has no occupation matrix at istep 0 / iter 1 unless it was // loaded from file; record an explicit placeholder instead of a silent @@ -429,28 +197,23 @@ void write_snapshot_section(std::ostream& os, // intentionally rejected by the strict restart reader. if (!occmat_ready) { - os << "\n N/A\n"; + os << "\n N/A\n\n"; return; } // Full occupation matrices together with per-atom magnetism. - write_occup_m(dftu, ucell, os, true, nspin, npol, - OCMAT_FMT_READABLE, soc_layout); + write_occup_m(dftu, ucell, os, soc_layout); + os << "\n"; } void append_ion_step_snapshot(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - bool conv_esolver, - bool occmat_ready, - double scf_thr, - double drho, + const OccmatStepInfo& step, + const bool conv_esolver, + const bool occmat_ready, const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) + const OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "append_ion_step_snapshot"); @@ -458,7 +221,7 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, // (Plus_U_Base::l_channel stays empty); writing here would dereference // a null l_channel.data() in write_occup_m. Match the documented contract // that out_occ_mat only takes effect for DFT+U calculations. Check this - // before validating iter/istep/nspin: when DFT+U is off this function is + // before validating the step: when DFT+U is off this function is // a no-op, so an invalid iter (e.g. from an upstream bug) must not crash // here. Mirrors the guard-first pattern in write_latest_occmat. if (!cfg.out_occ_mat || cfg.dft_plus_u <= 0) @@ -466,51 +229,19 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, return; } - if (nspin != 1 && nspin != 2 && nspin != 4) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "nspin must be 1, 2 or 4"); - } - if (istep < 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "istep must be >= 0"); - } - if (iter < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "iter must be >= 1"); - } - // out_freq_ion == 0 is the valid default (no numbered file) and is - // handled by the gates below; only negative values are an error. - if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", - "invalid occupation-matrix output frequency configuration"); - } + validate_output_configuration("DFTU_BASE::append_ion_step_snapshot", + step.istep, + step.iter, + cfg); - const bool ion_step_output = is_ion_step_output_step(istep, cfg); - const bool elec_trigger = is_elec_snapshot_trigger(iter, conv_esolver, cfg); + const bool ion_step_output = is_ion_step_output_step(step.istep, cfg); + const bool elec_trigger = is_elec_snapshot_trigger(step.iter, conv_esolver, cfg); if (!ion_step_output || !elec_trigger || GlobalV::MY_RANK != 0) { return; } - const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, istep); - - // The provenance header written by prepare_ion_step_file() must not be - // counted as a recorded section: scan the file for an existing section - // marker instead of testing the stream position. - bool section_exists = false; - { - std::ifstream checker(ion_step_fn.c_str()); - std::string line; - while (std::getline(checker, line)) - { - if (line.rfind("# Electronic step", 0) == 0) - { - section_exists = true; - break; - } - } - } + const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, step.istep); std::ofstream ofs_ion_step; ofs_ion_step.open(ion_step_fn, std::ios::app); @@ -520,34 +251,25 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, "Can't open per-ionic-step occupation-matrix file"); } - // Separate consecutive sections by a blank line; the first section - // directly follows the provenance header. - if (section_exists) + // A direct caller may not have prepared the file at iter 1. Add the + // provenance header when append mode opens an empty file. + if (ofs_ion_step.tellp() == std::streampos(0)) { - ofs_ion_step << "\n"; + write_provenance_header(ofs_ion_step, step.istep); } - write_snapshot_section(ofs_ion_step, dftu, ucell, iter, scf_thr, drho, - occmat_ready, nspin, npol, soc_layout); - - ofs_ion_step.close(); + write_snapshot_section(ofs_ion_step, dftu, ucell, step, occmat_ready, soc_layout); - GlobalV::ofs_running << " Append local occupation matrices of electronic step " << iter + GlobalV::ofs_running << " Append local occupation matrices of electronic step " << step.iter << " to file: " << ion_step_fn << std::endl; - return; } void write_latest_occmat(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - double scf_thr, - double drho, + const OccmatStepInfo& step, const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) + const OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "write_latest_occmat"); @@ -566,429 +288,11 @@ void write_latest_occmat(const Plus_U_Base& dftu, // occ_mat.txt is a single-section snapshot file: the provenance // header shared with the g files followed by the same section writer. - write_provenance_header(ofdftu, istep); + write_provenance_header(ofdftu, step.istep); // Callers of write_latest_occmat() guarantee the matrix exists: the PW // esolver gates on latest_ready, the LCAO esolver computes it every // electronic iteration. - write_snapshot_section(ofdftu, dftu, ucell, iter, scf_thr, drho, - true, nspin, npol, soc_layout); - - ofdftu.close(); - return; -} - - -void write_occup_m(const Plus_U_Base& dftu, - const UnitCell& ucell, - std::ostream& ofs, - bool diag, - int nspin, - int npol, - OccmatTextFormat fmt, - OccmatSocLayout soc_layout) -{ - ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); - - if (GlobalV::MY_RANK != 0) - { - return; - } - - for (int T = 0; T < ucell.ntype; T++) - { - if (!dftu.has_l_channel(T)) - { - continue; - } - const int NL = ucell.atoms[T].nwl + 1; - const int LC = dftu.get_l_channel(T); - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - - for (int l = 0; l < NL; l++) - { - if (l != dftu.get_l_channel(T)) - { - continue; - } - - const std::string& elem_label = ucell.atoms[T].label; - const int index_in_type = I + 1; - const ModuleBase::Vector3& tau = ucell.atoms[T].tau[I]; - const bool readable = (fmt == OCMAT_FMT_READABLE); - - // In the compact collinear/SOC layout the header carries the - // per-atom magnetism, which is only known after the Pauli - // block traces; that header is written in the spin branches. - const bool defer_header = readable - && (nspin == 1 || nspin == 2 || nspin == 4); - if (!defer_header) - { - if (readable) - { - // Readable header, e.g. " Fe Atom 1 L 2". - // The position is not repeated: it is recorded in - // STRU and constant within one ionic (g) step. - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l << std::endl; - } - else - { - // Legacy diagnostic layout. The compact-only restart - // reader intentionally rejects this representation. - ofs << "\n" << elem_label << " " << index_in_type; - ofs << " Atom=" << iat + 1; - ofs << " L=" << l; - ofs << std::setprecision(8) << std::fixed - << std::setw(14) << tau.x - << std::setw(14) << tau.y - << std::setw(14) << tau.z << std::endl; - } - } - - if (nspin == 1 || nspin == 2) - { - const int nm = 2 * l + 1; - double sum0[2] = {0.0, 0.0}; - - if (readable) - { - // Compact snapshot layout, e.g. - // " Fe Atom 1 L 2 mag 0.84215429" - // " spin 1 nelec 1.76153075" - // followed by the matrix rows. - if (diag) - { - for (int is = 0; is < 2; is++) - { - for (int m0 = 0; m0 < nm; m0++) - { - sum0[is] += dftu.occmat().get(iat, l, is, m0, m0); - } - } - } - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l; - if (diag) - { - ofs << " mag " << std::fixed << std::setprecision(8) - << sum0[0] - sum0[1]; - } - ofs << std::endl; - - ofs << std::fixed << std::setprecision(8); - for (int is = 0; is < 2; is++) - { - ofs << " spin " << is + 1; - if (diag) - { - ofs << " nelec " << sum0[is]; - } - ofs << std::endl; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - // One separating space before every - // value; positive values are 10 chars, - // negative ones 11. - ofs << " " << std::setw(10) - << dftu.occmat().get(iat, l, is, m0, m1); - } - ofs << std::endl; - } - } - } - else - { - const std::string trace_label = " Trace (electrons) "; - const std::string matrix_label = " spin= "; - for (int is = 0; is < 2; is++) - { - if (diag) - { - std::vector> A(nm, - std::vector(nm)); - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - A[m0][m1] = dftu.occmat().get(iat, l, is, m0, m1); - } - } - std::vector eigenvalues = CalculateEigenvalues(A, nm); - sum0[is] = 0.0; - ofs << " Eigenvalues for spin=" << is + 1 << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < nm; i++) - { - ofs << std::setw(12) << eigenvalues[i]; - sum0[is] += eigenvalues[i]; - } - ofs << std::endl; - ofs << std::fixed << std::setprecision(8); - ofs << trace_label << sum0[is] << std::endl; - } - ofs << matrix_label << is + 1 << std::endl; - // Matrix elements use fixed-point notation with 8 - // fractional digits (~1e-8 absolute precision). - ofs << std::fixed << std::setprecision(8); - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - // Legacy keeps its historical - // 12-character column alignment. - ofs << " " << std::setw(11) - << dftu.occmat().get(iat, l, is, m0, m1); - } - ofs << std::endl; - } - } - if (diag) - { - ofs << " mag " << std::fixed << std::setprecision(8) - << sum0[0] - sum0[1] << std::endl; - } - } - } - else if (nspin == 4) // SOC - { - const int m = 2 * l + 1; - std::vector> blocks; - extract_soc_pauli_blocks(dftu, iat, l, soc_layout, blocks); - - if (readable) - { - // Compact snapshot layout aligned with the - // nspin == 2 style, e.g. - // " Fe Atom 1 L 2 mag mx my mz" - // " spin 1 nelec ..." - // followed by the four spin blocks. - double mag_x = 0.0; - double mag_y = 0.0; - double mag_z = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - mag_x += blocks[1][k]; - mag_y += blocks[2][k]; - mag_z += blocks[3][k]; - } - // Spin-resolved electron counts: traces of the - // up/up and down/down blocks. - double n_up = 0.0; - double n_down = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - n_up += 0.5 * (blocks[0][k] + blocks[3][k]); - n_down += 0.5 * (blocks[0][k] - blocks[3][k]); - } - - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l; - if (diag) - { - ofs << " mag " << std::fixed << std::setprecision(8) - << mag_x << " " << mag_y << " " << mag_z; - } - ofs << std::endl; - - ofs << std::fixed << std::setprecision(8); - - // spin 1: spin up/up block (Re). - ofs << " spin 1"; - if (diag) - { - ofs << " nelec " << n_up; - } - ofs << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] + blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // Off-diagonal Re block. - ofs << " spin 12 re" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[1][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // Off-diagonal Im block. - ofs << " spin 12 im" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[2][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // spin 2: spin down/down block (Re). - ofs << " spin 2"; - if (diag) - { - ofs << " nelec " << n_down; - } - ofs << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] - blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - } - else if (diag) - { - // Occupation numbers of the correlated spinor shell. - std::vector eigenvalues = calculate_spinor_eigenvalues(blocks, m); - ofs << " Eigenvalues in spinor space:" << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < 2 * m; ++i) - { - ofs << std::setw(12) << eigenvalues[i]; - if ((i + 1) % m == 0) - { - ofs << std::endl; - } - } - - // Per-atom charge and spin-moment traces of the shell. - double n_tot = 0.0; - double mag_x = 0.0; - double mag_y = 0.0; - double mag_z = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - n_tot += blocks[0][k]; - mag_x += blocks[1][k]; - mag_y += blocks[2][k]; - mag_z += blocks[3][k]; - } - ofs << std::fixed << std::setprecision(8); - ofs << " Trace (electrons) " << n_tot << std::endl; - ofs << " Magnetism for atom " << iat + 1 << " (mx, my, mz):" - << std::setw(12) << mag_x << std::setw(12) << mag_y - << std::setw(12) << mag_z << std::endl; - - // Matrix elements use fixed-point notation with 8 - // fractional digits (~1e-8 absolute precision). - ofs << std::fixed << std::setprecision(8); - - // Hermitian occupation matrix shown spin block by spin block. - ofs << " Occupation matrix, spin up/up (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] + blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin up/down (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[1][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin up/down (Im):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[2][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin down/down (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] - blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - } - else - { - // Legacy real 2m x 2m matrix in spin basis; retained - // for diagnostic output only. - for (int m0 = 0; m0 < m; ++m0) - { - for (int ipol0 = 0; ipol0 < npol; ++ipol0) - { - for (int m1 = 0; m1 < m; ++m1) - { - for (int ipol1 = 0; ipol1 < npol; ++ipol1) - { - const int k = m0 * m + m1; - double val = 0.0; - if (ipol0 == 0 && ipol1 == 0) - { - val = 0.5 * (blocks[0][k] + blocks[3][k]); - } - else if (ipol0 == 1 && ipol1 == 1) - { - val = 0.5 * (blocks[0][k] - blocks[3][k]); - } - else - { - val = 0.5 * blocks[1][k]; - } - ofs << " " << std::setw(10) - << std::setprecision(8) - << std::fixed << val; - } - } - ofs << std::endl; - } - } - } - } - } // l - } // I - } // T - - return; + write_snapshot_section(ofdftu, dftu, ucell, step, true, soc_layout); } diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index babcc76908..4750cf0fbe 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -56,36 +56,24 @@ std::string find_first_existing_file(const std::string& dir, /// of the current electronic step, followed by the occupation matrices and /// per-atom magnetism. When occmat_ready is false, an "N/A" placeholder is /// recorded instead of the matrix body (the PW path has no matrix at -/// istep 0 / iter 1 unless it was loaded from file). The caller is -/// responsible for truncating the file at the first electronic step of the -/// ionic step (see prepare_ion_step_file()) and for invoking this function -/// only after drho of the current electronic step is computed. +/// istep 0 / iter 1 unless it was loaded from file). The caller invokes this +/// function only after drho of the current electronic step is computed. /// /// @param dftu DFT+U object holding the occupation matrices /// @param ucell unit cell /// @param global_out_dir output directory (including the trailing separator) -/// @param nspin number of spin components (1, 2 or 4) -/// @param npol number of polarizations -/// @param istep ionic-step index, starting from 0 -/// @param iter electronic-iteration index, starting from 1 +/// @param step electronic-step metadata and convergence values /// @param conv_esolver whether the electronic SCF is converged at this step /// @param occmat_ready whether the occupation matrix of this step exists; /// false records an "N/A" placeholder instead of the matrix body -/// @param scf_thr configured charge-density convergence threshold -/// @param drho actual charge-density residual of the current electronic step /// @param cfg frequency configuration /// @param soc_layout storage layout of the nspin == 4 occupation matrix void append_ion_step_snapshot(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, + const OccmatStepInfo& step, bool conv_esolver, bool occmat_ready, - double scf_thr, - double drho, const OccmatOutputCfg& cfg, OccmatSocLayout soc_layout); @@ -113,46 +101,18 @@ void read_occup_m(const UnitCell& ucell, /// Implemented in dftu_base_io.cpp (only available in MPI builds). void local_occup_bcast(const UnitCell& ucell, OccupationMatrix& occ, - const std::vector& l_channel, - int nspin, - int npol); + const std::vector& l_channel); /// Create (or truncate) the per-ionic-step file occ_matg{istep+1}.txt and /// write its header (rank 0 only). /// -/// Must be called once at the first electronic step (iter == 1) of an output -/// ionic step, before any append_ion_step_snapshot() section. Truncating here -/// also guarantees that a rerun in the same output directory can not append -/// snapshots of a previous calculation. -void prepare_ion_step_file(const std::string& global_out_dir, - const int istep, - const OccmatOutputCfg& cfg); - -/// Output DFT+U information (Hubbard U/J, local occupation matrices) to the -/// running log. -/// -/// When cfg.out_occ_mat is true, cfg.out_freq_ion is positive and istep is an -/// output ionic step, the per-ionic-step file occ_matg{istep+1}.txt is -/// created: at the first electronic step (iter == 1) it is truncated and -/// initialized with a provenance header, then append_ion_step_snapshot() -/// appends one section per recorded electronic step. -/// -/// Note: occ_mat.txt is NOT written here. It records the actual -/// charge-density residual drho, which is only known after the electronic -/// solve, so write_latest_occmat() writes it from the iter_finish stage. -/// -/// Extracted from Plus_U_Base::output as a free function so that IO logic is -/// decoupled from the Plus_U_Base class. The function only reads the -/// Plus_U_Base state via public accessors; no friend declaration needed. -void output(const Plus_U_Base& dftu, - const UnitCell& ucell, - const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout); +/// At iter == 1 of an enabled output ionic step, truncates the file and writes +/// its provenance header. Other calls are no-ops. Truncating here guarantees +/// that a rerun can not append snapshots from a previous calculation. +void prepare_ion_step_output(const std::string& global_out_dir, + int istep, + int iter, + const OccmatOutputCfg& cfg); /// Overwrite occ_mat.txt with the occupation matrix of the current /// electronic step (rank 0 only). @@ -167,41 +127,26 @@ void output(const Plus_U_Base& dftu, /// @param dftu DFT+U object holding the occupation matrices /// @param ucell unit cell /// @param global_out_dir output directory (including the trailing separator) -/// @param nspin number of spin components (1, 2 or 4) -/// @param npol number of polarizations -/// @param istep ionic-step index, starting from 0 -/// @param iter electronic-iteration index, starting from 1 -/// @param scf_thr configured charge-density convergence threshold -/// @param drho actual charge-density residual of the current electronic step +/// @param step electronic-step metadata and convergence values /// @param cfg output configuration; only cfg.out_occ_mat is consulted here /// @param soc_layout storage layout of the nspin == 4 occupation matrix void write_latest_occmat(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - double scf_thr, - double drho, + const OccmatStepInfo& step, const OccmatOutputCfg& cfg, OccmatSocLayout soc_layout); /// Write local occupation matrices to the given stream. /// -/// When diag is true, eigenvalues and magnetism are also printed; otherwise -/// only raw matrix elements. fmt selects between the legacy token layout -/// (used for running log) and the readable snapshot layout of -/// occ_mat.txt and occ_matg{#}.txt. soc_layout tells how the nspin == 4 -/// storage is arranged (PW Pauli blocks or LCAO real spin-basis matrix). +/// Writes the compact snapshot layout accepted by read_occup_m(). The output +/// includes per-atom magnetism, electron counts and the full matrix blocks. +/// soc_layout tells how the nspin == 4 storage is arranged (PW Pauli blocks +/// or LCAO real spin-basis matrix). /// Caller is responsible for opening/closing the stream. void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, std::ostream& ofs, - bool diag, - int nspin, - int npol, - OccmatTextFormat fmt, OccmatSocLayout soc_layout); } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_base_io_write.cpp b/source/source_pw/module_pwdft/dftu_base_io_write.cpp new file mode 100644 index 0000000000..d10fe5038c --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_io_write.cpp @@ -0,0 +1,329 @@ +#include "source_pw/module_pwdft/dftu_base_io.h" + +#include "source_base/global_function.h" +#include "source_base/module_external/lapack_connector.h" +#include "source_cell/unitcell.h" +#include "source_pw/module_pwdft/dftu_base.h" + +#include +#include +#include +#include +#include +#include + +namespace DFTU_BASE +{ + +struct SocCompactBlocks +{ + std::vector up_up; + std::vector up_down_real; + std::vector up_down_imag; + std::vector down_down; +}; + +// Retained for the follow-up occupation-eigensystem output refactor. The +// compact writer does not currently call these helpers: this PR only removes +// the legacy text format and separates the writer from the file orchestration. +namespace retained_eigenvalue_helpers +{ + +inline void jacobi_rotate(std::vector>& matrix, + const int p, + const int q, + const int size) +{ + if (std::abs(matrix[p][q]) <= 1.0e-10) + { + return; + } + + const double ratio = (matrix[q][q] - matrix[p][p]) / (2.0 * matrix[p][q]); + const double tangent = ratio >= 0.0 + ? 1.0 / (ratio + std::sqrt(1.0 + ratio * ratio)) + : -1.0 / (-ratio + std::sqrt(1.0 + ratio * ratio)); + const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent); + const double sine = tangent * cosine; + + matrix[p][p] -= tangent * matrix[p][q]; + matrix[q][q] += tangent * matrix[p][q]; + matrix[p][q] = matrix[q][p] = 0.0; + + for (int k = 0; k < size; ++k) + { + if (k == p || k == q) + { + continue; + } + const double matrix_kp = cosine * matrix[k][p] - sine * matrix[k][q]; + const double matrix_kq = sine * matrix[k][p] + cosine * matrix[k][q]; + matrix[k][p] = matrix[p][k] = matrix_kp; + matrix[k][q] = matrix[q][k] = matrix_kq; + } +} + +inline std::vector calculate_eigenvalues(std::vector>& matrix, + const int size) +{ + std::vector eigenvalues(size); + while (true) + { + int p = 0; + int q = 1; + for (int row = 0; row < size; ++row) + { + for (int column = row + 1; column < size; ++column) + { + if (std::abs(matrix[row][column]) > std::abs(matrix[p][q])) + { + p = row; + q = column; + } + } + } + + if (std::abs(matrix[p][q]) < 1.0e-10) + { + for (int i = 0; i < size; ++i) + { + eigenvalues[i] = matrix[i][i]; + } + break; + } + jacobi_rotate(matrix, p, q, size); + } + return eigenvalues; +} + +inline std::vector calculate_spinor_eigenvalues(const SocCompactBlocks& blocks, + const int orbital_size) +{ + const int spinor_size = 2 * orbital_size; + ModuleBase::ComplexMatrix spinor_matrix(spinor_size, spinor_size); + for (int row = 0; row < orbital_size; ++row) + { + for (int column = 0; column < orbital_size; ++column) + { + const int index = row * orbital_size + column; + spinor_matrix(row, column) = blocks.up_up[index]; + spinor_matrix(row, column + orbital_size) + = std::complex(blocks.up_down_real[index], + blocks.up_down_imag[index]); + spinor_matrix(row + orbital_size, column) + = std::complex(blocks.up_down_real[index], + -blocks.up_down_imag[index]); + spinor_matrix(row + orbital_size, column + orbital_size) + = blocks.down_down[index]; + } + } + + std::vector eigenvalues(spinor_size, 0.0); + std::vector> work(1, std::complex(0.0, 0.0)); + std::vector rwork(std::max(1, 3 * spinor_size - 2), 0.0); + int info = 0; + int lwork = -1; + LapackConnector::zheev('N', + 'U', + spinor_size, + spinor_matrix, + spinor_size, + eigenvalues.data(), + work.data(), + lwork, + rwork.data(), + &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev workspace query failed"); + } + + lwork = std::max(1, static_cast(work[0].real())); + work.resize(lwork); + LapackConnector::zheev('N', + 'U', + spinor_size, + spinor_matrix, + spinor_size, + eigenvalues.data(), + work.data(), + lwork, + rwork.data(), + &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev diagonalization failed"); + } + return eigenvalues; +} + +} // namespace retained_eigenvalue_helpers + +static SocCompactBlocks extract_soc_compact_blocks(const Plus_U_Base& dftu, + const int iat, + const int l, + const OccmatSocLayout layout) +{ + const int size = 2 * l + 1; + const int block_size = size * size; + SocCompactBlocks blocks; + blocks.up_up.resize(block_size, 0.0); + blocks.up_down_real.resize(block_size, 0.0); + blocks.up_down_imag.resize(block_size, 0.0); + blocks.down_down.resize(block_size, 0.0); + + const OccupationMatrix& occmat = dftu.occmat(); + if (layout == SOC_LAYOUT_PAULI) + { + const ModuleBase::matrix& occ = occmat.mat(iat, l, 0); + for (int index = 0; index < block_size; ++index) + { + const double charge = occ.c[0 * block_size + index]; + const double spin_x = occ.c[1 * block_size + index]; + const double spin_y = occ.c[2 * block_size + index]; + const double spin_z = occ.c[3 * block_size + index]; + blocks.up_up[index] = 0.5 * (charge + spin_z); + blocks.up_down_real[index] = 0.5 * spin_x; + blocks.up_down_imag[index] = 0.5 * spin_y; + blocks.down_down[index] = 0.5 * (charge - spin_z); + } + } + else + { + for (int row = 0; row < size; ++row) + { + for (int column = 0; column < size; ++column) + { + const int index = row * size + column; + const double up_up = occmat.get(iat, l, 0, row, column); + const double up_down = occmat.get(iat, l, 0, row, column + size); + const double down_up = occmat.get(iat, l, 0, row + size, column); + const double down_down = occmat.get(iat, l, 0, row + size, column + size); + blocks.up_up[index] = up_up; + blocks.up_down_real[index] = 0.5 * (up_down + down_up); + blocks.down_down[index] = down_down; + } + } + } + return blocks; +} + +static void write_matrix(std::ostream& os, const double* matrix, const int size) +{ + for (int row = 0; row < size; ++row) + { + for (int column = 0; column < size; ++column) + { + os << " " << std::setw(10) << matrix[row * size + column]; + } + os << std::endl; + } +} + +static void write_collinear_atom(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const int type, + const int iat, + const int l) +{ + const int size = 2 * l + 1; + double electron_count[2] = {0.0, 0.0}; + for (int spin = 0; spin < 2; ++spin) + { + for (int m = 0; m < size; ++m) + { + electron_count[spin] += dftu.occmat().get(iat, l, spin, m, m); + } + } + + os << "\n " << ucell.atoms[type].label; + os << " Atom " << iat + 1; + os << " L " << l; + os << " mag " << std::fixed << std::setprecision(8) + << electron_count[0] - electron_count[1] << std::endl; + + os << std::fixed << std::setprecision(8); + for (int spin = 0; spin < 2; ++spin) + { + os << " spin " << spin + 1; + os << " nelec " << electron_count[spin] << std::endl; + write_matrix(os, dftu.occmat().mat(iat, l, spin).c, size); + } +} + +static void write_soc_atom(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const int type, + const int iat, + const int l, + const OccmatSocLayout layout) +{ + const int size = 2 * l + 1; + const SocCompactBlocks blocks = extract_soc_compact_blocks(dftu, iat, l, layout); + double electron_up = 0.0; + double electron_down = 0.0; + double magnetization_x = 0.0; + double magnetization_y = 0.0; + double magnetization_z = 0.0; + for (int m = 0; m < size; ++m) + { + const int index = m * size + m; + electron_up += blocks.up_up[index]; + electron_down += blocks.down_down[index]; + magnetization_x += 2.0 * blocks.up_down_real[index]; + magnetization_y += 2.0 * blocks.up_down_imag[index]; + magnetization_z += blocks.up_up[index] - blocks.down_down[index]; + } + + os << "\n " << ucell.atoms[type].label; + os << " Atom " << iat + 1; + os << " L " << l; + os << " mag " << std::fixed << std::setprecision(8) + << magnetization_x << " " << magnetization_y << " " << magnetization_z + << std::endl; + os << std::fixed << std::setprecision(8); + + os << " spin 1 nelec " << electron_up << std::endl; + write_matrix(os, blocks.up_up.data(), size); + os << " spin 12 re" << std::endl; + write_matrix(os, blocks.up_down_real.data(), size); + os << " spin 12 im" << std::endl; + write_matrix(os, blocks.up_down_imag.data(), size); + os << " spin 2 nelec " << electron_down << std::endl; + write_matrix(os, blocks.down_down.data(), size); +} + +void write_occup_m(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); + const int nspin = dftu.occmat().nspin(); + for (int type = 0; type < ucell.ntype; ++type) + { + if (!dftu.has_l_channel(type)) + { + continue; + } + const int l = dftu.get_l_channel(type); + for (int atom = 0; atom < ucell.atoms[type].na; ++atom) + { + const int iat = ucell.itia2iat(type, atom); + if (nspin == 4) + { + write_soc_atom(dftu, ucell, os, type, iat, l, soc_layout); + } + else + { + write_collinear_atom(dftu, ucell, os, type, iat, l); + } + } + } +} + +} // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_types.h b/source/source_pw/module_pwdft/dftu_types.h index e1f1bee13d..a2e35b9bd8 100644 --- a/source/source_pw/module_pwdft/dftu_types.h +++ b/source/source_pw/module_pwdft/dftu_types.h @@ -23,11 +23,13 @@ struct OccmatOutputCfg int dft_plus_u; }; -/// Text format used by write_occup_m(). -enum OccmatTextFormat +/// Electronic-step metadata recorded in an occupation-matrix snapshot. +struct OccmatStepInfo { - OCMAT_FMT_LEGACY, ///< legacy diagnostic layout; not accepted on input - OCMAT_FMT_READABLE ///< compact human-readable layout for snapshot files + int istep; + int iter; + double scf_thr; + double drho; }; /// Storage layout of the real (2l+1)*npol matrix when nspin == 4. diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index 6e0332eac8..877445b88b 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -18,24 +18,15 @@ void iter_init_dftu_pw(const int iter, const OccmatOutputCfg& occmat_cfg, const int* isk, const int nspin, - const int kpar, - const int npol) + const int kpar) { if (!p_chgmix || occmat_cfg.dft_plus_u <= 0) { return; } - // Prepare the per-ionic-step file at the first electronic step, before - // the first occupation matrix exists. output() is skipped below at - // istep 0 / iter 1; without this call the g1 file would miss its header - // and could inherit stale sections from a previous run. - if (iter == 1 && occmat_cfg.out_occ_mat - && DFTU_BASE::is_ion_step_output_step(istep, occmat_cfg)) - { - DFTU_BASE::prepare_ion_step_file(global_out_dir, - istep, occmat_cfg); - } + // Prepare the per-ionic-step file before the first snapshot is appended. + DFTU_BASE::prepare_ion_step_output(global_out_dir, istep, iter, occmat_cfg); const bool first_scf_step = (iter == 1 && istep == 0); @@ -75,16 +66,6 @@ void iter_init_dftu_pw(const int iter, dftu.get_uterm_mat(), dftu.energy_ref()); } - // Keep the historical output cadence: output() is skipped at the very - // first electronic step. - if (first_scf_step) - { - return; - } - - DFTU_BASE::output(dftu, ucell, global_out_dir, - nspin, npol, istep, iter, occmat_cfg, - DFTU_BASE::SOC_LAYOUT_PAULI); } } diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.h b/source/source_pw/module_pwdft/setup_dftu_pw.h index 05c448e216..0b242ce852 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.h +++ b/source/source_pw/module_pwdft/setup_dftu_pw.h @@ -26,8 +26,7 @@ void iter_init_dftu_pw(const int iter, const OccmatOutputCfg& occmat_cfg, const int* isk, const int nspin, - const int kpar, - const int npol); + const int kpar); } diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 5fc6a0b16d..91e0ce07bc 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -39,6 +39,7 @@ addtest( ../dftu_base.cpp ../dftu_base_io.cpp ../dftu_base_io_read.cpp + ../dftu_base_io_write.cpp ../dftu_pw_tools.cpp ../yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 56f80a8f04..596da0b7ca 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -86,12 +86,25 @@ class DFTUBaseTest : public testing::Test /// Call init_base with the given Yukawa switch on a fresh d orbital void init_dftu(Plus_U_Base& dftu, const bool yukawa_potential) + { + init_dftu_for_spin(dftu, + yukawa_potential, + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + void init_dftu_for_spin(Plus_U_Base& dftu, + const bool yukawa_potential, + const int nspin, + const int npol, + const DFTU_BASE::OccmatSocLayout soc_layout) { const std::vector l_channel = {2}; const std::vector hubbard_u = {0.0}; dftu.init_base(ucell, - 1, // npol - 2, // nspin + npol, + nspin, l_channel, yukawa_potential, 0.5, // yukawa_lambda @@ -103,7 +116,7 @@ class DFTUBaseTest : public testing::Test 0.0, // uramping 0, // init_occ_mat 0, // mixing_dftu - DFTU_BASE::SOC_LAYOUT_PAULI); + soc_layout); } }; @@ -449,20 +462,17 @@ TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 1}; const std::string out_dir = testing::TempDir(); + const DFTU_BASE::OccmatStepInfo first_step = {0, 1, 1e-6, 0.5}; + const DFTU_BASE::OccmatStepInfo second_step = {1, 2, 1e-6, 0.5}; // fresh run: no occupation matrix has been computed or loaded ASSERT_FALSE(dftu.is_occmat_ready()); DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep -> occ_matg1.txt - 1, // iter + first_step, false, false, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -483,14 +493,9 @@ TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 1, // istep -> occ_matg2.txt - 2, // iter + second_step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -519,6 +524,7 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, false, 1}; const std::string out_dir = testing::TempDir(); + const DFTU_BASE::OccmatStepInfo step = {0, 1, 1e-6, 0.5}; // Remove stale files from previous test binary runs so the // "file not written" assertions are meaningful. @@ -528,14 +534,9 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep (an output ionic step) - 1, // iter + step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -546,12 +547,7 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) DFTU_BASE::write_latest_occmat(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 2, // iter - 1e-6, - 0.5, + step, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -576,18 +572,14 @@ TEST_F(DFTUBaseTest, DftPlusUDisabledWritesNothing) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; const std::string out_dir = testing::TempDir() + "guard_test/"; + const DFTU_BASE::OccmatStepInfo step = {0, 1, 1e-6, 0.5}; DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep (an output ionic step) - 1, // iter + step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -598,12 +590,7 @@ TEST_F(DFTUBaseTest, DftPlusUDisabledWritesNothing) DFTU_BASE::write_latest_occmat(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 2, // iter - 1e-6, - 0.5, + step, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -829,6 +816,239 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatSocPauliRoundtrip) } } +TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI); + for (int spin = 0; spin < 2; ++spin) + { + for (int row = 0; row < 5; ++row) + { + for (int column = 0; column < 5; ++column) + { + const double value = 0.1 * (spin + 1) + 0.01 * row + 0.001 * column; + source.occmat().set(0, 2, spin, row, column, value); + } + } + } + + const std::string filename = testing::TempDir() + "compact_writer_collinear.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + OccupationMatrix restored; + const std::vector l_channel = {2}; + restored.init(ucell, l_channel, 2, 1); + DFTU_BASE::read_occup_m(ucell, + restored, + l_channel, + filename, + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI); + + for (int spin = 0; spin < 2; ++spin) + { + for (int row = 0; row < 5; ++row) + { + for (int column = 0; column < 5; ++column) + { + EXPECT_NEAR(restored.get(0, 2, spin, row, column), + source.occmat().get(0, 2, spin, row, column), + 1e-8); + } + } + } + + // Assert the diagnostic header values that read_occup_m only validates + // for finiteness: per-spin nelec (trace) and collinear mag (up - down). + double expected_nelec[2] = {0.0, 0.0}; + for (int spin = 0; spin < 2; ++spin) + { + for (int m = 0; m < 5; ++m) + { + expected_nelec[spin] += source.occmat().get(0, 2, spin, m, m); + } + } + const double expected_mag = expected_nelec[0] - expected_nelec[1]; + + std::ifstream input(filename.c_str()); + ASSERT_TRUE(input.is_open()); + std::string token; + double parsed_mag = 0.0; + double parsed_nelec[2] = {0.0, 0.0}; + int spin_index = -1; + while (input >> token) + { + if (token == "mag") + { + ASSERT_TRUE(input >> parsed_mag); + } + else if (token == "spin") + { + ASSERT_TRUE(input >> spin_index); + ASSERT_GE(spin_index, 1); + ASSERT_LE(spin_index, 2); + } + else if (token == "nelec") + { + ASSERT_EQ(spin_index == 1 || spin_index == 2, true); + ASSERT_TRUE(input >> parsed_nelec[spin_index - 1]); + } + } + EXPECT_NEAR(parsed_mag, expected_mag, 1e-7); + EXPECT_NEAR(parsed_nelec[0], expected_nelec[0], 1e-7); + EXPECT_NEAR(parsed_nelec[1], expected_nelec[1], 1e-7); +} + +TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 4, 2, DFTU_BASE::SOC_LAYOUT_PAULI); + const int block_size = 25; + for (int block = 0; block < 4; ++block) + { + for (int index = 0; index < block_size; ++index) + { + source.occmat().mat(0, 2, 0).c[block * block_size + index] + = 0.1 * (block + 1) + 0.001 * index; + } + } + + const std::string filename = testing::TempDir() + "compact_writer_soc_pauli.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + OccupationMatrix restored; + const std::vector l_channel = {2}; + restored.init(ucell, l_channel, 4, 2); + DFTU_BASE::read_occup_m(ucell, + restored, + l_channel, + filename, + 4, + 2, + DFTU_BASE::SOC_LAYOUT_PAULI); + + for (int index = 0; index < 4 * block_size; ++index) + { + EXPECT_NEAR(restored.mat(0, 2, 0).c[index], + source.occmat().mat(0, 2, 0).c[index], + 1e-8); + } + + // Assert the SOC diagnostic header: mag x/y/z and per-spin nelec. + // Pauli blocks are charge/spin_x/spin_y/spin_z, so the writer reports + // mag_x = sum(spin_x diag), mag_y = sum(spin_y diag), mag_z = sum(spin_z diag) + // and nelec = 0.5 * sum(charge +/- spin_z) on the diagonal. + const ModuleBase::matrix& occ = source.occmat().mat(0, 2, 0); + double expected_mag[3] = {0.0, 0.0, 0.0}; + double expected_nelec_soc[2] = {0.0, 0.0}; + for (int m = 0; m < 5; ++m) + { + const int diag = m * 5 + m; + expected_mag[0] += occ.c[1 * block_size + diag]; + expected_mag[1] += occ.c[2 * block_size + diag]; + expected_mag[2] += occ.c[3 * block_size + diag]; + expected_nelec_soc[0] += 0.5 * (occ.c[0 * block_size + diag] + occ.c[3 * block_size + diag]); + expected_nelec_soc[1] += 0.5 * (occ.c[0 * block_size + diag] - occ.c[3 * block_size + diag]); + } + // The fixture fills every block entry, including a nonzero spin_y block, + // so mag_y must be nonzero for this assertion to be meaningful. + EXPECT_NE(expected_mag[1], 0.0); + + std::ifstream soc_input(filename.c_str()); + ASSERT_TRUE(soc_input.is_open()); + std::string token; + double parsed_mag_soc[3] = {0.0, 0.0, 0.0}; + double parsed_nelec_soc[2] = {0.0, 0.0}; + int spin_index = -1; + while (soc_input >> token) + { + if (token == "mag") + { + ASSERT_TRUE(soc_input >> parsed_mag_soc[0] >> parsed_mag_soc[1] >> parsed_mag_soc[2]); + } + else if (token == "spin") + { + ASSERT_TRUE(soc_input >> spin_index); + } + else if (token == "nelec") + { + ASSERT_EQ(spin_index == 1 || spin_index == 2, true); + ASSERT_TRUE(soc_input >> parsed_nelec_soc[spin_index - 1]); + } + } + EXPECT_NEAR(parsed_mag_soc[0], expected_mag[0], 1e-7); + EXPECT_NEAR(parsed_mag_soc[1], expected_mag[1], 1e-7); + EXPECT_NEAR(parsed_mag_soc[2], expected_mag[2], 1e-7); + EXPECT_NEAR(parsed_nelec_soc[0], expected_nelec_soc[0], 1e-7); + EXPECT_NEAR(parsed_nelec_soc[1], expected_nelec_soc[1], 1e-7); +} + +TEST_F(DFTUBaseTest, CompactWriterReaderSocRealSpinBasisRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 4, 2, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + const int orbital_size = 5; + for (int row = 0; row < orbital_size; ++row) + { + for (int column = 0; column < orbital_size; ++column) + { + const double up_up = 0.1 + 0.01 * row + 0.001 * column; + const double up_down = 0.2 + 0.01 * row + 0.001 * column; + const double down_down = 0.3 + 0.01 * row + 0.001 * column; + source.occmat().set(0, 2, 0, row, column, up_up); + source.occmat().set(0, 2, 0, row, column + orbital_size, up_down); + source.occmat().set(0, 2, 0, row + orbital_size, column, up_down); + source.occmat().set(0, 2, 0, row + orbital_size, column + orbital_size, down_down); + } + } + + const std::string filename = testing::TempDir() + "compact_writer_soc_real.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + } + + OccupationMatrix restored; + const std::vector l_channel = {2}; + restored.init(ucell, l_channel, 4, 2); + DFTU_BASE::read_occup_m(ucell, + restored, + l_channel, + filename, + 4, + 2, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + + for (int row = 0; row < 2 * orbital_size; ++row) + { + for (int column = 0; column < 2 * orbital_size; ++column) + { + EXPECT_NEAR(restored.get(0, 2, 0, row, column), + source.occmat().get(0, 2, 0, row, column), + 1e-8); + } + } +} + TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) { const std::string fn = testing::TempDir() + "compact_long_comment.txt"; @@ -1297,6 +1517,7 @@ TEST_F(DFTUBaseTest, AppendSnapshotGuardDftuOffBeforeIterValidation) // out_occ_mat=true but dft_plus_u=0: the guard must catch dft_plus_u // and return before the iter < 1 check fires. const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; + const DFTU_BASE::OccmatStepInfo step = {0, 0, 1e-6, 0.5}; // Use a unique subdirectory so files from other tests in the same // binary do not cause false positives. @@ -1307,14 +1528,9 @@ TEST_F(DFTUBaseTest, AppendSnapshotGuardDftuOffBeforeIterValidation) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 0, // iter — invalid, but DFT+U is off + step, false, false, - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI);