Feature: periodic self-consistent continuum solvation (SCCS) Step3: lowpass and core fill

source/Makefile.Objects

@@ -635,6 +635,7 @@ OBJS_SURCHEM=surchem.o\
 OBJS_SCCS=sccs_cavity.o\
     sccs_cav_deriv.o\
     sccs_functional.o\
+    sccs_lowpass.o\
     sccs_parameters.o\
     sccs_ion_charge.o\
     sccs_ion_force.o\
@@ -773,6 +774,7 @@ OBJS_IO=module_parameter/input_conv.o\
     module_parameter/read_inp_deepks.o\
     module_parameter/read_inp_model.o\
     module_parameter/read_inp_sccs.o\
+    module_parameter/sccs_inp_check.o\
     module_parameter/read_inp_postproc.o\
     module_parameter/read_inp_exx_dftu.o\
     module_parameter/read_inp_dfpt.o\

source/source_cell/cell_tools.cpp

@@ -3,9 +3,31 @@
  * @brief Implementation of cell tool free functions.
  */
 #include "cell_tools.h"
+#include "source_base/element_name.h"
+
+#include <algorithm>
 
 namespace unitcell
 {
+    int pseudo_atom_z(const Atom& atom)
+    {
+        const std::string& label = atom.ncpp.psd;
+        const std::size_t first = label.find_first_not_of(" \t");
+        if (first == std::string::npos)
+        {
+            return 0;
+        }
+        const std::size_t last = label.find_last_not_of(" \t");
+        const std::size_t symbol_length = last - first + 1;
+        const std::string symbol = label.substr(first, symbol_length);
+        const std::vector<std::string>::const_iterator element = std::find(ModuleBase::element_name.begin(), ModuleBase::element_name.end(), symbol);
+        if (element == ModuleBase::element_name.end())
+        {
+            return 0;
+        }
+        return element - ModuleBase::element_name.begin() + 1;
+    }
+
     std::vector<AtomData> get_atom_data(const Atom* atoms, const int ntype, const double lattice_scale)
     {
         int nat = 0;
@@ -17,12 +39,14 @@ namespace unitcell
         data.reserve(nat);
         for (int it = 0; it < ntype; ++it)
         {
+            const int atomic_number = pseudo_atom_z(atoms[it]);
             for (int ia = 0; ia < atoms[it].na; ++ia)
             {
                 AtomData atom;
                 atom.position = atoms[it].tau[ia] * lattice_scale;
                 atom.mass = atoms[it].mass;
                 atom.valence_charge = atoms[it].ncpp.zv;
+                atom.atomic_number = atomic_number;
                 data.push_back(atom);
             }
         }

source/source_cell/cell_tools.h

@@ -22,8 +22,13 @@ namespace unitcell
         ModuleBase::Vector3<double> position; ///< Bohr
         double mass = 0.0;
         double valence_charge = 0.0; ///< Positive ionic charge in units of e
+        int atomic_number = 0; ///< Nuclear charge from the pseudopotential element; 0 if unknown
     };
 
+    /// Atomic number of the pseudopotential element of one atom type; 0 if the
+    /// element symbol is not recognized.
+    int pseudo_atom_z(const Atom& atom);
+
     /// Extract positions, masses and pseudopotential valence charges.
     /// Requires initialized Atom storage and a positive lattice scale in Bohr.
     std::vector<AtomData> get_atom_data(const Atom* atoms, int ntype, double lattice_scale);

source/source_cell/test/test_cell_tools.cpp

@@ -128,3 +128,43 @@ TEST(CellTools, AcceptsNoAtomTypes)
     const std::vector<unitcell::AtomData> data = unitcell::get_atom_data(nullptr, 0, 10.0);
     EXPECT_TRUE(data.empty());
 }
