Fix(dft+u): stabilize occupation-matrix mixing and frozen potentials

Author: liu-687

source/source_base/module_mixing/broyden_mixing.cpp

@@ -4,6 +4,9 @@
 #include "source_base/module_container/base/third_party/blas.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
+
+#include <algorithm>
+
 namespace Base_Mixing
 {
 template void Broyden_Mixing::tem_push_data(Mixing_Data& mdata,
@@ -170,16 +173,29 @@ void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function<double
 
 		if (info != 0)
 		{
-			ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYSV.");
+			// A rank-deficient history can occur when successive SCF residuals
+			// become linearly dependent.  Aborting loses a perfectly usable
+			// latest fixed-point step; fall back to that step and let the next
+			// iteration rebuild the history.
+			ModuleBase::WARNING("Broyden_Mixing",
+			                    "DSYSV failed in Broyden mixing; using the latest linear step.");
+			std::fill(coef.begin(), coef.end(), 0.0);
+			coef[mdata.start] = 1.0;
+			// Discard the singular secant history before the next iteration.
+			ndim_cal_dF = 0;
+			start_dF = -1;
 		}
+		else
+		{
 
-        // after solving, gamma store the coeficients for mixing
-        coef[mdata.start] = 1 + gamma[dFindex_move(0)];
-        for (int i = 1; i < ndim_cal_dF; ++i)
-        {
-            coef[mdata.index_move(-i)] = gamma[dFindex_move(-i)] - gamma[dFindex_move(-i + 1)];
+            // after solving, gamma store the coeficients for mixing
+            coef[mdata.start] = 1 + gamma[dFindex_move(0)];
+            for (int i = 1; i < ndim_cal_dF; ++i)
+            {
+                coef[mdata.index_move(-i)] = gamma[dFindex_move(-i)] - gamma[dFindex_move(-i + 1)];
+            }
+            coef[mdata.index_move(-ndim_cal_dF)] = -gamma[dFindex_move(-ndim_cal_dF + 1)];
         }
-        coef[mdata.index_move(-ndim_cal_dF)] = -gamma[dFindex_move(-ndim_cal_dF + 1)];
 
         delete[] work;
         delete[] iwork;

source/source_base/module_mixing/test/mixing_test.cpp

@@ -198,6 +198,33 @@ TEST_F(Mixing_Test, BroydenSolveLinearEq)
     clear();
 }
 
+TEST_F(Mixing_Test, BroydenFallsBackForSingularHistory)
+{
+    init_method("broyden");
+
+    Base_Mixing::Mixing_Data mdata;
+    this->mixing->init_mixing_data(mdata, 2, sizeof(double));
+
+    // A rank-deficient residual history makes the symmetric linear system
+    // passed to DSYSV singular.  The mixer should retain the latest linear
+    // step instead of aborting the SCF calculation.
+    const auto singular_inner_product = [](double*, double*) { return 1.0; };
+    const std::vector<double> data_in = {0.0, 0.0};
+    std::vector<double> data_out(2);
+    for (int scale = 1; scale <= 3; ++scale)
+    {
+        data_out[0] = static_cast<double>(scale);
+        data_out[1] = static_cast<double>(2 * scale);
+        this->mixing->push_data(mdata, data_in.data(), data_out.data(), nullptr, true);
+        this->mixing->cal_coef(mdata, singular_inner_product);
+    }
+
+    std::vector<double> mixed(2);
+    this->mixing->mix_data(mdata, mixed.data());
+    EXPECT_DOUBLE_EQ(mixed[0], 1.8);
+    EXPECT_DOUBLE_EQ(mixed[1], 3.6);
+}
+
 TEST_F(Mixing_Test, PulaySolveLinearEq)
 {
 #ifdef _OPENMP
@@ -306,4 +333,4 @@ TEST_F(Mixing_Test, OtherCover)
     EXPECT_EQ(nodata.length, 0);
 
     clear();
-}
+}

source/source_estate/module_charge/charge_mixing.cpp

@@ -1,15 +1,19 @@
 #include "charge_mixing.h"
 
+#include <functional>
+
 #include "source_io/module_parameter/parameter.h"
 #include "source_base/module_mixing/broyden_mixing.h"
 #include "source_base/module_mixing/pulay_mixing.h"
 #include "source_base/parallel_common.h"
+#include "source_base/parallel_reduce.h"
 #include "source_base/timer.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 
 Charge_Mixing::Charge_Mixing()
 {
     this->mixing = nullptr;
+    this->mixing_uom = nullptr;
     this->mixing_highf = nullptr;
 }
 
@@ -26,6 +30,12 @@ Charge_Mixing::~Charge_Mixing()
 		delete this->mixing_highf;
         this->mixing_highf = nullptr;
 	}
