Feature: periodic self-consistent continuum solvation (SCCS) Step4: delayed start, relax, diagnostics and out_sol

source/Makefile.Objects

@@ -288,6 +288,10 @@ OBJS_ELECSTAT=elecstate.o\
     pcc_2d.o\
     pot_pcc.o\
     pot_sccs.o\
+    potential_sccs.o\
+    sccs_fixed_src.o\
+    sccs_output.o\
+    sccs_collect.o\
     sccs_pcc_0d.o\
     sccs_pcc_2d.o\
     elecstate_op.o\
@@ -651,6 +655,7 @@ OBJS_SURCHEM=surchem.o\
 
 OBJS_SCCS=sccs_cavity.o\
     sccs_cav_deriv.o\
+    sccs_diag.o\
     sccs_functional.o\
     sccs_lowpass.o\
     sccs_parameters.o\

source/source_esolver/esolver_ks.cpp

@@ -314,9 +314,12 @@ bool ESolver_KS::defer_iter_output(const int istep) const
 void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver)
 {
     const int iter_limit = this->scf_iter_limit(istep);
+    const bool two_fermi = PARAM.globalv.two_fermi;
+    const bool model_changed
+        = module_charge::update_pot_model(iter, this->drho, this->kv, *this->pelec->pot, this->p_chgmix);
 
     // 1.1) print out band gap 
-    if (!PARAM.globalv.two_fermi)
+    if (!two_fermi)
     {
         this->pelec->cal_bandgap();
     }
@@ -336,7 +339,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
 
     // 2.1) compute magnetization, only for spin==2
     ucell.magnet.compute_mag(ucell.omega, this->chr.nrxx, this->chr.nxyz, this->chr.rho,
-                                       this->inp_->nspin, PARAM.globalv.two_fermi, this->inp_->nelec,
+                                       this->inp_->nspin, two_fermi, this->inp_->nelec,
                                        this->pelec->nelec_spin.data());
 
     // 2.2) charge mixing
@@ -354,6 +357,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
     ctx.iter_limit = iter_limit;
     ctx.hsolver_error = hsolver_error;
     ctx.scf_thr = this->scf_thr;
+    ctx.model_changed = model_changed;
     ctx.scf_ene_thr = this->scf_ene_thr;
     ctx.converged_u = converged_u;
     ctx.ks_run = PARAM.globalv.ks_run;
@@ -435,7 +439,9 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &
 
     // print energies
     elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, iter_limit, drho,
-    dkin, duration, *this->inp_, PARAM.globalv.two_fermi, diag_ethr, 0, true, this->ds_rms_);
+    dkin, duration, *this->inp_, two_fermi, diag_ethr, 0, true, this->ds_rms_);
+
+    module_charge::print_corr_iter(this->kv, this->inp_->sccs_debug, this->drho, *this->pelec->pot);
 
 
     if (!this->defer_iter_output(istep))

source/source_estate/CMakeLists.txt