+
+TEST(CellTools, AtomDataCarriesPseudopotentialAtomicNumber)
+{
+    Atom atoms[2];
+    atoms[0].na = 2;
+    atoms[0].ncpp.psd = " C ";
+    atoms[0].tau = {ModuleBase::Vector3<double>(), ModuleBase::Vector3<double>()};
+    atoms[1].na = 1;
+    atoms[1].ncpp.psd = "H";
+    atoms[1].tau = {ModuleBase::Vector3<double>()};
+
+    const std::vector<unitcell::AtomData> data = unitcell::get_atom_data(atoms, 2, 10.0);
+    ASSERT_EQ(data.size(), 3u);
+    EXPECT_EQ(data[0].atomic_number, 6);
+    EXPECT_EQ(data[1].atomic_number, 6);
+    EXPECT_EQ(data[2].atomic_number, 1);
+}
+
+TEST(CellTools, UnknownPseudopotentialElementDoesNotGuessFromAtomLabel)
+{
+    Atom atom;
+    atom.na = 1;
+    atom.label = "H";
+    atom.ncpp.psd = "unknown";
+    atom.tau = {ModuleBase::Vector3<double>()};
+    const std::vector<unitcell::AtomData> data = unitcell::get_atom_data(&atom, 1, 1.0);
+    ASSERT_EQ(data.size(), 1u);
+    EXPECT_EQ(data[0].atomic_number, 0);
+}
+
+TEST(CellTools, PseudoAtomicNumberTrimsTheElementSymbol)
+{
+    Atom atom;
+    atom.ncpp.psd = " O ";
+    EXPECT_EQ(unitcell::pseudo_atom_z(atom), 8);
+    atom.ncpp.psd = "";
+    EXPECT_EQ(unitcell::pseudo_atom_z(atom), 0);
+    atom.ncpp.psd = "Oxygen";
+    EXPECT_EQ(unitcell::pseudo_atom_z(atom), 0);
+}

source/source_estate/module_pot/pot_sccs.cpp

@@ -44,7 +44,11 @@ void sccs_cfg_input(const Input_para& input,
     {
         config = ModuleSccs::preset_config(preset);
     }
+    config.core_electrons = input.sccs_solvent_mode == "full";
+    config.core_spreads = input.sccs_corespread;
     config.surface_regularization = input.sccs_surface_eta;
+    config.cavity.lowpass_p1 = input.sccs_lowpass_p1;
+    config.cavity.lowpass_p2 = input.sccs_lowpass_p2;
     if (input.assume_isolated == "pcc_0d")
     {
         config.boundary = ModuleSccs::Boundary::Pcc0d;
@@ -63,6 +67,38 @@ void sccs_cfg_input(const Input_para& input,
     solver.tolerance_max = input.sccs_tol_max;
 }
 
+std::string check_sccs_cell(const ModuleSccs::SccsConfig& config, const UnitCell& cell)
+{
+    if (!config.core_electrons)
+    {
+        return std::string();
+    }
+    const std::size_t width_count = config.core_spreads.size();
+    const std::size_t atom_count = cell.nat;
+    if (width_count != 1 && width_count != atom_count)
+    {
+        ModuleBase::WARNING_QUIT("check_sccs_cell", "sccs_corespread must contain one value or exactly nat values");
+    }
+    if (width_count != 1)
+    {
+        return std::string();
+    }
+    std::ostringstream message;
+    for (int it = 0; it < cell.ntype; ++it)
+    {
+        const int atomic_number = unitcell::pseudo_atom_z(cell.atoms[it]);
+        if (atomic_number != 0)
+        {
+            continue;
+        }
+        message << "sccs_solvent_mode full: the element of atom type " << cell.atoms[it].label
+                << " (pseudopotential element '" << cell.atoms[it].ncpp.psd
+                << "') is unknown, so the single sccs_corespread value is applied to it even if its"
+                << " pseudopotential has no core electrons; give one value per atom to control it.\n";
+    }
+    return message.str();
+}
+
 std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count)
 {
     double ionic_charge = 0.0;
@@ -140,6 +176,15 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::
     {
         coulomb.reset(new ModuleSccs::PeriodicCoulombOperator(basis, cell->tpiba));
     }
+    if (config_.core_electrons)
+    {
+        std::vector<double> core_density;
+        ModuleSccs::gauss_core_rho(atoms, basis, cell->tpiba, config_.core_spreads, core_density);
+        for (int ir = 0; ir < basis.nrxx; ++ir)
+        {
+            density[ir] += core_density[ir];
+        }
+    }
     ModuleSccs::SccsResponse response;
     ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis,
                                cell->tpiba, *coulomb, response);
@@ -154,6 +199,7 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::
         const double value = 2.0 * functional.electron_potential[ir];
         for (int spin = 0; spin < charge->nspin; ++spin) { potential(spin, ir) += value; }
     }
+    cavity_potential_ = std::move(functional.cavity_potential);
     restart_potential_ = std::move(response.restart_potential);
     ModuleBase::timer::end("PotSccs", "cal_v_eff");
 }