+
+    if (this->mixing_uom != nullptr)
+    {
+        delete this->mixing_uom;
+        this->mixing_uom = nullptr;
+    }
 }
 
 void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
@@ -112,7 +122,9 @@ void Charge_Mixing::init_mixing()
     ModuleBase::TITLE("Charge_Mixing", "init_mixing");
     ModuleBase::timer::start("Charge_Mixing", "init_mixing");
 
-    // (re)construct mixing object
+    // (re)construct the charge mixer.  The independent UOM mixer is
+    // initialized explicitly by the DFT+U call path.
+
     if (this->mixing_mode == "broyden")
     {
         delete this->mixing;
@@ -133,6 +145,9 @@ void Charge_Mixing::init_mixing()
         ModuleBase::WARNING_QUIT("Charge_Mixing", "This Mixing mode is not implemended yet,coming soon.");
     }
 
+    delete this->mixing_uom;
+    this->mixing_uom = nullptr;
+
     if ( PARAM.globalv.double_grid)
     {
         // ONLY smooth part of charge density is mixed by specific mixing method
@@ -190,6 +205,20 @@ void Charge_Mixing::init_mixing()
     return;
 }
 
+void Charge_Mixing::init_mixing_uom()
+{
+    ModuleBase::TITLE("Charge_Mixing", "init_mixing_uom");
+    ModuleBase::timer::start("Charge_Mixing", "init_mixing_uom");
+
+    delete this->mixing_uom;
+    // The occupation matrix is a small, strongly coupled nonlinear variable.
+    // A conservative plain step avoids binding its history to the much
+    // larger charge-density mixer.
+    this->mixing_uom = new Base_Mixing::Plain_Mixing(0.5 * this->mixing_beta);
+
+    ModuleBase::timer::end("Charge_Mixing", "init_mixing_uom");
+}
+
 void Charge_Mixing::set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in)
 {
     this->rhopw = rhopw_in;
@@ -205,6 +234,11 @@ void Charge_Mixing::mix_reset()
     {
         this->tau_mdata.reset();
     }
+    if (this->mixing_uom != nullptr)
+    {
+        this->mixing_uom->reset();
+        this->uom_mdata.reset();
+    }
 }
 
 bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold)
@@ -265,8 +299,12 @@ void Charge_Mixing::allocate_mixing_uom(int uom_size)
     ModuleBase::timer::start("Charge_Mixing", "allocate_mixing_uom");
     // For nspin=2, uom_size already includes both spin channels
     // (pot_uterm_pw.size() = pot_index * 2 for nspin=2)
-    // So uom_fold should always be 1
-    this->mixing->init_mixing_data(this->uom_mdata, uom_size, sizeof(double));
+    // So uom_fold should always be 1.  UOM has an independent history.
+    if (this->mixing_uom == nullptr)
+    {
+        ModuleBase::WARNING_QUIT("Charge_Mixing", "UOM mixing object is not initialized.");
+    }
+    this->mixing_uom->init_mixing_data(this->uom_mdata, uom_size, sizeof(double));
     this->uom_mdata.reset();
     ModuleBase::timer::end("Charge_Mixing", "allocate_mixing_uom");
     return;
