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_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 09e1b4b1177..ee3f22c8dc5 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_io/module_hs/write_vxc.hpp b/source/source_io/module_hs/write_vxc.hpp index d4719a17874..d6757f3cf0b 100644 --- a/source/source_io/module_hs/write_vxc.hpp +++ b/source/source_io/module_hs/write_vxc.hpp @@ -207,7 +207,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_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/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index d24c563e568..43e78b82a01 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/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp index 9f76ac2ba98..6cce705817d 100644 --- a/source/source_lcao/LCAO_set.cpp +++ b/source/source_lcao/LCAO_set.cpp @@ -78,10 +78,22 @@ 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, kv.get_nks(), &orb); + 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); } //! 4) init exact exchange calculations diff --git a/source/source_lcao/dftu_lcao.cpp b/source/source_lcao/dftu_lcao.cpp index d8b8421d6e7..5abde6a584b 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,30 @@ 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..023a9a5ae2d 100644 --- a/source/source_lcao/dftu_lcao.h +++ b/source/source_lcao/dftu_lcao.h @@ -1,9 +1,10 @@ #ifndef DFTU_LCAO_H #define DFTU_LCAO_H +#include +#include #include "source_cell/unitcell.h" #include "source_cell/klist.h" -#include "source_io/module_parameter/input_parameter.h" namespace ModuleESolver { @@ -16,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 @@ -26,7 +27,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,24 +42,32 @@ 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 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 dm_onsite.txt + * @param nspin Number of spin channels (1, 2, or 4) + * @param npol Number of polarizations */ 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/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.cpp b/source/source_lcao/module_dftu/dftu.cpp index ff06003c43d..28d49abb43c 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -43,9 +43,22 @@ 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, + 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 @@ -69,12 +82,22 @@ 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); + 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, // unitcell class if (omc != 0) { std::stringstream sst; - sst << "initial_onsite.dm"; - this->read_occup_m(cell,sst.str()); + 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); + this->local_occup_bcast(cell, nspin, npol); #endif mark_locale_initialized(); @@ -248,13 +271,13 @@ void Plus_U::init(UnitCell& cell, // unitcell class } else { - if (PARAM.inp.init_chg == "file") + if (this->init_chg == "file") { std::stringstream sst; - sst << PARAM.globalv.global_readin_dir << "onsite.dm"; - this->read_occup_m(cell,sst.str()); + 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); + 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 ff494d4fa3f..c41b3fe632a 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -25,13 +25,26 @@ class Plus_U public: // allocate relevant data strcutures - void init(UnitCell& cell, // unitcell class - const Parallel_Orbitals* pv, - const int nks + void init(UnitCell& cell, + 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); @@ -90,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>>>> @@ -113,18 +136,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 //============================================================= @@ -261,8 +287,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, @@ -318,13 +344,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, @@ -369,15 +396,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_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 3c2f2608abe..1f58d094d9f 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_force.cpp @@ -24,29 +24,30 @@ 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"); - 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'; @@ -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, @@ -75,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, @@ -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 @@ -122,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); @@ -182,6 +183,9 @@ 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 = this->nlocal; + assert(nlocal>0); + std::vector> dm_VU_dSm(pv.nloc); std::vector> dSm_k(pv.nloc); @@ -192,9 +196,9 @@ void Plus_U::cal_force_k(const UnitCell& ucell, #ifdef __MPI ScalapackConnector::gemm(transN, transC, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + nlocal, + nlocal, + nlocal, one, &dSm_k[0], one_int, @@ -230,9 +234,9 @@ void Plus_U::cal_force_k(const UnitCell& ucell, #ifdef __MPI ScalapackConnector::gemm(transN, transN, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, - PARAM.globalv.nlocal, + nlocal, + nlocal, + nlocal, one, &dSm_k[0], one_int, @@ -273,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); @@ -306,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; @@ -375,9 +379,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 = this->nlocal; + assert(nlocal>0); std::vector dm_VU_dSm(pv.nloc); @@ -400,9 +409,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 +447,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, @@ -483,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); @@ -527,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 e2c37039960..19c20a1c095 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,11 +33,11 @@ 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, - 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, @@ -50,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, @@ -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,11 +83,11 @@ 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, - 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, @@ -98,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, @@ -109,18 +110,18 @@ 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, - 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, @@ -128,7 +129,7 @@ void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) 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, @@ -139,18 +140,18 @@ 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, - 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, SR, one_int, one_int, this->paraV->desc, @@ -158,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_io.cpp b/source/source_lcao/module_dftu/dftu_io.cpp index d44113d1be9..9be899e244f 100644 --- a/source/source_lcao/module_dftu/dftu_io.cpp +++ b/source/source_lcao/module_dftu/dftu_io.cpp @@ -2,13 +2,20 @@ #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, + bool out_chg, + const std::string& global_out_dir, + int nspin, + int npol) { 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++) { @@ -20,28 +27,28 @@ 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) { - 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 { 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; - 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" + 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 << " 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; } } @@ -49,23 +56,25 @@ void Plus_U::output(const UnitCell &ucell) } } - GlobalV::ofs_running << "Local occupation matrices" << std::endl; - this->write_occup_m(ucell,GlobalV::ofs_running, true); - GlobalV::ofs_running << "//=======================================================//" << std::endl; - - //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"); - } + GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; + this->write_occup_m(ucell, GlobalV::ofs_running, true, nspin, npol); + + // Write dm_onsite.txt + if (out_chg && GlobalV::MY_RANK == 0) + { + std::ofstream ofdftu; + ofdftu.open(global_out_dir + "dm_onsite.txt"); + if (!ofdftu) + { + ModuleBase::WARNING_QUIT("Plus_U::output", "Can't create file dm_onsite.txt"); + } + 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); - ofdftu.close(); + + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + GlobalV::ofs_running << " | # END DFT+U INFO |" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl << std::endl; return; } @@ -74,59 +83,58 @@ 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, + int nspin, + int npol) { 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++) { const int iat = ucell.itia2iat(T, I); - ofs << "atoms" - << " " << iat << std::endl; 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" - << " " << l << std::endl; 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; + ofs << "\n Atom=" << iat+1; + 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++) { - 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,43 +146,43 @@ 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; + ofs << std::setprecision(8) << std::fixed; 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; } } - 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,60 +194,62 @@ 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 << " 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 < 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; + ofs << std::setw(12) << std::setprecision(8) << std::fixed + << locale[iat][l][n][0](m0_all, m1_all); + } + } + ofs << std::endl; + } } - ofs << std::setw(12) << std::setprecision(8) - << std::fixed<< "atomic mag: "< 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 dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); } 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; + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not find the file dm_onsite.txt. Please do scf calculation first"); } } - exit(0); + ModuleBase::WARNING_QUIT("Plus_U::read_occup_m", "Can not open dm_onsite.txt file"); } ifdftu.clear(); ifdftu.seekg(0); - char word[10]; + char word[20]; - int T=0; - int iat=0; - int spin=0; - int L=0; - int zeta=0; + int T = 0; + int iat = 0; + int spin = 0; + int L = 0; + int zeta = 0; ifdftu.rdstate(); while (ifdftu.good()) { ifdftu >> word; - if (ifdftu.eof()) - { - break; - } + if (ifdftu.eof()) + { + break; + } - if (strcmp("atoms", word) == 0) + if (strcmp("Atom=", word) == 0) { ifdftu >> iat; + iat -= 1; + ifdftu >> word; + + if (strcmp("L=", word) != 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) + { + 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'); T = ucell.iat2it[iat]; @@ -294,103 +316,66 @@ void Plus_U::read_occup_m(const UnitCell& ucell, for (int l = 0; l < NL; l++) { - if (l != get_orbital_corr(T)) - { - continue; - } - - ifdftu >> word; - - if (strcmp("L", word) == 0) + if (l != get_orbital_corr(T)) { - ifdftu >> L; - ifdftu.ignore(150, '\n'); + continue; + } - const int N = ucell.atoms[T].l_nchi[L]; - for (int n = 0; n < N; n++) + if (nspin == 1 || nspin == 2) + { + 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 (PARAM.inp.nspin == 1 || PARAM.inp.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 (PARAM.inp.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 < PARAM.globalv.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++) - { - 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 { - 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"); } } } - 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'); + } + } } } } 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(); @@ -404,16 +389,18 @@ 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"); 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,20 +409,20 @@ 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) + if (nspin == 1 || nspin == 2) { for (int spin = 0; spin < 2; spin++) { @@ -450,17 +437,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 @@ -482,15 +469,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 +490,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 +503,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; diff --git a/source/source_lcao/module_dftu/dftu_occup.cpp b/source/source_lcao/module_dftu/dftu_occup.cpp index 311861e0c75..13ba48c5385 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], @@ -289,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]; @@ -346,7 +323,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 +331,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 +396,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 +426,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 +435,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 +480,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 +537,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 +545,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 +579,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_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_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_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/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) { 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); } } diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index e6a5603546c..d81b71e0d61 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -121,7 +121,20 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- if (PARAM.inp.dft_plus_u) { - dftu.init(ucell, nullptr, kv.get_nks()); + const int nlocal_dftu = 0; + dftu.init(ucell, nullptr, + PARAM.globalv.npol, + 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; 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