@@ -171,6 +217,16 @@ void PotSccs::add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) con
     std::vector<ModuleBase::Vector3<double>> ionic_force;
     ModuleSccs::gauss_ion_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread,
                                 ionic_force);
+    if (config_.core_electrons)
+    {
+        std::vector<ModuleBase::Vector3<double>> core_force;
+        ModuleSccs::gauss_core_frc(atoms, cavity_potential_, basis, cell.tpiba, config_.core_spreads,
+                                   core_force);
+        for (int ia = 0; ia < cell.nat; ++ia)
+        {
+            ionic_force[ia] += core_force[ia];
+        }
+    }
     for (int ia = 0; ia < cell.nat; ++ia)
     {
         for (int axis = 0; axis < 3; ++axis)

source/source_estate/module_pot/pot_sccs.h

@@ -14,6 +14,12 @@ void sccs_cfg_input(const Input_para& input,
                     ModuleSccs::SccsConfig& config,
                     ModuleSccs::PolarizationSolverParameters& solver);
 
+// Structure-dependent checks that the INPUT reader cannot make: a
+// sccs_corespread count other than one or nat stops the run. Returns a warning
+// (empty if none) when a single width meets an unknown pseudopotential
+// element, whose automatic core exclusion is then not possible.
+std::string check_sccs_cell(const ModuleSccs::SccsConfig& config, const UnitCell& cell);
+
 // 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
@@ -44,6 +50,7 @@ class PotSccs : public PotBase
     double non_electrostatic_rydberg_ = 0.0;
     std::vector<double> electrostatic_potential_;
     std::vector<double> restart_potential_;
+    std::vector<double> cavity_potential_;
 };
 } // namespace elecstate
 #endif

source/source_estate/module_pot/potential_types.cpp

@@ -29,8 +29,16 @@ PotBase* Potential::make_sccs_pot(const Input_para& input)
     sccs_cfg_input(input, config, solver);
     if (this->ucell_ != nullptr)
     {
-        const std::string warning = check_charge(config, *this->ucell_, input.nelec);
-        if (!warning.empty()) { ModuleBase::WARNING("Potential::make_sccs_pot", warning); }
+        const std::string structure_warning = check_sccs_cell(config, *this->ucell_);
+        if (!structure_warning.empty())
+        {
+            ModuleBase::WARNING("Potential::make_sccs_pot", structure_warning);
+        }
+        const std::string charge_warning = check_charge(config, *this->ucell_, input.nelec);
+        if (!charge_warning.empty())
+        {
+            ModuleBase::WARNING("Potential::make_sccs_pot", charge_warning);
+        }
     }
     return new PotSccs(this->rho_basis_, config, solver, input.nelec);
 }

source/source_estate/test/test_pot_sccs.cpp

@@ -187,6 +187,30 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig)
     }
 }
 
+TEST_F(PotSccsTest, InputMapsCavityModelKeys)
+{
+    Input_para input;
+    ModuleSccs::SccsConfig config;
+    ModuleSccs::PolarizationSolverParameters solver;
+    elecstate::sccs_cfg_input(input, config, solver);
+    EXPECT_DOUBLE_EQ(config.cavity.lowpass_p1, -1.0);
+    EXPECT_DOUBLE_EQ(config.cavity.lowpass_p2, -1.0);
+    EXPECT_FALSE(config.core_electrons);
+    input.sccs_preset = "water-neutral";
+    input.assume_isolated = "pcc_0d";
+    input.sccs_lowpass_p1 = 10.0;
+    input.sccs_lowpass_p2 = 5.0;
+    input.sccs_solvent_mode = "full";
+    input.sccs_corespread = {0.8, 0.0, -1.0};
+    elecstate::sccs_cfg_input(input, config, solver);
+    EXPECT_DOUBLE_EQ(config.cavity.lowpass_p1, 10.0);
+    EXPECT_DOUBLE_EQ(config.cavity.lowpass_p2, 5.0);
+    EXPECT_TRUE(config.core_electrons);
+    EXPECT_EQ(config.core_spreads, input.sccs_corespread);
+    // The cavity model keys leave the preset dielectric unchanged.
+    EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3);
+}
+
 TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce)
 {
     Input_para input;
@@ -199,6 +223,57 @@ TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce)
     check_solvation_force(config, solver);
 }
 