@@ -276,16 +314,53 @@ void Charge_Mixing::mix_uom(std::vector<double>& uom_in, std::vector<double>& uo
 {
     ModuleBase::TITLE("Charge_Mixing", "mix_uom");
     ModuleBase::timer::start("Charge_Mixing", "mix_uom");
+
+    if (uom_in.empty())
+    {
+        ModuleBase::timer::end("Charge_Mixing", "mix_uom");
+        return;
+    }
+    if (this->mixing_uom == nullptr || this->uom_mdata.length != uom_in.size()
+        || uom_save_in.size() != uom_in.size())
+    {
+        ModuleBase::WARNING_QUIT("Charge_Mixing", "UOM mixing data is not initialized consistently.");
+    }
+
     double* uom_value_out = uom_in.data();
     double* uom_value_in = uom_save_in.data();
-    // For all nspin cases, uom_array layout is already fully sized
-    // and mixing operates on the entire array
-    this->mixing->push_data(this->uom_mdata, uom_value_in, uom_value_out, nullptr, false);
-    this->mixing->mix_data(this->uom_mdata, uom_value_out);
+    const std::size_t length = this->uom_mdata.length;
+
+    // The UOM vector has its own history and its own contractive step.  It
+    // must not reuse the charge-density mixer: the two vectors have different
+    // lengths and represent different nonlinear variables.
+    this->mixing_uom->push_data(this->uom_mdata,
+                                uom_value_in,
+                                uom_value_out,
+                                nullptr,
+                                true);
+
+    auto inner_product_uom = [length](double* first, double* second) {
+        double value = 0.0;
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+ : value)
+#endif
+        for (std::size_t i = 0; i < length; ++i)
+        {
+            value += first[i] * second[i];
+        }
+#ifdef __MPI
+        Parallel_Reduce::reduce_pool(value);
+#endif
+        return value;
+    };
+
+    this->mixing_uom->cal_coef(this->uom_mdata, inner_product_uom);
+    this->mixing_uom->mix_data(this->uom_mdata, uom_value_out);
+
     ModuleBase::timer::end("Charge_Mixing", "mix_uom");
 #ifdef __MPI
     // Synchronize mixed uom across all ranks to prevent divergence
-    // after multiple Pulay steps (same pattern as mix_dmr)
+    // after multiple mixing steps.
     Parallel_Common::bcast_double(uom_in.data(), uom_in.size());
 #endif
     return;

source/source_estate/module_charge/charge_mixing.h

@@ -57,6 +57,11 @@ class Charge_Mixing
      */
     void init_mixing();
 
+    /**
+     * @brief initialize the independent DFT+U occupation-matrix mixer
+     */
+    void init_mixing_uom();
+
     /**
      * @brief allocate memory of dmr_mdata
      * @param nnr size of real-space density matrix
@@ -128,6 +133,7 @@ class Charge_Mixing
 
     // mixing_data
     Base_Mixing::Mixing* mixing = nullptr; ///< Mixing object to mix charge density, kinetic energy density and compensation density
+    Base_Mixing::Mixing* mixing_uom = nullptr; ///< Independent mixer for the DFT+U occupation matrix
     Base_Mixing::Mixing_Data rho_mdata;    ///< Mixing data for charge density
     Base_Mixing::Mixing_Data tau_mdata;    ///< Mixing data for kinetic energy density
     Base_Mixing::Mixing_Data nhat_mdata;   ///< Mixing data for compensation density

source/source_estate/module_charge/chgmixing.cpp

@@ -130,6 +130,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
         p_chgmix->mixing_restart_step = inp.scf_nmax + 1;
         if (inp.dft_plus_u && inp.mixing_dftu)
         {
+            p_chgmix->init_mixing_uom();
             // enable mixing_dftu for DFT+U occupation mixing
             dftu.enable_mixing();
             // allocate memory for uom_mdata
@@ -143,6 +144,13 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number
         p_chgmix->init_mixing();
         p_chgmix->mixing_restart_count++;
 
+        if (inp.dft_plus_u && inp.mixing_dftu)
+        {
+            p_chgmix->init_mixing_uom();
+            dftu.enable_mixing();
+            p_chgmix->allocate_mixing_uom(dftu.get_size_pot_uterm_pw());
+        }
+
         if (inp.dft_plus_u)
         {
             if (dftu.get_uramping() > 0.01 && !dftu.u_converged())

source/source_estate/test/charge_mixing_test.cpp

@@ -285,6 +285,36 @@ TEST_F(ChargeMixingTest, InitMixingTest)
     EXPECT_EQ(CMtest.rho_mdata.length, 2 * pw_basis.nrxx);
 }
 
+TEST_F(ChargeMixingTest, DftuOccupationMatrixUsesIndependentMixer)
+{
+    Charge_Mixing CMtest;
+    CMtest.set_rhopw(&pw_basis, &pw_basis);
+    CMtest.set_mixing("broyden",
+                      0.8,
+                      3,
+                      0.0,
+                      false,
+                      0.0,
+                      0.0,
+                      0.1,
+                      0.0,
+                      false,
+                      ucell.omega,
+                      ucell.tpiba);
+    CMtest.init_mixing();
+    CMtest.init_mixing_uom();
+
+    EXPECT_NE(CMtest.mixing, CMtest.mixing_uom);
+    CMtest.allocate_mixing_uom(2);
+    std::vector<double> uom = {1.0, 2.0};
+    std::vector<double> uom_save = {0.0, 0.0};
+    CMtest.mix_uom(uom, uom_save);
+
+    // UOM uses 0.5 * mixing_beta, independently of the charge mixer.
+    EXPECT_DOUBLE_EQ(uom[0], 0.4);
+    EXPECT_DOUBLE_EQ(uom[1], 0.8);
+}
+
 TEST_F(ChargeMixingTest, InnerDotRealTest)
 {
     Charge_Mixing CMtest;

source/source_lcao/module_dftu/dftu_lcao_op.cpp

@@ -440,7 +440,8 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
     // Purpose: Ensure occ_mat is recomputed from updated DMR in next SCF iteration,
     // rather than using stale pre-read data from file.
     // TODO: This logic is confusing. Consider explicit variable like `is_last_spin_channel`.
-    if (this->current_spin == this->nspin - 1 || this->nspin == 4) 
+    if ((this->current_spin == this->nspin - 1 || this->nspin == 4)
+        && this->dftu->get_occ_mat_ctrl() != 2)
     {
         this->dftu->mark_occ_mat_dirty();
     }