From 61fc660721e0a046c96df327b92f576c19e5d16e Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Tue, 20 Oct 2026 15:11:13 +0800 Subject: [PATCH 1/4] Feature: SCCS derivative lowpass filter and full (core-electron) cavity --- source/Makefile.Objects | 2 + source/source_cell/cell_tools.cpp | 24 ++ source/source_cell/cell_tools.h | 5 + source/source_cell/test/test_cell_tools.cpp | 28 +++ source/source_estate/module_pot/pot_sccs.cpp | 56 +++++ source/source_estate/module_pot/pot_sccs.h | 7 + .../module_pot/potential_types.cpp | 12 +- source/source_estate/test/test_pot_sccs.cpp | 63 ++++- .../source_hamilt/module_sccs/CMakeLists.txt | 1 + .../module_sccs/sccs_cav_deriv.cpp | 19 +- .../module_sccs/sccs_cav_deriv.h | 11 +- .../source_hamilt/module_sccs/sccs_cavity.cpp | 5 + .../source_hamilt/module_sccs/sccs_cavity.h | 5 + .../module_sccs/sccs_functional.cpp | 2 + .../module_sccs/sccs_functional.h | 1 + .../module_sccs/sccs_ion_charge.cpp | 53 ++++- .../module_sccs/sccs_ion_charge.h | 26 ++ .../module_sccs/sccs_ion_force.cpp | 26 ++ .../module_sccs/sccs_ion_force.h | 17 ++ .../module_sccs/sccs_lowpass.cpp | 81 +++++++ .../source_hamilt/module_sccs/sccs_lowpass.h | 32 +++ .../module_sccs/sccs_parameters.h | 3 + .../module_sccs/sccs_response.cpp | 42 ++-- .../source_hamilt/module_sccs/sccs_response.h | 6 +- .../module_sccs/test/CMakeLists.txt | 3 +- .../module_sccs/test/test_sccs_cav_deriv.cpp | 27 ++- .../module_sccs/test/test_sccs_coul_0d.cpp | 5 +- .../module_sccs/test/test_sccs_coul_2d.cpp | 5 +- .../module_sccs/test/test_sccs_ion_charge.cpp | 59 +++++ .../module_sccs/test/test_sccs_ion_force.cpp | 133 +++++++++++ .../module_sccs/test/test_sccs_lowpass.cpp | 146 ++++++++++++ source/source_io/CMakeLists.txt | 1 + .../module_parameter/input_parameter.h | 4 + .../module_parameter/read_inp_model.cpp | 3 +- .../module_parameter/read_inp_sccs.cpp | 187 +++++---------- .../module_parameter/read_inp_sccs.h | 4 + .../module_parameter/sccs_inp_check.cpp | 224 ++++++++++++++++++ source/source_io/test_serial/CMakeLists.txt | 1 + .../test_serial/test_read_inp_sccs.cpp | 100 ++++++++ tests/01_PW/CASES_CPU.txt | 1 + tests/01_PW/scf_sccs_full/INPUT | 22 ++ tests/01_PW/scf_sccs_full/KPT | 4 + tests/01_PW/scf_sccs_full/README | 1 + tests/01_PW/scf_sccs_full/STRU | 24 ++ tests/01_PW/scf_sccs_full/result.ref | 6 + tests/01_PW/scf_sccs_full/threshold | 3 + tests/01_PW/scf_sccs_pcc_0d/INPUT | 2 + tests/01_PW/scf_sccs_pcc_0d/README | 2 +- tests/01_PW/scf_sccs_pcc_0d/result.ref | 14 +- 49 files changed, 1337 insertions(+), 171 deletions(-) create mode 100644 source/source_hamilt/module_sccs/sccs_lowpass.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_lowpass.h create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp create mode 100644 source/source_io/module_parameter/sccs_inp_check.cpp create mode 100644 tests/01_PW/scf_sccs_full/INPUT create mode 100644 tests/01_PW/scf_sccs_full/KPT create mode 100644 tests/01_PW/scf_sccs_full/README create mode 100644 tests/01_PW/scf_sccs_full/STRU create mode 100644 tests/01_PW/scf_sccs_full/result.ref create mode 100644 tests/01_PW/scf_sccs_full/threshold diff --git a/source/Makefile.Objects b/source/Makefile.Objects index ced2b16b823..b2361a72e26 100644 --- a/source/Makefile.Objects +++ b/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\ diff --git a/source/source_cell/cell_tools.cpp b/source/source_cell/cell_tools.cpp index 624f04145fd..74d26a63f99 100644 --- a/source/source_cell/cell_tools.cpp +++ b/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 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::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 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); } } diff --git a/source/source_cell/cell_tools.h b/source/source_cell/cell_tools.h index 037d152a282..09e7fbe8059 100644 --- a/source/source_cell/cell_tools.h +++ b/source/source_cell/cell_tools.h @@ -22,8 +22,13 @@ namespace unitcell ModuleBase::Vector3 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 get_atom_data(const Atom* atoms, int ntype, double lattice_scale); diff --git a/source/source_cell/test/test_cell_tools.cpp b/source/source_cell/test/test_cell_tools.cpp index f6f3cbde4dd..bd724de9dfd 100644 --- a/source/source_cell/test/test_cell_tools.cpp +++ b/source/source_cell/test/test_cell_tools.cpp @@ -96,11 +96,13 @@ TEST(CellTools, ExtractsPositionsMassesAndValenceChargesInAtomOrder) atoms[0].na = 2; atoms[0].mass = 12.0; atoms[0].ncpp.zv = 4.0; + atoms[0].ncpp.psd = " C "; atoms[0].tau = {ModuleBase::Vector3(0.1, 0.2, 0.3), ModuleBase::Vector3(0.4, 0.5, 0.6)}; atoms[1].na = 1; atoms[1].mass = 1.0; atoms[1].ncpp.zv = 1.0; + atoms[1].ncpp.psd = "H"; atoms[1].tau = {ModuleBase::Vector3(0.7, 0.8, 0.9)}; const std::vector data = unitcell::get_atom_data(atoms, 2, 10.0); @@ -120,6 +122,9 @@ TEST(CellTools, ExtractsPositionsMassesAndValenceChargesInAtomOrder) EXPECT_DOUBLE_EQ(data[0].valence_charge, 4.0); EXPECT_DOUBLE_EQ(data[1].valence_charge, 4.0); EXPECT_DOUBLE_EQ(data[2].valence_charge, 1.0); + EXPECT_EQ(data[0].atomic_number, 6); + EXPECT_EQ(data[1].atomic_number, 6); + EXPECT_EQ(data[2].atomic_number, 1); EXPECT_DOUBLE_EQ(atoms[0].tau[0].x, 0.1); } @@ -128,3 +133,26 @@ TEST(CellTools, AcceptsNoAtomTypes) const std::vector data = unitcell::get_atom_data(nullptr, 0, 10.0); EXPECT_TRUE(data.empty()); } + +TEST(CellTools, UnknownPseudopotentialElementDoesNotGuessFromAtomLabel) +{ + Atom atom; + atom.na = 1; + atom.label = "H"; + atom.ncpp.psd = "unknown"; + atom.tau = {ModuleBase::Vector3()}; + const std::vector 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); +} diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp index cfdf978761b..c021964caf5 100644 --- a/source/source_estate/module_pot/pot_sccs.cpp +++ b/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 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> ionic_force; ModuleSccs::gauss_ion_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, ionic_force); + if (config_.core_electrons) + { + std::vector> 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) diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h index ab68bb85010..77962f5c54a 100644 --- a/source/source_estate/module_pot/pot_sccs.h +++ b/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 electrostatic_potential_; std::vector restart_potential_; + std::vector cavity_potential_; }; } // namespace elecstate #endif diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 0cd6048e667..488bec38ea1 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/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); } diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp index 523d28d16f1..10ef2a00a56 100644 --- a/source/source_estate/test/test_pot_sccs.cpp +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -185,9 +185,23 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) EXPECT_EQ(config.boundary, ModuleSccs::Boundary::Pcc2d); EXPECT_EQ(config.pcc_2d_axis, axis); } + EXPECT_DOUBLE_EQ(config.cavity.lowpass_p1, -1.0); + EXPECT_FALSE(config.core_electrons); + 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); + EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); } -TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) +// Electrostatic forces in a dielectric, and the cavity terms of core +// Gaussians with surface tension and pressure. +TEST_F(PotSccsTest, SolvationForceMatchesTheEnergyDerivative) { Input_para input; input.sccs_epsilon = 5.0; @@ -197,6 +211,53 @@ TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) ModuleSccs::PolarizationSolverParameters solver; elecstate::sccs_cfg_input(input, config, solver); check_solvation_force(config, solver); + + Input_para full_input; + full_input.sccs_epsilon = 1.0; + full_input.sccs_tol_rms = 1e-13; + full_input.sccs_tol_max = 1e-12; + ModuleSccs::SccsConfig full; + ModuleSccs::PolarizationSolverParameters full_solver; + elecstate::sccs_cfg_input(full_input, full, full_solver); + full.core_electrons = true; + full.core_spreads = {2.0}; + full.surface_tension = 1e-5; + full.pressure = 1e-6; + check_solvation_force(full, full_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) diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt index d6f210a273a..5760b9f5b5d 100644 --- a/source/source_hamilt/module_sccs/CMakeLists.txt +++ b/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 diff --git a/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp index ddfbad02463..54c920bf8ce 100644 --- a/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp +++ b/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& density, double tpiba, bool open_boundary, SccsResponse& response, - std::vector& coefficient) + CavityDerivatives& derivatives) { + derivatives.filter = switch_filter(cavity, basis); + const bool lowpass = uses_lowpass(cavity); + derivatives.gradient.clear(); + std::vector& 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& density, // grid, since its chain coefficient is not converged at the cavity edge. std::vector> 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& density, } } eps_of_boundary(cavity, gradient, laplacian, response, coefficient); + if (open_boundary && lowpass) + { + derivatives.gradient.swap(gradient); + } } + } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_cav_deriv.h b/source/source_hamilt/module_sccs/sccs_cav_deriv.h index eb2d86a7521..d06fd8a8c90 100644 --- a/source/source_hamilt/module_sccs/sccs_cav_deriv.h +++ b/source/source_hamilt/module_sccs/sccs_cav_deriv.h @@ -1,6 +1,8 @@ #ifndef SCCS_CAV_DERIV_H #define SCCS_CAV_DERIV_H +#include "source_base/vector3.h" + #include namespace ModulePW @@ -13,6 +15,13 @@ namespace ModuleSccs struct CavityParameters; struct SccsResponse; +struct CavityDerivatives +{ + std::vector coefficient; + std::vector filter; + std::vector> gradient; +}; + // Collective over the PW pool. Periodic cells use density chain derivatives; // open boundaries differentiate solute on the FFT grid. void make_cav_derivs(const std::vector& density, @@ -21,7 +30,7 @@ void make_cav_derivs(const std::vector& density, double tpiba, bool open_boundary, SccsResponse& response, - std::vector& coefficient); + CavityDerivatives& derivatives); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_cavity.cpp b/source/source_hamilt/module_sccs/sccs_cavity.cpp index 08aeb097989..182ee672182 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.cpp +++ b/source/source_hamilt/module_sccs/sccs_cavity.cpp @@ -5,6 +5,11 @@ namespace ModuleSccs { +bool uses_lowpass(const CavityParameters& parameters) +{ + return parameters.lowpass_p1 > 0.0 && parameters.lowpass_p2 > 0.0; +} + CavityPoint evaluate_cavity(double density, const CavityParameters& parameters) { CavityPoint point; diff --git a/source/source_hamilt/module_sccs/sccs_cavity.h b/source/source_hamilt/module_sccs/sccs_cavity.h index ae4d987eda8..9c36bd2a50e 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.h +++ b/source/source_hamilt/module_sccs/sccs_cavity.h @@ -9,6 +9,9 @@ struct CavityParameters double density_min = 0.0; double density_max = 0.0; double epsilon_bulk = 0.0; + // Both positive enable the PCC switching-function derivative filter. + double lowpass_p1 = -1.0; + double lowpass_p2 = -1.0; }; struct CavityPoint @@ -20,6 +23,8 @@ struct CavityPoint double depsilon_drho = 0.0; }; +bool uses_lowpass(const CavityParameters& parameters); + // Parameters are validated by the INPUT reader. Negative Fourier ringing // uses the bulk limit. CavityPoint evaluate_cavity(double density, const CavityParameters& parameters); diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp index 19636b61ac8..d56867b5247 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -39,6 +39,7 @@ void eval_functional(const std::vector& charge, FunctionalResult candidate; candidate.reaction_potential.resize(size); candidate.electron_potential.resize(size); + candidate.cavity_potential = response.cavity_potential; for (std::size_t i = 0; i < size; ++i) { const double reaction = response.polarization.potential[i] - vacuum[i]; @@ -73,6 +74,7 @@ void eval_functional(const std::vector& charge, const double nonel = (config.pressure - config.surface_tension * divergence[i]) * response.dsolute_drho[i]; candidate.electron_potential[i] += nonel; + candidate.cavity_potential[i] += nonel; } result = std::move(candidate); } diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h index 9b9111abc86..249379dd913 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.h +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -24,6 +24,7 @@ struct FunctionalResult double volume = 0.0; // Bohr^3 std::vector reaction_potential; // dielectric minus vacuum, Hartree std::vector electron_potential; // derivative of all three energy terms + std::vector cavity_potential; // derivative w.r.t. cavity density, Hartree }; // Collective over the PW pool. diff --git a/source/source_hamilt/module_sccs/sccs_ion_charge.cpp b/source/source_hamilt/module_sccs/sccs_ion_charge.cpp index 0876ce5e0c2..0ccebd2dd25 100644 --- a/source/source_hamilt/module_sccs/sccs_ion_charge.cpp +++ b/source/source_hamilt/module_sccs/sccs_ion_charge.cpp @@ -14,11 +14,23 @@ void gauss_ion_rho(const std::vector& atoms, double tpiba, double spread, std::vector& density) +{ + const std::vector spreads(atoms.size(), spread); + gauss_ion_rho(atoms, basis, tpiba, spreads, density); +} + +void gauss_ion_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& spreads, + std::vector& density) { std::vector> charge(basis.npw, 0.0); const double tpiba2 = tpiba * tpiba; - for (const unitcell::AtomData& atom : atoms) + for (std::size_t ia = 0; ia < atoms.size(); ++ia) { + const unitcell::AtomData& atom = atoms[ia]; + const double spread = spreads[ia]; for (int ig = 0; ig < basis.npw; ++ig) { const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; @@ -30,4 +42,43 @@ void gauss_ion_rho(const std::vector& atoms, density.resize(basis.nrxx); basis.recip2real(charge.data(), density.data()); } + +void prep_core_gauss(const std::vector& atoms, + const std::vector& atom_spreads, + std::vector& core_atoms, + std::vector& widths) +{ + core_atoms = atoms; + widths.assign(atoms.size(), gaussian_ion_spread); + // A valence charge equal to the atomic number within round-off means no core. + const double charge_tolerance = 1e-8; + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + const std::size_t index = atom_spreads.size() == 1 ? 0 : ia; + const double width = atom_spreads[index]; + const double charge_difference = atoms[ia].atomic_number - atoms[ia].valence_charge; + const bool no_core = atoms[ia].atomic_number > 0 && std::abs(charge_difference) < charge_tolerance; + const bool automatic_skip = atom_spreads.size() == 1 && no_core; + if (width <= 0.0 || automatic_skip) + { + core_atoms[ia].valence_charge = 0.0; + } + else + { + widths[ia] = width; + } + } +} + +void gauss_core_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& atom_spreads, + std::vector& density) +{ + std::vector core_atoms; + std::vector widths; + prep_core_gauss(atoms, atom_spreads, core_atoms, widths); + gauss_ion_rho(core_atoms, basis, tpiba, widths, density); +} } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ion_charge.h b/source/source_hamilt/module_sccs/sccs_ion_charge.h index 4cbebd72af8..cbbd52e0459 100644 --- a/source/source_hamilt/module_sccs/sccs_ion_charge.h +++ b/source/source_hamilt/module_sccs/sccs_ion_charge.h @@ -24,6 +24,32 @@ void gauss_ion_rho(const std::vector& atoms, double tpiba, double spread, std::vector& density); + +// One positive width per atom. +void gauss_ion_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& spreads, + std::vector& density); + +// Resolve full-cavity widths without changing the physical ionic charge. +// atom_spreads holds one value or one per atom (checked when SCCS is set up). +void prep_core_gauss(const std::vector& atoms, + const std::vector& atom_spreads, + std::vector& core_atoms, + std::vector& widths); + +// Full cavity, as Environ 3.1.1 solvent_mode 'full' (core_electrons in +// ions.f90, used by boundary_electronic.f90): one Gaussian per atom with its +// valence charge and width corespread, added to the electron density that +// defines the cavity. +// A singleton skips known Z == zv atoms; per-atom lists override +// that choice. Nonpositive widths always disable the atom. +void gauss_core_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& atom_spreads, + std::vector& density); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_ion_force.cpp b/source/source_hamilt/module_sccs/sccs_ion_force.cpp index d76daa99d20..1af380e6cce 100644 --- a/source/source_hamilt/module_sccs/sccs_ion_force.cpp +++ b/source/source_hamilt/module_sccs/sccs_ion_force.cpp @@ -1,4 +1,5 @@ #include "sccs_ion_force.h" +#include "sccs_ion_charge.h" #include "source_base/constants.h" #include "source_base/parallel_reduce.h" @@ -16,6 +17,17 @@ void gauss_ion_force(const std::vector& atoms, double tpiba, double spread, std::vector>& forces) +{ + const std::vector spreads(atoms.size(), spread); + gauss_ion_force(atoms, reaction_potential, basis, tpiba, spreads, forces); +} + +void gauss_ion_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& spreads, + std::vector>& forces) { std::vector> potential_g(basis.npw); basis.real2recip(reaction_potential.data(), potential_g.data()); @@ -26,6 +38,7 @@ void gauss_ion_force(const std::vector& atoms, for (std::size_t ia = 0; ia < atoms.size(); ++ia) { const unitcell::AtomData& atom = atoms[ia]; + const double spread = spreads[ia]; for (int ig = 0; ig < basis.npw; ++ig) { const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; @@ -42,4 +55,17 @@ void gauss_ion_force(const std::vector& atoms, Parallel_Reduce::reduce_pool(forces[ia].z); } } + +void gauss_core_frc(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& atom_spreads, + std::vector>& forces) +{ + std::vector core_atoms; + std::vector widths; + prep_core_gauss(atoms, atom_spreads, core_atoms, widths); + gauss_ion_force(core_atoms, reaction_potential, basis, tpiba, widths, forces); +} } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ion_force.h b/source/source_hamilt/module_sccs/sccs_ion_force.h index a405ac4e551..e69c15a124c 100644 --- a/source/source_hamilt/module_sccs/sccs_ion_force.h +++ b/source/source_hamilt/module_sccs/sccs_ion_force.h @@ -26,5 +26,22 @@ void gauss_ion_force(const std::vector& atoms, double tpiba, double spread, std::vector>& forces); + +// One positive width per atom. +void gauss_ion_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& spreads, + std::vector>& forces); + +// Full cavity: a singleton skips known Z == zv atoms; per-atom lists override +// that choice. Nonpositive widths always disable the atom. +void gauss_core_frc(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector& atom_spreads, + std::vector>& forces); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_lowpass.cpp b/source/source_hamilt/module_sccs/sccs_lowpass.cpp new file mode 100644 index 00000000000..8f0758e318e --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_lowpass.cpp @@ -0,0 +1,81 @@ +#include "sccs_lowpass.h" +#include "sccs_cavity.h" +#include "sccs_cav_deriv.h" +#include "sccs_response.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" + +#include +#include +#include + +namespace ModuleSccs +{ +std::vector switch_filter(const CavityParameters& cavity, const ModulePW::PW_Basis& basis) +{ + const bool lowpass = uses_lowpass(cavity); + std::vector filter; + if (lowpass) + { + filter.resize(basis.npw); + for (int ig = 0; ig < basis.npw; ++ig) + { + // Environ deriv_lowpass: Gcut is the density cutoff, not ecutwfc. + const double argument = cavity.lowpass_p1 * basis.gg[ig] / basis.ggecut - cavity.lowpass_p2; + filter[ig] = 0.5 * std::erfc(argument); + } + } + return filter; +} + +void lowpass_cav_pot(const std::vector& charge, + const CavityParameters& cavity, + const CavityDerivatives& derivatives, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& response) +{ + const std::size_t size = charge.size(); + const double log_bulk = std::log(cavity.epsilon_bulk); + const std::vector& potential = response.polarization.potential; + std::vector weight(size); + for (std::size_t i = 0; i < size; ++i) + { + weight[i] = response.epsilon[i] * potential[i] * potential[i] / (8.0 * ModuleBase::PI); + } + // A = sqrt(eps) C^-1 sqrt(eps) + F is symmetric. For E=q^T A^-1 q/2, + // dE/ds = L/2 (q v + lapl b + L div(b grad s)), b=eps v^2/(8 pi). + // The filtered derivative transposes share one inverse FFT. + std::vector> transpose_g(basis.npw); + basis.real2recip(weight.data(), transpose_g.data()); + for (int ig = 0; ig < basis.npw; ++ig) + { + transpose_g[ig] *= -tpiba * tpiba * basis.gg[ig] * derivatives.filter[ig]; + } + std::vector component(size); + std::vector> component_g(basis.npw); + for (int d = 0; d < 3; ++d) + { + for (std::size_t i = 0; i < size; ++i) + { + component[i] = weight[i] * derivatives.gradient[i][d]; + } + basis.real2recip(component.data(), component_g.data()); + for (int ig = 0; ig < basis.npw; ++ig) + { + transpose_g[ig] += log_bulk * ModuleBase::IMAG_UNIT * tpiba * basis.gcar[ig][d] + * component_g[ig] * derivatives.filter[ig]; + } + } + std::vector transpose(size); + basis.recip2real(transpose_g.data(), transpose.data()); + std::vector candidate(size); + for (std::size_t i = 0; i < size; ++i) + { + const double boundary_potential = 0.5 * log_bulk * (charge[i] * potential[i] + transpose[i]); + candidate[i] = boundary_potential * response.dsolute_drho[i]; + } + response.cavity_potential = std::move(candidate); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_lowpass.h b/source/source_hamilt/module_sccs/sccs_lowpass.h new file mode 100644 index 00000000000..99cc5938ebc --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_lowpass.h @@ -0,0 +1,32 @@ +#ifndef SCCS_LOWPASS_H +#define SCCS_LOWPASS_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct CavityParameters; +struct CavityDerivatives; +struct SccsResponse; + +// Reciprocal-space filter on the local G vectors; empty when disabled. As +// Environ 3.1.1 core_fft_lowpass.f90 (filter_lowpass): f(G) = +// erfc(p1 G^2 / Gcut^2 - p2) / 2, applied before every FFT derivative. +std::vector switch_filter(const CavityParameters& cavity, const ModulePW::PW_Basis& basis); + +// Exact derivative of the filtered, discrete sqrt-CG reaction energy through +// the electronic-density cavity; collective over the PW pool. +void lowpass_cav_pot(const std::vector& charge, + const CavityParameters& cavity, + const CavityDerivatives& derivatives, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& response); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h index bcc805af4ac..145546b1999 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.h +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -4,6 +4,7 @@ #include "sccs_cavity.h" #include +#include namespace ModuleSccs { @@ -16,6 +17,8 @@ struct SccsConfig CavityParameters cavity; Boundary boundary = Boundary::Periodic; int pcc_2d_axis = 2; + bool core_electrons = false; + std::vector core_spreads = {0.5}; // Bohr: singleton skips Z == zv; per-atom list overrides; <= 0 disables double surface_tension = 0.0; // Hartree/Bohr^2 double pressure = 0.0; // Hartree/Bohr^3 double surface_regularization = 0.0; diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index e7e5cc6b40d..9904db6d915 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -1,5 +1,6 @@ #include "sccs_response.h" #include "sccs_cav_deriv.h" +#include "sccs_lowpass.h" #include "sccs_parameters.h" #include "sccs_pw_coulomb.h" @@ -216,6 +217,7 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, void finish_response(const SqrtCgOperator& system, const std::vector& charge, const CavityParameters& cavity, + const CavityDerivatives& derivatives, double tpiba, std::vector& potential, SccsResponse& response) @@ -242,23 +244,34 @@ void finish_response(const SqrtCgOperator& system, value -= mean; } } - std::vector> potential_g(basis.npw); - basis.real2recip(potential.data(), potential_g.data()); - polarization.gradient.resize(size); - XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); + const bool lowpass = uses_lowpass(cavity); + if (!lowpass) + { + std::vector> potential_g(basis.npw); + basis.real2recip(potential.data(), potential_g.data()); + polarization.gradient.resize(size); + XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); + } polarization.potential.swap(potential); if (open_boundary) { finish_open_bc(charge, coefficient, invsqrt, cavity, basis, response); } - // Continuum derivative, as in Environ: exact for the discrete energy only - // when |grad v|^2 is resolved by the FFT grid. - response.cavity_potential.resize(size); - for (std::size_t i = 0; i < size; ++i) + if (lowpass) + { + lowpass_cav_pot(charge, cavity, derivatives, basis, tpiba, response); + } + else { - const ModuleBase::Vector3& gradient = polarization.gradient[i]; - const double gradient_square = gradient * gradient; - response.cavity_potential[i] = -response.depsilon_drho[i] * gradient_square / (8.0 * ModuleBase::PI); + // Continuum derivative, as in Environ: exact for the discrete energy only + // when |grad v|^2 is resolved by the FFT grid. + response.cavity_potential.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const ModuleBase::Vector3& gradient = polarization.gradient[i]; + const double gradient_square = gradient * gradient; + response.cavity_potential[i] = -response.depsilon_drho[i] * gradient_square / (8.0 * ModuleBase::PI); + } } } } // namespace @@ -293,8 +306,9 @@ void solve_response(const std::vector& density, const bool open_boundary = coulomb.has_correction(); SccsResponse candidate; - std::vector coefficient; - make_cav_derivs(density, cavity, basis, tpiba, open_boundary, candidate, coefficient); + CavityDerivatives derivatives; + make_cav_derivs(density, cavity, basis, tpiba, open_boundary, candidate, derivatives); + const std::vector& coefficient = derivatives.coefficient; const std::size_t size = density.size(); std::vector invsqrt(size); for (std::size_t i = 0; i < size; ++i) @@ -313,7 +327,7 @@ void solve_response(const std::vector& density, } iterate_sqrt_cg(system, solver, state, polarization); std::vector& potential = state.potential; - finish_response(system, charge, cavity, tpiba, potential, candidate); + finish_response(system, charge, cavity, derivatives, tpiba, potential, candidate); result = std::move(candidate); } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index 076ffbb6459..eda88069a08 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -24,7 +24,7 @@ struct PolarizationResult double residual_max = 0.0; bool warm_started = false; std::vector potential; // Hartree; zero mean only for periodic cells - std::vector> gradient; // Hartree/Bohr + std::vector> gradient; // Hartree/Bohr; empty with lowpass }; struct SccsResponse @@ -35,7 +35,9 @@ struct SccsResponse std::vector depsilon_drho; std::vector> grad_log_epsilon; PolarizationResult polarization; - std::vector cavity_potential; // -epsilon' |grad v|^2/(8 pi), Hartree + // Hartree: discrete cavity derivative with lowpass, otherwise the continuum + // derivative -epsilon' |grad v|^2/(8 pi). + std::vector cavity_potential; std::vector restart_potential; // unshifted solution for explicit warm starts double far_field_polarization_charge = 0.0; // open-boundary screened charge minus solute charge }; diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt index deb7efe8302..3a8e88d9dfe 100644 --- a/source/source_hamilt/module_sccs/test/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -6,7 +6,7 @@ AddTest(TARGET MODULE_HAMILT_sccs_parameters AddTest(TARGET MODULE_HAMILT_sccs_poisson LIBS parameter planewave device base - SOURCES test_sccs_cav_deriv.cpp test_sccs_coul_0d.cpp test_sccs_coul_2d.cpp test_sccs_pw_coulomb.cpp test_sccs_ion_charge.cpp test_sccs_ion_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + SOURCES test_sccs_lowpass.cpp test_sccs_cav_deriv.cpp test_sccs_coul_0d.cpp test_sccs_coul_2d.cpp test_sccs_pw_coulomb.cpp test_sccs_ion_charge.cpp test_sccs_ion_force.cpp test_sccs_response.cpp test_sccs_functional.cpp ../../../source_estate/test/test_pot_sccs.cpp ../../../source_estate/test/test_sccs_pcc_0d.cpp ../../../source_estate/test/test_sccs_pcc_2d.cpp @@ -21,6 +21,7 @@ AddTest(TARGET MODULE_HAMILT_sccs_poisson ../sccs_coul_0d.cpp ../sccs_coul_2d.cpp ../sccs_moments.cpp + ../sccs_lowpass.cpp ../sccs_cav_deriv.cpp ../sccs_functional.cpp ../sccs_parameters.cpp diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp index 0c3cfd6eba8..bb941d08b90 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp @@ -40,9 +40,10 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf } const bool open_boundary = true; ModuleSccs::SccsResponse response; - std::vector coefficient; + ModuleSccs::CavityDerivatives derivatives; + const std::vector& coefficient = derivatives.coefficient; ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, - response, coefficient); + response, derivatives); const double log_bulk = std::log(cavity.epsilon_bulk); for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -65,9 +66,10 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf } const bool periodic_boundary = false; ModuleSccs::SccsResponse periodic; - std::vector periodic_coefficient; + ModuleSccs::CavityDerivatives periodic_derivatives; + const std::vector& periodic_coefficient = periodic_derivatives.coefficient; ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, periodic_boundary, - periodic, periodic_coefficient); + periodic, periodic_derivatives); double difference = 0.0; for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -81,4 +83,21 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf Parallel_Reduce::reduce_max_pool(basis.poolnproc, difference); EXPECT_GT(difference, 1e-8); + cavity.lowpass_p1 = basis.ggecut; + cavity.lowpass_p2 = 0.5; + ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + response, derivatives); + const double weight = 0.5 * std::erfc(0.5); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const double angle = ModuleBase::TWO_PI * ix / basis.nx; + const double gradient = -amplitude * tpiba * std::sin(angle) * weight; + const double laplacian = -amplitude * tpiba * tpiba * mode[ir] * weight; + const double expected = response.epsilon[ir] * (-0.5 * log_bulk * laplacian + + 0.25 * log_bulk * log_bulk * gradient * gradient) + / ModuleBase::FOUR_PI; + EXPECT_NEAR(coefficient[ir], expected, 1e-12); + EXPECT_NEAR(derivatives.gradient[ir].x, gradient, 1e-12); + } } diff --git a/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp b/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp index 9e3d9b1cf34..4acaf4e3a81 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp @@ -186,9 +186,10 @@ TEST_F(SccsPcc0dCoulombTest, NonuniformResponseSatisfiesFixedPointAndFarFieldCha EXPECT_GT(response.polarization.iterations, 1); const bool open_boundary = true; ModuleSccs::SccsResponse derivatives; - std::vector coefficient; + ModuleSccs::CavityDerivatives cavity_derivatives; + const std::vector& coefficient = cavity_derivatives.coefficient; ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, - derivatives, coefficient); + derivatives, cavity_derivatives); std::vector source(basis.nrxx); for (int ir = 0; ir < basis.nrxx; ++ir) { diff --git a/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp b/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp index da3ed193790..c4ced8cb89d 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp @@ -173,9 +173,10 @@ TEST_F(SccsPcc2dCoulombTest, NonuniformResponseSatisfiesFixedPointAndWarmStart) basis, tpiba, coulomb, response); const bool open_boundary = true; ModuleSccs::SccsResponse derivatives; - std::vector coefficient; + ModuleSccs::CavityDerivatives cavity_derivatives; + const std::vector& coefficient = cavity_derivatives.coefficient; ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, - derivatives, coefficient); + derivatives, cavity_derivatives); std::vector source(basis.nrxx); for (int ir = 0; ir < basis.nrxx; ++ir) { diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp index 19fedc4cc11..ca4ea6ac9fe 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp @@ -51,3 +51,62 @@ TEST_F(SccsIonicChargeTest, ValenceNormalizationAndFourierPhase) EXPECT_NEAR(translated[ir], density[ir], 1e-13); } } + +// Atoms without a known atomic number keep the singleton width; a known Z == +// zv atom is skipped by a singleton but not by a per-atom list; widths <= 0 +// disable an atom; the atom data are never modified. +TEST_F(SccsIonicChargeTest, CoreGaussianSelection) +{ + std::vector atoms(3); + atoms[0].valence_charge = 6.0; + atoms[1].valence_charge = 1.0; + atoms[2].valence_charge = 7.0; + atoms[0].position.x = length / 4.0; + std::vector widths = {0.5}; + std::vector density; + auto integral = [&]() { + double sum = 0.0; + for (double value : density) + { + sum += value; + } + Parallel_Reduce::reduce_pool(sum); + return sum * basis.omega / basis.nxyz; + }; + ModuleSccs::gauss_core_rho(atoms, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 14.0, 1e-12); + widths = {0.8}; + ModuleSccs::gauss_core_rho(atoms, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 14.0, 1e-12); + widths = {0.7, 0.0, -1.0}; + ModuleSccs::gauss_core_rho(atoms, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 6.0, 1e-12); + std::vector selected; + std::vector resolved; + ModuleSccs::prep_core_gauss(atoms, widths, selected, resolved); + EXPECT_DOUBLE_EQ(resolved[0], 0.7); + EXPECT_DOUBLE_EQ(selected[1].valence_charge, 0.0); + EXPECT_DOUBLE_EQ(atoms[1].valence_charge, 1.0); + widths = {0.0}; + ModuleSccs::gauss_core_rho(atoms, basis, tpiba, widths, density); + EXPECT_DOUBLE_EQ(integral(), 0.0); + + std::vector known(4); + known[0].atomic_number = 8; + known[0].valence_charge = 6.0; + known[1].atomic_number = 1; + known[1].valence_charge = 1.0; + known[2].atomic_number = 2; + known[2].valence_charge = 2.0; + known[3].valence_charge = 3.0; // Unknown element: do not silently exclude. + widths = {0.8}; + ModuleSccs::gauss_core_rho(known, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 9.0, 1e-12); + widths = {0.8, 0.8, 0.8, 0.8}; + ModuleSccs::gauss_core_rho(known, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 12.0, 1e-12); + widths = {0.8, 0.0, -1.0, 0.8}; + ModuleSccs::gauss_core_rho(known, basis, tpiba, widths, density); + EXPECT_NEAR(integral(), 9.0, 1e-12); + EXPECT_DOUBLE_EQ(known[1].valence_charge, 1.0); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp index 851a1ce74f4..0fb5ff96d3e 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp @@ -6,6 +6,8 @@ #include "../sccs_functional.h" #include "source_cell/cell_tools.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "../sccs_coul_0d.h" using SccsIonicForceTest = SccsTest::PwTest; @@ -95,3 +97,134 @@ TEST_F(SccsIonicForceTest, ReactionEnergyFiniteDifferenceAtFixedNonuniformCavity EXPECT_NEAR(force.z, 0.0, 1e-12); } } + +// The core force follows the same width selection as the core density. +TEST_F(SccsIonicForceTest, CoreForceSelection) +{ + std::vector atoms(2); + atoms[0].valence_charge = 6.0; + atoms[1].valence_charge = 1.0; + for (unitcell::AtomData& atom : atoms) + { + atom.position.x = length / 4.0; + } + const std::vector potential = cosine_mode(0); + std::vector> force; + std::vector widths = {0.8, 0.0}; + ModuleSccs::gauss_core_frc(atoms, potential, basis, tpiba, widths, force); + const double gaussian_exponent = -0.25 * 0.8 * 0.8 * tpiba * tpiba; + const double expected = 6.0 * tpiba * std::exp(gaussian_exponent); + EXPECT_NEAR(force[0].x, expected, 1e-12); + EXPECT_DOUBLE_EQ(force[1].x, 0.0); + widths = {0.8}; + ModuleSccs::gauss_core_frc(atoms, potential, basis, tpiba, widths, force); + const double expected_hydrogen_force = expected / 6.0; + EXPECT_NEAR(force[1].x, expected_hydrogen_force, 1e-12); + widths = {0.0}; + ModuleSccs::gauss_core_frc(atoms, potential, basis, tpiba, widths, force); + EXPECT_DOUBLE_EQ(force[0].x, 0.0); + EXPECT_DOUBLE_EQ(force[1].x, 0.0); + + std::vector known(2); + known[0].atomic_number = 8; + known[0].valence_charge = 6.0; + known[1].atomic_number = 1; + known[1].valence_charge = 1.0; + for (unitcell::AtomData& atom : known) + { + atom.position.x = length / 4.0; + } + std::vector> forces; + widths = {0.8}; + ModuleSccs::gauss_core_frc(known, potential, basis, tpiba, widths, forces); + EXPECT_DOUBLE_EQ(forces[1].x, 0.0); + const double oxygen_force = forces[0].x; + widths = {0.8, 0.8}; + ModuleSccs::gauss_core_frc(known, potential, basis, tpiba, widths, forces); + EXPECT_DOUBLE_EQ(forces[0].x, oxygen_force); + const double hydrogen_force = oxygen_force / 6.0; + EXPECT_NEAR(forces[1].x, hydrogen_force, 1e-12); +} + +TEST_F(SccsIonicForceTest, FullCavityPccLowpassTotalCorrectionForceFiniteDifference) +{ + std::vector atoms(2); + atoms[0].valence_charge = 1.0; + atoms[0].position = ModuleBase::Vector3(3.1, 3.2, 4.3); + atoms[1].valence_charge = 1.0; + atoms[1].position = ModuleBase::Vector3(6.4, 5.1, 5.7); + const ModuleBase::Matrix3 lattice; + std::vector> positions; + ModulePW::grid_positions(basis, lattice, length, positions); + const double half = 0.5 * length; + const ModuleBase::Vector3 center(half, half, half); + for (ModuleBase::Vector3& position : positions) + { + position -= center; + } + elecstate::Pcc0dParameters parameters; + parameters.length = length; + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, parameters); + ModuleSccs::SccsConfig config; + config.boundary = ModuleSccs::Boundary::Pcc0d; + config.core_electrons = true; + config.core_spreads = {2.0, 0.0}; + config.cavity.epsilon_bulk = 5.0; + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + config.cavity.lowpass_p1 = basis.ggecut; + config.cavity.lowpass_p2 = 0.5; + config.surface_regularization = 1e-8; + config.surface_tension = 1e-5; + config.pressure = 1e-6; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-13; + solver.tolerance_max = 1e-13; + const std::vector cold; + auto evaluate = [&](const std::vector& displaced, + ModuleSccs::FunctionalResult& functional) { + std::vector ions; + std::vector core; + ModuleSccs::gauss_ion_rho(displaced, basis, tpiba, ModuleSccs::gaussian_ion_spread, ions); + ModuleSccs::gauss_core_rho(displaced, basis, tpiba, config.core_spreads, core); + std::vector cavity_density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double electron_density = 2.0 / basis.omega; + cavity_density[ir] = electron_density + core[ir]; + ions[ir] -= electron_density; + } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(cavity_density, ions, config.cavity, solver, cold, basis, tpiba, coulomb, + response); + ModuleSccs::eval_functional(ions, response, config, basis, tpiba, coulomb, functional); + }; + ModuleSccs::FunctionalResult baseline; + evaluate(atoms, baseline); + std::vector> ionic; + std::vector> core; + ModuleSccs::gauss_ion_force(atoms, baseline.reaction_potential, basis, + tpiba, ModuleSccs::gaussian_ion_spread, ionic); + ModuleSccs::gauss_core_frc(atoms, baseline.cavity_potential, basis, + tpiba, config.core_spreads, core); + const double step = 1e-4; + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + std::vector plus = atoms; + std::vector minus = atoms; + plus[ia].position[axis] += step; + minus[ia].position[axis] -= step; + ModuleSccs::FunctionalResult positive; + ModuleSccs::FunctionalResult negative; + evaluate(plus, positive); + evaluate(minus, negative); + const double energy_plus = positive.reaction_energy + positive.surface_energy + positive.volume_energy; + const double energy_minus = negative.reaction_energy + negative.surface_energy + negative.volume_energy; + const double finite_difference = -(energy_plus - energy_minus) / (2.0 * step); + const double predicted = ionic[ia][axis] + core[ia][axis]; + EXPECT_NEAR(predicted, finite_difference, 1e-8); + } + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp new file mode 100644 index 00000000000..dd202dabf41 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp @@ -0,0 +1,146 @@ +#include "sccs_test.h" +#include "../sccs_lowpass.h" +#include "../sccs_cav_deriv.h" +#include "../sccs_functional.h" +#include "../sccs_parameters.h" +#include "../sccs_coul_0d.h" +#include "../sccs_coul_2d.h" +#include "../sccs_response.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_grid_geometry.h" + + +class SccsLowpassTest : public SccsTest::PwTest +{ +protected: + void check_electron_derivative(ModuleSccs::CoulombOperator& coulomb, ModuleSccs::Boundary boundary) + { + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + config.boundary = boundary; + config.cavity.epsilon_bulk = 5.0; + // Strong filtering on the first shell makes omitting either transpose + // detectable with a small FFT grid. + config.cavity.lowpass_p1 = basis.ggecut; + config.cavity.lowpass_p2 = 0.5; + config.surface_tension = 0.0; + config.pressure = 0.0; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-13; + solver.tolerance_max = 1e-13; + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + std::vector direction(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 + 1e-4 * mode[ir]; + charge[ir] = 5e-4 + 2e-4 * mode[ir]; + direction[ir] = 1e-4 * mode[ir]; + } + const std::vector cold; + ModuleSccs::SccsResponse center_response; + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + basis, tpiba, coulomb, center_response); + EXPECT_TRUE(center_response.polarization.gradient.empty()); + ModuleSccs::FunctionalResult center; + ModuleSccs::eval_functional(charge, center_response, config, basis, + tpiba, coulomb, center); + double predicted = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + predicted += center.electron_potential[ir] * direction[ir]; + } + Parallel_Reduce::reduce_pool(predicted); + predicted *= basis.omega / basis.nxyz; + const double steps[] = {1e-3, 5e-4}; + for (double step : steps) + { + double energies[2]; + for (int side = 0; side < 2; ++side) + { + const double sign = side == 0 ? 1.0 : -1.0; + std::vector displaced_density(basis.nrxx); + std::vector displaced_charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double change = sign * step * direction[ir]; + displaced_density[ir] = density[ir] + change; + displaced_charge[ir] = charge[ir] - change; + } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(displaced_density, displaced_charge, + config.cavity, solver, cold, basis, tpiba, coulomb, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::eval_functional(displaced_charge, response, config, + basis, tpiba, coulomb, functional); + energies[side] = functional.reaction_energy; + } + const double finite_difference = (energies[0] - energies[1]) / (2.0 * step); + EXPECT_GT(std::abs(predicted), 1e-5); + EXPECT_NEAR(predicted, finite_difference, 1e-8); + } + } + + std::vector> centered_positions() + { + const ModuleBase::Matrix3 lattice; + std::vector> positions; + ModulePW::grid_positions(basis, lattice, length, positions); + const double half_length = 0.5 * length; + const ModuleBase::Vector3 center(half_length, half_length, half_length); + for (ModuleBase::Vector3& position : positions) + { + position -= center; + } + return positions; + } +}; + +TEST_F(SccsLowpassTest, PccElectronPotentialDifferentiatesDiscreteEnergy) +{ + elecstate::Pcc0dParameters molecule; + molecule.length = length; + const std::vector> positions = centered_positions(); + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, molecule); + check_electron_derivative(coulomb, ModuleSccs::Boundary::Pcc0d); + elecstate::Pcc2dParameters slab; + slab.length = length; + slab.area = length * length; + std::vector> projected = centered_positions(); + for (ModuleBase::Vector3& position : projected) + { + position.x = 0.0; + position.y = 0.0; + } + const ModuleBase::Vector3 normal(0.0, 0.0, 1.0); + ModuleSccs::Pcc2dCoulombOperator slab_coulomb(basis, tpiba, projected, normal, slab); + check_electron_derivative(slab_coulomb, ModuleSccs::Boundary::Pcc2d); +} + +TEST_F(SccsLowpassTest, FilterUsesDensityCutoff) +{ + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + std::vector filter = ModuleSccs::switch_filter(config.cavity, basis); + EXPECT_TRUE(filter.empty()); + config.cavity.lowpass_p1 = 10.0; + config.cavity.lowpass_p2 = 5.0; + filter = ModuleSccs::switch_filter(config.cavity, basis); + ASSERT_EQ(filter.size(), static_cast(basis.npw)); + const double half_cutoff = basis.ggecut / 2.0; + for (int ig = 0; ig < basis.npw; ++ig) + { + if (basis.gg[ig] == 0.0) + { + EXPECT_NEAR(filter[ig], 1.0, 1e-12); + } + if (basis.gg[ig] > half_cutoff) + { + EXPECT_LT(filter[ig], 0.5); + } + else + { + EXPECT_GE(filter[ig], 0.5); + } + } +} diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 32a97fbdb37..6444267178b 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -133,6 +133,7 @@ add_library( module_parameter/read_inp_deepks.cpp module_parameter/read_inp_model.cpp module_parameter/read_inp_sccs.cpp + module_parameter/sccs_inp_check.cpp module_parameter/read_inp_postproc.cpp module_parameter/read_inp_exx_dftu.cpp module_parameter/read_inp_other.cpp diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1c624aa5a02..75be1fba4b8 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -598,6 +598,10 @@ struct Input_para double sccs_tol_rms = 1e-10; double sccs_tol_max = 1e-8; double sccs_surface_eta = 1e-8; // Bohr^-1 + double sccs_lowpass_p1 = -1.0; // PCC switching-derivative filter slope + double sccs_lowpass_p2 = -1.0; // PCC switching-derivative filter offset + std::string sccs_solvent_mode = "electronic"; + std::vector sccs_corespread = {0.5}; // Bohr: one value or one per atom // ============== #Parameters (14.vdW Correction) =========================== // ========================================================== diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index 961002946cf..37aa6d915a3 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -197,7 +197,8 @@ void ReadInput::item_model() "pcc_2d; no electric/gate field, DFT-1/2, DeePKS or dm_to_rho. Use assume_isolated pcc_0d " "or pcc_2d for charged systems: with assume_isolated none a charged cell only prints a " "warning and runs, and its energy then depends on the cell size because the periodic " - "Poisson solver drops the G = 0 net charge.\n\n" + "Poisson solver drops the G = 0 net charge. With pcc_0d or pcc_2d, see sccs_lowpass_p1 " + "for the accuracy of forces.\n\n" "Legacy Boolean values (true/false) are read as 1/0."; item.default_value = "0"; item.read_value = [](const Input_Item& item, Parameter& para) { diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 816ea7c6568..06e6a1dcf1a 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -3,139 +3,8 @@ #include "read_input_tool.h" #include "source_base/tool_quit.h" -#include -#include -#include -#include - namespace ModuleIO { -bool parse_imp_sol(const std::string& value, int& model, std::string& error) -{ - error.clear(); - std::string normalized = value; - for (char& c : normalized) - { - const unsigned char byte = static_cast(c); - const int lower = std::tolower(byte); - c = static_cast(lower); - } - const std::vector yes = {"true", "1", "t", "yes", "y", "on", ".true."}; - const std::vector no = {"false", "0", "f", "no", "n", "off", ".false."}; - if (normalized == "2") { model = 2; } - else if (std::find(yes.begin(), yes.end(), normalized) != yes.end()) { model = 1; } - else if (std::find(no.begin(), no.end(), normalized) != no.end()) { model = 0; } - else - { - error = "imp_sol must be 0, 1, 2 or a legacy Boolean value"; - return false; - } - return true; -} - -namespace -{ -// Calculation settings that SCCS supports; empty when all are supported. -std::string context_error(const Input_para& input) -{ - std::string error; - if (input.assume_isolated != "none" && input.assume_isolated != "pcc_0d" - && input.assume_isolated != "pcc_2d") - { - error = "SCCS currently supports assume_isolated none, pcc_0d or pcc_2d"; - } - else if (input.device != "cpu" || input.esolver_type != "ksdft") - { - error = "SCCS requires CPU KS-DFT"; - } - else if (input.basis_type != "pw" && input.basis_type != "lcao") - { - error = "SCCS requires basis_type pw or lcao"; - } - else if (input.calculation != "scf" || input.cal_stress) - { - error = "SCCS currently supports SCF energies and forces without stress"; - } - else if (input.nspin != 1 && input.nspin != 2) - { - error = "SCCS requires nspin 1 or 2"; - } - else if (input.efield_flag || input.gate_flag || input.dfthalf_type != 0 - || input.deepks_scf || input.deepks_out_labels || input.deepks_bandgap - || input.deepks_v_delta || input.dm_to_rho) - { - error = "SCCS does not support electric/gate fields, DFT-1/2, DeePKS or dm_to_rho"; - } - else if (!input.vl_in_h || !input.vion_in_h || !input.vh_in_h) - { - error = "SCCS requires the local ionic and Hartree potentials in the Hamiltonian"; - } - return error; -} - -// The shared real parser (std::stod) accepts nan and inf, which pass the -// ordered comparisons below; signed values stay allowed for gamma/pressure. -std::string nonfinite_error(const Input_para& input) -{ - const std::vector> values = { - {"sccs_epsilon", input.sccs_epsilon}, - {"sccs_rho_min", input.sccs_rho_min}, - {"sccs_rho_max", input.sccs_rho_max}, - {"sccs_gamma", input.sccs_gamma}, - {"sccs_pressure", input.sccs_pressure}, - {"sccs_tol_rms", input.sccs_tol_rms}, - {"sccs_tol_max", input.sccs_tol_max}, - {"sccs_surface_eta", input.sccs_surface_eta}}; - for (const std::pair& value : values) - { - if (!std::isfinite(value.second)) { return std::string(value.first) + " must be finite"; } - } - return std::string(); -} - -// Cavity and surface model parameters, starting with the preset. -std::string model_error(const Input_para& input) -{ - const std::vector presets = {"custom", "vacuum", "water-neutral", "water-cation", "water-anion"}; - if (std::find(presets.begin(), presets.end(), input.sccs_preset) == presets.end()) - { - return "Unknown sccs_preset"; - } - std::string error; - if (input.sccs_epsilon < 1.0) { error = "sccs_epsilon must be at least 1"; } - else if (input.sccs_rho_min <= 0.0 || input.sccs_rho_max <= input.sccs_rho_min) - { error = "sccs_rho_min and sccs_rho_max must satisfy 0 < sccs_rho_min < sccs_rho_max"; } - else if (input.sccs_surface_eta <= 0.0) { error = "sccs_surface_eta must be positive"; } - return error; -} - -// sqrt-CG solver controls. -std::string solver_error(const Input_para& input) -{ - std::string error; - if (input.sccs_maxiter <= 0) { error = "sccs_maxiter must be positive"; } - else if (input.sccs_tol_rms <= 0.0) { error = "sccs_tol_rms must be positive"; } - else if (input.sccs_tol_max <= 0.0) { error = "sccs_tol_max must be positive"; } - return error; -} -} // namespace - -bool check_sccs_inp(const Input_para& input, std::string& error) -{ - error.clear(); - if (input.imp_sol < 0 || input.imp_sol > 2) - { - error = "imp_sol must be 0, 1 or 2"; - return false; - } - if (input.imp_sol != 2) { return true; } - error = context_error(input); - if (error.empty()) { error = nonfinite_error(input); } - if (error.empty()) { error = model_error(input); } - if (error.empty()) { error = solver_error(input); } - return error.empty(); -} - void ReadInput::item_sccs() { { @@ -258,5 +127,61 @@ void ReadInput::item_sccs() read_sync_double(input.sccs_surface_eta); this->add_item(item); } + { + Input_Item item("sccs_lowpass_p1"); + item.annotation = "SCCS switching-derivative low-pass slope"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Low-pass filter of the SCCS switching-function derivatives, as Environ deriv_lowpass_p1 with deriv_method fft. Only with assume_isolated pcc_0d or pcc_2d. The default -1 turns it off and reproduces Environ deriv_method fft (continuum cavity potential). With lowpass disabled (the default), analytical forces may differ from finite differences of the self-consistent energy. For geometry optimization with PCC, consider enabling lowpass and check force accuracy against finite differences. With lowpass disabled, the cavity potential uses the FFT gradient of the PCC-corrected potential, which oscillates around the potential step at the cell boundary half a cell from the system center; keep the dielectric transition region several bohr away from that boundary. Typical values are sccs_lowpass_p1 10 and sccs_lowpass_p2 5. The filter changes the model energy, so compare only results computed with the same values."; + item.default_value = "-1"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_lowpass_p1); + this->add_item(item); + } + { + Input_Item item("sccs_lowpass_p2"); + item.annotation = "SCCS switching-derivative low-pass offset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Offset of the SCCS switching-function low-pass filter, as Environ deriv_lowpass_p2; see sccs_lowpass_p1. Both must be positive or both non-positive. Default -1 (off)."; + item.default_value = "-1"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_lowpass_p2); + this->add_item(item); + } + { + Input_Item item("sccs_solvent_mode"); + item.annotation = "density that defines the SCCS cavity"; + item.category = "Implicit solvation model"; + item.type = "String"; + item.description = "Allowed values: electronic (default) and full, as Environ solvent_mode. electronic builds the dielectric cavity from the valence electron density. full adds valence-charge Gaussians selected by sccs_corespread."; + item.default_value = "electronic"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_string(input.sccs_solvent_mode); + this->add_item(item); + } + { + Input_Item item("sccs_corespread"); + item.annotation = "Widths of the full-cavity Gaussians"; + item.category = "Implicit solvation model"; + item.type = "Vector of Real"; + item.description = "Widths of the valence-charge Gaussians used only with sccs_solvent_mode full. One value applies to every atom except those whose known atomic number equals the pseudopotential valence charge (within 1e-8). More than one value requires exactly nat values in STRU atom order (grouped by atom type), and overrides this automatic exclusion. A value <= 0 disables the cavity Gaussian on that atom."; + item.default_value = "0.5"; + item.unit = "bohr"; + item.set_availability("imp_sol==2"); + item.read_value = [](const Input_Item& item, Parameter& para) { + std::string error; + const bool valid = parse_spreads(item.str_values, para.input.sccs_corespread, error); + if (!valid) + { + ModuleBase::WARNING_QUIT("ReadInput", error); + } + }; + sync_doublevec(input.sccs_corespread, para.input.sccs_corespread.size(), 0.0); + this->add_item(item); + } } } // namespace ModuleIO diff --git a/source/source_io/module_parameter/read_inp_sccs.h b/source/source_io/module_parameter/read_inp_sccs.h index cb7c926c52c..38e7af025c2 100644 --- a/source/source_io/module_parameter/read_inp_sccs.h +++ b/source/source_io/module_parameter/read_inp_sccs.h @@ -1,10 +1,14 @@ #ifndef READ_INP_SCCS_H #define READ_INP_SCCS_H #include +#include struct Input_para; namespace ModuleIO { bool parse_imp_sol(const std::string& value, int& model, std::string& error); bool check_sccs_inp(const Input_para& input, std::string& error); +// Parse sccs_corespread words as finite real numbers; on failure spreads is +// unchanged and error names the offending word. +bool parse_spreads(const std::vector& words, std::vector& spreads, std::string& error); } #endif diff --git a/source/source_io/module_parameter/sccs_inp_check.cpp b/source/source_io/module_parameter/sccs_inp_check.cpp new file mode 100644 index 00000000000..d0567753214 --- /dev/null +++ b/source/source_io/module_parameter/sccs_inp_check.cpp @@ -0,0 +1,224 @@ +// Parsing and validation of the SCCS INPUT keys; the keys themselves are +// registered in read_inp_sccs.cpp. +#include "read_inp_sccs.h" +#include "input_parameter.h" + +#include +#include +#include +#include +#include + +namespace ModuleIO +{ +bool parse_imp_sol(const std::string& value, int& model, std::string& error) +{ + error.clear(); + std::string normalized = value; + for (char& c : normalized) + { + const unsigned char byte = static_cast(c); + const int lower = std::tolower(byte); + c = static_cast(lower); + } + const std::vector yes = {"true", "1", "t", "yes", "y", "on", ".true."}; + const std::vector no = {"false", "0", "f", "no", "n", "off", ".false."}; + if (normalized == "2") + { + model = 2; + } + else if (std::find(yes.begin(), yes.end(), normalized) != yes.end()) + { + model = 1; + } + else if (std::find(no.begin(), no.end(), normalized) != no.end()) + { + model = 0; + } + else + { + error = "imp_sol must be 0, 1, 2 or a legacy Boolean value"; + return false; + } + return true; +} + +bool parse_spreads(const std::vector& words, std::vector& spreads, std::string& error) +{ + error.clear(); + std::vector values; + for (const std::string& word : words) + { + const char* begin = word.c_str(); + char* end = nullptr; + const double value = std::strtod(begin, &end); + const bool complete = end != begin && *end == '\0'; + if (!complete || !std::isfinite(value)) + { + error = "sccs_corespread values must be finite real numbers, got '" + word + "'"; + return false; + } + values.push_back(value); + } + spreads.swap(values); + return true; +} + +namespace +{ +// Calculation settings that SCCS supports; empty when all are supported. +std::string context_error(const Input_para& input) +{ + std::string error; + if (input.assume_isolated != "none" && input.assume_isolated != "pcc_0d" + && input.assume_isolated != "pcc_2d") + { + error = "SCCS currently supports assume_isolated none, pcc_0d or pcc_2d"; + } + else if (input.device != "cpu" || input.esolver_type != "ksdft") + { + error = "SCCS requires CPU KS-DFT"; + } + else if (input.basis_type != "pw" && input.basis_type != "lcao") + { + error = "SCCS requires basis_type pw or lcao"; + } + else if (input.calculation != "scf" || input.cal_stress) + { + error = "SCCS currently supports SCF energies and forces without stress"; + } + else if (input.nspin != 1 && input.nspin != 2) + { + error = "SCCS requires nspin 1 or 2"; + } + else if (input.efield_flag || input.gate_flag || input.dfthalf_type != 0 + || input.deepks_scf || input.deepks_out_labels || input.deepks_bandgap + || input.deepks_v_delta || input.dm_to_rho) + { + error = "SCCS does not support electric/gate fields, DFT-1/2, DeePKS or dm_to_rho"; + } + else if (!input.vl_in_h || !input.vion_in_h || !input.vh_in_h) + { + error = "SCCS requires the local ionic and Hartree potentials in the Hamiltonian"; + } + return error; +} + +// nan and inf pass the ordered comparisons below. The INPUT reader already +// rejects them; this also covers an Input_para filled in code. Signed values +// stay allowed for gamma/pressure and for the lowpass pair, where non-positive +// values turn the filter off. +std::string nonfinite_error(const Input_para& input) +{ + const std::vector> values = { + {"sccs_epsilon", input.sccs_epsilon}, + {"sccs_rho_min", input.sccs_rho_min}, + {"sccs_rho_max", input.sccs_rho_max}, + {"sccs_gamma", input.sccs_gamma}, + {"sccs_pressure", input.sccs_pressure}, + {"sccs_tol_rms", input.sccs_tol_rms}, + {"sccs_tol_max", input.sccs_tol_max}, + {"sccs_surface_eta", input.sccs_surface_eta}, + {"sccs_lowpass_p1", input.sccs_lowpass_p1}, + {"sccs_lowpass_p2", input.sccs_lowpass_p2}}; + for (const std::pair& value : values) + { + if (!std::isfinite(value.second)) + { + return std::string(value.first) + " must be finite"; + } + } + return std::string(); +} + +// Cavity and surface model parameters, starting with the preset. +std::string model_error(const Input_para& input) +{ + const std::vector presets = {"custom", "vacuum", "water-neutral", "water-cation", "water-anion"}; + if (std::find(presets.begin(), presets.end(), input.sccs_preset) == presets.end()) + { + return "Unknown sccs_preset"; + } + const bool p1_positive = input.sccs_lowpass_p1 > 0.0; + const bool p2_positive = input.sccs_lowpass_p2 > 0.0; + const bool pcc = input.assume_isolated == "pcc_0d" || input.assume_isolated == "pcc_2d"; + std::string error; + if (input.sccs_epsilon < 1.0) + { + error = "sccs_epsilon must be at least 1"; + } + else if (input.sccs_rho_min <= 0.0 || input.sccs_rho_max <= input.sccs_rho_min) + { + error = "sccs_rho_min and sccs_rho_max must satisfy 0 < sccs_rho_min < sccs_rho_max"; + } + else if (input.sccs_surface_eta <= 0.0) + { + error = "sccs_surface_eta must be positive"; + } + else if (p1_positive != p2_positive) + { + error = "sccs_lowpass_p1 and sccs_lowpass_p2 must be both positive or both non-positive"; + } + else if (p1_positive && !pcc) + { + error = "sccs_lowpass_p1 and sccs_lowpass_p2 require assume_isolated pcc_0d or pcc_2d"; + } + else if (input.sccs_solvent_mode != "electronic" && input.sccs_solvent_mode != "full") + { + error = "sccs_solvent_mode must be electronic or full"; + } + else if (input.sccs_corespread.empty()) + { + error = "sccs_corespread must contain one value or exactly nat values"; + } + return error; +} + +// sqrt-CG solver controls. +std::string solver_error(const Input_para& input) +{ + std::string error; + if (input.sccs_maxiter <= 0) + { + error = "sccs_maxiter must be positive"; + } + else if (input.sccs_tol_rms <= 0.0) + { + error = "sccs_tol_rms must be positive"; + } + else if (input.sccs_tol_max <= 0.0) + { + error = "sccs_tol_max must be positive"; + } + return error; +} +} // namespace + +bool check_sccs_inp(const Input_para& input, std::string& error) +{ + error.clear(); + if (input.imp_sol < 0 || input.imp_sol > 2) + { + error = "imp_sol must be 0, 1 or 2"; + return false; + } + if (input.imp_sol != 2) + { + return true; + } + error = context_error(input); + if (error.empty()) + { + error = nonfinite_error(input); + } + if (error.empty()) + { + error = model_error(input); + } + if (error.empty()) + { + error = solver_error(input); + } + return error.empty(); +} +} // namespace ModuleIO diff --git a/source/source_io/test_serial/CMakeLists.txt b/source/source_io/test_serial/CMakeLists.txt index 03226d0b02e..991c2d77821 100644 --- a/source/source_io/test_serial/CMakeLists.txt +++ b/source/source_io/test_serial/CMakeLists.txt @@ -18,6 +18,7 @@ add_library( ../module_parameter/read_inp_deepks.cpp ../module_parameter/read_inp_model.cpp ../module_parameter/read_inp_sccs.cpp + ../module_parameter/sccs_inp_check.cpp ../module_parameter/read_inp_postproc.cpp ../module_parameter/read_inp_exx_dftu.cpp ../module_parameter/read_inp_other.cpp diff --git a/source/source_io/test_serial/test_read_inp_sccs.cpp b/source/source_io/test_serial/test_read_inp_sccs.cpp index ea18435b482..37dcec861b2 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -2,6 +2,8 @@ #include "source_io/module_parameter/input_parameter.h" #include +#include +#include #include TEST(ReadInpSccs, LegacyBooleansAndModelSelection) @@ -67,3 +69,101 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) input.sccs_preset = "unused"; EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); } + +TEST(ReadInpSccs, LowpassRequiresPairedValuesAndPcc) +{ + Input_para input; + input.imp_sol = 2; + std::string error; + EXPECT_DOUBLE_EQ(input.sccs_lowpass_p1, -1.0); + EXPECT_DOUBLE_EQ(input.sccs_lowpass_p2, -1.0); + input.sccs_lowpass_p1 = 10.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_lowpass_p2 = 5.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.assume_isolated = "pcc_0d"; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.assume_isolated = "pcc_2d"; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_lowpass_p1 = 0.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_lowpass_p2 = 0.0; + input.assume_isolated = "none"; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); +} + +TEST(ReadInpSccs, FullCavityWidthsAndModeValidation) +{ + Input_para input; + input.imp_sol = 2; + input.sccs_solvent_mode = "full"; + std::string error; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_corespread = {0.7, 0.0, -1.0}; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_corespread = {0.0}; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_corespread.clear(); + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_corespread = {0.5}; + input.sccs_solvent_mode = "invalid"; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); +} + +TEST(ReadInpSccs, CoreSpreadWordsMustBeFiniteNumbers) +{ + std::vector spreads = {0.5}; + std::string error; + ASSERT_TRUE(ModuleIO::parse_spreads({"0.7", "-1", "1e-1"}, spreads, error)) << error; + ASSERT_EQ(spreads.size(), 3u); + EXPECT_DOUBLE_EQ(spreads[0], 0.7); + EXPECT_DOUBLE_EQ(spreads[1], -1.0); + EXPECT_DOUBLE_EQ(spreads[2], 0.1); + const std::vector invalid[] = {{"abc"}, {"0.5abc"}, {"0.5", "x"}, {"nan"}, {"inf"}, {""}}; + for (const std::vector& words : invalid) + { + EXPECT_FALSE(ModuleIO::parse_spreads(words, spreads, error)); + EXPECT_NE(error.find("sccs_corespread"), std::string::npos); + EXPECT_EQ(spreads.size(), 3u); + } +} + +// nan and inf pass the ordered range checks, so check_sccs_inp tests every real +// SCCS key for finiteness. Signed gamma and pressure stay valid, and imp_sol 0 +// and 1 skip the SCCS checks. +TEST(ReadInpSccs, RejectsNonfiniteReals) +{ + Input_para input; + input.imp_sol = 2; + input.device = "cpu"; + input.basis_type = "pw"; + input.sccs_gamma = -1.0; + input.sccs_pressure = -0.36; + std::string error; + ASSERT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; + + double* const reals[] = {&input.sccs_epsilon, &input.sccs_rho_min, &input.sccs_rho_max, + &input.sccs_gamma, &input.sccs_pressure, &input.sccs_tol_rms, + &input.sccs_tol_max, &input.sccs_surface_eta, &input.sccs_lowpass_p1, + &input.sccs_lowpass_p2}; + const double infinity = std::numeric_limits::infinity(); + const double nonfinite[] = {std::nan(""), infinity, -infinity}; + for (double* value : reals) + { + const double saved = *value; + for (double bad : nonfinite) + { + *value = bad; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + const std::size_t position = error.find("must be finite"); + EXPECT_NE(position, std::string::npos) << error; + } + *value = saved; + } + + input.sccs_surface_eta = std::nan(""); + input.imp_sol = 0; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; + input.imp_sol = 1; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; +} diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index e3958064132..94fcad72296 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -130,3 +130,4 @@ scf_pcc_2d scf_sccs scf_sccs_pcc_2d scf_sccs_pcc_0d +scf_sccs_full diff --git a/tests/01_PW/scf_sccs_full/INPUT b/tests/01_PW/scf_sccs_full/INPUT new file mode 100644 index 00000000000..80dff139de0 --- /dev/null +++ b/tests/01_PW/scf_sccs_full/INPUT @@ -0,0 +1,22 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated none +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +sccs_solvent_mode full diff --git a/tests/01_PW/scf_sccs_full/KPT b/tests/01_PW/scf_sccs_full/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/scf_sccs_full/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_full/README b/tests/01_PW/scf_sccs_full/README new file mode 100644 index 00000000000..8492a26b809 --- /dev/null +++ b/tests/01_PW/scf_sccs_full/README @@ -0,0 +1 @@ +PW neutral water with periodic SCCS and the full cavity (sccs_solvent_mode full, default corespread 0.5): the O valence Gaussian is added to the cavity density, H is skipped (Z = zv). ecutwfc 20 and ecutrho 240. diff --git a/tests/01_PW/scf_sccs_full/STRU b/tests/01_PW/scf_sccs_full/STRU new file mode 100644 index 00000000000..ed974dd43a5 --- /dev/null +++ b/tests/01_PW/scf_sccs_full/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +12.0 0.0 0.0 +0.0 12.0 0.0 +0.0 0.0 12.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +7.433 7.111 6.0 0 0 0 +4.567 7.111 6.0 0 0 0 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/scf_sccs_full/result.ref b/tests/01_PW/scf_sccs_full/result.ref new file mode 100644 index 00000000000..ee13708521d --- /dev/null +++ b/tests/01_PW/scf_sccs_full/result.ref @@ -0,0 +1,6 @@ +etotref -441.6484408867916613 +etotperatomref -147.2161469623 +totalforceref 16.644048 +esolelref -0.2882959320 +esolcavref 0.0807339547 +totaltimeref 1.17 diff --git a/tests/01_PW/scf_sccs_full/threshold b/tests/01_PW/scf_sccs_full/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/scf_sccs_full/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/scf_sccs_pcc_0d/INPUT b/tests/01_PW/scf_sccs_pcc_0d/INPUT index f287686f377..c50c86ca36f 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/INPUT +++ b/tests/01_PW/scf_sccs_pcc_0d/INPUT @@ -19,3 +19,5 @@ sccs_maxiter 200 sccs_tol_rms 1.0e-10 sccs_tol_max 1.0e-8 sccs_surface_eta 1.0e-8 +sccs_lowpass_p1 10 +sccs_lowpass_p2 5 diff --git a/tests/01_PW/scf_sccs_pcc_0d/README b/tests/01_PW/scf_sccs_pcc_0d/README index ebc67d35015..44a134bed92 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/README +++ b/tests/01_PW/scf_sccs_pcc_0d/README @@ -1 +1 @@ -PW neutral water with SCCS (water-neutral preset) and 0D PCC, ecutwfc 20 and ecutrho 240 +PW neutral water with SCCS (water-neutral preset), 0D PCC and the cavity low-pass filter (p1 10, p2 5), ecutwfc 20 and ecutrho 240 diff --git a/tests/01_PW/scf_sccs_pcc_0d/result.ref b/tests/01_PW/scf_sccs_pcc_0d/result.ref index 6b87966cdf4..a0057579af3 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/result.ref +++ b/tests/01_PW/scf_sccs_pcc_0d/result.ref @@ -1,7 +1,7 @@ -etotref -441.6509526968036994 -etotperatomref -147.2169842323 -totalforceref 16.638619 -esolelref -0.3079260694 -esolcavref 0.0807450240 -epccref 0.0166182392 -totaltimeref 4.31 +etotref -441.6282836551653759 +etotperatomref -147.2094278851 +totalforceref 16.655714 +esolelref -0.2903851652 +esolcavref 0.0807683679 +epccref 0.0170866295 +totaltimeref 3.56 From 5f305c8f7feda99e5eecb03f739658a7362c68dc Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Mon, 26 Oct 2026 09:51:46 +0800 Subject: [PATCH 2/4] Keep existing SCCS tests unchanged and add lowpass/full-cavity tests separately - Restore tests/01_PW/scf_sccs_pcc_0d to its Step 2 input and reference, so the default PCC0D path keeps its integration coverage. The lowpass run moves to a new case, scf_sccs_pcc_0d_lowpass (same input as before plus sccs_lowpass_p1 10 / sccs_lowpass_p2 5, same reference values). - Restore PotSccsTest.IonicForceAddsRydbergDerivativeOnce and PotSccsTest.InputMapsOntoSccsConfig; the full-cavity force check and the new key mapping move to FullCavityForceMatchesTheEnergyDerivative and InputMapsCavityModelKeys. - Restore CellTools.ExtractsPositionsMassesAndValenceChargesInAtomOrder; the atomic-number check moves to AtomDataCarriesPseudopotentialAtomicNumber. - SccsCavityDerivativesTest.OpenBoundaryCoefficientUsesSwitchingFunctionFft keeps only the adaptation to the CavityDerivatives interface; the lowpass check moves to OpenBoundaryLowpassFiltersSwitchingFunctionDerivatives. No change to source code or to reference values. --- source/source_cell/test/test_cell_tools.cpp | 22 ++++++-- source/source_estate/test/test_pot_sccs.cpp | 44 ++++++++++------ .../module_sccs/test/test_sccs_cav_deriv.cpp | 52 ++++++++++++++++++- tests/01_PW/CASES_CPU.txt | 1 + tests/01_PW/scf_sccs_pcc_0d/INPUT | 2 - tests/01_PW/scf_sccs_pcc_0d/README | 2 +- tests/01_PW/scf_sccs_pcc_0d/result.ref | 14 ++--- tests/01_PW/scf_sccs_pcc_0d_lowpass/INPUT | 23 ++++++++ tests/01_PW/scf_sccs_pcc_0d_lowpass/KPT | 4 ++ tests/01_PW/scf_sccs_pcc_0d_lowpass/README | 1 + tests/01_PW/scf_sccs_pcc_0d_lowpass/STRU | 24 +++++++++ .../01_PW/scf_sccs_pcc_0d_lowpass/result.ref | 7 +++ tests/01_PW/scf_sccs_pcc_0d_lowpass/threshold | 3 ++ 13 files changed, 168 insertions(+), 31 deletions(-) create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/INPUT create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/KPT create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/README create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/STRU create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/result.ref create mode 100644 tests/01_PW/scf_sccs_pcc_0d_lowpass/threshold diff --git a/source/source_cell/test/test_cell_tools.cpp b/source/source_cell/test/test_cell_tools.cpp index bd724de9dfd..75405c51a7f 100644 --- a/source/source_cell/test/test_cell_tools.cpp +++ b/source/source_cell/test/test_cell_tools.cpp @@ -96,13 +96,11 @@ TEST(CellTools, ExtractsPositionsMassesAndValenceChargesInAtomOrder) atoms[0].na = 2; atoms[0].mass = 12.0; atoms[0].ncpp.zv = 4.0; - atoms[0].ncpp.psd = " C "; atoms[0].tau = {ModuleBase::Vector3(0.1, 0.2, 0.3), ModuleBase::Vector3(0.4, 0.5, 0.6)}; atoms[1].na = 1; atoms[1].mass = 1.0; atoms[1].ncpp.zv = 1.0; - atoms[1].ncpp.psd = "H"; atoms[1].tau = {ModuleBase::Vector3(0.7, 0.8, 0.9)}; const std::vector data = unitcell::get_atom_data(atoms, 2, 10.0); @@ -122,9 +120,6 @@ TEST(CellTools, ExtractsPositionsMassesAndValenceChargesInAtomOrder) EXPECT_DOUBLE_EQ(data[0].valence_charge, 4.0); EXPECT_DOUBLE_EQ(data[1].valence_charge, 4.0); EXPECT_DOUBLE_EQ(data[2].valence_charge, 1.0); - EXPECT_EQ(data[0].atomic_number, 6); - EXPECT_EQ(data[1].atomic_number, 6); - EXPECT_EQ(data[2].atomic_number, 1); EXPECT_DOUBLE_EQ(atoms[0].tau[0].x, 0.1); } @@ -134,6 +129,23 @@ TEST(CellTools, AcceptsNoAtomTypes) EXPECT_TRUE(data.empty()); } +TEST(CellTools, AtomDataCarriesPseudopotentialAtomicNumber) +{ + Atom atoms[2]; + atoms[0].na = 2; + atoms[0].ncpp.psd = " C "; + atoms[0].tau = {ModuleBase::Vector3(), ModuleBase::Vector3()}; + atoms[1].na = 1; + atoms[1].ncpp.psd = "H"; + atoms[1].tau = {ModuleBase::Vector3()}; + + const std::vector 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; diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp index 10ef2a00a56..2e3ab0ac001 100644 --- a/source/source_estate/test/test_pot_sccs.cpp +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -185,8 +185,19 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) EXPECT_EQ(config.boundary, ModuleSccs::Boundary::Pcc2d); EXPECT_EQ(config.pcc_2d_axis, axis); } +} + +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"; @@ -196,12 +207,11 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) 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); } -// Electrostatic forces in a dielectric, and the cavity terms of core -// Gaussians with surface tension and pressure. -TEST_F(PotSccsTest, SolvationForceMatchesTheEnergyDerivative) +TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) { Input_para input; input.sccs_epsilon = 5.0; @@ -211,19 +221,23 @@ TEST_F(PotSccsTest, SolvationForceMatchesTheEnergyDerivative) ModuleSccs::PolarizationSolverParameters solver; elecstate::sccs_cfg_input(input, config, solver); check_solvation_force(config, solver); +} - Input_para full_input; - full_input.sccs_epsilon = 1.0; - full_input.sccs_tol_rms = 1e-13; - full_input.sccs_tol_max = 1e-12; - ModuleSccs::SccsConfig full; - ModuleSccs::PolarizationSolverParameters full_solver; - elecstate::sccs_cfg_input(full_input, full, full_solver); - full.core_electrons = true; - full.core_spreads = {2.0}; - full.surface_tension = 1e-5; - full.pressure = 1e-6; - check_solvation_force(full, full_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) diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp index bb941d08b90..d23ef2d0374 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp @@ -83,10 +83,60 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf Parallel_Reduce::reduce_max_pool(basis.poolnproc, difference); EXPECT_GT(difference, 1e-8); +} + +namespace +{ +// Density whose switching function is 0.5 + amplitude * mode, found by +// bisection of the cavity function at every grid point. +std::vector shell_density(const ModuleSccs::CavityParameters& cavity, + const std::vector& mode, + double amplitude) +{ + std::vector density(mode.size()); + for (std::size_t ir = 0; ir < mode.size(); ++ir) + { + const double solute = 0.5 + amplitude * mode[ir]; + double lower = cavity.density_min; + double upper = cavity.density_max; + for (int iteration = 0; iteration < 64; ++iteration) + { + const double middle = 0.5 * (lower + upper); + const ModuleSccs::CavityPoint point = ModuleSccs::evaluate_cavity(middle, cavity); + if (point.solute < solute) + { + lower = middle; + } + else + { + upper = middle; + } + } + density[ir] = 0.5 * (lower + upper); + } + return density; +} +} // namespace + +// The lowpass multiplies the single cosine shell by f(G) = erfc(G^2 - 0.5)/2 +// = erfc(0.5)/2 before the FFT derivatives of the switching function. +TEST_F(SccsCavityDerivativesTest, OpenBoundaryLowpassFiltersSwitchingFunctionDerivatives) +{ + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; cavity.lowpass_p1 = basis.ggecut; cavity.lowpass_p2 = 0.5; + const std::vector mode = cosine_mode(0); + const double amplitude = 0.45; + const std::vector density = shell_density(cavity, mode, amplitude); + const bool open_boundary = true; + ModuleSccs::SccsResponse response; + ModuleSccs::CavityDerivatives derivatives; ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, response, derivatives); + const double log_bulk = std::log(cavity.epsilon_bulk); const double weight = 0.5 * std::erfc(0.5); for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -97,7 +147,7 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf const double expected = response.epsilon[ir] * (-0.5 * log_bulk * laplacian + 0.25 * log_bulk * log_bulk * gradient * gradient) / ModuleBase::FOUR_PI; - EXPECT_NEAR(coefficient[ir], expected, 1e-12); + EXPECT_NEAR(derivatives.coefficient[ir], expected, 1e-12); EXPECT_NEAR(derivatives.gradient[ir].x, gradient, 1e-12); } } diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 94fcad72296..0b97a367812 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -130,4 +130,5 @@ scf_pcc_2d scf_sccs scf_sccs_pcc_2d scf_sccs_pcc_0d +scf_sccs_pcc_0d_lowpass scf_sccs_full diff --git a/tests/01_PW/scf_sccs_pcc_0d/INPUT b/tests/01_PW/scf_sccs_pcc_0d/INPUT index c50c86ca36f..f287686f377 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/INPUT +++ b/tests/01_PW/scf_sccs_pcc_0d/INPUT @@ -19,5 +19,3 @@ sccs_maxiter 200 sccs_tol_rms 1.0e-10 sccs_tol_max 1.0e-8 sccs_surface_eta 1.0e-8 -sccs_lowpass_p1 10 -sccs_lowpass_p2 5 diff --git a/tests/01_PW/scf_sccs_pcc_0d/README b/tests/01_PW/scf_sccs_pcc_0d/README index 44a134bed92..ebc67d35015 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/README +++ b/tests/01_PW/scf_sccs_pcc_0d/README @@ -1 +1 @@ -PW neutral water with SCCS (water-neutral preset), 0D PCC and the cavity low-pass filter (p1 10, p2 5), ecutwfc 20 and ecutrho 240 +PW neutral water with SCCS (water-neutral preset) and 0D PCC, ecutwfc 20 and ecutrho 240 diff --git a/tests/01_PW/scf_sccs_pcc_0d/result.ref b/tests/01_PW/scf_sccs_pcc_0d/result.ref index a0057579af3..6b87966cdf4 100644 --- a/tests/01_PW/scf_sccs_pcc_0d/result.ref +++ b/tests/01_PW/scf_sccs_pcc_0d/result.ref @@ -1,7 +1,7 @@ -etotref -441.6282836551653759 -etotperatomref -147.2094278851 -totalforceref 16.655714 -esolelref -0.2903851652 -esolcavref 0.0807683679 -epccref 0.0170866295 -totaltimeref 3.56 +etotref -441.6509526968036994 +etotperatomref -147.2169842323 +totalforceref 16.638619 +esolelref -0.3079260694 +esolcavref 0.0807450240 +epccref 0.0166182392 +totaltimeref 4.31 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/INPUT b/tests/01_PW/scf_sccs_pcc_0d_lowpass/INPUT new file mode 100644 index 00000000000..c50c86ca36f --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/INPUT @@ -0,0 +1,23 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated pcc_0d +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +sccs_lowpass_p1 10 +sccs_lowpass_p2 5 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/KPT b/tests/01_PW/scf_sccs_pcc_0d_lowpass/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/README b/tests/01_PW/scf_sccs_pcc_0d_lowpass/README new file mode 100644 index 00000000000..44a134bed92 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/README @@ -0,0 +1 @@ +PW neutral water with SCCS (water-neutral preset), 0D PCC and the cavity low-pass filter (p1 10, p2 5), ecutwfc 20 and ecutrho 240 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/STRU b/tests/01_PW/scf_sccs_pcc_0d_lowpass/STRU new file mode 100644 index 00000000000..ed974dd43a5 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +12.0 0.0 0.0 +0.0 12.0 0.0 +0.0 0.0 12.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +7.433 7.111 6.0 0 0 0 +4.567 7.111 6.0 0 0 0 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/result.ref b/tests/01_PW/scf_sccs_pcc_0d_lowpass/result.ref new file mode 100644 index 00000000000..a0057579af3 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/result.ref @@ -0,0 +1,7 @@ +etotref -441.6282836551653759 +etotperatomref -147.2094278851 +totalforceref 16.655714 +esolelref -0.2903851652 +esolcavref 0.0807683679 +epccref 0.0170866295 +totaltimeref 3.56 diff --git a/tests/01_PW/scf_sccs_pcc_0d_lowpass/threshold b/tests/01_PW/scf_sccs_pcc_0d_lowpass/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d_lowpass/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 From 8c35761fc729a414382f41fe2bb12ab71b01f3fd Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Mon, 26 Oct 2026 14:05:51 +0800 Subject: [PATCH 3/4] Disable spatial symmetry for per-atom SCCS cavity widths; keep the cavity-density derivative only for the full cavity - A list of nat values for sccs_corespread can give symmetry-equivalent atoms different core Gaussians, which the symmetry analysis does not see. The symmetry reset now lowers symmetry 1 to 0 for imp_sol 2 with sccs_solvent_mode full and more than one width, as it already does for efield. A single width is chosen per element and stays compatible with symmetry. The sccs_corespread help states this. - eval_functional copies and accumulates the cavity-density derivative only when the cavity has core Gaussians; its only reader is the core-Gaussian force of the full cavity. The electronic mode no longer keeps this grid. - Tests: InputTest.PerAtomSccsCavityWidthsDisableSpatialSymmetry and SccsFunctionalTest.CavityDensityDerivativeOnlyForTheFullCavity. --- .../module_sccs/sccs_functional.cpp | 13 ++++++- .../module_sccs/sccs_functional.h | 2 +- .../module_sccs/test/test_sccs_functional.cpp | 39 +++++++++++++++++++ .../module_parameter/read_inp_sccs.cpp | 2 +- .../module_parameter/read_inp_sys.cpp | 7 ++++ .../test_serial/read_input_item_test.cpp | 27 +++++++++++++ 6 files changed, 86 insertions(+), 4 deletions(-) diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp index d56867b5247..a2b28a516b1 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -39,7 +39,13 @@ void eval_functional(const std::vector& charge, FunctionalResult candidate; candidate.reaction_potential.resize(size); candidate.electron_potential.resize(size); - candidate.cavity_potential = response.cavity_potential; + // The cavity-density derivative is needed only for the forces on the core + // Gaussians of the full cavity. + const bool core_electrons = config.core_electrons; + if (core_electrons) + { + candidate.cavity_potential = response.cavity_potential; + } for (std::size_t i = 0; i < size; ++i) { const double reaction = response.polarization.potential[i] - vacuum[i]; @@ -74,7 +80,10 @@ void eval_functional(const std::vector& charge, const double nonel = (config.pressure - config.surface_tension * divergence[i]) * response.dsolute_drho[i]; candidate.electron_potential[i] += nonel; - candidate.cavity_potential[i] += nonel; + if (core_electrons) + { + candidate.cavity_potential[i] += nonel; + } } result = std::move(candidate); } diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h index 249379dd913..8912ec57f68 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.h +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -24,7 +24,7 @@ struct FunctionalResult double volume = 0.0; // Bohr^3 std::vector reaction_potential; // dielectric minus vacuum, Hartree std::vector electron_potential; // derivative of all three energy terms - std::vector cavity_potential; // derivative w.r.t. cavity density, Hartree + std::vector cavity_potential; // derivative w.r.t. cavity density, Hartree; full cavity only }; // Collective over the PW pool. diff --git a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp index a6c91db1d84..1d335c5d0ec 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -121,3 +121,42 @@ TEST_F(SccsFunctionalTest, ExplicitPeriodicOperatorPreservesResponseAndFunctiona EXPECT_EQ(legacy_result.reaction_potential, explicit_result.reaction_potential); EXPECT_EQ(legacy_result.electron_potential, explicit_result.electron_potential); } + +// The cavity-density derivative is produced only for the core-Gaussian forces +// of the full cavity; the electron potential does not depend on it. +TEST_F(SccsFunctionalTest, CavityDensityDerivativeOnlyForTheFullCavity) +{ + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 * (1.0 + 0.3 * mode[ir]); + charge[ir] = -density[ir]; + } + ModuleSccs::SccsConfig config; + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + config.cavity.epsilon_bulk = 5.0; + config.surface_regularization = 1e-8; + config.surface_tension = 1e-5; + config.pressure = 1e-6; + ModuleSccs::PolarizationSolverParameters solver; + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + basis, tpiba, response); + ModuleSccs::FunctionalResult electronic; + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, electronic); + EXPECT_TRUE(electronic.cavity_potential.empty()); + config.core_electrons = true; + ModuleSccs::FunctionalResult full; + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, full); + ASSERT_EQ(full.cavity_potential.size(), static_cast(basis.nrxx)); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_EQ(full.electron_potential[ir], electronic.electron_potential[ir]); + const double expected = full.electron_potential[ir] + full.reaction_potential[ir]; + EXPECT_NEAR(full.cavity_potential[ir], expected, 1e-14); + } +} diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 06e6a1dcf1a..2c8761c9a55 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -168,7 +168,7 @@ void ReadInput::item_sccs() item.annotation = "Widths of the full-cavity Gaussians"; item.category = "Implicit solvation model"; item.type = "Vector of Real"; - item.description = "Widths of the valence-charge Gaussians used only with sccs_solvent_mode full. One value applies to every atom except those whose known atomic number equals the pseudopotential valence charge (within 1e-8). More than one value requires exactly nat values in STRU atom order (grouped by atom type), and overrides this automatic exclusion. A value <= 0 disables the cavity Gaussian on that atom."; + item.description = "Widths of the valence-charge Gaussians used only with sccs_solvent_mode full. One value applies to every atom except those whose known atomic number equals the pseudopotential valence charge (within 1e-8). More than one value requires exactly nat values in STRU atom order (grouped by atom type), and overrides this automatic exclusion. A value <= 0 disables the cavity Gaussian on that atom. A per-atom list sets symmetry 1 to 0, because the symmetry analysis does not distinguish equivalent atoms by their widths."; item.default_value = "0.5"; item.unit = "bohr"; item.set_availability("imp_sol==2"); diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 70e69d9736f..43adbfc8169 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -371,6 +371,13 @@ Socket mode always computes energy. Force and stress extraction follows cal_forc { para.input.symmetry = "-1"; // disable kpoint reduce } + // Per-atom SCCS cavity widths can differ between equivalent atoms, + // which the symmetry analysis does not see. + const bool sccs_full = para.input.imp_sol == 2 && para.input.sccs_solvent_mode == "full"; + if (sccs_full && para.input.sccs_corespread.size() > 1 && para.input.symmetry == "1") + { + para.input.symmetry = "0"; + } }; this->add_item(item); } diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 3b4365dbff3..b1221ff0b77 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -2789,3 +2789,30 @@ TEST_F(InputTest, LRFixedFieldRejectsForcesAndRelaxation) check_item(it->second, param); } } + +// Per-atom full-cavity widths can differ between equivalent atoms, which the +// symmetry analysis does not see, so symmetry 1 is reduced to 0. +TEST_F(InputTest, PerAtomSccsCavityWidthsDisableSpatialSymmetry) +{ + ModuleIO::ReadInput readinput(0); + Parameter param; + auto it = find_label("symmetry", readinput.input_lists); + ASSERT_NE(it, readinput.input_lists.end()); + Input_para& input = TestParameters::input(param); + input.imp_sol = 2; + input.sccs_solvent_mode = "full"; + input.sccs_corespread = {0.5}; + input.symmetry = "1"; + it->second.reset_value(it->second, param); + EXPECT_EQ(input.symmetry, "1"); + input.sccs_corespread = {0.5, 0.0}; + it->second.reset_value(it->second, param); + EXPECT_EQ(input.symmetry, "0"); + input.symmetry = "-1"; + it->second.reset_value(it->second, param); + EXPECT_EQ(input.symmetry, "-1"); + input.symmetry = "1"; + input.sccs_solvent_mode = "electronic"; + it->second.reset_value(it->second, param); + EXPECT_EQ(input.symmetry, "1"); +} From 81cfdc530f3c76dd33996fe38ae2f95066ab2b69 Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Mon, 26 Oct 2026 16:32:47 +0800 Subject: [PATCH 4/4] Select the lowpass cavity derivative only with open boundaries finish_response chose lowpass_cav_pot from uses_lowpass(cavity) alone, but make_cav_derivs keeps the filtered gradient that it reads only when the boundary is open as well. Use open_boundary && uses_lowpass(cavity) in both. Input validation already rejects lowpass without PCC, so results do not change; the path no longer depends on that check to stay in bounds. Test: SccsLowpassTest.PeriodicResponseIgnoresTheLowpass sets the lowpass parameters on a periodic cavity in code and checks that the response is unchanged; without the fix it reads the empty filtered gradient and crashes. --- .../module_sccs/sccs_response.cpp | 3 +- .../module_sccs/test/test_sccs_lowpass.cpp | 33 +++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index 9904db6d915..6a94f5593bf 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -244,7 +244,8 @@ void finish_response(const SqrtCgOperator& system, value -= mean; } } - const bool lowpass = uses_lowpass(cavity); + // The same condition under which make_cav_derivs keeps the filtered gradient. + const bool lowpass = open_boundary && uses_lowpass(cavity); if (!lowpass) { std::vector> potential_g(basis.npw); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp index dd202dabf41..cc69a79a966 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp @@ -144,3 +144,36 @@ TEST_F(SccsLowpassTest, FilterUsesDensityCutoff) } } } + +// The periodic chain path takes no FFT derivative of s, so lowpass parameters +// set in code (the INPUT reader rejects them without PCC) leave the periodic +// response unchanged. +TEST_F(SccsLowpassTest, PeriodicResponseIgnoresTheLowpass) +{ + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + ModuleSccs::PolarizationSolverParameters solver; + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 + 1e-4 * mode[ir]; + charge[ir] = 2e-4 * mode[ir]; + } + const std::vector cold; + ModuleSccs::SccsResponse plain; + ModuleSccs::solve_response(density, charge, cavity, solver, cold, basis, tpiba, plain); + cavity.lowpass_p1 = basis.ggecut; + cavity.lowpass_p2 = 0.5; + ModuleSccs::SccsResponse filtered; + ModuleSccs::solve_response(density, charge, cavity, solver, cold, basis, tpiba, filtered); + EXPECT_FALSE(filtered.polarization.gradient.empty()); + ASSERT_EQ(filtered.cavity_potential.size(), plain.cavity_potential.size()); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_EQ(filtered.cavity_potential[ir], plain.cavity_potential[ir]); + } +}