From 0eb07961f068a93ea6513ab78f874f459a300b53 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 30 Jun 2026 12:54:17 +0800 Subject: [PATCH 1/3] 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 2/3] 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 3/3] 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); } }