@@ -22,6 +22,10 @@ list(APPEND objects
     module_pot/pot_local.cpp
     module_pot/pot_pcc.cpp
     module_pot/pot_sccs.cpp
+    module_pot/potential_sccs.cpp
+    module_pot/sccs_fixed_src.cpp
+    module_pot/sccs_output.cpp
+    module_pot/sccs_collect.cpp
     module_pot/sccs_pcc_0d.cpp
     module_pot/sccs_pcc_2d.cpp
     module_pot/potential_new.cpp

source/source_estate/module_charge/chg_mix.cpp

@@ -194,6 +194,10 @@ void Charge_Mixing::mix_reset()
 {
     this->mixing->reset();
     this->rho_mdata.reset();
+    // mix_dmr and mix_uom reuse the Pulay coefficients of the rho history, so
+    // the density-matrix and DFT+U occupation histories must restart with it.
+    this->dmr_mdata.reset();
+    this->uom_mdata.reset();
     // initailize tau_mdata
     if (cfg_.mixing_tau)
     {
@@ -209,6 +213,12 @@ void Charge_Mixing::reset_scf_history(int iter_limit)
     this->mixing_restart_count = 0;
 }
 
+void Charge_Mixing::restart_for_new_model(const int iter)
+{
+    this->mix_reset();
+    this->mixing_restart_last = iter;
+}
+
 bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho,
                                      const int iternum_used, const double threshold)
 {

source/source_estate/module_charge/chg_mix.h

@@ -77,7 +77,15 @@ class Charge_Mixing
 
     /** @brief Reset restart counters and residual history for an independently limited SCF trial. */
     void reset_scf_history(int iter_limit);
-
+
+    /**
+     * @brief restart mixing after the potential model changed during the SCF (SCCS delayed start)
+     *
+     * The drho of earlier iterations belongs to the old model, so the window of
+     * if_scf_oscillate also restarts at iteration iter, as after a mixing restart.
+     */
+    void restart_for_new_model(const int iter);
+
     // extracting parameters normally these parameters will not be used outside charge mixing
     // while Exx is using them as well as some other places
     const std::string& get_mixing_mode() const {return cfg_.mixing_mode;}

source/source_estate/module_charge/chg_routine.cpp

@@ -2,8 +2,70 @@
 #include "source_estate/module_charge/chg_dmr.h"
 
 #include "source_base/parallel_comm.h"
+#include "source_base/parallel_common.h"
 #include "source_pw/module_pwdft/dftu_base.h" // Plus_U_Base members used below
 #include "source_estate/update_pot.h"
+#include "source_cell/klist.h"
+
+#include <iostream>
+
+namespace
+{
+// SCCS reports and the start notice come from the first rank of the first pool.
+bool reports_here(const K_Vectors& kv)
+{
+    return kv.para_k.my_pool == 0 && kv.para_k.rank_in_pool == 0;
+}
+
+// sccs_debug reports go to the running log of that rank.
+std::ostream* corr_log(const K_Vectors& kv, const int sccs_debug)
+{
+    if (sccs_debug > 0 && reports_here(kv))
+    {
+        return &GlobalV::ofs_running;
+    }
+    return nullptr;
+}
+} // namespace
+
+bool module_charge::update_pot_model(const int iter,
+        const double drho,
+        const K_Vectors& kv,
+        elecstate::Potential& pot,
+        Charge_Mixing* p_chgmix)
+{
+    // Only while a potential component waits to change its model: ordinary
+    // runs skip the broadcast. The state is synchronized, so every rank takes
+    // the same branch.
+    if (!pot.model_pending())
+    {
+        return false;
+    }
+    // Band-parallel groups may differ slightly in drho until mixing broadcasts
+    // it; broadcast first so every rank changes the model in the same iteration.
+    double density_residual = drho;
+    Parallel_Common::bcast_double(density_residual);
+    const bool reports = reports_here(kv);
+    std::ostream* const notice = reports ? &std::cout : nullptr;
+    const bool model_changed = pot.update_model(iter, density_residual, notice);
+    if (model_changed)
+    {
+        p_chgmix->restart_for_new_model(iter);
+    }
+    return model_changed;
+}
+
+void module_charge::print_corr_iter(const K_Vectors& kv,
+        const int sccs_debug,
+        const double drho,
+        const elecstate::Potential& pot)
+{
+    std::ostream* const log = corr_log(kv, sccs_debug);
+    if (log != nullptr)
+    {
+        pot.write_corr_iter(*log, sccs_debug, drho);
+    }
+}
 
 void module_charge::chgmixing_ks(const int iter,
         UnitCell& ucell,
@@ -42,7 +104,7 @@ void module_charge::chgmixing_ks(const int iter,
         // not ground state
         bool not_restart_step = !(iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0);
 
-        conv_esolver = (drho < scf_thr && not_restart_step && converged_u);
+        conv_esolver = (drho < scf_thr && not_restart_step && converged_u && !ctx.model_changed);
 
         // add energy threshold for SCF convergence
         if (scf_ene_thr > 0.0)

source/source_estate/module_charge/chg_routine.h

@@ -10,6 +10,7 @@
 
 // Plus_U_Base forward declaration, full definition in source_pw/module_pwdft/dftu_base.h
 class Plus_U_Base;
+class K_Vectors;
 
 namespace module_charge
 {
@@ -23,11 +24,28 @@ struct ScfMixingCtx
     double scf_ene_thr;    ///< energy convergence threshold
     bool converged_u;      ///< whether DFT+U has converged
     bool ks_run;           ///< whether the current run is a KS calculation (PARAM.globalv.ks_run)
+    bool model_changed = false; ///< a potential component changed its model in this iteration (SCCS delayed start); the SCF is not converged
     double drho;            ///< charge density deviation (in/out)
     bool oscillate_esolver; ///< whether esolver oscillates (out)
     bool conv_esolver;      ///< whether esolver converged (out)
 };
 
+/// Lets a potential component change its model before mixing (SCCS delayed
+/// start); true only in the iteration in which it does, which then restarts the
+/// mixing. The first rank of the first pool announces the change on the screen.
+bool update_pot_model(const int iter,
+        const double drho,
+        const K_Vectors& kv,
+        elecstate::Potential& pot,
+        Charge_Mixing* p_chgmix);
+
+/// sccs_debug report of the SCCS correction in this SCF iteration, in the
+/// running log of the first rank of the first pool.
+void print_corr_iter(const K_Vectors& kv,
+        const int sccs_debug,
+        const double drho,
+        const elecstate::Potential& pot);
+
 void chgmixing_ks(const int iter,
         UnitCell& ucell,
         elecstate::ElecState* pelec,

source/source_estate/module_charge/unittests/test_chg_mix.cpp

@@ -1,6 +1,7 @@
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 #include "../chg_mix.h"
+#include "../chg_dmr.h"
 #include "../chg_drho.h"
 #include "../chg_drho_detail.h"
 #include "../chg_precond.h"
@@ -193,6 +194,97 @@ class ChargeMixingTest : public ::testing::Test
     }
 };
 
+// A model change (SCCS delayed start) raises drho; the oscillation check must
+// not compare it with the drho of the old model.
+TEST_F(ChargeMixingTest, ModelChangeRestartsTheOscillationWindow)
+{
+    const int scf_nmax = 20;
+    cfg.scf_nmax = scf_nmax;
+    Charge_Mixing plain;
+    Charge_Mixing restarted;
+    plain.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba);
+    restarted.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba);
+    plain.init_mixing();
+    restarted.init_mixing();
+    const int scf_os_ndim = 3;
+    const double scf_os_thr = -0.05;
+    const int model_change = 6;
+    // Steady decrease, a jump when the new model starts, then a steady decrease.
+    const double drho[9] = {1e-2, 3e-3, 9e-4, 2.7e-4, 8.1e-5, 2.43e-5, 1e-2, 5e-3, 2.5e-3};
+    for (int iter = 1; iter <= 9; ++iter)
+    {
+        if (iter == model_change)
+        {
+            restarted.restart_for_new_model(iter);
+        }
+        const double residual = drho[iter - 1];
+        const bool plain_oscillates = plain.if_scf_oscillate(iter, residual, scf_os_ndim, scf_os_thr);
+        const bool restarted_oscillates = restarted.if_scf_oscillate(iter, residual, scf_os_ndim, scf_os_thr);
+        EXPECT_FALSE(restarted_oscillates) << "iteration " << iter;
+        if (iter == model_change + 1)
+        {
+            EXPECT_TRUE(plain_oscillates);
+        }
+    }
+    EXPECT_EQ(restarted.mixing_restart_last, model_change);
+}
+
+// mix_uom reuses the coefficients of the rho history, so after a model change
+// (SCCS delayed start) the first occupation update must be a single simple
+// mixing step, not a combination of the old occupations.
+TEST_F(ChargeMixingTest, ModelChangeRestartsUomHistory)
+{
+    cfg.mixing_mode = "pulay";
+    Charge_Mixing charge_mixing;
+    charge_mixing.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba);
+    charge_mixing.init_mixing();
+    const int uom_size = 4;
+    charge_mixing.allocate_mixing_uom(uom_size);
+    for (int step = 0; step < 2; ++step)
+    {
+        const double uom_value = 1.0 + step;
+        const double uom_save_value = 0.5 + step;
+        std::vector<double> uom(uom_size, uom_value);
+        std::vector<double> uom_save(uom_size, uom_save_value);
+        charge_mixing.mix_uom(uom, uom_save);
+    }
+
+    const int model_change = 3;
+    charge_mixing.restart_for_new_model(model_change);
+    std::vector<double> uom = {2.0, 3.0, 4.0, 5.0};
+    std::vector<double> uom_save = {1.0, 1.0, 1.0, 1.0};
+    charge_mixing.mix_uom(uom, uom_save);
+
+    const double beta = cfg.mixing_beta;
+    for (int i = 0; i < uom_size; ++i)
+    {
+        const double expected = 1.0 + beta * (2.0 + i - 1.0);
+        EXPECT_NEAR(uom[i], expected, 1e-12) << "element " << i;
+    }
+}
+
+TEST_F(ChargeMixingTest, MixResetClearsDmrHistory)
+{
+    cfg.mixing_dmr = true;
+    cfg.scf_thr_type = 2;
+    Charge_Mixing charge_mixing;
+    charge_mixing.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba);
+    charge_mixing.init_mixing();
+
+    Base_Mixing::Mixing_Data& dmr_mdata = charge_mixing.get_dmr_mdata();
+    module_charge::init_mixing_dmr(charge_mixing.get_mixing(), dmr_mdata, 4,
+                                   charge_mixing.get_mixing_config());
+    const double dmr[4] = {1.0, 2.0, 3.0, 4.0};
+    dmr_mdata.push(dmr);
+    ASSERT_EQ(dmr_mdata.ndim_history, 1);
+
+    charge_mixing.mix_reset();
+
+    EXPECT_EQ(dmr_mdata.ndim_use, 0);
+    EXPECT_EQ(dmr_mdata.ndim_history, 0);
+    EXPECT_EQ(dmr_mdata.start, -1);
+}
+
 TEST_F(ChargeMixingTest, SetMixingTest)
 {
 #ifdef _OPENMP

source/source_estate/module_pot/pot_base.h

@@ -9,6 +9,8 @@
 
 namespace elecstate
 {
+struct SolventGridField;
+
 /** This class is the base class of Potential module
  1. Main class Potential is derived from it.
  2. components of potentials on real space grids can derived from it and will be registered into Potential.
@@ -40,6 +42,8 @@ class PotBase
     virtual void add_sol_force(const UnitCell&, ModuleBase::matrix&) const {}
     // Electronic electrostatic correction only; excludes cavity/surface derivatives.
     virtual const std::vector<double>* sol_elec_pot() const { return nullptr; }
+    // Implicit-solvent fields of the last correction for out_sol; none by default.
+    virtual void add_sol_fields(std::vector<SolventGridField>&) const {}
 
     bool fixed_mode = 0;
     bool dynamic_mode = 0;