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

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`.
+- **Default**: -1.0.
+
 ### omc
 
 - **Type**: Integer

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -796,6 +796,13 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
 
 void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<double, double>* DM)
 {
+    if (GlobalV::MIXING_RESTART <= 0 
+        || this->mixing_restart_count <= 0 
+        || !GlobalV::MIXING_DMR )
+    {
+        return;
+    }
+
     // Notice that DensityMatrix object is a Template class
     ModuleBase::TITLE("Charge_Mixing", "mix_dmr");
     ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");
@@ -895,6 +902,13 @@ void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<double, double>* DM)
 
 void Charge_Mixing::mix_dmr(elecstate::DensityMatrix<std::complex<double>, double>* DM)
 {
+    if (GlobalV::MIXING_RESTART <= 0 
+        || this->mixing_restart_count <= 0 
+        || !GlobalV::MIXING_DMR )
+    {
+        return;
+    }
+
     // Notice that DensityMatrix object is a Template class
     ModuleBase::TITLE("Charge_Mixing", "mix_dmr");
     ModuleBase::timer::tick("Charge_Mixing", "mix_dmr");

source/module_elecstate/module_charge/charge_mixing.h

@@ -102,6 +102,7 @@ class Charge_Mixing
 
     // for mixing restart
     int mixing_restart = 0; //which step to restart mixing during SCF
+    int mixing_restart_count = 0; // the number of restart mixing during SCF
 
   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"
@@ -449,9 +450,17 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 				}
 
 				// 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 && 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)
@@ -532,7 +541,7 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
              && iter == this->p_chgmix->mixing_restart - 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

@@ -553,12 +553,42 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
     {
         this->p_chgmix->init_mixing(); // init mixing
         this->p_chgmix->mixing_restart = 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 
         && 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 = 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();
@@ -896,12 +928,9 @@ 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 )
-    {
-        elecstate::DensityMatrix<TK, double>* dm
-                    = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
-        this->p_chgmix->mix_dmr(dm);
-    }
+    elecstate::DensityMatrix<TK, double>* dm
+                = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
+    this->p_chgmix->mix_dmr(dm);
 
     //-----------------------------------
     // save charge density

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