Feature: periodic self-consistent continuum solvation (SCCS) Step1: SCCS basic function

source/CMakeLists.txt

@@ -751,6 +751,7 @@ target_link_libraries(
     io_basic
     io_input
     surchem
+    sccs
     neighbor_search
     domain_decomposition
     neighbor

source/Makefile.Objects

@@ -49,6 +49,7 @@ VPATH=./src_global:\
 ./source_hamilt:\
 ./source_hamilt/module_ewald:\
 ./source_hamilt/module_surchem:\
+./source_hamilt/module_sccs:\
 ./source_hamilt/module_xc:\
 ./source_hamilt/module_xc/kernels:\
 ./source_pw/module_proj:\
@@ -105,6 +106,7 @@ ${OBJS_PSI_INITIALIZER}\
 ${OBJS_PW}\
 ${OBJS_RELAXATION}\
 ${OBJS_SURCHEM}\
+${OBJS_SCCS}\
 ${OBJS_SYMMETRY}\
 ${OBJS_XC}\
 ${OBJS_IO}\
@@ -280,6 +282,7 @@ OBJS_ELECSTAT=elecstate.o\
     pcc_0d.o\
     pcc_2d.o\
     pot_pcc.o\
+    pot_sccs.o\
     elecstate_op.o\
     efield.o\
     gatefield.o\
@@ -623,6 +626,14 @@ OBJS_SURCHEM=surchem.o\
     minimize_cg.o\
     sol_force.o\
 
+OBJS_SCCS=sccs_cavity.o\
+    sccs_functional.o\
+    sccs_parameters.o\
+    sccs_ionic_charge.o\
+    sccs_ionic_force.o\
+    sccs_pw_coulomb.o\
+    sccs_response.o
+
 OBJS_SYMMETRY=symm_other.o\
     symm_rot_spin.o\
     symm_basic.o\
@@ -750,6 +761,7 @@ OBJS_IO=module_parameter/input_conv.o\
     module_parameter/read_inp_bse.o\
     module_parameter/read_inp_deepks.o\
     module_parameter/read_inp_model.o\
+    module_parameter/read_inp_sccs.o\
     module_parameter/read_inp_postproc.o\
     module_parameter/read_inp_exx_dftu.o\
     module_parameter/read_inp_dfpt.o\

source/source_esolver/esolver_double_xc.cpp

@@ -427,7 +427,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
     const FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h,
                         this->inp_->sc_mag_switch, this->inp_->device,
                         PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad,
-                        PARAM.globalv.gamma_only_pw};
+                        PARAM.globalv.gamma_only_pw, this->inp_->imp_sol};
 
     fsl.getForceStress(ucell,
                        this->get_vdw_result(),

source/source_esolver/esolver_ks_lcao.cpp

@@ -362,7 +362,8 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
                               PARAM.globalv.domag,
                               PARAM.globalv.domag_z,
                               this->inp_->gga_grad,
-                              PARAM.globalv.gamma_only_pw};
+                              PARAM.globalv.gamma_only_pw,
+                              this->inp_->imp_sol};
 
     fsl.getForceStress(ucell,
                        this->get_vdw_result(),

source/source_estate/CMakeLists.txt

@@ -21,6 +21,7 @@ list(APPEND objects
     module_pot/pot_xc.cpp
     module_pot/pot_local.cpp
     module_pot/pot_pcc.cpp
+    module_pot/pot_sccs.cpp
     module_pot/potential_new.cpp
     module_pot/potential_types.cpp
     module_pot/pot_surchem.cpp

source/source_estate/elecstate.h

@@ -157,7 +157,7 @@ class ElecState
     /// @param dft_plus_u 0 disables the DFT+U energy term
     /// @param assume_isolated "makov-payne" adds the isolated-system correction
     void cal_energies(const int type,
-                      const bool imp_sol,
+                      const int imp_sol,
                       const bool sc_mag_switch,
                       const int dft_plus_u,
                       const std::string& assume_isolated);

source/source_estate/elecstate_energy.cpp

@@ -316,7 +316,7 @@ void ElecState::cal_converged()
  * @param type: 2 means Kohn-Sham functional;
  */
 void ElecState::cal_energies(const int type,
-                             const bool imp_sol,
+                             const int imp_sol,
                              const bool sc_mag_switch,
                              const int dft_plus_u,
                              const std::string& assume_isolated)
@@ -331,11 +331,20 @@ void ElecState::cal_energies(const int type,
     this->f_en.gatefield = get_etot_gatefield();
 
     //! energy from implicit solvation model
-    if (imp_sol)
+    if (imp_sol == 1)
     {
         this->f_en.esol_el = get_solvent_model_Ael();
         this->f_en.esol_cav = get_solvent_model_Acav();
     }
+    else if (imp_sol == 2)
+    {
+        this->pot->get_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav);
+    }
+    else
+    {
+        this->f_en.esol_el = 0.0;
+        this->f_en.esol_cav = 0.0;
+    }
 
     //! spin constrained energy
     if (sc_mag_switch)

source/source_estate/module_pot/pot_base.h

@@ -28,7 +28,18 @@ class PotBase
 
     virtual void cal_fixed_v(double* vl_pseudo){}
 
-    virtual double get_energy() const { return 0.0; } 
+    virtual double get_energy() const { return 0.0; }
+
+    // Solvation contributions in Ry, added once after the band-energy potential subtraction.
+    virtual void get_solvation_energy(double& electrostatic, double& non_electrostatic) const
+    {
+        electrostatic = 0.0;
+        non_electrostatic = 0.0;
+    }
+    // Explicit ionic solvation force in Ry/Bohr, added to existing contributions.
+    virtual void add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {}
+    // Electronic electrostatic correction only; excludes cavity/surface derivatives.
+    virtual const std::vector<double>* solvent_electrostatic_potential() const { return nullptr; }
 
     bool fixed_mode = 0;
     bool dynamic_mode = 0;

source/source_estate/module_pot/pot_sccs.cpp

@@ -0,0 +1,157 @@
+#include "pot_sccs.h"
+
+#include "source_base/parallel_reduce.h"
+#include "source_base/timer.h"
+#include "source_base/tool_quit.h"
+#include "source_cell/cell_tools.h"
+#include "source_hamilt/module_sccs/sccs_functional.h"
+#include "source_hamilt/module_sccs/sccs_ionic_charge.h"
+#include "source_hamilt/module_sccs/sccs_ionic_force.h"
+#include "source_hamilt/module_sccs/sccs_response.h"
+#include "source_io/module_parameter/input_parameter.h"
+
+#include <cmath>
+#include <sstream>
+#include <utility>
+
+namespace elecstate
+{
+void make_sccs_config_from_input(const Input_para& input,
+                                 ModuleSccs::SccsConfig& config,
+                                 ModuleSccs::PolarizationSolverParameters& solver)
+{
+    const ModuleSccs::Preset preset = ModuleSccs::parse_preset(input.sccs_preset);
+    if (preset == ModuleSccs::Preset::Custom)
+    {
+        config = ModuleSccs::SccsConfig();
+        config.cavity.density_min = input.sccs_rho_min;
+        config.cavity.density_max = input.sccs_rho_max;
+        config.cavity.epsilon_bulk = input.sccs_epsilon;
+        config.surface_tension = ModuleSccs::dyn_per_cm_to_hartree_per_bohr2(input.sccs_gamma);
+        config.pressure = ModuleSccs::gpa_to_hartree_per_bohr3(input.sccs_pressure);
+    }
+    else
+    {
+        config = ModuleSccs::make_sccs_config(preset);
+    }
+    config.surface_regularization = input.sccs_surface_eta;
+    solver.max_iterations = input.sccs_maxiter;
+    solver.tolerance_rms = input.sccs_tol_rms;
+    solver.tolerance_max = input.sccs_tol_max;
+}
+
+std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count)
+{
+    double ionic_charge = 0.0;
+    for (int it = 0; it < cell.ntype; ++it)
+    {
+        ionic_charge += cell.atoms[it].ncpp.zv * cell.atoms[it].na;
+    }
+    const double net_charge = ionic_charge - electron_count;
+    if (config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= 1e-6) { return std::string(); }
+    return "charged SCCS with periodic boundaries: the periodic Poisson solver drops the G = 0 component "
+           "of the net charge, so the energy depends on the cell size.";
+}
+
+PotSccs::PotSccs(const ModulePW::PW_Basis* basis,
+                 const ModuleSccs::SccsConfig& config,
+                 const ModuleSccs::PolarizationSolverParameters& solver,
+                 double electron_count)
+    : config_(config), solver_(solver), electron_count_(electron_count)
+{
+    this->rho_basis_ = basis;
+    this->dynamic_mode = true;
+}
+
+void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential)
+{
+    ModuleBase::timer::start("PotSccs", "cal_v_eff");
+    const ModulePW::PW_Basis& basis = *this->rho_basis_;
+    const std::vector<unitcell::AtomData> atoms = unitcell::get_atom_data(cell->atoms, cell->ntype, cell->lat0);
+    std::vector<double> ions;
+    ModuleSccs::gaussian_ionic_density(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions);
+    std::vector<double> density(basis.nrxx, 0.0);
+    std::vector<double> solute_charge(basis.nrxx);
+    double grid_charges[2] = {0.0, 0.0}; // electrons, ions
+    for (int ir = 0; ir < basis.nrxx; ++ir)
+    {
+        for (int spin = 0; spin < charge->nspin; ++spin) { density[ir] += charge->rho[spin][ir]; }
+        solute_charge[ir] = ions[ir] - density[ir];
+        grid_charges[0] += density[ir];
+        grid_charges[1] += ions[ir];
+    }
+    // Pool-reduced, so every rank takes the same decision.
+    Parallel_Reduce::reduce_pool(grid_charges, 2);
+    const double dv = basis.omega / basis.nxyz;
+    const double electron_error = grid_charges[0] * dv - electron_count_;
+    double valence_charge = 0.0;
+    for (const unitcell::AtomData& atom : atoms) { valence_charge += atom.valence_charge; }
+    const double ionic_error = grid_charges[1] * dv - valence_charge;
+    if (std::abs(ionic_error) > 1e-6)
+    {
+        ModuleBase::WARNING_QUIT("PotSccs::cal_v_eff", "SCCS ionic density normalization does not match the valence charge");
+    }
+    if (std::abs(electron_error) > 1e-6)
+    {
+        std::ostringstream message;
+        message << "SCCS grid electron count differs from the expected value by " << electron_error
+                << " e; the solute charge uses the grid density";
+        ModuleBase::WARNING("PotSccs::cal_v_eff", message.str());
+    }
+    ModuleSccs::SccsResponse response;
+    ModuleSccs::solve_sccs_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis,
+                                    cell->tpiba, response);
+    ModuleSccs::FunctionalResult functional;
+    ModuleSccs::evaluate_functional(solute_charge, response, config_, basis, cell->tpiba, functional);
+    electrostatic_rydberg_ = 2.0 * functional.reaction_energy;
+    non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy);
+    electrostatic_potential_.resize(basis.nrxx);
+    for (int ir = 0; ir < basis.nrxx; ++ir)
+    {
+        electrostatic_potential_[ir] = -2.0 * functional.reaction_potential[ir];
+        const double value = 2.0 * functional.electron_potential[ir];
+        for (int spin = 0; spin < charge->nspin; ++spin) { potential(spin, ir) += value; }
+    }
+    restart_potential_ = std::move(response.restart_potential);
+    ModuleBase::timer::end("PotSccs", "cal_v_eff");
+}
+
+void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const
+{
+    ModuleBase::timer::start("PotSccs", "add_solvation_force");
+    const ModulePW::PW_Basis& basis = *this->rho_basis_;
+    const std::vector<unitcell::AtomData> atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0);
+    std::vector<double> reaction(basis.nrxx);
+    for (int ir = 0; ir < basis.nrxx; ++ir)
+    {
+        reaction[ir] = -0.5 * electrostatic_potential_[ir];
+    }
+    std::vector<ModuleBase::Vector3<double>> ionic_force;
+    ModuleSccs::gaussian_ionic_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread,
+                                     ionic_force);
+    for (int ia = 0; ia < cell.nat; ++ia)
+    {
+        for (int axis = 0; axis < 3; ++axis)
+        {
+            force(ia, axis) += 2.0 * ionic_force[ia][axis];
+        }
+    }
+    ModuleBase::timer::end("PotSccs", "add_solvation_force");
+}
+
+double PotSccs::get_energy() const
+{
+    return electrostatic_rydberg_ + non_electrostatic_rydberg_;
+}
+
+void PotSccs::get_solvation_energy(double& electrostatic, double& non_electrostatic) const
+{
+    electrostatic = electrostatic_rydberg_;
+    non_electrostatic = non_electrostatic_rydberg_;
+}
+
+const std::vector<double>* PotSccs::solvent_electrostatic_potential() const
+{
+    return &electrostatic_potential_;
+}
+} // namespace elecstate

