Feature: add so-called U-Ramping method to improve DFTU convergence (#…

…3890)

* add a new input parameter uramping, default is -1.0

* add a new member uramping in dftu.h and value it in Input_conv::Convert()

* add a new member U0 in dftu.h and value it in Input_conv::Convert()

* set U zero if uramping > 0.1 in Input_conv::Convert()

* add a new member function DFTU::uramping_update() to change U during SCF

* add a new member function DFTU::U_converged() to check if U=U0

* add uramping in esolver_lcao::iter_init() and change conv_elec

* make some minor changes

* fix bug which only restart once

* add a new member mixing_restart_count in Charge_Mixing

* add U output at iter 1

* fix build error without LCAO

* add some output of U value during U-Raming

* add some docs about uramping

* modify under some comments

* add a space in output

* add docs abouht DFTU in converge.md

* add recommendations of mixing_restart in input-main.md

* add some comment to explain why mixing_restart_count is int instead of bool

* add ref & in mix_dmr()

* add a new CI case 260_NO_DJ_PK_PU_AFM_URAMPING

* rename chgmix->mixing_restart as mixing_restart_step to avoid misunderstanding between GlobalV::MIXING_RESTART

* keep verbosity of judge in esolver to remind others mixing_dmr is only called in this case.

---------

Co-authored-by: Mohan Chen <mohan.chen.chen.mohan@gmail.com>

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 `drho<mixing_restart`, U value will increase by `uramping` eV. SCF will repeat above calcuations until U values reach target defined in `hubbard_u`. As for `uramping=1.0 eV`, the recommendations of `mixing_restart` is around `5e-4`.  
+- **Default**: -1.0.
+
 ### omc
 
 - **Type**: Integer

docs/advanced/scf/converge.md

@@ -16,6 +16,8 @@ For magnetic calculations, `mixing_beta_mag` and `mixing_gg0_mag` are activated.
 
 An example showcasing different charge mixing methods can be found in our [repository](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/charge_mixing/pw_Al). Four INPUT files are provided, with description given in README.
 
+As for DFT+U calculations, where the hamiltonian is not only dependent on charge density, but also dependent on density matrix. You can try `mixing_restart>0` 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.

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -801,7 +801,7 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<double, double>* DM)
     ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
     //
     std::vector<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
-    std::vector<std::vector<double>> dmr_save = DM->get_DMR_save();
+    std::vector<std::vector<double>>& 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<std::complex<double>, doubl
     ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
     //
     std::vector<hamilt::HContainer<double>*> dmr = DM->get_DMR_vector();
-    std::vector<std::vector<double>> dmr_save = DM->get_DMR_save();
+    std::vector<std::vector<double>>& dmr_save = DM->get_DMR_save();
     //
     //const int dmr_nspin = (GlobalV::NSPIN == 2) ? 2 : 1;
     double* dmr_in;

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:
 

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<T, Device>::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<T, Device>::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<T, Device>::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!"<<std::endl;
+				std::cout<<" SCF restart after this step!"<<std::endl;
 			}
 		}
 #ifdef __RAPIDJSON

source/module_esolver/esolver_ks_lcao.cpp

@@ -552,13 +552,43 @@ void ESolver_KS_LCAO<TK, TR>::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<TK, TR>::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<TK, double>* dm
             = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
@@ -895,11 +927,13 @@ void ESolver_KS_LCAO<TK, TR>::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<TK, double>* dm
-                    = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
+            = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
         this->p_chgmix->mix_dmr(dm);
     }
 

source/module_esolver/esolver_ks_pw.cpp

@@ -587,10 +587,10 @@ void ESolver_KS_PW<T, Device>::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();
     }

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<std::complex<double>, double>* dmr)
 {
     this->dm_in_dftu_cd = dmr;

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<double> 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
 

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)

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

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<double>(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

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);

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)"));

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 ";

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