diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e2a9504c0e..3c3d488a8b 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -276,6 +276,7 @@ - [hubbard\_u](#hubbard_u) - [yukawa\_potential](#yukawa_potential) - [yukawa\_lambda](#yukawa_lambda) + - [uramping](#uramping) - [omc](#omc) - [onsite\_radius](#onsite_radius) - [vdW correction](#vdw-correction) @@ -2610,6 +2611,14 @@ These variables are used to control DFT+U correlated parameters - **Description**: The screen length of Yukawa potential. If left to default, the screen length will be calculated as an average of the entire system. It's better to stick to the default setting unless there is a very good reason. - **Default**: Calculated on the fly. +### uramping + +- **Type**: Real +- **Unit**: eV +- **Availability**: DFT+U calculations with `mixing_restart > 0`. +- **Description**: Once `uramping` > 0.15 eV. DFT+U calculations will start SCF with U = 0 eV, namely normal LDA/PBE calculations. Once SCF restarts when `drho0` and `mixing_dmr=1` to improve convergence. For case extremely hard to converge, you can use so-called U-Ramping method by setting a finite positive `uramping` with `mixing_restart>0` and `mixing_dmr=1`. + ## Smearing Thermal smearing is an efficient tool for accelerating SCF convergence by allowing fractional occupation of molecular orbitals near the band edge. It is important for metallic systems. diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index a44c050a1a..57083301b2 100755 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -801,7 +801,7 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix* DM) ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); // std::vector*> dmr = DM->get_DMR_vector(); - std::vector> dmr_save = DM->get_DMR_save(); + std::vector>& dmr_save = DM->get_DMR_save(); // //const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; double* dmr_in; @@ -900,7 +900,7 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix, doubl ModuleBase::timer::tick("Charge_Mixing", "mix_dmr"); // std::vector*> dmr = DM->get_DMR_vector(); - std::vector> dmr_save = DM->get_DMR_save(); + std::vector>& dmr_save = DM->get_DMR_save(); // //const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1; double* dmr_in; diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 3eb0127ec0..ee812b0138 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -101,7 +101,8 @@ class Charge_Mixing Base_Mixing::Mixing* get_mixing() const {return mixing;} // for mixing restart - int mixing_restart = 0; //which step to restart mixing during SCF + int mixing_restart_step = 0; //which step to restart mixing during SCF + int mixing_restart_count = 0; // the number of restart mixing during SCF. Do not set mixing_restart_count as bool since I want to keep some flexibility in the future private: diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 02b21e71e1..59db2849fd 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -13,6 +13,7 @@ //--------------Temporary---------------- #include "module_base/global_variable.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_hamilt_lcao/module_dftu/dftu.h" //--------------------------------------- #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" @@ -443,15 +444,23 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) // mixing will restart at this->p_chgmix->mixing_restart steps if (drho <= GlobalV::MIXING_RESTART && GlobalV::MIXING_RESTART > 0.0 - && this->p_chgmix->mixing_restart > iter) + && this->p_chgmix->mixing_restart_step > iter) { - this->p_chgmix->mixing_restart = iter + 1; + this->p_chgmix->mixing_restart_step = iter + 1; } // drho will be 0 at this->p_chgmix->mixing_restart step, which is not ground state + bool not_restart_step = !(iter==this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0); + // SCF will continue if U is not converged for uramping calculation + bool is_U_converged = true; + // to avoid unnecessary dependence on dft+u, refactor is needed +#ifdef __LCAO + if (GlobalV::dft_plus_u) is_U_converged = GlobalC::dftu.u_converged(); +#endif + // this->conv_elec = (drho < this->scf_thr - && !(iter==this->p_chgmix->mixing_restart - && GlobalV::MIXING_RESTART > 0.0)); + && not_restart_step + && is_U_converged); // If drho < hsolver_error in the first iter or drho < scf_thr, we do not change rho. if (drho < hsolver_error || this->conv_elec) @@ -466,7 +475,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) //----------charge mixing--------------- // mixing will restart after this->p_chgmix->mixing_restart steps if (GlobalV::MIXING_RESTART > 0 - && iter == this->p_chgmix->mixing_restart - 1) + && iter == this->p_chgmix->mixing_restart_step - 1) { // do not mix charge density } @@ -529,10 +538,10 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) // notice for restart if (GlobalV::MIXING_RESTART > 0 - && iter == this->p_chgmix->mixing_restart - 1 + && iter == this->p_chgmix->mixing_restart_step - 1 && iter != GlobalV::SCF_NMAX) { - std::cout<<"SCF restart after this step!"<::iter_init(const int istep, const int iter) if (iter == 1) { this->p_chgmix->init_mixing(); // init mixing - this->p_chgmix->mixing_restart = GlobalV::SCF_NMAX + 1; + this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; + this->p_chgmix->mixing_restart_count = 0; + // this output will be removed once the feeature is stable + if(GlobalC::dftu.uramping > 0.01) + { + std::cout << " U-Ramping! Current U = " ; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + } + std::cout << " eV " << std::endl; + } } // for mixing restart - if (iter == this->p_chgmix->mixing_restart + if (iter == this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); + this->p_chgmix->mixing_restart_count++; + if (GlobalV::dft_plus_u) + { + GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 + if(GlobalC::dftu.uramping > 0.01) + { + std::cout << " U-Ramping! Current U = " ; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + } + std::cout << " eV " << std::endl; + } + if(GlobalC::dftu.uramping > 0.01 + && !GlobalC::dftu.u_converged()) + { + this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; + } + } if (GlobalV::MIXING_DMR) // for mixing_dmr { // allocate memory for dmr_mdata @@ -687,7 +717,9 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) // save input rho this->pelec->charge->save_rho_before_sum_band(); // save density matrix for mixing - if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= this->p_chgmix->mixing_restart) + if (GlobalV::MIXING_RESTART > 0 + && GlobalV::MIXING_DMR + && this->p_chgmix->mixing_restart_count > 0) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); @@ -895,11 +927,13 @@ void ESolver_KS_LCAO::iter_finish(int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); - // mix density matrix - if (GlobalV::MIXING_RESTART > 0 && iter >= this->p_chgmix->mixing_restart && GlobalV::MIXING_DMR ) + // mix density matrix if mixing_restart + mixing_dmr + not first mixing_restart at every iter + if (GlobalV::MIXING_RESTART > 0 + && this->p_chgmix->mixing_restart_count > 0 + && GlobalV::MIXING_DMR) { elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec)->get_DM(); this->p_chgmix->mix_dmr(dm); } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 8e1054f6b0..d5e7b9a195 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -587,10 +587,10 @@ void ESolver_KS_PW::iter_init(const int istep, const int iter) if (iter == 1) { this->p_chgmix->init_mixing(); - this->p_chgmix->mixing_restart = GlobalV::SCF_NMAX + 1; + this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; } // for mixing restart - if (iter == this->p_chgmix->mixing_restart && GlobalV::MIXING_RESTART > 0.0) + if (iter == this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); } diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp index d69e755b90..627a1a0d86 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu.cpp @@ -355,6 +355,36 @@ void DFTU::cal_energy_correction(const int istep) return; } +void DFTU::uramping_update() +{ + // if uramping < 0.1, use the original U + if(this->uramping < 0.01) return; + // loop to change U + for(int i = 0; i < this->U0.size(); i++) + { + if (this->U[i] + this->uramping < this->U0[i] ) + { + this->U[i] += this->uramping; + } + else + { + this->U[i] = this->U0[i]; + } + } +} + +bool DFTU::u_converged() +{ + for(int i = 0; i < this->U0.size(); i++) + { + if (this->U[i] != this->U0[i]) + { + return false; + } + } + return true; +} + void DFTU::set_dmr(const elecstate::DensityMatrix, double>* dmr) { this->dm_in_dftu_cd = dmr; diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index e493f47205..29c19f3c51 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -44,9 +44,13 @@ class DFTU // calculate the energy correction void cal_energy_correction(const int istep); double get_energy(){return EU;} + void uramping_update(); // update U by uramping + bool u_converged(); // check if U is converged double* U; // U (Hubbard parameter U) + std::vector U0; // U0 (target Hubbard parameter U0) int* orbital_corr; // + double uramping; // increase U by uramping, default is -1.0 int omc; // occupation matrix control int mixing_dftu; //whether to mix locale diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 9588132f25..8e65755cad 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -557,6 +557,7 @@ void Input::Default(void) yukawa_potential = false; yukawa_lambda = -1.0; omc = 0; + uramping = -1.0; // -1.0 means no ramping //========================================================== // DFT+DMFT Xin Qu added on 2020-08 @@ -2111,14 +2112,12 @@ bool Input::Read(const std::string& fn) { read_bool(ifs, test_skip_ewald); } - //-------------- - //---------------------------------------------------------------------------------- - // Xin Qu added on 2020-10-29 for DFT+U - //---------------------------------------------------------------------------------- + //---------------------------------------------------------- else if (strcmp("dft_plus_u", word) == 0) { read_value(ifs, dft_plus_u); } + // ignore to avoid error else if (strcmp("yukawa_potential", word) == 0) ifs.ignore(150, '\n'); else if (strcmp("hubbard_u", word) == 0) @@ -2129,6 +2128,10 @@ bool Input::Read(const std::string& fn) ifs.ignore(150, '\n'); else if (strcmp("yukawa_lambda", word) == 0) ifs.ignore(150, '\n'); + else if (strcmp("uramping", word) == 0) + { + ifs.ignore(150, '\n'); + } //---------------------------------------------------------------------------------- // Xin Qu added on 2020-08 for DFT+DMFT //---------------------------------------------------------------------------------- @@ -2471,6 +2474,11 @@ bool Input::Read(const std::string& fn) { ifs >> yukawa_lambda; } + else if (strcmp("uramping", word) == 0) + { + read_value(ifs, uramping); + uramping /= ModuleBase::Ry_to_eV; + } else if (strcmp("hubbard_u", word) == 0) { for (int i = 0; i < ntype; i++) @@ -3633,6 +3641,7 @@ void Input::Bcast() //----------------------------------------------------------------------------------- Parallel_Common::bcast_int(dft_plus_u); Parallel_Common::bcast_bool(yukawa_potential); + Parallel_Common::bcast_double(uramping); Parallel_Common::bcast_int(omc); Parallel_Common::bcast_double(yukawa_lambda); if (GlobalV::MY_RANK != 0) diff --git a/source/module_io/input.h b/source/module_io/input.h index 91f5347c20..7929b8a682 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -492,6 +492,7 @@ class Input int omc; ///< whether turn on occupation matrix control method or not bool yukawa_potential; ///< default:false double yukawa_lambda; ///< default:-1.0, which means we calculate lambda + double uramping; ///< default:-1.0, which means we do not use U-Ramping method //========================================================== // DFT+DMFT Xin Qu added on 2021-08 diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index e18661e7f0..6b77f06b3f 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -413,6 +413,7 @@ void Input_Conv::Convert(void) GlobalC::dftu.Yukawa = INPUT.yukawa_potential; GlobalC::dftu.omc = INPUT.omc; GlobalC::dftu.orbital_corr = INPUT.orbital_corr; + GlobalC::dftu.uramping = INPUT.uramping; GlobalC::dftu.mixing_dftu = INPUT.mixing_dftu; if (INPUT.yukawa_potential && INPUT.hubbard_u == nullptr) { @@ -422,6 +423,11 @@ void Input_Conv::Convert(void) INPUT.hubbard_u = new double[GlobalC::ucell.ntype]; } GlobalC::dftu.U = INPUT.hubbard_u; + GlobalC::dftu.U0 = std::vector(INPUT.hubbard_u, INPUT.hubbard_u + GlobalC::ucell.ntype); + if (INPUT.uramping > 0.01) + { + ModuleBase::GlobalFunc::ZEROS(GlobalC::dftu.U, GlobalC::ucell.ntype); + } } GlobalV::onsite_radius = INPUT.onsite_radius; #endif diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index 4a661bf582..304a1bf253 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -317,6 +317,7 @@ TEST_F(InputParaTest, Bcast) EXPECT_EQ(INPUT.dft_plus_u, 0); EXPECT_FALSE(INPUT.yukawa_potential); EXPECT_DOUBLE_EQ(INPUT.yukawa_lambda, -1.0); + EXPECT_DOUBLE_EQ(INPUT.uramping, -1.0); EXPECT_EQ(INPUT.omc, 0); EXPECT_FALSE(INPUT.dft_plus_dmft); EXPECT_FALSE(INPUT.rpa); diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index 64c48d0a9c..ccdb67cbc2 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -839,6 +839,7 @@ TEST_F(write_input, DFTU20) testing::HasSubstr( "dft_plus_u 0 #1/2:new/old DFT+U correction method; 0: standard DFT calcullation(default)")); EXPECT_THAT(output, testing::HasSubstr("yukawa_lambda -1 #default:0.0")); + EXPECT_THAT(output, testing::HasSubstr("uramping -1 #increasing U values during SCF")); EXPECT_THAT(output, testing::HasSubstr("yukawa_potential 0 #default: false")); EXPECT_THAT(output, testing::HasSubstr("omc 0 #the mode of occupation matrix control")); EXPECT_THAT(output, testing::HasSubstr("hubbard_u 0 #Hubbard Coulomb interaction parameter U(ev)")); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 5974a54003..ff4e523ada 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -454,6 +454,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs, "dft_plus_u", dft_plus_u, "1/2:new/old DFT+U correction method; 0: standard DFT calcullation(default)"); ModuleBase::GlobalFunc::OUTP(ofs, "yukawa_lambda", yukawa_lambda, "default:0.0"); ModuleBase::GlobalFunc::OUTP(ofs, "yukawa_potential", yukawa_potential, "default: false"); + ModuleBase::GlobalFunc::OUTP(ofs, "uramping", uramping, "increasing U values during SCF"); ModuleBase::GlobalFunc::OUTP(ofs, "omc", omc, "the mode of occupation matrix control"); ModuleBase::GlobalFunc::OUTP(ofs, "onsite_radius", onsite_radius, "radius of the sphere for onsite projection (Bohr)"); ofs << std::setw(20) << "hubbard_u "; diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/INPUT b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/INPUT new file mode 100644 index 0000000000..3ede777db9 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/INPUT @@ -0,0 +1,35 @@ +INPUT_PARAMETERS +suffix autotest +nbands 40 + +calculation scf +ecutwfc 20 +scf_thr 1.0e-4 +scf_nmax 500 +out_chg 0 + +smearing_method gaussian +smearing_sigma 0.01 + +cal_force 1 +cal_stress 1 + +mixing_type broyden +mixing_beta 0.4 +mixing_restart 5e-3 +mixing_dmr 1 +mixing_ndim 15 + +ks_solver genelpa +basis_type lcao +gamma_only 0 +symmetry 0 +nspin 2 + +#Parameter DFT+U +dft_plus_u 1 +orbital_corr 2 -1 +hubbard_u 5.0 0.0 +uramping 2.5 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/KPT b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/README b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/README new file mode 100644 index 0000000000..d5afa236a8 --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/README @@ -0,0 +1,7 @@ +FeO Anti-Ferromagnetic+LDA +U-Ramping method to find the ground state +NSPIN 2 +lcao +init magnet: up + down +SG15 + diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/STRU b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/STRU new file mode 100644 index 0000000000..65d8f7c98f --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/STRU @@ -0,0 +1,29 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf +O 1.000 O_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_4s2p2d1f.orb +O_gga_7au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 +0.51 0.51 0.51 mag -1.0 + +O +0.0 +2 +0.25 0.25 0.25 1 1 1 +0.75 0.75 0.75 1 1 1 diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/jd b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/jd new file mode 100644 index 0000000000..e407617cfa --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/jd @@ -0,0 +1 @@ +DFTU + U-Ramping + NSPIN2 AFM, Fe2O2, multi-k case diff --git a/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/result.ref b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/result.ref new file mode 100644 index 0000000000..407a2fed4e --- /dev/null +++ b/tests/integrate/260_NO_DJ_PK_PU_AFM_URAMPING/result.ref @@ -0,0 +1,5 @@ +etotref -7645.9759390680846991 +etotperatomref -1911.4939847670 +totalforceref 5.988604 +totalstressref 2367.340815 +totaltimeref 29.88 diff --git a/tests/integrate/CASES b/tests/integrate/CASES index f82151cab6..24d8025003 100644 --- a/tests/integrate/CASES +++ b/tests/integrate/CASES @@ -182,6 +182,7 @@ 260_NO_DJ_PK_PU_FM 260_NO_DJ_PK_PU_SO 260_NO_DJ_PK_PU_S1 +260_NO_DJ_PK_PU_AFM_URAMPING 270_NO_MD_1O 270_NO_MD_2O 281_NO_KP_HSE