From 0eb07961f068a93ea6513ab78f874f459a300b53 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 30 Jun 2026 12:54:17 +0800 Subject: [PATCH 01/10] format dftu_io.cpp --- source/source_lcao/module_dftu/dftu_io.cpp | 272 +++++++++++---------- 1 file changed, 144 insertions(+), 128 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index d44113d1be9..61366127927 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -2,9 +2,10 @@ #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" #include +#include -void Plus_U::output(const UnitCell &ucell) +void Plus_U::output(const UnitCell& ucell) { ModuleBase::TITLE("Plus_U", "output"); @@ -20,10 +21,10 @@ void Plus_U::output(const UnitCell &ucell) if (L >= get_orbital_corr(T) && has_correlated_orbital(T)) { - if (L != get_orbital_corr(T)) - { - continue; - } + if (L != get_orbital_corr(T)) + { + continue; + } if (!Yukawa) { @@ -34,11 +35,11 @@ void Plus_U::output(const UnitCell &ucell) { for (int n = 0; n < N; n++) { - if (n != 0) - { - continue; - } - double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV; + if (n != 0) + { + continue; + } + double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV; GlobalV::ofs_running << "atom_type=" << T << " L=" << L << " chi=" << n << " U=" << this->U_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV " << "J=" << this->J_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV" @@ -50,21 +51,24 @@ void Plus_U::output(const UnitCell &ucell) } GlobalV::ofs_running << "Local occupation matrices" << std::endl; - this->write_occup_m(ucell,GlobalV::ofs_running, true); + this->write_occup_m(ucell, GlobalV::ofs_running, true); GlobalV::ofs_running << "//=======================================================//" << std::endl; - - //Write onsite.dm + + // Write onsite.dm std::ofstream ofdftu; - if(PARAM.inp.out_chg[0]){ - if(GlobalV::MY_RANK == 0){ - ofdftu.open(PARAM.globalv.global_out_dir + "onsite.dm"); - } + if (PARAM.inp.out_chg[0]) + { + if (GlobalV::MY_RANK == 0) + { + ofdftu.open(PARAM.globalv.global_out_dir + "onsite.dm"); + } } - if(!ofdftu){ - std::cout << "Plus_U::write_occup_m. Can't create file onsite.dm!" << std::endl; - exit(0); - } - this->write_occup_m(ucell,ofdftu); + if (!ofdftu) + { + std::cout << "Plus_U::write_occup_m. Can't create file onsite.dm!" << std::endl; + exit(0); + } + this->write_occup_m(ucell, ofdftu); ofdftu.close(); return; @@ -74,23 +78,23 @@ void Plus_U::output(const UnitCell &ucell) std::vector CalculateEigenvalues(std::vector>& A, int n); void Plus_U::write_occup_m(const UnitCell& ucell, - std::ofstream &ofs, - bool diag) + std::ofstream& ofs, + bool diag) { ModuleBase::TITLE("Plus_U", "write_occup_m"); - if(GlobalV::MY_RANK != 0) - { - return; - } + if (GlobalV::MY_RANK != 0) + { + return; + } for (int T = 0; T < ucell.ntype; T++) { - if (!has_correlated_orbital(T)) - { - continue; - } - const int NL = ucell.atoms[T].nwl + 1; + if (!has_correlated_orbital(T)) + { + continue; + } + const int NL = ucell.atoms[T].nwl + 1; const int LC = get_orbital_corr(T); for (int I = 0; I < ucell.atoms[T].na; I++) @@ -101,10 +105,10 @@ void Plus_U::write_occup_m(const UnitCell& ucell, for (int l = 0; l < NL; l++) { - if (l != get_orbital_corr(T)) - { - continue; - } + if (l != get_orbital_corr(T)) + { + continue; + } const int N = ucell.atoms[T].l_nchi[l]; ofs << "L" @@ -113,10 +117,10 @@ void Plus_U::write_occup_m(const UnitCell& ucell, for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) - { - continue; - } + if (n != 0) + { + continue; + } ofs << "zeta" << " " << n << std::endl; @@ -126,7 +130,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, double sum0[2]; for (int is = 0; is < 2; is++) { - if(diag)// diagonalization for local occupation matrix and print the eigenvalues + if (diag) // diagonalization for local occupation matrix and print the eigenvalues { std::vector> A(2 * l + 1, std::vector(2 * l + 1)); for (int m0 = 0; m0 < 2 * l + 1; m0++) @@ -138,7 +142,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); sum0[is] = 0.0; - ofs<< "eigenvalues" + ofs << "eigenvalues" << " " << is << std::endl; for (int i = 0; i < 2 * l + 1; i++) { @@ -161,20 +165,20 @@ void Plus_U::write_occup_m(const UnitCell& ucell, ofs << std::endl; } } - if(diag) + if (diag) { - ofs << std::setw(12) << std::setprecision(8) - << std::fixed<< "atomic mag: "<> A(2 * l + 1, std::vector(2 * l + 1)); int index = 0; - for(int is=0;is<4;is++) + for (int is = 0; is < 4; is++) { for (int m0 = 0; m0 < 2 * l + 1; m0++) { @@ -186,7 +190,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); sum0[is] = 0.0; - ofs<< "eigenvalues" + ofs << "eigenvalues" << " " << is << std::endl; for (int i = 0; i < 2 * l + 1; i++) { @@ -197,49 +201,49 @@ void Plus_U::write_occup_m(const UnitCell& ucell, ofs << std::setw(12) << std::setprecision(8) << std::fixed << sum0[is] << std::endl; } - ofs << std::setw(12) << std::setprecision(8) - << std::fixed<< "atomic mag: "<> word; - if (ifdftu.eof()) - { - break; - } + if (ifdftu.eof()) + { + break; + } if (strcmp("atoms", word) == 0) { @@ -294,10 +298,10 @@ void Plus_U::read_occup_m(const UnitCell& ucell, for (int l = 0; l < NL; l++) { - if (l != get_orbital_corr(T)) - { - continue; - } + if (l != get_orbital_corr(T)) + { + continue; + } ifdftu >> word; @@ -310,10 +314,10 @@ void Plus_U::read_occup_m(const UnitCell& ucell, for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) - { - continue; - } + if (n != 0) + { + continue; + } ifdftu >> word; if (strcmp("zeta", word) == 0) @@ -375,7 +379,8 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" + << std::endl; exit(0); } } @@ -410,10 +415,10 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell) for (int T = 0; T < ucell.ntype; T++) { - if (!has_correlated_orbital(T)) - { - continue; - } + if (!has_correlated_orbital(T)) + { + continue; + } for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -422,18 +427,18 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell) for (int l = 0; l <= ucell.atoms[T].nwl; l++) { - if (l != get_orbital_corr(T)) - { - continue; - } + if (l != get_orbital_corr(T)) + { + continue; + } for (int n = 0; n < ucell.atoms[T].l_nchi[l]; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) - { - continue; - } + if (n != 0) + { + continue; + } if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { @@ -482,15 +487,18 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell) return; } -inline void JacobiRotate(std::vector>& A, int p, int q, int n) +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]); + 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) { + if (r >= 0) + { t = 1.0 / (r + sqrt(1.0 + r * r)); - } else { + } + else + { t = -1.0 / (-r + sqrt(1.0 + r * r)); } double c = 1.0 / sqrt(1.0 + t * t); @@ -500,8 +508,10 @@ inline void JacobiRotate(std::vector>& A, int p, int q, int 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) { + 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; @@ -511,22 +521,28 @@ inline void JacobiRotate(std::vector>& A, int p, int q, int } } -inline std::vector CalculateEigenvalues(std::vector>& A, int n) +inline std::vector CalculateEigenvalues(std::vector>& A, int n) { std::vector eigenvalues(n); - while (true) { + 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])) { + 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++) { + if (std::abs(A[p][q]) < 1e-10) + { + for (int i = 0; i < n; i++) + { eigenvalues[i] = A[i][i]; } break; From 0e95a1096d4649982d886a23029e5f494ce41128 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 30 Jun 2026 13:39:00 +0800 Subject: [PATCH 02/10] update output formats of DFT+U --- source/source_lcao/module_dftu/dftu_io.cpp | 51 +++++++++++----------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index 61366127927..454dd8a22b3 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -9,7 +9,9 @@ void Plus_U::output(const UnitCell& ucell) { ModuleBase::TITLE("Plus_U", "output"); - GlobalV::ofs_running << "//=========================L(S)DA+U===========================//" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + GlobalV::ofs_running << " | #DFT+U INFORMATION# |" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; for (int T = 0; T < ucell.ntype; T++) { @@ -28,8 +30,8 @@ void Plus_U::output(const UnitCell& ucell) if (!Yukawa) { - GlobalV::ofs_running << "atom_type=" << T << " L=" << L << " chi=" << 0 - << " U=" << this->U[T] * ModuleBase::Ry_to_eV << "eV" << std::endl; + GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << 0 + << " U=" << this->U[T] * ModuleBase::Ry_to_eV << " eV" << std::endl; } else { @@ -40,9 +42,9 @@ void Plus_U::output(const UnitCell& ucell) continue; } double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV; - GlobalV::ofs_running << "atom_type=" << T << " L=" << L << " chi=" << n - << " U=" << this->U_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV " - << "J=" << this->J_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV" + GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << n + << " U=" << this->U_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << " eV" + << " J=" << this->J_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << " eV" << std::endl; } } @@ -50,9 +52,8 @@ void Plus_U::output(const UnitCell& ucell) } } - GlobalV::ofs_running << "Local occupation matrices" << std::endl; + GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; this->write_occup_m(ucell, GlobalV::ofs_running, true); - GlobalV::ofs_running << "//=======================================================//" << std::endl; // Write onsite.dm std::ofstream ofdftu; @@ -65,12 +66,16 @@ void Plus_U::output(const UnitCell& ucell) } if (!ofdftu) { - std::cout << "Plus_U::write_occup_m. Can't create file onsite.dm!" << std::endl; + std::cout << " Plus_U::write_occup_m. Can't create file onsite.dm" << std::endl; exit(0); } this->write_occup_m(ucell, ofdftu); ofdftu.close(); + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + GlobalV::ofs_running << " | # END DFT+U INFO |" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl << std::endl; + return; } @@ -100,8 +105,6 @@ void Plus_U::write_occup_m(const UnitCell& ucell, for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); - ofs << "atoms" - << " " << iat << std::endl; for (int l = 0; l < NL; l++) { @@ -111,8 +114,6 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } const int N = ucell.atoms[T].l_nchi[l]; - ofs << "L" - << " " << l << std::endl; for (int n = 0; n < N; n++) { @@ -122,8 +123,9 @@ void Plus_U::write_occup_m(const UnitCell& ucell, continue; } - ofs << "zeta" - << " " << n << std::endl; + ofs << "\n Atom=" << iat+1; + ofs << " L=" << l; + ofs << " ORBITAL=" << n << std::endl; if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { @@ -142,24 +144,22 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); sum0[is] = 0.0; - ofs << "eigenvalues" - << " " << is << std::endl; + ofs << " Eigenvalues for spin=" << is+1 << std::endl; + ofs << std::setprecision(8) << std::fixed; for (int i = 0; i < 2 * l + 1; i++) { - ofs << std::setw(12) << std::setprecision(8) << std::fixed - << eigenvalues[i]; + ofs << std::setw(12) << eigenvalues[i]; sum0[is] += eigenvalues[i]; } - ofs << std::setw(12) << std::setprecision(8) << std::fixed - << sum0[is] << std::endl; + ofs << std::endl; + ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } - ofs << "spin" - << " " << is << std::endl; + ofs << " spin=" << is+1 << std::endl; for (int m0 = 0; m0 < 2 * l + 1; m0++) { for (int m1 = 0; m1 < 2 * l + 1; m1++) { - ofs << std::setw(12) << std::setprecision(8) << std::fixed + ofs << std::setw(12) << locale[iat][l][n][is](m0, m1); } ofs << std::endl; @@ -168,7 +168,8 @@ void Plus_U::write_occup_m(const UnitCell& ucell, if (diag) { ofs << std::setw(12) << std::setprecision(8) - << std::fixed << "atomic mag: " << iat << " " << sum0[0] - sum0[1] << std::endl; + << std::fixed << " Magnetism for atom " << iat+1 << ": " << sum0[0] - sum0[1] + << std::endl; } } else if (PARAM.inp.nspin == 4) // SOC From 87b9a0962474e0f451425b917982d848755f5c73 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 30 Jun 2026 17:33:11 +0800 Subject: [PATCH 03/10] update dftu, remove PARAM --- source/source_esolver/esolver_ks_lcao.cpp | 4 +- source/source_lcao/dftu_lcao.cpp | 39 ++++++++---- source/source_lcao/dftu_lcao.h | 14 ++-- source/source_lcao/module_dftu/dftu.cpp | 8 +-- source/source_lcao/module_dftu/dftu.h | 21 ++++-- source/source_lcao/module_dftu/dftu_io.cpp | 74 ++++++++++++---------- source/source_pw/module_pwdft/dftu_pw.cpp | 2 +- 7 files changed, 100 insertions(+), 62 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 232f75ab543..25e902cc4f6 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -370,7 +370,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const } #endif - init_dftu_lcao(istep, iter, PARAM.inp, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx); + init_dftu_lcao(istep, iter, PARAM.inp.dft_plus_u, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx); #ifdef __MLALGO // the density matrixes of DeePKS have been updated in each iter @@ -481,7 +481,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& const std::vector>& dm_vec = this->dmat.dm->get_DMK_vector(); // 1) calculate the local occupation number matrix and energy correction in DFT+U - finish_dftu_lcao(iter, conv_esolver, PARAM.inp, &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao); + finish_dftu_lcao(iter, conv_esolver, PARAM.inp.dft_plus_u, PARAM.inp.out_chg[0], &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol); // 2) for deepks, calculate delta_e, output labels during electronic steps this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp); diff --git a/source/source_lcao/dftu_lcao.cpp b/source/source_lcao/dftu_lcao.cpp index d8b8421d6e7..de9c6db943e 100644 --- a/source/source_lcao/dftu_lcao.cpp +++ b/source/source_lcao/dftu_lcao.cpp @@ -9,14 +9,14 @@ namespace ModuleESolver template void init_dftu_lcao(const int istep, const int iter, - const Input_para& inp, + int dft_plus_u, void* dftu, void* dm, const UnitCell& ucell, double** rho, const int nrxx) { - if (!inp.dft_plus_u) + if (!dft_plus_u) { return; } @@ -36,15 +36,19 @@ void init_dftu_lcao(const int istep, template void finish_dftu_lcao(const int iter, const bool conv_esolver, - const Input_para& inp, + int dft_plus_u, + bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, const K_Vectors& kv, const double mixing_beta, - void* hamilt_lcao) + void* hamilt_lcao, + const std::string& global_out_dir, + int nspin, + int npol) { - if (!inp.dft_plus_u) + if (!dft_plus_u) { return; } @@ -54,7 +58,7 @@ void finish_dftu_lcao(const int iter, /// old DFT+U method calculates energy correction in esolver, /// new DFT+U method calculates energy in Hamiltonian - if (inp.dft_plus_u == 2) + if (dft_plus_u == 2) { if (dftu_ptr->omc != 2) { @@ -63,7 +67,7 @@ void finish_dftu_lcao(const int iter, } dftu_ptr->cal_energy_correction(ucell, iter); } - dftu_ptr->output(ucell); + dftu_ptr->output(ucell, out_chg, global_out_dir, nspin, npol); /// use the converged occupation matrix for next MD/Relax SCF calculation if (conv_esolver) @@ -75,15 +79,16 @@ void finish_dftu_lcao(const int iter, /// Template instantiation template void init_dftu_lcao(const int istep, const int iter, - const Input_para& inp, + int dft_plus_u, void* dftu, void* dm, const UnitCell& ucell, double** rho, const int nrxx); + template void init_dftu_lcao>(const int istep, const int iter, - const Input_para& inp, + int dft_plus_u, void* dftu, void* dm, const UnitCell& ucell, @@ -92,21 +97,29 @@ template void init_dftu_lcao>(const int istep, template void finish_dftu_lcao(const int iter, const bool conv_esolver, - const Input_para& inp, + int dft_plus_u, + bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, const K_Vectors& kv, const double mixing_beta, - void* hamilt_lcao); + void* hamilt_lcao, + const std::string& global_out_dir, + int nspin, + int npol); template void finish_dftu_lcao>(const int iter, const bool conv_esolver, - const Input_para& inp, + int dft_plus_u, + bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>>& dm_vec, const K_Vectors& kv, const double mixing_beta, - void* hamilt_lcao); + void* hamilt_lcao, + const std::string& global_out_dir, + int nspin, + int npol); } // namespace ModuleESolver diff --git a/source/source_lcao/dftu_lcao.h b/source/source_lcao/dftu_lcao.h index 5138b66256f..73e5c0448f0 100644 --- a/source/source_lcao/dftu_lcao.h +++ b/source/source_lcao/dftu_lcao.h @@ -3,7 +3,6 @@ #include "source_cell/unitcell.h" #include "source_cell/klist.h" -#include "source_io/module_parameter/input_parameter.h" namespace ModuleESolver { @@ -26,7 +25,7 @@ namespace ModuleESolver template void init_dftu_lcao(const int istep, const int iter, - const Input_para& inp, + int dft_plus_u, void* dftu, void* dm, const UnitCell& ucell, @@ -41,7 +40,8 @@ void init_dftu_lcao(const int istep, * * @param iter Current SCF iteration * @param conv_esolver Whether ESolver has converged - * @param inp Input parameters + * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new) + * @param out_chg Whether to output onsite.dm * @param dftu DFT+U object * @param ucell Unit cell * @param dm_vec Density matrix vector @@ -52,13 +52,17 @@ void init_dftu_lcao(const int istep, template void finish_dftu_lcao(const int iter, const bool conv_esolver, - const Input_para& inp, + int dft_plus_u, + bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, const K_Vectors& kv, const double mixing_beta, - void* hamilt_lcao); + void* hamilt_lcao, + const std::string& global_out_dir, + int nspin, + int npol); } // namespace ModuleESolver diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp index ff06003c43d..f3f967f86f2 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -238,9 +238,9 @@ void Plus_U::init(UnitCell& cell, // unitcell class { std::stringstream sst; sst << "initial_onsite.dm"; - this->read_occup_m(cell,sst.str()); + this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); #ifdef __MPI - this->local_occup_bcast(cell); + this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); #endif mark_locale_initialized(); @@ -252,9 +252,9 @@ void Plus_U::init(UnitCell& cell, // unitcell class { std::stringstream sst; sst << PARAM.globalv.global_readin_dir << "onsite.dm"; - this->read_occup_m(cell,sst.str()); + this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); #ifdef __MPI - this->local_occup_bcast(cell); + this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); #endif mark_locale_initialized(); } diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h index ff494d4fa3f..5b3ca2bcc08 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -369,15 +369,26 @@ class Plus_U // For reading/writing/broadcasting/copying relevant data structures //============================================================= public: - void output(const UnitCell& ucell); + void output(const UnitCell& ucell, + bool out_chg, + const std::string& global_out_dir, + int nspin, + int npol); private: void write_occup_m(const UnitCell& ucell, - std::ofstream& ofs, - bool diag=false); + std::ofstream& ofs, + bool diag, + int nspin, + int npol); void read_occup_m(const UnitCell& ucell, - const std::string& fn); - void local_occup_bcast(const UnitCell& ucell); + const std::string& fn, + const std::string& init_chg, + int nspin, + int npol); + void local_occup_bcast(const UnitCell& ucell, + int nspin, + int npol); //============================================================= // In dftu_yukawa.cpp diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index 454dd8a22b3..82928f6c880 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -5,7 +5,11 @@ #include -void Plus_U::output(const UnitCell& ucell) +void Plus_U::output(const UnitCell& ucell, + bool out_chg, + const std::string& global_out_dir, + int nspin, + int npol) { ModuleBase::TITLE("Plus_U", "output"); @@ -53,15 +57,15 @@ void Plus_U::output(const UnitCell& ucell) } GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; - this->write_occup_m(ucell, GlobalV::ofs_running, true); + this->write_occup_m(ucell, GlobalV::ofs_running, true, nspin, npol); // Write onsite.dm std::ofstream ofdftu; - if (PARAM.inp.out_chg[0]) + if (out_chg) { if (GlobalV::MY_RANK == 0) { - ofdftu.open(PARAM.globalv.global_out_dir + "onsite.dm"); + ofdftu.open(global_out_dir + "onsite.dm"); } } if (!ofdftu) @@ -69,7 +73,7 @@ void Plus_U::output(const UnitCell& ucell) std::cout << " Plus_U::write_occup_m. Can't create file onsite.dm" << std::endl; exit(0); } - this->write_occup_m(ucell, ofdftu); + this->write_occup_m(ucell, ofdftu, false, nspin, npol); ofdftu.close(); GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; @@ -84,7 +88,9 @@ std::vector CalculateEigenvalues(std::vector>& A, in void Plus_U::write_occup_m(const UnitCell& ucell, std::ofstream& ofs, - bool diag) + bool diag, + int nspin, + int npol) { ModuleBase::TITLE("Plus_U", "write_occup_m"); @@ -127,7 +133,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, ofs << " L=" << l; ofs << " ORBITAL=" << n << std::endl; - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (nspin == 1 || nspin == 2) { double sum0[2]; for (int is = 0; is < 2; is++) @@ -145,13 +151,13 @@ void Plus_U::write_occup_m(const UnitCell& ucell, std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); sum0[is] = 0.0; ofs << " Eigenvalues for spin=" << is+1 << std::endl; - ofs << std::setprecision(8) << std::fixed; + ofs << std::setprecision(8) << std::fixed; for (int i = 0; i < 2 * l + 1; i++) { ofs << std::setw(12) << eigenvalues[i]; sum0[is] += eigenvalues[i]; } - ofs << std::endl; + ofs << std::endl; ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } ofs << " spin=" << is+1 << std::endl; @@ -159,7 +165,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, { for (int m1 = 0; m1 < 2 * l + 1; m1++) { - ofs << std::setw(12) + ofs << std::setw(12) << locale[iat][l][n][is](m0, m1); } ofs << std::endl; @@ -169,10 +175,10 @@ void Plus_U::write_occup_m(const UnitCell& ucell, { ofs << std::setw(12) << std::setprecision(8) << std::fixed << " Magnetism for atom " << iat+1 << ": " << sum0[0] - sum0[1] - << std::endl; + << std::endl; } } - else if (PARAM.inp.nspin == 4) // SOC + else if (nspin == 4) // SOC { if (diag) // diagonalization for local occupation matrix and print the eigenvalues { // output the eigenvalues for rho , mag_x, mag_y, mag_z @@ -191,32 +197,31 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); sum0[is] = 0.0; - ofs << "eigenvalues" - << " " << is << std::endl; + ofs << " Eigenvalues for is=" << is << std::endl; + ofs << std::setprecision(8) << std::fixed; for (int i = 0; i < 2 * l + 1; i++) { - ofs << std::setw(12) << std::setprecision(8) << std::fixed - << eigenvalues[i]; + ofs << std::setw(12) << eigenvalues[i]; sum0[is] += eigenvalues[i]; } - ofs << std::setw(12) << std::setprecision(8) << std::fixed - << sum0[is] << std::endl; + ofs << std::endl; + ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } ofs << std::setw(12) << std::setprecision(8) - << std::fixed << "atomic mag: " << iat << " " + << std::fixed << " Magnetism for atom " << iat + 1 << ": " << sum0[1] << " " << sum0[2] << " " << sum0[3] << std::endl; } else { for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < npol; ipol0++) { const int m0_all = m0 + (2 * l + 1) * ipol0; for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ipol1++) { int m1_all = m1 + (2 * l + 1) * ipol1; ofs << std::setw(12) << std::setprecision(8) << std::fixed @@ -237,7 +242,10 @@ void Plus_U::write_occup_m(const UnitCell& ucell, } void Plus_U::read_occup_m(const UnitCell& ucell, - const std::string& fn) + const std::string& fn, + const std::string& init_chg, + int nspin, + int npol) { ModuleBase::TITLE("Plus_U", "read_occup_m"); @@ -258,7 +266,7 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } else { - if (PARAM.inp.init_chg == "file") + if (init_chg == "file") { std::cout << "Plus_U::read_occup_m. Can not find the file onsite.dm . Please do scf calculation first" << std::endl; @@ -326,7 +334,7 @@ void Plus_U::read_occup_m(const UnitCell& ucell, ifdftu >> zeta; ifdftu.ignore(150, '\n'); - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (nspin == 1 || nspin == 2) { for (int is = 0; is < 2; is++) { @@ -355,18 +363,18 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } } } - else if (PARAM.inp.nspin == 4) // SOC + else if (nspin == 4) // SOC { double value = 0.0; for (int m0 = 0; m0 < 2 * L + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < npol; ipol0++) { const int m0_all = m0 + (2 * L + 1) * ipol0; for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ipol1++) { int m1_all = m1 + (2 * L + 1) * ipol1; ifdftu >> value; @@ -410,7 +418,9 @@ void Plus_U::read_occup_m(const UnitCell& ucell, return; } -void Plus_U::local_occup_bcast(const UnitCell& ucell) +void Plus_U::local_occup_bcast(const UnitCell& ucell, + int nspin, + int npol) { ModuleBase::TITLE("Plus_U", "local_occup_bcast"); @@ -441,7 +451,7 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell) continue; } - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (nspin == 1 || nspin == 2) { for (int spin = 0; spin < 2; spin++) { @@ -456,17 +466,17 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell) } } } - else if (PARAM.inp.nspin == 4) // SOC + else if (nspin == 4) // SOC { for (int m0 = 0; m0 < 2 * L + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < npol; ipol0++) { const int m0_all = m0 + (2 * L + 1) * ipol0; for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ipol1++) { int m1_all = m1 + (2 * L + 1) * ipol1; #ifdef __MPI diff --git a/source/source_pw/module_pwdft/dftu_pw.cpp b/source/source_pw/module_pwdft/dftu_pw.cpp index 97cb8d5ccce..667612e23bd 100644 --- a/source/source_pw/module_pwdft/dftu_pw.cpp +++ b/source/source_pw/module_pwdft/dftu_pw.cpp @@ -28,7 +28,7 @@ void iter_init_dftu_pw(const int iter, { dftu.cal_occ_pw(iter, psi, wg, ucell, p_chgmix, isk); } - dftu.output(ucell); + dftu.output(ucell, PARAM.inp.out_chg[0], PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol); } } From a5a86f3e9c9cd333d07ec6c72e56c9dfaca4c778 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 11:47:32 +0800 Subject: [PATCH 04/10] fix read_occup_m function --- source/source_lcao/module_dftu/dftu_io.cpp | 118 +++++++++------------ 1 file changed, 50 insertions(+), 68 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index 82928f6c880..89458a0c344 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -278,7 +278,7 @@ void Plus_U::read_occup_m(const UnitCell& ucell, ifdftu.clear(); ifdftu.seekg(0); - char word[10]; + char word[20]; int T = 0; int iat = 0; @@ -296,9 +296,26 @@ void Plus_U::read_occup_m(const UnitCell& ucell, break; } - if (strcmp("atoms", word) == 0) + if (strcmp("Atom=", word) == 0) { ifdftu >> iat; + iat -= 1; + ifdftu >> word; + + if (strcmp("L=", word) != 0) + { + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; + exit(0); + } + ifdftu >> L; + ifdftu >> word; + + if (strcmp("ORBITAL=", word) != 0) + { + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; + exit(0); + } + ifdftu >> zeta; ifdftu.ignore(150, '\n'); T = ucell.iat2it[iat]; @@ -312,78 +329,26 @@ void Plus_U::read_occup_m(const UnitCell& ucell, continue; } - ifdftu >> word; - - if (strcmp("L", word) == 0) + if (nspin == 1 || nspin == 2) { - ifdftu >> L; - ifdftu.ignore(150, '\n'); - - const int N = ucell.atoms[T].l_nchi[L]; - for (int n = 0; n < N; n++) + for (int is = 0; is < 2; is++) { - // if(!Yukawa && n!=0) continue; - if (n != 0) - { - continue; - } - ifdftu >> word; - if (strcmp("zeta", word) == 0) + if (strcmp("spin=", word) == 0) { - ifdftu >> zeta; + ifdftu >> spin; + spin -= 1; ifdftu.ignore(150, '\n'); - if (nspin == 1 || nspin == 2) - { - for (int is = 0; is < 2; is++) - { - ifdftu >> word; - if (strcmp("spin", word) == 0) - { - ifdftu >> spin; - ifdftu.ignore(150, '\n'); - - double value = 0.0; - for (int m0 = 0; m0 < 2 * L + 1; m0++) - { - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - ifdftu >> value; - locale[iat][L][zeta][spin](m0, m1) = value; - } - ifdftu.ignore(150, '\n'); - } - } - else - { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" - << std::endl; - exit(0); - } - } - } - else if (nspin == 4) // SOC + double value = 0.0; + for (int m0 = 0; m0 < 2 * L + 1; m0++) { - double value = 0.0; - for (int m0 = 0; m0 < 2 * L + 1; m0++) + for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol0 = 0; ipol0 < npol; ipol0++) - { - const int m0_all = m0 + (2 * L + 1) * ipol0; - - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - for (int ipol1 = 0; ipol1 < npol; ipol1++) - { - int m1_all = m1 + (2 * L + 1) * ipol1; - ifdftu >> value; - locale[iat][L][zeta][0](m0_all, m1_all) = value; - } - } - ifdftu.ignore(150, '\n'); - } + ifdftu >> value; + locale[iat][L][zeta][spin](m0, m1) = value; } + ifdftu.ignore(150, '\n'); } } else @@ -394,10 +359,27 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } } } - else + else if (nspin == 4) // SOC { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; - exit(0); + double value = 0.0; + for (int m0 = 0; m0 < 2 * L + 1; m0++) + { + for (int ipol0 = 0; ipol0 < npol; ipol0++) + { + const int m0_all = m0 + (2 * L + 1) * ipol0; + + for (int m1 = 0; m1 < 2 * L + 1; m1++) + { + for (int ipol1 = 0; ipol1 < npol; ipol1++) + { + int m1_all = m1 + (2 * L + 1) * ipol1; + ifdftu >> value; + locale[iat][L][zeta][0](m0_all, m1_all) = value; + } + } + ifdftu.ignore(150, '\n'); + } + } } } } From 1f031bf0ef164322e7a8af09f833942db00ccf54 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 11:54:56 +0800 Subject: [PATCH 05/10] DFT+U I/O: Fix multiple issues with onsite.dm format, precision, and MPI This commit addresses several issues in the DFT+U I/O code related to the onsite.dm file handling: 1. Format compatibility fix: Updated read_occup_m() to parse the new output format with labels Atom=, L=, ORBITAL=, and spin= instead of the old tokens atoms, L, zeta, and spin. The writer was already using the new format but the reader was not updated, causing new onsite.dm files to be unreadable. 2. Off-by-one fix: The writer outputs 1-based indices (iat+1, is+1) for human readability, but read_occup_m() was using these values directly as array indices. Added iat -= 1 and spin -= 1 to convert back to 0-based indices. 3. Precision fix: Added std::setprecision(8) << std::fixed for the collinear (nspin=1,2) diag=false output path in write_occup_m(). The eigenvalues path and SOC path already had this, but the matrix values path was missing it. This ensures stable numeric precision for restart data. 4. out_chg logic fix: The out_chg parameter is supposed to control whether onsite.dm is written, but the implementation was still checking if(!ofdftu) even when out_chg == false. Moved the file-open check and write operations inside the out_chg && MY_RANK == 0 block. 5. MPI fix: Only rank 0 opens the onsite.dm file, but all ranks were executing the if(!ofdftu) check and write_occup_m() call. Non-root ranks would see an unopened stream and potentially fail. Now all file operations are rank-0-only. 6. Header fix: Added missing #include and #include to dftu_lcao.h. The header was using std::string and std::vector in function signatures but relying on transitive includes. 7. Documentation fix: Updated stale comments in dftu_lcao.h: - Changed @param inp to @param dft_plus_u in init_dftu_lcao() - Added documentation for global_out_dir, nspin, and npol parameters in finish_dftu_lcao() 8. Error handling: Replaced all exit(0) calls with ModuleBase::WARNING_QUIT() for consistent error handling across the codebase. Files modified: - source/source_lcao/module_dftu/dftu_io.cpp - source/source_lcao/dftu_lcao.h - source/source_lcao/dftu_lcao.cpp --- source/source_lcao/dftu_lcao.cpp | 1 + source/source_lcao/dftu_lcao.h | 7 +++- source/source_lcao/module_dftu/dftu_io.cpp | 41 ++++++++-------------- 3 files changed, 22 insertions(+), 27 deletions(-) diff --git a/source/source_lcao/dftu_lcao.cpp b/source/source_lcao/dftu_lcao.cpp index de9c6db943e..5abde6a584b 100644 --- a/source/source_lcao/dftu_lcao.cpp +++ b/source/source_lcao/dftu_lcao.cpp @@ -108,6 +108,7 @@ template void finish_dftu_lcao(const int iter, const std::string& global_out_dir, int nspin, int npol); + template void finish_dftu_lcao>(const int iter, const bool conv_esolver, int dft_plus_u, diff --git a/source/source_lcao/dftu_lcao.h b/source/source_lcao/dftu_lcao.h index 73e5c0448f0..245b6d3f415 100644 --- a/source/source_lcao/dftu_lcao.h +++ b/source/source_lcao/dftu_lcao.h @@ -1,6 +1,8 @@ #ifndef DFTU_LCAO_H #define DFTU_LCAO_H +#include +#include #include "source_cell/unitcell.h" #include "source_cell/klist.h" @@ -15,7 +17,7 @@ namespace ModuleESolver * * @param istep Current ionic step * @param iter Current SCF iteration - * @param inp Input parameters + * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new) * @param dftu DFT+U object * @param dm Density matrix * @param ucell Unit cell @@ -48,6 +50,9 @@ void init_dftu_lcao(const int istep, * @param kv K-vectors * @param mixing_beta Mixing beta parameter * @param hamilt_lcao Hamiltonian LCAO object + * @param global_out_dir Output directory for onsite.dm + * @param nspin Number of spin channels (1, 2, or 4) + * @param npol Number of polarizations */ template void finish_dftu_lcao(const int iter, diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index 89458a0c344..1f04684c103 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -60,21 +60,17 @@ void Plus_U::output(const UnitCell& ucell, this->write_occup_m(ucell, GlobalV::ofs_running, true, nspin, npol); // Write onsite.dm - std::ofstream ofdftu; - if (out_chg) + if (out_chg && GlobalV::MY_RANK == 0) { - if (GlobalV::MY_RANK == 0) + std::ofstream ofdftu; + ofdftu.open(global_out_dir + "onsite.dm"); + if (!ofdftu) { - ofdftu.open(global_out_dir + "onsite.dm"); + ModuleBase::WARNING_QUIT("Plus_U::output", "Can't create file onsite.dm"); } + this->write_occup_m(ucell, ofdftu, false, nspin, npol); + ofdftu.close(); } - if (!ofdftu) - { - std::cout << " Plus_U::write_occup_m. Can't create file onsite.dm" << std::endl; - exit(0); - } - this->write_occup_m(ucell, ofdftu, false, nspin, npol); - ofdftu.close(); GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; GlobalV::ofs_running << " | # END DFT+U INFO |" << std::endl; @@ -161,6 +157,7 @@ void Plus_U::write_occup_m(const UnitCell& ucell, ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } ofs << " spin=" << is+1 << std::endl; + ofs << std::setprecision(8) << std::fixed; for (int m0 = 0; m0 < 2 * l + 1; m0++) { for (int m1 = 0; m1 < 2 * l + 1; m1++) @@ -260,19 +257,16 @@ void Plus_U::read_occup_m(const UnitCell& ucell, { if (omc > 0) { - std::cout - << "Plus_U::read_occup_m. Can not find the file initial_onsite.dm . Please check your initial_onsite.dm" - << std::endl; + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file initial_onsite.dm. Please check your initial_onsite.dm"); } else { if (init_chg == "file") { - std::cout << "Plus_U::read_occup_m. Can not find the file onsite.dm . Please do scf calculation first" - << std::endl; + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file onsite.dm. Please do scf calculation first"); } } - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not open onsite.dm file"); } ifdftu.clear(); @@ -304,16 +298,14 @@ void Plus_U::read_occup_m(const UnitCell& ucell, if (strcmp("L=", word) != 0) { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); } ifdftu >> L; ifdftu >> word; if (strcmp("ORBITAL=", word) != 0) { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); } ifdftu >> zeta; ifdftu.ignore(150, '\n'); @@ -353,9 +345,7 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" - << std::endl; - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); } } } @@ -385,8 +375,7 @@ void Plus_U::read_occup_m(const UnitCell& ucell, } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); } ifdftu.rdstate(); From b5572e43085843a023197d3e42acdb71e90a0fc1 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 12:13:19 +0800 Subject: [PATCH 06/10] DFT+U I/O: Rename onsite.dm to dm_onsite.txt and adjust file paths MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This commit renames the DFT+U occupation matrix files and adjusts their read/write paths: 1. File name changes: - onsite.dm → dm_onsite.txt (output from DFT+U calculations) - initial_onsite.dm → dm_onsite_ini.txt (user-provided initial occupation matrix) 2. Path adjustments: - dm_onsite.txt: Both written and read from global_out_dir (OUT.prefix) - Previously read from global_readin_dir - This ensures the file is always in the output directory - dm_onsite_ini.txt: Read from global_readin_dir (set via read_file_dir parameter) - Previously read from global_out_dir - This allows users to place initial files in any directory 3. Code changes: - dftu_io.cpp: Updated output filename and error messages - dftu.cpp: Updated read paths for both files - dftu_lcao.h and module_operator_lcao/dftu_lcao.cpp: Updated comments 4. Documentation updates: - read_input_item_exx_dftu.cpp: Updated omc parameter description - parameters.yaml: Updated omc parameter description - input-main.md: Updated omc parameter description - band.md and dos.md: Updated file references - tests/17_DS_DFTU/README.md: Updated test documentation - tests/17_DS_DFTU/run_scf_nscf.sh: Updated copy logic for dm_onsite.txt 5. Backward compatibility: - The new format with labels Atom=, L=, ORBITAL=, spin= was already in the writer - This commit updates the reader to parse the new format - 1-based indices in output are converted to 0-based when reading Files modified: - source/source_lcao/module_dftu/dftu_io.cpp - source/source_lcao/module_dftu/dftu.cpp - source/source_lcao/dftu_lcao.h - source/source_lcao/module_operator_lcao/dftu_lcao.cpp - source/source_io/module_parameter/read_input_item_exx_dftu.cpp - docs/parameters.yaml - docs/advanced/input_files/input-main.md - docs/advanced/elec_properties/band.md - docs/advanced/elec_properties/dos.md - tests/17_DS_DFTU/README.md - tests/17_DS_DFTU/run_scf_nscf.sh --- docs/advanced/elec_properties/band.md | 2 +- docs/advanced/elec_properties/dos.md | 2 +- docs/advanced/input_files/input-main.md | 6 +++--- docs/parameters.yaml | 6 +++--- .../read_input_item_exx_dftu.cpp | 6 +++--- source/source_lcao/dftu_lcao.h | 4 ++-- source/source_lcao/module_dftu/dftu.cpp | 4 ++-- source/source_lcao/module_dftu/dftu_io.cpp | 12 +++++------ .../module_operator_lcao/dftu_lcao.cpp | 6 +++--- tests/17_DS_DFTU/README.md | 4 ++-- tests/17_DS_DFTU/run_scf_nscf.sh | 20 +++++++++++-------- 11 files changed, 38 insertions(+), 34 deletions(-) diff --git a/docs/advanced/elec_properties/band.md b/docs/advanced/elec_properties/band.md index 44fb88e1c51..d6b66fa2c2a 100644 --- a/docs/advanced/elec_properties/band.md +++ b/docs/advanced/elec_properties/band.md @@ -8,7 +8,7 @@ out_chg 1 ``` With this input parameter, the converged charge density will be output in the files such as `chgs1.cube`, `chgs2.cube`, etc. -Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file onsite.dm if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: +Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: ``` INPUT_PARAMETERS diff --git a/docs/advanced/elec_properties/dos.md b/docs/advanced/elec_properties/dos.md index c89e9d47d8d..82777abf0f9 100644 --- a/docs/advanced/elec_properties/dos.md +++ b/docs/advanced/elec_properties/dos.md @@ -10,7 +10,7 @@ out_chg 1 ``` this will produce the converged charge density, which is contained in the file SPIN1_CHG.cube. -Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file onsite.dm if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: +Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: ``` INPUT_PARAMETERS diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 504c4df731a..6d5efef43be 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3736,10 +3736,10 @@ - **Type**: Integer - **Description**: The parameter controls the form of occupation matrix control used. - 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step. - - 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated. - - 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation. + - 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated. + - 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation. - > Note: The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward. + > Note: The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward. - **Default**: 0 ### onsite_radius diff --git a/docs/parameters.yaml b/docs/parameters.yaml index bcb562fa693..d39db4c0417 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -4387,10 +4387,10 @@ parameters: description: | The parameter controls the form of occupation matrix control used. * 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step. - * 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated. - * 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation. + * 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated. + * 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation. - [NOTE] The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward. + [NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward. default_value: "0" unit: "" availability: "" diff --git a/source/source_io/module_parameter/read_input_item_exx_dftu.cpp b/source/source_io/module_parameter/read_input_item_exx_dftu.cpp index 8f5f926a174..e4e2d4f4ec5 100644 --- a/source/source_io/module_parameter/read_input_item_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_input_item_exx_dftu.cpp @@ -834,10 +834,10 @@ void ReadInput::item_dftu() item.type = "Integer"; item.description = R"(The parameter controls the form of occupation matrix control used. * 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step. -* 1: The first SCF step will use an initial density matrix read from a file named initial_onsite.dm, but for later steps, the onsite density matrix will be updated. -* 2: The same onsite density matrix from initial_onsite.dm will be used throughout the entire calculation. +* 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated. +* 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation. -[NOTE] The easiest way to create initial_onsite.dm is to run a DFT+U calculation, look for a file named onsite.dm in the OUT.prefix directory, and make replacements there. The format of the file is rather straight-forward.)"; +[NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.)"; item.default_value = "0"; item.unit = ""; item.availability = ""; diff --git a/source/source_lcao/dftu_lcao.h b/source/source_lcao/dftu_lcao.h index 245b6d3f415..023a9a5ae2d 100644 --- a/source/source_lcao/dftu_lcao.h +++ b/source/source_lcao/dftu_lcao.h @@ -43,14 +43,14 @@ void init_dftu_lcao(const int istep, * @param iter Current SCF iteration * @param conv_esolver Whether ESolver has converged * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new) - * @param out_chg Whether to output onsite.dm + * @param out_chg Whether to output dm_onsite.txt * @param dftu DFT+U object * @param ucell Unit cell * @param dm_vec Density matrix vector * @param kv K-vectors * @param mixing_beta Mixing beta parameter * @param hamilt_lcao Hamiltonian LCAO object - * @param global_out_dir Output directory for onsite.dm + * @param global_out_dir Output directory for dm_onsite.txt * @param nspin Number of spin channels (1, 2, or 4) * @param npol Number of polarizations */ diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp index f3f967f86f2..5671c879eeb 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -237,7 +237,7 @@ void Plus_U::init(UnitCell& cell, // unitcell class if (omc != 0) { std::stringstream sst; - sst << "initial_onsite.dm"; + sst << PARAM.globalv.global_readin_dir << "dm_onsite_ini.txt"; this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); #ifdef __MPI this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); @@ -251,7 +251,7 @@ void Plus_U::init(UnitCell& cell, // unitcell class if (PARAM.inp.init_chg == "file") { std::stringstream sst; - sst << PARAM.globalv.global_readin_dir << "onsite.dm"; + sst << PARAM.globalv.global_out_dir << "dm_onsite.txt"; this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); #ifdef __MPI this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); diff --git a/source/source_lcao/module_dftu/dftu_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index 1f04684c103..9be899e244f 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -59,14 +59,14 @@ void Plus_U::output(const UnitCell& ucell, GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; this->write_occup_m(ucell, GlobalV::ofs_running, true, nspin, npol); - // Write onsite.dm + // Write dm_onsite.txt if (out_chg && GlobalV::MY_RANK == 0) { std::ofstream ofdftu; - ofdftu.open(global_out_dir + "onsite.dm"); + ofdftu.open(global_out_dir + "dm_onsite.txt"); if (!ofdftu) { - ModuleBase::WARNING_QUIT("Plus_U::output", "Can't create file onsite.dm"); + ModuleBase::WARNING_QUIT("Plus_U::output", "Can't create file dm_onsite.txt"); } this->write_occup_m(ucell, ofdftu, false, nspin, npol); ofdftu.close(); @@ -257,16 +257,16 @@ void Plus_U::read_occup_m(const UnitCell& ucell, { if (omc > 0) { - ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file initial_onsite.dm. Please check your initial_onsite.dm"); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); } else { if (init_chg == "file") { - ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file onsite.dm. Please do scf calculation first"); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file dm_onsite.txt. Please do scf calculation first"); } } - ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not open onsite.dm file"); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not open dm_onsite.txt file"); } ifdftu.clear(); diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp index d6e7bda9ba1..9885878c8cb 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp @@ -190,7 +190,7 @@ void hamilt::DFTU>::cal_nlm_all(const Parallel_Orbi * * For nspin=1: occ is scaled by 0.5 (since only one spin channel computed) * - Subsequent iterations: locale is computed fresh each iteration from updated DMR * - * Case 2: Locale IS initialized (is_locale_initialized, i.e., read from onsite.dm file) + * Case 2: Locale IS initialized (is_locale_initialized, i.e., read from dm_onsite.txt file) * - First electronic iteration: uses pre-read locale directly without DMR calculation * * Skips DMR-based occ calculation entirely * * Reads locale from stored data via get_locale() @@ -334,8 +334,8 @@ void hamilt::DFTU>::contributeHR() // BRANCH 2: Locale IS initialized (use pre-read data) // ============================================================ // This branch is taken when: - // - is_locale_initialized() == true (locale read from onsite.dm file) - // - OR omc != 0 (occupation matrix control with initial_onsite.dm) + // - is_locale_initialized() == true (locale read from dm_onsite.txt file) + // - OR omc != 0 (occupation matrix control with dm_onsite_ini.txt) // Typical scenario: first SCF iteration with file input, or restart calculation else { diff --git a/tests/17_DS_DFTU/README.md b/tests/17_DS_DFTU/README.md index c9a6e6a7abf..a3b87cd3d01 100644 --- a/tests/17_DS_DFTU/README.md +++ b/tests/17_DS_DFTU/README.md @@ -175,9 +175,9 @@ The following test cases are disabled in `CASES_CPU.txt` (commented out with `#` - 09 (PW DFT+U + noncollinear): Only supports **2-process MPI** execution, `result.ref` reference files provided - The following test cases set `kpar=2` in INPUT and require at least **2 MPI processes** to run: 11, 12, 14, 15, 16, 18, 19, 21, 37, 39, 41, 43, 45 -- 62 (LCAO_DFTU_NSCF_Band_XY): Single-thread and multi-thread results are inconsistent; investigation shows HR, HK, and SK are consistent across threads, but eigenvalues from genelpa differ; switching to scalapack_gvx produces consistent results across thread counts. Note: this test is named "NSCF" but actually runs with `calculation = scf` (`scf_nmax = 1`), using pre-shipped charge density and onsite.dm files as initial guess +- 62 (LCAO_DFTU_NSCF_Band_XY): Single-thread and multi-thread results are inconsistent; investigation shows HR, HK, and SK are consistent across threads, but eigenvalues from genelpa differ; switching to scalapack_gvx produces consistent results across thread counts. Note: this test is named "NSCF" but actually runs with `calculation = scf` (`scf_nmax = 1`), using pre-shipped charge density and dm_onsite.txt files as initial guess - All NSCF tests (55, 60, 61, 62, 63, 64) have been **converted to SCF+NSCF workflow**: - - Pre-converged `autotest-CHARGE-DENSITY.restart` and `onsite.dm` files have been removed + - Pre-converged `autotest-CHARGE-DENSITY.restart` and `dm_onsite.txt` files have been removed - Each test directory contains a `scf/` subdirectory with SCF input files - Run with: `bash ../run_scf_nscf.sh [mpi_np]` - These tests are **disabled in CI** (commented out in CASES_CPU.txt) diff --git a/tests/17_DS_DFTU/run_scf_nscf.sh b/tests/17_DS_DFTU/run_scf_nscf.sh index 483e6856af4..3a73506cdc7 100755 --- a/tests/17_DS_DFTU/run_scf_nscf.sh +++ b/tests/17_DS_DFTU/run_scf_nscf.sh @@ -7,7 +7,7 @@ # # This script: # 1. Runs SCF calculation in scf/ subdirectory -# 2. Copies charge density (and onsite.dm for DFT+U) from SCF output +# 2. Copies charge density (and dm_onsite.txt for DFT+U) from SCF output # 3. Runs NSCF calculation in current directory # 4. Compares output with reference # @@ -56,11 +56,11 @@ if [ ! -f "${TEST_DIR}/STRU" ]; then exit 1 fi -# Check if this is a DFT+U test (needs onsite.dm handling) +# Check if this is a DFT+U test (needs dm_onsite.txt handling) IS_DFTU=false if grep -q "dft_plus_u" "${TEST_DIR}/INPUT" 2>/dev/null; then IS_DFTU=true - echo "DFT+U test detected: will handle onsite.dm files" + echo "DFT+U test detected: will handle dm_onsite.txt files" fi # ------------------------------------------------------- @@ -135,14 +135,18 @@ NSCF_CHG_FILE="${NSCF_SUFFIX}-${CHG_BASENAME#*-}" cp "${CHG_FILE}" "${TEST_DIR}/${NSCF_CHG_FILE}" echo " Copied to: ${TEST_DIR}/${NSCF_CHG_FILE}" -# For DFT+U tests, also copy onsite.dm if it exists +# For DFT+U tests, also copy dm_onsite.txt if it exists. +# dm_onsite.txt is read from OUT.prefix (global_out_dir), so we need to copy it +# to the NSCF's OUT directory (TEST_DIR/OUT.suffix/) if [ "${IS_DFTU}" = true ]; then - ONSITE_FILE=$(find "${SCF_OUT}" -name "onsite.dm" 2>/dev/null | head -1) + ONSITE_FILE=$(find "${SCF_OUT}" -name "dm_onsite.txt" 2>/dev/null | head -1) if [ -n "${ONSITE_FILE}" ]; then - cp "${ONSITE_FILE}" "${TEST_DIR}/onsite.dm" - echo " Copied onsite.dm" + NSCF_OUT_DIR="${TEST_DIR}/${NSCF_SUFFIX}" + mkdir -p "${NSCF_OUT_DIR}" + cp "${ONSITE_FILE}" "${NSCF_OUT_DIR}/dm_onsite.txt" + echo " Copied dm_onsite.txt to: ${NSCF_OUT_DIR}/dm_onsite.txt" else - echo " WARNING: onsite.dm not found in SCF output" + echo " WARNING: dm_onsite.txt not found in SCF output" fi fi From e9491f321cac05cd2e67f9a49d7d57ca7fe23569 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 13:31:28 +0800 Subject: [PATCH 07/10] remove some usage of PARAM --- source/source_lcao/LCAO_set.cpp | 5 ++++- source/source_lcao/module_dftu/dftu.cpp | 14 +++++++------- source/source_lcao/module_dftu/dftu.h | 7 +++++-- source/source_pw/module_pwdft/setup_pot.cpp | 4 +++- 4 files changed, 19 insertions(+), 11 deletions(-) diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp index 9f76ac2ba98..ea8ba028390 100644 --- a/source/source_lcao/LCAO_set.cpp +++ b/source/source_lcao/LCAO_set.cpp @@ -81,7 +81,10 @@ void LCAO_domain::set_pot( //! 3) initialize DFT+U if (inp.dft_plus_u) { - dftu.init(ucell, &pv, kv.get_nks(), &orb); + dftu.init(ucell, &pv, + PARAM.globalv.npol, + inp.nspin, inp.orbital_corr, inp.yukawa_lambda, + &orb); } //! 4) init exact exchange calculations diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp index 5671c879eeb..1b4fac4a457 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -43,9 +43,12 @@ Plus_U::Plus_U() Plus_U::~Plus_U() {} -void Plus_U::init(UnitCell& cell, // unitcell class +void Plus_U::init(UnitCell& cell, const Parallel_Orbitals* pv, - const int nks + const int npol, + const int nspin, + const std::vector& orbital_corr, + const double yukawa_lambda #ifdef __LCAO , const LCAO_Orbitals* orb #endif @@ -69,12 +72,9 @@ void Plus_U::init(UnitCell& cell, // unitcell class ucell = &cell; #endif - // needs reconstructions in future - // global parameters, need to be removed in future - const int npol = PARAM.globalv.npol; // number of polarization directions - const int nlocal = PARAM.globalv.nlocal; // number of total local orbitals - const int nspin = PARAM.inp.nspin; // number of spins Plus_U::nspin = nspin; + Plus_U::orbital_corr = orbital_corr; + Plus_U::Yukawa = (yukawa_lambda > 0); // mohan update 2025-11-06 Plus_U::energy_u = 0.0; diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h index 5b3ca2bcc08..1a98c266e9d 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -25,9 +25,12 @@ class Plus_U public: // allocate relevant data strcutures - void init(UnitCell& cell, // unitcell class + void init(UnitCell& cell, const Parallel_Orbitals* pv, - const int nks + const int npol, + const int nspin, + const std::vector& orbital_corr, + const double yukawa_lambda #ifdef __LCAO , const LCAO_Orbitals* orb = nullptr #endif diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index e6a5603546c..c57ddcc7a86 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -121,7 +121,9 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- if (PARAM.inp.dft_plus_u) { - dftu.init(ucell, nullptr, kv.get_nks()); + dftu.init(ucell, nullptr, + PARAM.globalv.npol, + PARAM.inp.nspin, PARAM.inp.orbital_corr, PARAM.inp.yukawa_lambda); } return; From 8cb95c8244bfaf1b582255b96a65ae0361aaba5f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 14:07:26 +0800 Subject: [PATCH 08/10] remove C++17 codes, reduce PARAM usage in dftu --- .../source_io/module_ctrl/ctrl_scf_lcao.cpp | 41 ++++++++++++------- source/source_io/module_hs/write_vxc.hpp | 2 +- source/source_lcao/FORCE_STRESS.cpp | 2 +- source/source_lcao/hamilt_lcao.cpp | 14 ++++--- source/source_lcao/module_dftu/dftu.h | 22 ++++++---- source/source_lcao/module_dftu/dftu_force.cpp | 34 ++++++++------- .../source_lcao/module_dftu/dftu_hamilt.cpp | 17 ++++---- source/source_lcao/module_dftu/dftu_tools.cpp | 20 ++++----- .../module_operator_lcao/op_dftu_lcao.cpp | 6 +-- .../module_operator_lcao/op_dftu_lcao.h | 13 +++--- source/source_lcao/spar_u.cpp | 4 +- 11 files changed, 101 insertions(+), 74 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index 91125f64637..7aafa470678 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -38,6 +38,28 @@ #include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24 #include "source_lcao/module_operator_lcao/overlap.h" // use hamilt::Overlap for NAMD +#ifdef __EXX +template +void set_exx_interface(ModuleIO::WriteDHParams& params, Exx_NAO& exx_nao) {} + +template <> +void set_exx_interface(ModuleIO::WriteDHParams& params, Exx_NAO& exx_nao) +{ + if (exx_nao.exd) { params.exd = exx_nao.exd.get(); } + if (exx_nao.exc) { params.exc = exx_nao.exc.get(); } +} + +template +void set_exx_interface(ModuleIO::WriteHParams& params, Exx_NAO& exx_nao) {} + +template <> +void set_exx_interface(ModuleIO::WriteHParams& params, Exx_NAO& exx_nao) +{ + if (exx_nao.exd) { params.exd = exx_nao.exd.get(); } + if (exx_nao.exc) { params.exc = exx_nao.exc.get(); } +} +#endif + template void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, const Input_para& inp, @@ -359,13 +381,9 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, #ifdef __EXX // dV^EXX/dR output is wired for the gamma (TK==double) exx interfaces. exd/exc are // mutually exclusive (real vs complex Hexx); write_dH_exx picks by info_ri.real_number. - if constexpr (std::is_same::value) + if (GlobalC::exx_info.info_global.cal_exx) { - if (GlobalC::exx_info.info_global.cal_exx) - { - if (exx_nao.exd) { dh_params.exd = exx_nao.exd.get(); } - if (exx_nao.exc) { dh_params.exc = exx_nao.exc.get(); } - } + set_exx_interface(dh_params, exx_nao); } #endif ModuleIO::write_dH_components(dh_params); @@ -419,15 +437,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, if (inp.out_mat_h_exx[0] && GlobalC::exx_info.info_global.cal_exx) { // V^EXX(R) output is wired for the gamma (TK==double) exx interfaces. - if constexpr (std::is_same::value) - { - if (GlobalC::exx_info.info_global.cal_exx) - { - if (exx_nao.exd) { h_params.exd = exx_nao.exd.get(); } - if (exx_nao.exc) { h_params.exc = exx_nao.exc.get(); } - ModuleIO::write_h_exx(h_params); - } - } + set_exx_interface(h_params, exx_nao); + ModuleIO::write_h_exx(h_params); } #endif } diff --git a/source/source_io/module_hs/write_vxc.hpp b/source/source_io/module_hs/write_vxc.hpp index 881f267e43e..01f364d26cc 100644 --- a/source/source_io/module_hs/write_vxc.hpp +++ b/source/source_io/module_hs/write_vxc.hpp @@ -206,7 +206,7 @@ void write_Vxc(const int nspin, &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); std::vector> e_orb_exx; // orbital energy (EXX) #endif - hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk); + hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk, PARAM.globalv.npol); // 4. calculate and write the MO-matrix Exc Parallel_2D p2d; diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index c5a0ca6e6b2..a88e2ac25ef 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -420,7 +420,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, std::vector>* dmk_d = nullptr; std::vector>>* dmk_c = nullptr; assign_dmk_ptr(dmat.dm, dmk_d, dmk_c, PARAM.globalv.gamma_only_local); - dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv); + dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr_dftu, force_u, stress_u, kv, PARAM.globalv.npol); } else { diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index 1d54414c590..20b81e86622 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -225,9 +225,10 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, { plus_u = new OperatorDFTU>(this->hsk, this->kv->kvec_d, - this->hR, // no explicit call yet - p_dftu, // mohan add 2025-11-07 - this->kv->isk); + this->hR, + p_dftu, + this->kv->isk, + PARAM.globalv.npol); } else { @@ -381,9 +382,10 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, { plus_u = new OperatorDFTU>(this->hsk, this->kv->kvec_d, - this->hR, // no explicit call yet - p_dftu, // mohan add 2025-11-07 - this->kv->isk); + this->hR, + p_dftu, + this->kv->isk, + PARAM.globalv.npol); } else { diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h index 1a98c266e9d..f5789b6f85a 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -116,18 +116,21 @@ class Plus_U void cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::vector& isk, - const std::complex* sk); + const std::complex* sk, + const int npol); void cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, - const double* sk); + const double* sk, + const int npol); - void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR); + void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, const int npol); void cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, - std::complex* HR); + std::complex* HR, + const int npol); #endif //============================================================= @@ -264,8 +267,8 @@ class Plus_U // for both Hamiltonian and force/stress //============================================================= - void cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU); - void cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU); + void cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU, const int npol); + void cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU, const int npol); double get_onebody_eff_pot(const int T, const int iat, @@ -321,13 +324,14 @@ class Plus_U public: void force_stress(const UnitCell& ucell, const Grid_Driver& gd, - std::vector>* dmk_d, // mohan modify 2025-11-02 - std::vector>>* dmk_c, // dmat.get_dm()->get_DMK_vector(); + std::vector>* dmk_d, + std::vector>>* dmk_c, const Parallel_Orbitals& pv, ForceStressArrays& fsr, ModuleBase::matrix& force_dftu, ModuleBase::matrix& stress_dftu, - const K_Vectors& kv); + const K_Vectors& kv, + const int npol); private: void cal_force_k(const UnitCell& ucell, diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_force.cpp index 3c2f2608abe..8e2ff284895 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_force.cpp @@ -24,13 +24,14 @@ void Plus_U::force_stress(const UnitCell& ucell, const Grid_Driver& gd, - std::vector>* dmk_d, // mohan modify 2025-11-02 - std::vector>>* dmk_c, // dmat.get_dm()->get_DMK_vector(); + std::vector>* dmk_d, + std::vector>>* dmk_c, const Parallel_Orbitals& pv, - ForceStressArrays& fsr, // mohan add 2024-06-16 + ForceStressArrays& fsr, ModuleBase::matrix& force_dftu, ModuleBase::matrix& stress_dftu, - const K_Vectors& kv) + const K_Vectors& kv, + const int npol) { ModuleBase::TITLE("Plus_U", "force_stress"); ModuleBase::timer::start("Plus_U", "force_stress"); @@ -63,7 +64,7 @@ void Plus_U::force_stress(const UnitCell& ucell, double* VU = new double[pv.nloc]; - this->cal_VU_pot_mat_real(spin, false, VU); + this->cal_VU_pot_mat_real(spin, false, VU, npol); #ifdef __MPI ScalapackConnector::gemm(transT, transN, nlocal, nlocal, nlocal, @@ -110,7 +111,7 @@ void Plus_U::force_stress(const UnitCell& ucell, std::complex* VU = new std::complex[pv.nloc]; - this->cal_VU_pot_mat_complex(spin, false, VU); + this->cal_VU_pot_mat_complex(spin, false, VU, npol); #ifdef __MPI @@ -375,9 +376,14 @@ void Plus_U::cal_force_gamma(const UnitCell& ucell, { ModuleBase::TITLE("Plus_U", "cal_force_gamma"); ModuleBase::timer::start("Plus_U", "cal_force_gamma"); - const char transN = 'N', transT = 'T'; + const char transN = 'N'; + const char transT = 'T'; const int one_int = 1; - const double one = 1.0, zero = 0.0, minus_one = -1.0; + const double one = 1.0; + const double zero = 0.0; + const double minus_one = -1.0; + const int nlocal = PARAM.globalv.nlocal; + assert(nlocal>0); std::vector dm_VU_dSm(pv.nloc); @@ -400,9 +406,9 @@ void Plus_U::cal_force_gamma(const UnitCell& ucell, #ifdef __MPI ScalapackConnector::gemm(transN, transT, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + nlocal, + nlocal, + nlocal, one, tmp_ptr, 1, @@ -438,9 +444,9 @@ void Plus_U::cal_force_gamma(const UnitCell& ucell, #ifdef __MPI ScalapackConnector::gemm(transN, transT, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + nlocal, + nlocal, + nlocal, one, tmp_ptr, 1, diff --git a/source/source_lcao/module_dftu/dftu_hamilt.cpp b/source/source_lcao/module_dftu/dftu_hamilt.cpp index e2c37039960..019366979ff 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/source_lcao/module_dftu/dftu_hamilt.cpp @@ -8,7 +8,8 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::vector& isk, - const std::complex* sk) + const std::complex* sk, + const int npol) { ModuleBase::TITLE("Plus_U", "cal_eff_pot_c"); if (!is_locale_initialized()) @@ -32,7 +33,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik, const std::complex zero = 0.0; std::vector> VU(this->paraV->nloc); - this->cal_VU_pot_mat_complex(spin, true, &VU[0]); + this->cal_VU_pot_mat_complex(spin, true, &VU[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, @@ -61,7 +62,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik, return; } -void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk) +void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk, const int npol) { ModuleBase::TITLE("Plus_U", "cal_eff_pot_r"); if (!is_locale_initialized()) @@ -82,7 +83,7 @@ void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vect double alpha = 1.0, beta = 0.0, half = 0.5, one = 1.0; std::vector VU(this->paraV->nloc); - this->cal_VU_pot_mat_real(spin, 1, &VU[0]); + this->cal_VU_pot_mat_real(spin, 1, &VU[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, @@ -109,14 +110,14 @@ void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vect return; } -void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) +void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, const int npol) { const char transN = 'N', transT = 'T'; const int one_int = 1; const double alpha = 1.0, beta = 0.0, one = 1.0, half = 0.5; std::vector VU(this->paraV->nloc); - this->cal_VU_pot_mat_real(ispin, 1, &VU[0]); + this->cal_VU_pot_mat_real(ispin, 1, &VU[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, @@ -139,14 +140,14 @@ void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) return; } -void Plus_U::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR) +void Plus_U::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR, const int npol) { const char transN = 'N', transT = 'T'; const int one_int = 1; const std::complex zero = 0.0, one = 1.0, half = 0.5; std::vector> VU(this->paraV->nloc); - this->cal_VU_pot_mat_complex(ispin, 1, &VU[0]); + this->cal_VU_pot_mat_complex(ispin, 1, &VU[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, diff --git a/source/source_lcao/module_dftu/dftu_tools.cpp b/source/source_lcao/module_dftu/dftu_tools.cpp index 1ae144a8c28..d66a4b2d9a5 100644 --- a/source/source_lcao/module_dftu/dftu_tools.cpp +++ b/source/source_lcao/module_dftu/dftu_tools.cpp @@ -3,14 +3,14 @@ #include "source_io/module_parameter/parameter.h" #ifdef __LCAO -void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU) +void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU, const int npol) { ModuleBase::TITLE("Plus_U", "cal_VU_pot_mat_complex"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); for (int it = 0; it < this->ucell->ntype; ++it) { - if (PARAM.inp.orbital_corr[it] == -1) + if (Plus_U::orbital_corr[it] == -1) { continue; } @@ -19,7 +19,7 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::c const int iat = this->ucell->itia2iat(it, ia); for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != PARAM.inp.orbital_corr[it]) + if (L != Plus_U::orbital_corr[it]) { continue; } @@ -33,7 +33,7 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::c for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ipol1++) { const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) @@ -43,7 +43,7 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::c for (int m2 = 0; m2 < 2 * L + 1; m2++) { - for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) + for (int ipol2 = 0; ipol2 < npol; ipol2++) { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); @@ -67,14 +67,14 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::c return; } -void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) +void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU, const int npol) { ModuleBase::TITLE("Plus_U", "cal_VU_pot_mat_real"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); for (int it = 0; it < this->ucell->ntype; ++it) { - if (PARAM.inp.orbital_corr[it] == -1) + if (Plus_U::orbital_corr[it] == -1) { continue; } @@ -83,7 +83,7 @@ void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* V const int iat = this->ucell->itia2iat(it, ia); for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != PARAM.inp.orbital_corr[it]) + if (L != Plus_U::orbital_corr[it]) { continue; } @@ -96,7 +96,7 @@ void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* V } for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ipol1++) { const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) @@ -105,7 +105,7 @@ void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* V } for (int m2 = 0; m2 < 2 * L + 1; m2++) { - for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) + for (int ipol2 = 0; ipol2 < npol; ipol2++) { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); diff --git a/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp index cab3203a41c..65198b46934 100644 --- a/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp @@ -27,7 +27,7 @@ void OperatorDFTU>::contributeHk(int ik) // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector eff_pot(this->hsk->get_pv()->nloc); - this->dftu->cal_eff_pot_mat_real(ik, &eff_pot[0], isk, this->hsk->get_sk()); + this->dftu->cal_eff_pot_mat_real(ik, &eff_pot[0], isk, this->hsk->get_sk(), this->npol); double* hk = this->hsk->get_hk(); @@ -48,7 +48,7 @@ void OperatorDFTU, double>>::contributeHk(int // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk(), this->npol); std::complex* hk = this->hsk->get_hk(); @@ -68,7 +68,7 @@ void OperatorDFTU, std::complex>>::con // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk(), this->npol); std::complex* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) diff --git a/source/source_lcao/module_operator_lcao/op_dftu_lcao.h b/source/source_lcao/module_operator_lcao/op_dftu_lcao.h index a8b275f0450..d2ef6c9e795 100644 --- a/source/source_lcao/module_operator_lcao/op_dftu_lcao.h +++ b/source/source_lcao/module_operator_lcao/op_dftu_lcao.h @@ -25,12 +25,13 @@ class OperatorDFTU> : public OperatorLCAO OperatorDFTU>(HS_Matrix_K* hsk_in, const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, - Plus_U* dftu_in, // mohan add 2025-11-05 - const std::vector& isk_in) - : isk(isk_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) + Plus_U* dftu_in, + const std::vector& isk_in, + const int npol_in) + : isk(isk_in), npol(npol_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_dftu; - this->dftu = dftu_in; // mohan add 2025-11-07 + this->dftu = dftu_in; } virtual void contributeHR() override; @@ -39,11 +40,13 @@ class OperatorDFTU> : public OperatorLCAO private: - Plus_U *dftu; // mohan add 20251107 + Plus_U *dftu; bool HR_fixed_done = false; const std::vector& isk; + + const int npol; }; } // namespace hamilt #endif diff --git a/source/source_lcao/spar_u.cpp b/source/source_lcao/spar_u.cpp index 2276ed94718..462b43c8e6e 100644 --- a/source/source_lcao/spar_u.cpp +++ b/source/source_lcao/spar_u.cpp @@ -73,7 +73,7 @@ void sparse_format::cal_HR_dftu( } } - dftu.cal_eff_pot_mat_R_double(current_spin, SR_tmp, HR_tmp); + dftu.cal_eff_pot_mat_R_double(current_spin, SR_tmp, HR_tmp, PARAM.globalv.npol); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { @@ -192,7 +192,7 @@ void sparse_format::cal_HR_dftu_soc( } } - dftu.cal_eff_pot_mat_R_complex_double(current_spin, SR_soc_tmp, HR_soc_tmp); + dftu.cal_eff_pot_mat_R_complex_double(current_spin, SR_soc_tmp, HR_soc_tmp, PARAM.globalv.npol); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { From 2f2b12fd64edebf9d4b147f8ede520fd6c430d8d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 15:40:13 +0800 Subject: [PATCH 09/10] update --- source/source_lcao/LCAO_set.cpp | 11 +++- source/source_lcao/module_dftu/dftu.cpp | 56 ++++++++++++----- source/source_lcao/module_dftu/dftu.h | 42 +++++++++---- .../source_lcao/module_dftu/dftu_folding.cpp | 20 +++--- source/source_lcao/module_dftu/dftu_force.cpp | 32 +++++----- .../source_lcao/module_dftu/dftu_hamilt.cpp | 16 ++--- source/source_lcao/module_dftu/dftu_occup.cpp | 63 +++++++++---------- source/source_lcao/module_dftu/dftu_pw.cpp | 30 ++++----- .../source_lcao/module_dftu/dftu_yukawa.cpp | 12 ++-- source/source_pw/module_pwdft/setup_pot.cpp | 13 +++- 10 files changed, 176 insertions(+), 119 deletions(-) diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp index ea8ba028390..6cce705817d 100644 --- a/source/source_lcao/LCAO_set.cpp +++ b/source/source_lcao/LCAO_set.cpp @@ -78,12 +78,21 @@ void LCAO_domain::set_pot( &(pelec->f_en.etxc), &(pelec->f_en.vtxc)); } - //! 3) initialize DFT+U if (inp.dft_plus_u) { dftu.init(ucell, &pv, PARAM.globalv.npol, inp.nspin, inp.orbital_corr, inp.yukawa_lambda, + PARAM.globalv.global_readin_dir, + PARAM.globalv.global_out_dir, + inp.init_chg, + pv.nloc, + PARAM.globalv.gamma_only_local, + inp.ks_solver, + inp.cal_force, + inp.cal_stress, + inp.device, + inp.kpar, &orb); } diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp index 1b4fac4a457..28d49abb43c 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -48,7 +48,17 @@ void Plus_U::init(UnitCell& cell, const int npol, const int nspin, const std::vector& orbital_corr, - const double yukawa_lambda + const double yukawa_lambda, + const std::string& global_readin_dir, + const std::string& global_out_dir, + const std::string& init_chg, + const int nlocal, + const bool gamma_only_local, + const std::string& ks_solver, + const bool cal_force, + const bool cal_stress, + const std::string& device, + const int kpar #ifdef __LCAO , const LCAO_Orbitals* orb #endif @@ -75,6 +85,19 @@ void Plus_U::init(UnitCell& cell, Plus_U::nspin = nspin; Plus_U::orbital_corr = orbital_corr; Plus_U::Yukawa = (yukawa_lambda > 0); + this->yukawa_lambda = yukawa_lambda; + + this->global_readin_dir = global_readin_dir; + this->global_out_dir = global_out_dir; + this->init_chg = init_chg; + this->npol = npol; + this->nlocal = nlocal; + this->gamma_only_local = gamma_only_local; + this->ks_solver = ks_solver; + this->cal_force = cal_force; + this->cal_stress = cal_stress; + this->device = device; + this->kpar = kpar; // mohan update 2025-11-06 Plus_U::energy_u = 0.0; @@ -237,10 +260,10 @@ void Plus_U::init(UnitCell& cell, if (omc != 0) { std::stringstream sst; - sst << PARAM.globalv.global_readin_dir << "dm_onsite_ini.txt"; - this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); + sst << this->global_readin_dir << "dm_onsite_ini.txt"; + this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol); #ifdef __MPI - this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); + this->local_occup_bcast(cell, nspin, npol); #endif mark_locale_initialized(); @@ -248,13 +271,13 @@ void Plus_U::init(UnitCell& cell, } else { - if (PARAM.inp.init_chg == "file") + if (this->init_chg == "file") { std::stringstream sst; - sst << PARAM.globalv.global_out_dir << "dm_onsite.txt"; - this->read_occup_m(cell, sst.str(), PARAM.inp.init_chg, PARAM.inp.nspin, PARAM.globalv.npol); + sst << this->global_out_dir << "dm_onsite.txt"; + this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol); #ifdef __MPI - this->local_occup_bcast(cell, PARAM.inp.nspin, PARAM.globalv.npol); + this->local_occup_bcast(cell, nspin, npol); #endif mark_locale_initialized(); } @@ -319,7 +342,7 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell, continue; } - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { for (int spin = 0; spin < 2; spin++) { @@ -346,21 +369,21 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell, } } } - else if (PARAM.inp.nspin == 4) // SOC + else if (Plus_U::nspin == 4) { double nm_trace = 0.0; double nm2_trace = 0.0; for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < this->npol; ipol0++) { const int m0_all = m0 + (2 * l + 1) * ipol0; nm_trace += this->locale[iat][l][n][0](m0_all, m0_all); for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) { int m1_all = m1 + (2 * l + 1) * ipol1; @@ -381,19 +404,18 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell, } } - // calculate the double counting term included in eband for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) { const int m1_all = m1 + ipol1 * (2 * l + 1); for (int m2 = 0; m2 < 2 * l + 1; m2++) { - for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) + for (int ipol2 = 0; ipol2 < this->npol; ipol2++) { const int m2_all = m2 + ipol2 * (2 * l + 1); - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { for (int is = 0; is < 2; is++) { @@ -402,7 +424,7 @@ void Plus_U::cal_energy_correction(const UnitCell& ucell, energy_dc += VU * this->locale[iat][l][n][is](m1_all, m2_all); } } - else if (PARAM.inp.nspin == 4) // SOC + else if (Plus_U::nspin == 4) { double VU = 0.0; VU = get_onebody_eff_pot(T, iat, l, n, 0, m1_all, m2_all, false); diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h index f5789b6f85a..c41b3fe632a 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -26,15 +26,25 @@ class Plus_U public: // allocate relevant data strcutures void init(UnitCell& cell, - const Parallel_Orbitals* pv, - const int npol, - const int nspin, - const std::vector& orbital_corr, - const double yukawa_lambda + const Parallel_Orbitals* pv, + const int npol, + const int nspin, + const std::vector& orbital_corr, + const double yukawa_lambda, + const std::string& global_readin_dir, + const std::string& global_out_dir, + const std::string& init_chg, + const int nlocal, + const bool gamma_only_local, + const std::string& ks_solver, + const bool cal_force, + const bool cal_stress, + const std::string& device, + const int kpar #ifdef __LCAO - , const LCAO_Orbitals* orb = nullptr + , const LCAO_Orbitals* orb = nullptr #endif - ); + ); // calculate the energy correction void cal_energy_correction(const UnitCell& ucell, const int istep); @@ -93,15 +103,25 @@ class Plus_U static double energy_u; //+U energy, mohan update 2025-11-06, change this to private const Parallel_Orbitals* paraV = nullptr; - int cal_type = 3; // 1:dftu_tpye=1, dc=1; 2:dftu_type=1, dc=2; 3:dftu_tpye=2, dc=1; 4:dftu_tpye=2, dc=2; + int cal_type = 3; - // FIXME: the following variable does not have static lifetime; - // while the present class is used via a global variable. This has - // potential to cause dangling pointer issues. #ifdef __LCAO const LCAO_Orbitals* ptr_orb_ = nullptr; std::vector orb_cutoff_; #endif + + std::string global_readin_dir; + std::string global_out_dir; + double yukawa_lambda = 0.0; + std::string init_chg; + int npol = 1; + int nlocal = 0; + bool gamma_only_local = false; + std::string ks_solver; + bool cal_force = false; + bool cal_stress = false; + std::string device; + int kpar = 1; // transform between iwt index and it, ia, L, N and m index std::vector>>>> diff --git a/source/source_lcao/module_dftu/dftu_folding.cpp b/source/source_lcao/module_dftu/dftu_folding.cpp index acef2112bc5..1164ce8bc80 100644 --- a/source/source_lcao/module_dftu/dftu_folding.cpp +++ b/source/source_lcao/module_dftu/dftu_folding.cpp @@ -88,9 +88,9 @@ void Plus_U::fold_dSR_gamma(const UnitCell& ucell, if (adj) { - for (int jj = 0; jj < atom1->nw * PARAM.globalv.npol; ++jj) + for (int jj = 0; jj < atom1->nw * this->npol; ++jj) { - const int jj0 = jj / PARAM.globalv.npol; + const int jj0 = jj / this->npol; const int iw1_all = start1 + jj0; const int mu = pv.global2local_row(iw1_all); if (mu < 0) @@ -98,9 +98,9 @@ void Plus_U::fold_dSR_gamma(const UnitCell& ucell, continue; } - for (int kk = 0; kk < atom2->nw * PARAM.globalv.npol; ++kk) + for (int kk = 0; kk < atom2->nw * this->npol; ++kk) { - const int kk0 = kk / PARAM.globalv.npol; + const int kk0 = kk / this->npol; const int iw2_all = start2 + kk0; const int nu = pv.global2local_col(iw2_all); if (nu < 0) @@ -227,7 +227,7 @@ void Plus_U::folding_matrix_k(const UnitCell& ucell, // calculate how many matrix elements are in // this processor. //-------------------------------------------------- - for (int ii = 0; ii < atom1->nw * PARAM.globalv.npol; ii++) + for (int ii = 0; ii < atom1->nw * this->npol; ii++) { // the index of orbitals in this processor const int iw1_all = start1 + ii; @@ -237,7 +237,7 @@ void Plus_U::folding_matrix_k(const UnitCell& ucell, continue; } - for (int jj = 0; jj < atom2->nw * PARAM.globalv.npol; jj++) + for (int jj = 0; jj < atom2->nw * this->npol; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); @@ -247,7 +247,7 @@ void Plus_U::folding_matrix_k(const UnitCell& ucell, } int iic = 0; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(this->ks_solver)) { iic = mu + nu * pv.nrow; } @@ -285,20 +285,20 @@ void Plus_U::folding_matrix_k_new(const int ik, ModuleBase::timer::start("Plus_U", "folding_matrix_k_new"); int hk_type = 0; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(this->ks_solver)) { hk_type = 1; } // get SR and fold to mat_k - if(PARAM.globalv.gamma_only_local) + if(this->gamma_only_local) { dynamic_cast*>(p_ham) ->updateSk(ik, hk_type); } else { - if(PARAM.inp.nspin != 4) + if(Plus_U::nspin != 4) { dynamic_cast, double>*>(p_ham) ->updateSk(ik, hk_type); diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_force.cpp index bfc8f4d1d89..1f58d094d9f 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_force.cpp @@ -36,18 +36,18 @@ void Plus_U::force_stress(const UnitCell& ucell, ModuleBase::TITLE("Plus_U", "force_stress"); ModuleBase::timer::start("Plus_U", "force_stress"); - const int nlocal = PARAM.globalv.nlocal; + const int nlocal = this->nlocal; - if (PARAM.inp.cal_force) + if (this->cal_force) { force_dftu.zero_out(); } - if (PARAM.inp.cal_stress) + if (this->cal_stress) { stress_dftu.zero_out(); } - if (PARAM.globalv.gamma_only_local) + if (this->gamma_only_local) { const char transN = 'N'; const char transT = 'T'; @@ -76,12 +76,12 @@ void Plus_U::force_stress(const UnitCell& ucell, delete[] VU; - if (PARAM.inp.cal_force) + if (this->cal_force) { this->cal_force_gamma(ucell,&rho_VU[0], pv, fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z, force_dftu); } - if (PARAM.inp.cal_stress) + if (this->cal_stress) { this->cal_stress_gamma(ucell, pv, @@ -123,23 +123,23 @@ void Plus_U::force_stress(const UnitCell& ucell, delete[] VU; - if (PARAM.inp.cal_force) + if (this->cal_force) { cal_force_k(ucell, gd, fsr, pv, ik, &rho_VU[0], force_dftu, kv.kvec_d[ik]); } - if (PARAM.inp.cal_stress) + if (this->cal_stress) { cal_stress_k(ucell, gd, fsr, pv, ik, &rho_VU[0], stress_dftu, kv.kvec_d[ik]); } } // ik } - if (PARAM.inp.cal_force) + if (this->cal_force) { Parallel_Reduce::reduce_pool(force_dftu.c, force_dftu.nr * force_dftu.nc); } - if (PARAM.inp.cal_stress) + if (this->cal_stress) { Parallel_Reduce::reduce_pool(stress_dftu.c, stress_dftu.nr * stress_dftu.nc); @@ -183,7 +183,7 @@ void Plus_U::cal_force_k(const UnitCell& ucell, const std::complex zero(0.0, 0.0); const std::complex one(1.0, 0.0); - const int nlocal = PARAM.globalv.nlocal; + const int nlocal = this->nlocal; assert(nlocal>0); std::vector> dm_VU_dSm(pv.nloc); @@ -277,7 +277,7 @@ void Plus_U::cal_force_k(const UnitCell& ucell, for (int m = 0; m < 2 * l + 1; m++) { - for (int ipol = 0; ipol < PARAM.globalv.npol; ipol++) + for (int ipol = 0; ipol < this->npol; ipol++) { const int iwt = this->iatlnmipol2iwt[iat][l][n][m][ipol]; const int mu = pv.global2local_row(iwt); @@ -310,7 +310,7 @@ void Plus_U::cal_stress_k(const UnitCell& ucell, ModuleBase::TITLE("Plus_U", "cal_stress_k"); ModuleBase::timer::start("Plus_U", "cal_stress_k"); - const int nlocal = PARAM.globalv.nlocal; + const int nlocal = this->nlocal; const char transN = 'N'; const int one_int = 1; @@ -385,7 +385,7 @@ void Plus_U::cal_force_gamma(const UnitCell& ucell, const double one = 1.0; const double zero = 0.0; const double minus_one = -1.0; - const int nlocal = PARAM.globalv.nlocal; + const int nlocal = this->nlocal; assert(nlocal>0); std::vector dm_VU_dSm(pv.nloc); @@ -492,7 +492,7 @@ void Plus_U::cal_force_gamma(const UnitCell& ucell, // Calculate the local occupation number matrix for (int m = 0; m < 2 * l + 1; m++) { - for (int ipol = 0; ipol < PARAM.globalv.npol; ipol++) + for (int ipol = 0; ipol < this->npol; ipol++) { const int iwt = this->iatlnmipol2iwt[iat][l][n][m][ipol]; const int mu = pv.global2local_row(iwt); @@ -536,7 +536,7 @@ void Plus_U::cal_stress_gamma(const UnitCell& ucell, std::vector dSR_gamma(pv.nloc); std::vector dm_VU_sover(pv.nloc); - const int nlocal = PARAM.globalv.nlocal; + const int nlocal = this->nlocal; for (int dim1 = 0; dim1 < 3; dim1++) { diff --git a/source/source_lcao/module_dftu/dftu_hamilt.cpp b/source/source_lcao/module_dftu/dftu_hamilt.cpp index 019366979ff..19c20a1c095 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/source_lcao/module_dftu/dftu_hamilt.cpp @@ -37,7 +37,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik, #ifdef __MPI ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, PARAM.globalv.nlocal, PARAM.globalv.nlocal, + this->nlocal, this->nlocal, this->nlocal, half, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), one_int, one_int, this->paraV->desc, sk, one_int, one_int, this->paraV->desc, @@ -51,7 +51,7 @@ void Plus_U::cal_eff_pot_mat_complex(const int ik, } #ifdef __MPI - ScalapackConnector::tranu(PARAM.globalv.nlocal, PARAM.globalv.nlocal, + ScalapackConnector::tranu(this->nlocal, this->nlocal, one, &VU[0], one_int, one_int, this->paraV->desc, one, @@ -87,7 +87,7 @@ void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vect #ifdef __MPI ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, PARAM.globalv.nlocal, PARAM.globalv.nlocal, + this->nlocal, this->nlocal, this->nlocal, half, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), 1, 1, this->paraV->desc, sk, 1, 1, this->paraV->desc, @@ -99,7 +99,7 @@ void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vect VU[irc] = eff_pot[irc]; #ifdef __MPI - pdtran_(&PARAM.globalv.nlocal, &PARAM.globalv.nlocal, + pdtran_(&this->nlocal, &this->nlocal, &one, &VU[0], &one_int, &one_int, const_cast(this->paraV->desc), &one, @@ -121,7 +121,7 @@ void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, c #ifdef __MPI ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, PARAM.globalv.nlocal, PARAM.globalv.nlocal, + this->nlocal, this->nlocal, this->nlocal, half, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), 1, 1, this->paraV->desc, SR, 1, 1, this->paraV->desc, @@ -129,7 +129,7 @@ void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, c HR, 1, 1, this->paraV->desc); ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, PARAM.globalv.nlocal, PARAM.globalv.nlocal, + this->nlocal, this->nlocal, this->nlocal, half, SR, 1, 1, this->paraV->desc, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), 1, 1, this->paraV->desc, @@ -151,7 +151,7 @@ void Plus_U::cal_eff_pot_mat_R_complex_double(const int ispin, std::complexnlocal, this->nlocal, this->nlocal, half, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), one_int, one_int, this->paraV->desc, SR, one_int, one_int, this->paraV->desc, @@ -159,7 +159,7 @@ void Plus_U::cal_eff_pot_mat_R_complex_double(const int ispin, std::complexparaV->desc); ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, PARAM.globalv.nlocal, PARAM.globalv.nlocal, + this->nlocal, this->nlocal, this->nlocal, half, SR, one_int, one_int, this->paraV->desc, ModuleBase::GlobalFunc::VECTOR_TO_PTR(VU), one_int, one_int, this->paraV->desc, diff --git a/source/source_lcao/module_dftu/dftu_occup.cpp b/source/source_lcao/module_dftu/dftu_occup.cpp index 311861e0c75..e44c8203e90 100644 --- a/source/source_lcao/module_dftu/dftu_occup.cpp +++ b/source/source_lcao/module_dftu/dftu_occup.cpp @@ -27,10 +27,9 @@ void Plus_U::copy_locale(const UnitCell& ucell) { const int iat = ucell.itia2iat(T, I); - if (PARAM.inp.nspin == 4) + if (Plus_U::nspin == 4) { locale_save[iat][target_l][0][0] = locale[iat][target_l][0][0]; - // nspin=4 locale matrix already contains all spin components interleaved if(this->uom_save.size() != 0) { const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc; @@ -40,11 +39,10 @@ void Plus_U::copy_locale(const UnitCell& ucell) } } } - else if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + else if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { locale_save[iat][target_l][0][0] = locale[iat][target_l][0][0]; locale_save[iat][target_l][0][1] = locale[iat][target_l][0][1]; - // save locale matrix for spin=0,1 to uom_save if(this->uom_save.size() != 0) { const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc; @@ -83,11 +81,11 @@ void Plus_U::zero_locale(const UnitCell& ucell) for (int n = 0; n < N; n++) { - if (PARAM.inp.nspin == 4) + if (Plus_U::nspin == 4) { locale[iat][l][n][0].zero_out(); } - else if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + else if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { locale[iat][l][n][0].zero_out(); locale[iat][l][n][1].zero_out(); @@ -117,7 +115,7 @@ void Plus_U::mix_locale(const UnitCell& ucell, { const int iat = ucell.itia2iat(T, I); - if (PARAM.inp.nspin == 4) + if (Plus_U::nspin == 4) { const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc; for (int mm = 0; mm < size; mm++) @@ -132,7 +130,7 @@ void Plus_U::mix_locale(const UnitCell& ucell, } } } - else if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + else if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc; const int half_size = this->uom_save.size() / 2; @@ -173,18 +171,18 @@ void Plus_U::set_locale(const UnitCell& ucell) for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); - if (PARAM.inp.nspin == 4) + if (Plus_U::nspin == 4) { for(int mm = 0; mm < locale[iat][l][0][0].nr * locale[iat][l][0][0].nc; mm++) locale[iat][l][0][0].c[mm] = this->uom_array[eff_pot_pw_index[iat] + mm]; } - else if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + else if (Plus_U::nspin == 1 || Plus_U::nspin == 2) { const int half_size = this->uom_array.size() / 2; for(int mm = 0; mm < locale[iat][l][0][0].nr * locale[iat][l][0][0].nc; mm++) { locale[iat][l][0][0].c[mm] = this->uom_array[eff_pot_pw_index[iat] + mm]; - if (PARAM.inp.nspin == 2) + if (Plus_U::nspin == 2) { locale[iat][l][0][1].c[mm] = this->uom_array[half_size + eff_pot_pw_index[iat] + mm]; } @@ -259,7 +257,7 @@ void Plus_U::cal_occup_m_k(const int iter, std::complex* s_k_pointer = nullptr; - if(PARAM.inp.nspin != 4) + if(Plus_U::nspin != 4) { s_k_pointer = dynamic_cast, double>*>(p_ham)->getSk(); } @@ -271,16 +269,15 @@ void Plus_U::cal_occup_m_k(const int iter, #ifdef __MPI ScalapackConnector::gemm(transN, transT, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + this->nlocal, + this->nlocal, + this->nlocal, alpha, s_k_pointer, one_int, one_int, &this->paraV->desc[0], dm_k[ik].data(), - //dm_k[ik].c, one_int, one_int, &this->paraV->desc[0], @@ -346,7 +343,7 @@ void Plus_U::cal_occup_m_k(const int iter, // Calculate the local occupation number matrix for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < this->npol; ipol0++) { const int iwt0 = this->iatlnmipol2iwt[iat][l][n][m0][ipol0]; const int mu = this->paraV->global2local_row(iwt0); @@ -354,7 +351,7 @@ void Plus_U::cal_occup_m_k(const int iter, for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) { const int iwt1 = this->iatlnmipol2iwt[iat][l][n][m1][ipol1]; const int nu = this->paraV->global2local_col(iwt1); @@ -419,17 +416,17 @@ void Plus_U::cal_occup_m_k(const int iter, // set the local occupation mumber matrix of spin up and down zeros #ifdef __MPI - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 4) + if (Plus_U::nspin == 1 || Plus_U::nspin == 4) { ModuleBase::matrix temp(locale[iat][l][n][0]); MPI_Allreduce(&temp(0, 0), &locale[iat][l][n][0](0, 0), - (2 * l + 1) * PARAM.globalv.npol * (2 * l + 1) * PARAM.globalv.npol, + (2 * l + 1) * this->npol * (2 * l + 1) * this->npol, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } - else if (PARAM.inp.nspin == 2) + else if (Plus_U::nspin == 2) { ModuleBase::matrix temp0(locale[iat][l][n][0]); MPI_Allreduce(&temp0(0, 0), @@ -449,8 +446,7 @@ void Plus_U::cal_occup_m_k(const int iter, } #endif - // for the case spin independent calculation - switch (PARAM.inp.nspin) + switch (Plus_U::nspin) { case 1: locale[iat][l][n][0] += transpose(locale[iat][l][n][0]); @@ -459,11 +455,11 @@ void Plus_U::cal_occup_m_k(const int iter, break; case 2: - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < Plus_U::nspin; is++) locale[iat][l][n][is] += transpose(locale[iat][l][n][is]); break; - case 4: // SOC + case 4: locale[iat][l][n][0] += transpose(locale[iat][l][n][0]); break; @@ -504,17 +500,16 @@ void Plus_U::cal_occup_m_gamma(const int iter, const double alpha = 1.0, beta = 0.0; std::vector srho(this->paraV->nloc); - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < Plus_U::nspin; is++) { - // srho(mu,nu) = \sum_{iw} S(mu,iw)*dm_gamma(iw,nu) double* s_gamma_pointer = dynamic_cast*>(p_ham)->getSk(); #ifdef __MPI ScalapackConnector::gemm(transN, transT, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + this->nlocal, + this->nlocal, + this->nlocal, alpha, s_gamma_pointer, one_int, @@ -562,7 +557,7 @@ void Plus_U::cal_occup_m_gamma(const int iter, // Calculate the local occupation number matrix for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) + for (int ipol0 = 0; ipol0 < this->npol; ipol0++) { const int iwt0 = this->iatlnmipol2iwt[iat][l][n][m0][ipol0]; const int mu = this->paraV->global2local_row(iwt0); @@ -570,7 +565,7 @@ void Plus_U::cal_occup_m_gamma(const int iter, for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) { const int iwt1 = this->iatlnmipol2iwt[iat][l][n][m1][ipol1]; const int nu = this->paraV->global2local_col(iwt1); @@ -604,14 +599,14 @@ void Plus_U::cal_occup_m_gamma(const int iter, #ifdef __MPI MPI_Allreduce(&temp(0, 0), &locale[iat][l][n][is](0, 0), - (2 * l + 1) * PARAM.globalv.npol * (2 * l + 1) * PARAM.globalv.npol, + (2 * l + 1) * this->npol * (2 * l + 1) * this->npol, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif // for the case spin independent calculation - switch (PARAM.inp.nspin) + switch (Plus_U::nspin) { case 1: locale[iat][l][n][0] += transpose(locale[iat][l][n][0]); diff --git a/source/source_lcao/module_dftu/dftu_pw.cpp b/source/source_lcao/module_dftu/dftu_pw.cpp index c0386ffc36e..6ae77c0cd1a 100644 --- a/source/source_lcao/module_dftu/dftu_pw.cpp +++ b/source/source_lcao/module_dftu/dftu_pw.cpp @@ -29,7 +29,7 @@ void Plus_U::cal_occ_pw(const int iter, this->copy_locale(cell); this->zero_locale(cell); - if(PARAM.inp.device == "cpu") + if(this->device == "cpu") { auto* onsite_p = projectors::OnsiteProjector::get_instance(); const psi::Psi>* psi_p = (const psi::Psi>*)psi_in; @@ -37,7 +37,7 @@ void Plus_U::cal_occ_pw(const int iter, const int npol = psi_p->get_npol(); for(int ik = 0; ik < psi_p->get_nk(); ik++) { - int is = (PARAM.inp.nspin == 2) ? isk[ik] : 0; + int is = (Plus_U::nspin == 2) ? isk[ik] : 0; psi_p->fix_k(ik); onsite_p->tabulate_atomic(ik); @@ -59,7 +59,7 @@ void Plus_U::cal_occ_pw(const int iter, const int m_begin = target_l * target_l; const int tlp1 = 2 * target_l + 1; const int tlp1_2 = tlp1 * tlp1; - if(PARAM.inp.nspin == 4) + if(Plus_U::nspin == 4) { for(int ib = 0;ibget_npol(); for(int ik = 0; ik < psi_p->get_nk(); ik++) { - int is = (PARAM.inp.nspin == 2) ? isk[ik] : 0; + int is = (Plus_U::nspin == 2) ? isk[ik] : 0; psi_p->fix_k(ik); onsite_p->tabulate_atomic(ik); @@ -138,7 +138,7 @@ void Plus_U::cal_occ_pw(const int iter, const int m_begin = target_l * target_l; const int tlp1 = 2 * target_l + 1; const int tlp1_2 = tlp1 * tlp1; - if(PARAM.inp.nspin == 4) + if(Plus_U::nspin == 4) { for(int ib = 0;ibkpar, GlobalV::NPROC_IN_POOL, this->locale[iat][target_l][0][0].c, size); - if(PARAM.inp.nspin == 2) + if(Plus_U::nspin == 2) { - Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, + Parallel_Reduce::reduce_double_allpool(this->kpar, GlobalV::NPROC_IN_POOL, this->locale[iat][target_l][0][1].c, size); @@ -215,7 +215,7 @@ void Plus_U::cal_occ_pw(const int iter, } else { - Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, + Parallel_Reduce::reduce_double_allpool(this->kpar, GlobalV::NPROC_IN_POOL, this->locale[iat][target_l][0][0].c, size * 4); @@ -228,7 +228,7 @@ void Plus_U::cal_occ_pw(const int iter, { this->uom_array[eff_pot_pw_index[iat]+mm] = this->locale[iat][target_l][0][0].c[mm]; } - if(PARAM.inp.nspin == 2) + if(Plus_U::nspin == 2) { const int half_size = this->uom_array.size() / 2; for(int mm=0;mm* vu_iat = &(this->eff_pot_pw[this->eff_pot_pw_index[iat]]); const int m_size = 2 * target_l + 1; - if(PARAM.inp.nspin == 4) + if(Plus_U::nspin == 4) { for (int m1 = 0; m1 < m_size; m1++) { @@ -328,7 +328,7 @@ void Plus_U::cal_occ_pw(const int iter, } } // spin-down channel for nspin=2 - if(PARAM.inp.nspin == 2) + if(Plus_U::nspin == 2) { std::complex* vu_iat1 = &(this->eff_pot_pw[this->eff_pot_pw.size()/2 + this->eff_pot_pw_index[iat]]); for (int m1 = 0; m1 < m_size; m1++) diff --git a/source/source_lcao/module_dftu/dftu_yukawa.cpp b/source/source_lcao/module_dftu/dftu_yukawa.cpp index c05e8acfee5..a5909267b18 100644 --- a/source/source_lcao/module_dftu/dftu_yukawa.cpp +++ b/source/source_lcao/module_dftu/dftu_yukawa.cpp @@ -18,17 +18,17 @@ void Plus_U::cal_yukawa_lambda(double** rho, const int& nrxx) { ModuleBase::TITLE("Plus_U", "cal_yukawa_lambda"); - if (PARAM.inp.yukawa_lambda > 0) + if (this->yukawa_lambda > 0) { - this->lambda = PARAM.inp.yukawa_lambda; + this->lambda = this->yukawa_lambda; return; } double sum_rho = 0.0; double sum_rho_lambda = 0.0; - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < Plus_U::nspin; is++) { - if(PARAM.inp.nspin == 4 && is > 0) + if(Plus_U::nspin == 4 && is > 0) { continue;// for non-collinear spin case, first spin contains the charge density } @@ -140,9 +140,9 @@ void Plus_U::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) { const int N = ucell.atoms[T].l_nchi[L]; - if (L >= PARAM.inp.orbital_corr[T] && PARAM.inp.orbital_corr[T] != -1) + if (L >= Plus_U::get_orbital_corr(T) && Plus_U::get_orbital_corr(T) != -1) { - if (L != PARAM.inp.orbital_corr[T]) + if (L != Plus_U::get_orbital_corr(T)) { continue; } diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index c57ddcc7a86..d81b71e0d61 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -121,9 +121,20 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- if (PARAM.inp.dft_plus_u) { + const int nlocal_dftu = 0; dftu.init(ucell, nullptr, PARAM.globalv.npol, - PARAM.inp.nspin, PARAM.inp.orbital_corr, PARAM.inp.yukawa_lambda); + PARAM.inp.nspin, PARAM.inp.orbital_corr, PARAM.inp.yukawa_lambda, + PARAM.globalv.global_readin_dir, + PARAM.globalv.global_out_dir, + PARAM.inp.init_chg, + nlocal_dftu, + PARAM.globalv.gamma_only_local, + PARAM.inp.ks_solver, + PARAM.inp.cal_force, + PARAM.inp.cal_stress, + PARAM.inp.device, + PARAM.inp.kpar); } return; From cd86a431f14ba92ac78c0d269176bfe618d473bf Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sun, 12 Jul 2026 15:41:01 +0800 Subject: [PATCH 10/10] remove all PARAM in dftu modules --- source/source_lcao/module_dftu/dftu_occup.cpp | 20 ------------------- 1 file changed, 20 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_occup.cpp b/source/source_lcao/module_dftu/dftu_occup.cpp index e44c8203e90..13ba48c5385 100644 --- a/source/source_lcao/module_dftu/dftu_occup.cpp +++ b/source/source_lcao/module_dftu/dftu_occup.cpp @@ -286,26 +286,6 @@ void Plus_U::cal_occup_m_k(const int iter, one_int, one_int, &this->paraV->desc[0]); - /*pzgemm_(&transN, - &transT, - &PARAM.globalv.nlocal, - &PARAM.globalv.nlocal, - &PARAM.globalv.nlocal, - &alpha, - s_k_pointer, - &one_int, - &one_int, - this->paraV->desc, - dm_k[ik].data(), - //dm_k[ik].c, - &one_int, - &one_int, - this->paraV->desc, - &beta, - &srho[0], - &one_int, - &one_int, - this->paraV->desc);*/ #endif const int spin = kv.isk[ik];