source/source_estate/module_pot/pot_sccs.h

@@ -0,0 +1,49 @@
+#ifndef POT_SCCS_H
+#define POT_SCCS_H
+
+#include "pot_base.h"
+#include "source_hamilt/module_sccs/sccs_parameters.h"
+
+#include <string>
+
+struct Input_para;
+namespace elecstate
+{
+// Map INPUT values, already validated by the INPUT reader, onto the solver.
+void make_sccs_config_from_input(const Input_para& input,
+                                 ModuleSccs::SccsConfig& config,
+                                 ModuleSccs::PolarizationSolverParameters& solver);
+
+// Returns a warning (empty if none) for a charged cell in a dielectric: the
+// periodic Poisson solver drops the G = 0 component of the net charge, so the
+// energy depends on the cell size. electron_count is the resolved nelec (the
+// INPUT value, or the valence count when it is not set).
+std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count);
+
+class PotSccs : public PotBase
+{
+public:
+    // electron_count is the expected number of electrons (resolved nelec); a grid
+    // density that differs from it by more than 1e-6 e is reported in the
+    // warning log, and the solute charge uses the grid density.
+    PotSccs(const ModulePW::PW_Basis* basis,
+            const ModuleSccs::SccsConfig& config,
+            const ModuleSccs::PolarizationSolverParameters& solver,
+            double electron_count);
+    void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override;
+    double get_energy() const override;
+    void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override;
+    void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override;
+    const std::vector<double>* solvent_electrostatic_potential() const override;
+
+private:
+    const ModuleSccs::SccsConfig config_;
+    const ModuleSccs::PolarizationSolverParameters solver_;
+    const double electron_count_;
+    double electrostatic_rydberg_ = 0.0;
+    double non_electrostatic_rydberg_ = 0.0;
+    std::vector<double> electrostatic_potential_;
+    std::vector<double> restart_potential_;
+};
+} // namespace elecstate
+#endif