+// The cavity terms of core Gaussians, with surface tension and pressure.
+TEST_F(PotSccsTest, FullCavityForceMatchesTheEnergyDerivative)
+{
+    Input_para input;
+    input.sccs_epsilon = 1.0;
+    input.sccs_tol_rms = 1e-13;
+    input.sccs_tol_max = 1e-12;
+    ModuleSccs::SccsConfig config;
+    ModuleSccs::PolarizationSolverParameters solver;
+    elecstate::sccs_cfg_input(input, config, solver);
+    config.core_electrons = true;
+    config.core_spreads = {2.0};
+    config.surface_tension = 1e-5;
+    config.pressure = 1e-6;
+    check_solvation_force(config, solver);
+}
+
+TEST_F(PotSccsTest, FullCavityStructureChecksWarnOnceAndRejectWidthCount)
+{
+    UnitCell cell;
+    Atom atom;
+    cell.ntype = 1;
+    cell.nat = 2;
+    cell.atoms = &atom;
+    atom.na = 2;
+    atom.label = "X";
+    atom.ncpp.psd = "Xx";
+    ModuleSccs::SccsConfig config;
+    config.core_electrons = true;
+    config.core_spreads = {0.5};
+    auto warnings = [&]() { return elecstate::check_sccs_cell(config, cell); };
+    const std::string unknown = warnings();
+    EXPECT_NE(unknown.find("unknown"), std::string::npos);
+    atom.ncpp.psd = "O";
+    EXPECT_EQ(warnings(), "");
+    atom.ncpp.psd = "Xx";
+    config.core_spreads = {0.5, 0.6};
+    EXPECT_EQ(warnings(), "");
+    config.core_electrons = false;
+    config.core_spreads = {0.5};
+    EXPECT_EQ(warnings(), "");
+    config.core_electrons = true;
+    config.core_spreads = {0.5, 0.6, 0.7};
+    if (SccsTest::pool_size == 1)
+    {
+        testing::internal::CaptureStdout();
+        EXPECT_EXIT(elecstate::check_sccs_cell(config, cell), ::testing::ExitedWithCode(1), "");
+        testing::internal::GetCapturedStdout();
+    }
+}
+
 TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns)
 {
     UnitCell cell;

source/source_hamilt/module_sccs/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(sccs OBJECT
   sccs_parameters.cpp
   sccs_functional.cpp
   sccs_ion_charge.cpp
+  sccs_lowpass.cpp
   sccs_ion_force.cpp
   sccs_coul_0d.cpp
   sccs_coul_2d.cpp

source/source_hamilt/module_sccs/sccs_cav_deriv.cpp

@@ -1,5 +1,6 @@
 #include "sccs_cav_deriv.h"
 #include "sccs_cavity.h"
+#include "sccs_lowpass.h"
 #include "sccs_response.h"
 
 #include "source_base/constants.h"
@@ -45,8 +46,12 @@ void make_cav_derivs(const std::vector<double>& density,
                      double tpiba,
                      bool open_boundary,
                      SccsResponse& response,
-                     std::vector<double>& coefficient)
+                     CavityDerivatives& derivatives)
 {
+    derivatives.filter = switch_filter(cavity, basis);
+    const bool lowpass = uses_lowpass(cavity);
+    derivatives.gradient.clear();
+    std::vector<double>& coefficient = derivatives.coefficient;
     const std::size_t size = density.size();
     response.solute.resize(size);
     response.dsolute_drho.resize(size);
@@ -72,6 +77,13 @@ void make_cav_derivs(const std::vector<double>& density,
         // grid, since its chain coefficient is not converged at the cavity edge.
         std::vector<std::complex<double>> boundary_g(basis.npw);
         basis.real2recip(response.solute.data(), boundary_g.data());
+        if (lowpass)
+        {
+            for (int ig = 0; ig < basis.npw; ++ig)
+            {
+                boundary_g[ig] *= derivatives.filter[ig];
+            }
+        }
         XC_Functional::grad_rho(boundary_g.data(), gradient.data(), &basis, tpiba);
         XC_Functional::laplacian_rho(boundary_g.data(), laplacian.data(), &basis, tpiba);
     }
@@ -94,5 +106,10 @@ void make_cav_derivs(const std::vector<double>& density,
         }
     }
     eps_of_boundary(cavity, gradient, laplacian, response, coefficient);
+    if (open_boundary && lowpass)
+    {
+        derivatives.gradient.swap(gradient);
+    }
 }
+
 } // namespace ModuleSccs