From bb7cb74c237265d1578799659d3d719aff1e6f0d Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Fri, 9 Oct 2026 16:10:03 -0600 Subject: [PATCH 1/5] Feature: periodic self-consistent continuum solvation (SCCS) --- source/CMakeLists.txt | 1 + source/Makefile.Objects | 12 + source/source_estate/CMakeLists.txt | 1 + source/source_estate/elecstate.h | 2 +- source/source_estate/elecstate_energy.cpp | 13 +- source/source_estate/module_pot/pot_base.h | 13 +- source/source_estate/module_pot/pot_sccs.cpp | 157 +++++++++ source/source_estate/module_pot/pot_sccs.h | 49 +++ .../module_pot/potential_new.cpp | 39 ++- .../source_estate/module_pot/potential_new.h | 10 +- .../module_pot/potential_types.cpp | 25 +- .../test/elecstate_base_test.cpp | 4 + .../test/elecstate_energy_test.cpp | 32 +- .../source_estate/test/elecstate_pw_test.cpp | 4 + .../source_estate/test/potential_new_test.cpp | 34 +- source/source_estate/test/test_pot_sccs.cpp | 244 ++++++++++++++ source/source_estate/write_elecstat_pot.cpp | 14 +- source/source_estate/write_elecstat_pot.h | 1 + source/source_hamilt/CMakeLists.txt | 1 + .../source_hamilt/module_sccs/CMakeLists.txt | 16 + .../source_hamilt/module_sccs/sccs_cavity.cpp | 39 +++ .../source_hamilt/module_sccs/sccs_cavity.h | 27 ++ .../module_sccs/sccs_functional.cpp | 68 ++++ .../module_sccs/sccs_functional.h | 36 ++ .../module_sccs/sccs_ionic_charge.cpp | 33 ++ .../module_sccs/sccs_ionic_charge.h | 29 ++ .../module_sccs/sccs_ionic_force.cpp | 45 +++ .../module_sccs/sccs_ionic_force.h | 30 ++ .../module_sccs/sccs_parameters.cpp | 79 +++++ .../module_sccs/sccs_parameters.h | 36 ++ .../module_sccs/sccs_pw_coulomb.cpp | 32 ++ .../module_sccs/sccs_pw_coulomb.h | 31 ++ .../module_sccs/sccs_response.cpp | 316 ++++++++++++++++++ .../source_hamilt/module_sccs/sccs_response.h | 56 ++++ .../module_sccs/test/CMakeLists.txt | 27 ++ .../module_sccs/test/sccs_test.h | 60 ++++ .../module_sccs/test/test_sccs_cavity.cpp | 110 ++++++ .../module_sccs/test/test_sccs_functional.cpp | 84 +++++ .../test/test_sccs_ionic_charge.cpp | 53 +++ .../test/test_sccs_ionic_force.cpp | 97 ++++++ .../module_sccs/test/test_sccs_parameters.cpp | 43 +++ .../module_sccs/test/test_sccs_pw_coulomb.cpp | 68 ++++ .../module_sccs/test/test_sccs_response.cpp | 149 +++++++++ .../test/hsolver_supplementary_mock.h | 15 + source/source_io/CMakeLists.txt | 1 + .../source_io/module_ctrl/ctrl_output_fp.cpp | 2 + .../module_parameter/input_parameter.h | 13 +- .../module_parameter/read_inp_model.cpp | 27 +- .../module_parameter/read_inp_sccs.cpp | 220 ++++++++++++ .../module_parameter/read_inp_sccs.h | 10 + .../source_io/module_parameter/read_input.h | 1 + source/source_io/test_serial/CMakeLists.txt | 3 +- .../test_serial/test_read_inp_sccs.cpp | 56 ++++ source/source_lcao/force_stress_lcao.cpp | 2 + source/source_lcao/force_stress_terms.cpp | 2 +- source/source_lcao/hamilt_lcao.cpp | 3 +- source/source_lcao/hamilt_lcao_factory.cpp | 4 +- source/source_pw/module_pwdft/force_pw.cpp | 16 +- source/source_pw/module_pwdft/hamilt_pw.cpp | 5 +- .../module_stodft/test/test_sto_hamilt_pw.cpp | 3 + tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT | 21 ++ tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT | 4 + tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README | 2 + tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU | 24 ++ .../220_PW_SCCS_PERIODIC_WATER/result.ref | 6 + .../220_PW_SCCS_PERIODIC_WATER/threshold | 3 + tests/01_PW/CASES_CPU.txt | 1 + .../integrate/validation_tools/props_basic.sh | 2 +- 68 files changed, 2630 insertions(+), 36 deletions(-) create mode 100644 source/source_estate/module_pot/pot_sccs.cpp create mode 100644 source/source_estate/module_pot/pot_sccs.h create mode 100644 source/source_estate/test/test_pot_sccs.cpp create mode 100644 source/source_hamilt/module_sccs/CMakeLists.txt create mode 100644 source/source_hamilt/module_sccs/sccs_cavity.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_cavity.h create mode 100644 source/source_hamilt/module_sccs/sccs_functional.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_functional.h create mode 100644 source/source_hamilt/module_sccs/sccs_ionic_charge.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_ionic_charge.h create mode 100644 source/source_hamilt/module_sccs/sccs_ionic_force.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_ionic_force.h create mode 100644 source/source_hamilt/module_sccs/sccs_parameters.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_parameters.h create mode 100644 source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_pw_coulomb.h create mode 100644 source/source_hamilt/module_sccs/sccs_response.cpp create mode 100644 source/source_hamilt/module_sccs/sccs_response.h create mode 100644 source/source_hamilt/module_sccs/test/CMakeLists.txt create mode 100644 source/source_hamilt/module_sccs/test/sccs_test.h create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_functional.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp create mode 100644 source/source_hamilt/module_sccs/test/test_sccs_response.cpp create mode 100644 source/source_io/module_parameter/read_inp_sccs.cpp create mode 100644 source/source_io/module_parameter/read_inp_sccs.h create mode 100644 source/source_io/test_serial/test_read_inp_sccs.cpp create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref create mode 100644 tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 2322dda4976..910fbe15b4d 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -744,6 +744,7 @@ target_link_libraries( io_basic io_input surchem + sccs neighbor_search domain_decomposition neighbor diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 43a4072889c..81c1dcc20c4 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -49,6 +49,7 @@ VPATH=./src_global:\ ./source_hamilt:\ ./source_hamilt/module_ewald:\ ./source_hamilt/module_surchem:\ +./source_hamilt/module_sccs:\ ./source_hamilt/module_xc:\ ./source_hamilt/module_xc/kernels:\ ./source_pw/module_proj:\ @@ -103,6 +104,7 @@ ${OBJS_PSI_INITIALIZER}\ ${OBJS_PW}\ ${OBJS_RELAXATION}\ ${OBJS_SURCHEM}\ +${OBJS_SCCS}\ ${OBJS_SYMMETRY}\ ${OBJS_XC}\ ${OBJS_IO}\ @@ -278,6 +280,7 @@ OBJS_ELECSTAT=elecstate.o\ pcc_0d.o\ pcc_2d.o\ pot_pcc.o\ + pot_sccs.o\ elecstate_op.o\ efield.o\ gatefield.o\ @@ -620,6 +623,14 @@ OBJS_SURCHEM=surchem.o\ minimize_cg.o\ sol_force.o\ +OBJS_SCCS=sccs_cavity.o\ + sccs_functional.o\ + sccs_parameters.o\ + sccs_ionic_charge.o\ + sccs_ionic_force.o\ + sccs_pw_coulomb.o\ + sccs_response.o + OBJS_SYMMETRY=symm_other.o\ symm_rot_spin.o\ symm_basic.o\ @@ -746,6 +757,7 @@ OBJS_IO=module_parameter/input_conv.o\ module_parameter/read_inp_bse.o\ module_parameter/read_inp_deepks.o\ module_parameter/read_inp_model.o\ + module_parameter/read_inp_sccs.o\ module_parameter/read_inp_postproc.o\ module_parameter/read_inp_exx_dftu.o\ module_parameter/read_inp_dfpt.o\ diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index 15c09ac3ce1..6be568facfa 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -21,6 +21,7 @@ list(APPEND objects module_pot/pot_xc.cpp module_pot/pot_local.cpp module_pot/pot_pcc.cpp + module_pot/pot_sccs.cpp module_pot/potential_new.cpp module_pot/potential_types.cpp module_pot/pot_surchem.cpp diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 68d637a4d4f..ba96ff14576 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -157,7 +157,7 @@ class ElecState /// @param dft_plus_u 0 disables the DFT+U energy term /// @param assume_isolated "makov-payne" adds the isolated-system correction void cal_energies(const int type, - const bool imp_sol, + const int imp_sol, const bool sc_mag_switch, const int dft_plus_u, const std::string& assume_isolated); diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp index 043a7057d8c..af77832329c 100644 --- a/source/source_estate/elecstate_energy.cpp +++ b/source/source_estate/elecstate_energy.cpp @@ -316,7 +316,7 @@ void ElecState::cal_converged() * @param type: 2 means Kohn-Sham functional; */ void ElecState::cal_energies(const int type, - const bool imp_sol, + const int imp_sol, const bool sc_mag_switch, const int dft_plus_u, const std::string& assume_isolated) @@ -331,11 +331,20 @@ void ElecState::cal_energies(const int type, this->f_en.gatefield = get_etot_gatefield(); //! energy from implicit solvation model - if (imp_sol) + if (imp_sol == 1) { this->f_en.esol_el = get_solvent_model_Ael(); this->f_en.esol_cav = get_solvent_model_Acav(); } + else if (imp_sol == 2) + { + this->pot->get_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav); + } + else + { + this->f_en.esol_el = 0.0; + this->f_en.esol_cav = 0.0; + } //! spin constrained energy if (sc_mag_switch) diff --git a/source/source_estate/module_pot/pot_base.h b/source/source_estate/module_pot/pot_base.h index aaa578bc5bf..658962f276e 100644 --- a/source/source_estate/module_pot/pot_base.h +++ b/source/source_estate/module_pot/pot_base.h @@ -28,7 +28,18 @@ class PotBase virtual void cal_fixed_v(double* vl_pseudo){} - virtual double get_energy() const { return 0.0; } + virtual double get_energy() const { return 0.0; } + + // Solvation contributions in Ry, added once after the band-energy potential subtraction. + virtual void get_solvation_energy(double& electrostatic, double& non_electrostatic) const + { + electrostatic = 0.0; + non_electrostatic = 0.0; + } + // Explicit ionic solvation force in Ry/Bohr, added to existing contributions. + virtual void add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} + // Electronic electrostatic correction only; excludes cavity/surface derivatives. + virtual const std::vector* solvent_electrostatic_potential() const { return nullptr; } bool fixed_mode = 0; bool dynamic_mode = 0; diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp new file mode 100644 index 00000000000..7d5dd323b19 --- /dev/null +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -0,0 +1,157 @@ +#include "pot_sccs.h" + +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_base/tool_quit.h" +#include "source_cell/cell_tools.h" +#include "source_hamilt/module_sccs/sccs_functional.h" +#include "source_hamilt/module_sccs/sccs_ionic_charge.h" +#include "source_hamilt/module_sccs/sccs_ionic_force.h" +#include "source_hamilt/module_sccs/sccs_response.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include +#include + +namespace elecstate +{ +void make_sccs_config_from_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver) +{ + const ModuleSccs::Preset preset = ModuleSccs::parse_preset(input.sccs_preset); + if (preset == ModuleSccs::Preset::Custom) + { + config = ModuleSccs::SccsConfig(); + config.cavity.density_min = input.sccs_rho_min; + config.cavity.density_max = input.sccs_rho_max; + config.cavity.epsilon_bulk = input.sccs_epsilon; + config.surface_tension = ModuleSccs::dyn_per_cm_to_hartree_per_bohr2(input.sccs_gamma); + config.pressure = ModuleSccs::gpa_to_hartree_per_bohr3(input.sccs_pressure); + } + else + { + config = ModuleSccs::make_sccs_config(preset); + } + config.surface_regularization = input.sccs_surface_eta; + solver.max_iterations = input.sccs_maxiter; + solver.tolerance_rms = input.sccs_tol_rms; + solver.tolerance_max = input.sccs_tol_max; +} + +std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count) +{ + double ionic_charge = 0.0; + for (int it = 0; it < cell.ntype; ++it) + { + ionic_charge += cell.atoms[it].ncpp.zv * cell.atoms[it].na; + } + const double net_charge = ionic_charge - electron_count; + if (config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= 1e-6) { return std::string(); } + return "charged SCCS with periodic boundaries: the periodic Poisson solver drops the G = 0 component " + "of the net charge, so the energy depends on the cell size."; +} + +PotSccs::PotSccs(const ModulePW::PW_Basis* basis, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver, + double electron_count) + : config_(config), solver_(solver), electron_count_(electron_count) +{ + this->rho_basis_ = basis; + this->dynamic_mode = true; +} + +void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) +{ + ModuleBase::timer::start("PotSccs", "cal_v_eff"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + const std::vector atoms = unitcell::get_atom_data(cell->atoms, cell->ntype, cell->lat0); + std::vector ions; + ModuleSccs::gaussian_ionic_density(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions); + std::vector density(basis.nrxx, 0.0); + std::vector solute_charge(basis.nrxx); + double grid_charges[2] = {0.0, 0.0}; // electrons, ions + for (int ir = 0; ir < basis.nrxx; ++ir) + { + for (int spin = 0; spin < charge->nspin; ++spin) { density[ir] += charge->rho[spin][ir]; } + solute_charge[ir] = ions[ir] - density[ir]; + grid_charges[0] += density[ir]; + grid_charges[1] += ions[ir]; + } + // Pool-reduced, so every rank takes the same decision. + Parallel_Reduce::reduce_pool(grid_charges, 2); + const double dv = basis.omega / basis.nxyz; + const double electron_error = grid_charges[0] * dv - electron_count_; + double valence_charge = 0.0; + for (const unitcell::AtomData& atom : atoms) { valence_charge += atom.valence_charge; } + const double ionic_error = grid_charges[1] * dv - valence_charge; + if (std::abs(ionic_error) > 1e-6) + { + ModuleBase::WARNING_QUIT("PotSccs::cal_v_eff", "SCCS ionic density normalization does not match the valence charge"); + } + if (std::abs(electron_error) > 1e-6) + { + std::ostringstream message; + message << "SCCS grid electron count differs from the expected value by " << electron_error + << " e; the solute charge uses the grid density"; + ModuleBase::WARNING("PotSccs::cal_v_eff", message.str()); + } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, + cell->tpiba, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::evaluate_functional(solute_charge, response, config_, basis, cell->tpiba, functional); + electrostatic_rydberg_ = 2.0 * functional.reaction_energy; + non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); + electrostatic_potential_.resize(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + electrostatic_potential_[ir] = -2.0 * functional.reaction_potential[ir]; + const double value = 2.0 * functional.electron_potential[ir]; + for (int spin = 0; spin < charge->nspin; ++spin) { potential(spin, ir) += value; } + } + restart_potential_ = std::move(response.restart_potential); + ModuleBase::timer::end("PotSccs", "cal_v_eff"); +} + +void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + ModuleBase::timer::start("PotSccs", "add_solvation_force"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::vector reaction(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + reaction[ir] = -0.5 * electrostatic_potential_[ir]; + } + std::vector> ionic_force; + ModuleSccs::gaussian_ionic_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, + ionic_force); + for (int ia = 0; ia < cell.nat; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + force(ia, axis) += 2.0 * ionic_force[ia][axis]; + } + } + ModuleBase::timer::end("PotSccs", "add_solvation_force"); +} + +double PotSccs::get_energy() const +{ + return electrostatic_rydberg_ + non_electrostatic_rydberg_; +} + +void PotSccs::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = electrostatic_rydberg_; + non_electrostatic = non_electrostatic_rydberg_; +} + +const std::vector* PotSccs::solvent_electrostatic_potential() const +{ + return &electrostatic_potential_; +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h new file mode 100644 index 00000000000..00373447a83 --- /dev/null +++ b/source/source_estate/module_pot/pot_sccs.h @@ -0,0 +1,49 @@ +#ifndef POT_SCCS_H +#define POT_SCCS_H + +#include "pot_base.h" +#include "source_hamilt/module_sccs/sccs_parameters.h" + +#include + +struct Input_para; +namespace elecstate +{ +// Map INPUT values, already validated by the INPUT reader, onto the solver. +void make_sccs_config_from_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver); + +// Returns a warning (empty if none) for a charged cell in a dielectric: the +// periodic Poisson solver drops the G = 0 component of the net charge, so the +// energy depends on the cell size. electron_count is the resolved nelec (the +// INPUT value, or the valence count when it is not set). +std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count); + +class PotSccs : public PotBase +{ +public: + // electron_count is the expected number of electrons (resolved nelec); a grid + // density that differs from it by more than 1e-6 e is reported in the + // warning log, and the solute charge uses the grid density. + PotSccs(const ModulePW::PW_Basis* basis, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver, + double electron_count); + void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override; + double get_energy() const override; + void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; + void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; + const std::vector* solvent_electrostatic_potential() const override; + +private: + const ModuleSccs::SccsConfig config_; + const ModuleSccs::PolarizationSolverParameters solver_; + const double electron_count_; + double electrostatic_rydberg_ = 0.0; + double non_electrostatic_rydberg_ = 0.0; + std::vector electrostatic_potential_; + std::vector restart_potential_; +}; +} // namespace elecstate +#endif diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index b23b8710bb5..51dbced99f1 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -62,6 +62,11 @@ Potential::~Potential() } void Potential::pot_register(const std::vector& components_list) +{ + this->pot_register(components_list, PARAM.inp); +} + +void Potential::pot_register(const std::vector& components_list, const Input_para& input) { ModuleBase::TITLE("Potential", "pot_register"); // delete old components first. @@ -80,7 +85,7 @@ void Potential::pot_register(const std::vector& components_list) //--------------------------- for (auto comp: components_list) { - PotBase* tmp = this->get_pot_type(comp); + PotBase* tmp = this->get_pot_type(comp, input); this->components.push_back(tmp); } @@ -398,4 +403,36 @@ void Potential::add_pcc_force(const UnitCell& cell, ModuleBase::matrix& force) c } } +void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + for (const PotBase* component : this->components) + { + component->add_solvation_force(cell, force); + } +} + +void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = 0.0; + non_electrostatic = 0.0; + for (const PotBase* component : this->components) + { + double el = 0.0; + double nonel = 0.0; + component->get_solvation_energy(el, nonel); + electrostatic += el; + non_electrostatic += nonel; + } +} + +const std::vector* Potential::solvent_electrostatic_potential() const +{ + for (const PotBase* component : this->components) + { + const std::vector* value = component->solvent_electrostatic_potential(); + if (value != nullptr) { return value; } + } + return nullptr; +} + } // namespace elecstate diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index e1af8269816..9bf4683a139 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -10,6 +10,8 @@ #include #include +struct Input_para; + namespace elecstate { class PotPcc; @@ -73,12 +75,17 @@ class Potential : public PotBase void init_pot(const Charge* const chg); // initialize potential components before SCF void pot_register(const std::vector& components_list); + void pot_register(const std::vector& components_list, const Input_para& input); // update potential from current charge void update_from_charge(const Charge* const chg, const UnitCell* const ucell); // interface for SCF-converged, etxc vtxc for Energy, vnew for force_scc void get_vnew(const Charge* chg, ModuleBase::matrix& vnew); - PotBase* get_pot_type(const std::string& pot_type); + PotBase* get_pot_type(const std::string& pot_type, const Input_para& input); + + void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; + void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; + const std::vector* solvent_electrostatic_potential() const override; /** * @brief Inject the shared RT-TDDFT field state before potential setup. @@ -228,6 +235,7 @@ class Potential : public PotBase // the registered PCC component, or nullptr when PCC is off const PotPcc* pcc_component() const; + PotBase* make_sccs_potential(const Input_para& input); std::vector v_eff_fixed; ModuleBase::matrix v_eff; diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 4280bb8210f..6890d083bdc 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -3,6 +3,7 @@ #include "gatefield.h" #include "pot_local.h" #include "pot_pcc.h" +#include "pot_sccs.h" #include "pot_sep.h" #include "pot_surchem.h" #include "pot_xc.h" @@ -21,7 +22,20 @@ namespace elecstate { -PotBase* Potential::get_pot_type(const std::string& pot_type) +PotBase* Potential::make_sccs_potential(const Input_para& input) +{ + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + make_sccs_config_from_input(input, config, solver); + if (this->ucell_ != nullptr) + { + const std::string warning = check_sccs_charge(config, *this->ucell_, input.nelec); + if (!warning.empty()) { ModuleBase::WARNING("Potential::make_sccs_potential", warning); } + } + return new PotSccs(this->rho_basis_, config, solver, input.nelec); +} + +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& input) { ModuleBase::TITLE("Potential", "get_pot_type"); if (pot_type == "local") @@ -34,11 +48,10 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) } else if (pot_type == "xc") { - return new PotXC(PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); + return new PotXC(PARAM.globalv.domag, PARAM.globalv.domag_z, input.gga_grad, this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); } else if (pot_type == "pcc" || pot_type == "pcc_2d") { - const Input_para& input = PARAM.inp; const double electron_count = input.nelec; if (pot_type == "pcc") { @@ -47,13 +60,17 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) const int open_axis = input.pcc_2d_axis; return new PotPcc(this->rho_basis_, PotPcc::Dimension::slab, open_axis, electron_count); } + else if (pot_type == "sccs") + { + return this->make_sccs_potential(input); + } else if (pot_type == "surchem") { return new PotSurChem(this->rho_basis_, this->structure_factors_, this->v_eff_fixed.data(), this->solvent_); } else if (pot_type == "efield") { - return new PotEfield(this->rho_basis_, this->ucell_, this->solvent_, PARAM.inp.dip_cor_flag); + return new PotEfield(this->rho_basis_, this->ucell_, this->solvent_, input.dip_cor_flag); } else if (pot_type == "gatefield") { diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index 317ef2805a2..8bcae6b32de 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -37,6 +37,10 @@ int main(int argc, char** argv) // Mock functions for testing elecstate.cpp namespace elecstate { +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/elecstate_energy_test.cpp b/source/source_estate/test/elecstate_energy_test.cpp index 0bb1fed1b15..0ccadc08dbf 100644 --- a/source/source_estate/test/elecstate_energy_test.cpp +++ b/source/source_estate/test/elecstate_energy_test.cpp @@ -12,6 +12,18 @@ int XC_Functional::func_type = 1; bool XC_Functional::ked_flag = false; namespace elecstate { +Potential::~Potential() {} +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::cal_v_eff(const Charge*, const UnitCell*, ModuleBase::matrix&) {} +void Potential::cal_fixed_v(double*) {} +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + +void Potential::get_solvation_energy(double& el, double& cav) const +{ + el = 0.8; + cav = 0.9; +} + double Potential::pcc_energy_rydberg() const { return 0.0; @@ -80,7 +92,7 @@ class ElecStateEnergyTest : public ::testing::Test /// explicitly, so the fixture owns them instead of staging them in the /// global parameter singleton. The values mirror the Input_para defaults, /// except sc_mag_switch, which the original fixture turned on. - bool imp_sol = false; + int imp_sol = 0; bool sc_mag_switch = true; int dft_plus_u = 0; std::string assume_isolated = "none"; @@ -271,3 +283,21 @@ TEST_F(ElecStateEnergyTest, CalBandgapUpDwBoundaryConditions) EXPECT_DOUBLE_EQ(elecstate->bandgap_up, 5.0); EXPECT_DOUBLE_EQ(elecstate->bandgap_dw, 5.0); } + +TEST_F(ElecStateEnergyTest, AddsSccsEnergyOnceAndResetsVacuum) +{ + elecstate::Potential potential; + elecstate->pot = &potential; + elecstate->f_en.deband = 0.1; + imp_sol = 2; + elecstate->cal_energies(2, imp_sol, sc_mag_switch, dft_plus_u, assume_isolated); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_el, 0.8); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_cav, 0.9); + EXPECT_DOUBLE_EQ(elecstate->f_en.etot, 2.4); + imp_sol = 0; + elecstate->cal_energies(2, imp_sol, sc_mag_switch, dft_plus_u, assume_isolated); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_el, 0.0); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_cav, 0.0); + EXPECT_DOUBLE_EQ(elecstate->f_en.etot, 0.7); + elecstate->pot = nullptr; +} diff --git a/source/source_estate/test/elecstate_pw_test.cpp b/source/source_estate/test/elecstate_pw_test.cpp index 1c47c8356b8..fc74f035816 100644 --- a/source/source_estate/test/elecstate_pw_test.cpp +++ b/source/source_estate/test/elecstate_pw_test.cpp @@ -26,6 +26,10 @@ class TestParameters int XC_Functional::func_type = 1; namespace elecstate { +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/potential_new_test.cpp b/source/source_estate/test/potential_new_test.cpp index 6dc01b9f321..026b2222d70 100644 --- a/source/source_estate/test/potential_new_test.cpp +++ b/source/source_estate/test/potential_new_test.cpp @@ -130,6 +130,15 @@ class MockPotComponent : public PotBase } } + void add_solvation_force(const UnitCell&, ModuleBase::matrix& force) const override + { + if (type_ != "solvation") { return; } + for (int ia = 0; ia < force.nr; ++ia) + { + for (int axis = 0; axis < force.nc; ++axis) { force(ia, axis) += 1.25; } + } + } + static void reset() { created = 0; @@ -153,7 +162,7 @@ int MockPotComponent::destroyed = 0; int MockPotComponent::fixed_calls = 0; int MockPotComponent::dynamic_calls = 0; -PotBase* Potential::get_pot_type(const std::string& pot_type) +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para&) { const int grid_size = this->get_rho_basis() == nullptr ? 0 : this->get_rho_basis()->nrxx; return new MockPotComponent(pot_type, grid_size); @@ -411,3 +420,26 @@ TEST_F(PotentialNewTest, DifferentBasisObjectsRequireMatchingGammaOnly) EXPECT_EXIT(potential->update_from_charge(&charge, ucell.get()), ::testing::ExitedWithCode(1), ""); } + +TEST_F(PotentialNewTest, SolvationForcesUseRegisteredComponentsOnce) +{ + setup_smooth_basis(); + create_potential(smooth_basis.get(), smooth_basis.get()); + potential->pot_register({"fixed", "solvation"}); + ModuleBase::matrix force(2, 3); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { force(ia, axis) = 7.0; } + } + potential->add_solvation_force(*ucell, force); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } + } + potential->pot_register({"fixed"}); + potential->add_solvation_force(*ucell, force); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } + } +} diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp new file mode 100644 index 00000000000..4bf0aa26977 --- /dev/null +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -0,0 +1,244 @@ +#include "source_estate/module_pot/pot_sccs.h" +#include "source_io/module_parameter/input_parameter.h" +#include "source_hamilt/module_sccs/test/sccs_test.h" +#include "source_base/global_variable.h" + +#include +#include +#include +#include +#include + +// Public storage fixtures: no INPUT initialization or privately owned atom maps. +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +SepPot::SepPot() {} +SepPot::~SepPot() {} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} +Charge::Charge() {} +Charge::~Charge() {} + +class PotSccsTest : public SccsTest::PwTest +{ +protected: + // Two atoms in a uniform electron density: the solvation force is the + // finite difference of the Rydberg energy and is added once to an existing + // force. + void check_solvation_force(const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver) + { + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 2; + cell.atoms = &atom; + atom.na = 2; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.21, 0.32, 0.43), + ModuleBase::Vector3(0.64, 0.51, 0.27)}; + const double density_value = 2.0 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + elecstate::PotSccs component(&basis, config, solver, 2.0); + ModuleBase::matrix potential(1, basis.nrxx); + component.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix force(2, 3); + component.add_solvation_force(cell, force); + const double step = 1e-4; + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + const double original = atom.tau[ia][axis]; + atom.tau[ia][axis] = original + step / length; + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + const double positive = component.get_energy(); + atom.tau[ia][axis] = original - step / length; + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + const double negative = component.get_energy(); + atom.tau[ia][axis] = original; + const double finite_difference = -(positive - negative) / (2.0 * step); + EXPECT_NEAR(force(ia, axis), finite_difference, 1e-8); + } + } + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix accumulated(2, 3); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { accumulated(ia, axis) = 7.0; } + } + component.add_solvation_force(cell, accumulated); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + const double expected_force = 7.0 + force(ia, axis); + EXPECT_NEAR(accumulated(ia, axis), expected_force, 1e-10); + } + } + } +}; + +TEST_F(PotSccsTest, IndependentInstancesAndSpinChannels) +{ + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 1.0 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data(), density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + Input_para input; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + elecstate::PotSccs first(&basis, config, solver, 1.0); + config.cavity.epsilon_bulk = 1.0; + elecstate::PotSccs vacuum(&basis, config, solver, 1.0); + EXPECT_DOUBLE_EQ(first.get_energy(), 0.0); + ModuleBase::matrix potential(1, basis.nrxx); + first.cal_v_eff(&charge, &cell, potential); + const double first_energy = first.get_energy(); + EXPECT_LT(first_energy, 0.0); + double el = 0.0; + double cav = 0.0; + first.get_solvation_energy(el, cav); + EXPECT_DOUBLE_EQ(el, first_energy); + EXPECT_DOUBLE_EQ(cav, 0.0); + ModuleBase::matrix zero(1, basis.nrxx); + vacuum.cal_v_eff(&charge, &cell, zero); + EXPECT_NEAR(vacuum.get_energy(), 0.0, 1e-12); + EXPECT_DOUBLE_EQ(first.get_energy(), first_energy); + for (int ir = 0; ir < basis.nrxx; ++ir) { EXPECT_NEAR(zero(0, ir), 0.0, 1e-11); } + for (double& value : density) { value *= 0.5; } + charge.nspin = 2; + ModuleBase::matrix spin_potential(2, basis.nrxx); + first.cal_v_eff(&charge, &cell, spin_potential); + EXPECT_NEAR(first.get_energy(), first_energy, 1e-12); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(spin_potential(0, ir), potential(0, ir), 1e-11); + EXPECT_DOUBLE_EQ(spin_potential(0, ir), spin_potential(1, ir)); + } +} + +TEST_F(PotSccsTest, InputMapsOntoSccsConfig) +{ + // The custom defaults follow Environ environ_type input: + // env_static_permittivity 1, rhomin 1e-4, rhomax 5e-3, no surface tension + // or pressure. + Input_para input; + ModuleSccs::SccsConfig defaults; + ModuleSccs::PolarizationSolverParameters default_solver; + elecstate::make_sccs_config_from_input(input, defaults, default_solver); + EXPECT_EQ(input.sccs_preset, "custom"); + EXPECT_DOUBLE_EQ(defaults.cavity.epsilon_bulk, 1.0); + EXPECT_DOUBLE_EQ(defaults.cavity.density_min, 1e-4); + EXPECT_DOUBLE_EQ(defaults.cavity.density_max, 5e-3); + EXPECT_DOUBLE_EQ(defaults.surface_tension, 0.0); + EXPECT_DOUBLE_EQ(defaults.pressure, 0.0); + EXPECT_EQ(default_solver.max_iterations, 200); + input.sccs_preset = "water-neutral"; + input.sccs_epsilon = 2.0; // replaced by the preset + input.sccs_maxiter = 42; + input.sccs_surface_eta = 1e-6; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); + EXPECT_DOUBLE_EQ(config.cavity.density_min, 1e-4); + EXPECT_LT(config.pressure, 0.0); + EXPECT_DOUBLE_EQ(config.surface_regularization, 1e-6); + EXPECT_EQ(solver.max_iterations, 42); +} + +TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) +{ + Input_para input; + input.sccs_epsilon = 5.0; + input.sccs_tol_rms = 1e-13; + input.sccs_tol_max = 1e-12; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + check_solvation_force(config, solver); +} + +TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) +{ + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + Input_para input; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + EXPECT_TRUE(elecstate::check_sccs_charge(config, cell, 1.0).empty()); + EXPECT_NE(elecstate::check_sccs_charge(config, cell, 0.5).find("G = 0"), std::string::npos); + ModuleSccs::SccsConfig vacuum_config = config; + vacuum_config.cavity.epsilon_bulk = 1.0; + EXPECT_TRUE(elecstate::check_sccs_charge(vacuum_config, cell, 0.5).empty()); + + const double density_value = 0.5 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + // One log per rank: the MPI variant runs this test on several ranks at once. + const std::string log_name = "pot_sccs_charge_warning_" + std::to_string(basis.poolrank) + ".log"; + std::ofstream& warning_log = GlobalV::ofs_warning; + warning_log.open(log_name.c_str()); + elecstate::PotSccs matched(&basis, config, solver, 0.5); + ModuleBase::matrix matched_potential(1, basis.nrxx); + matched.cal_v_eff(&charge, &cell, matched_potential); + warning_log.flush(); + std::ifstream matched_log(log_name.c_str()); + std::stringstream matched_text; + matched_text << matched_log.rdbuf(); + EXPECT_EQ(matched_text.str().find("SCCS grid electron count"), std::string::npos); + EXPECT_TRUE(std::isfinite(matched.get_energy())); + + elecstate::PotSccs mismatched(&basis, config, solver, 1.0); + ModuleBase::matrix mismatched_potential(1, basis.nrxx); + mismatched.cal_v_eff(&charge, &cell, mismatched_potential); + warning_log.close(); + std::ifstream mismatched_log(log_name.c_str()); + std::stringstream mismatched_text; + mismatched_text << mismatched_log.rdbuf(); + EXPECT_NE(mismatched_text.str().find("SCCS grid electron count differs"), std::string::npos); + EXPECT_DOUBLE_EQ(mismatched.get_energy(), matched.get_energy()); + std::remove(log_name.c_str()); +} diff --git a/source/source_estate/write_elecstat_pot.cpp b/source/source_estate/write_elecstat_pot.cpp index 5f0d0d68c37..f530c3d3d0c 100644 --- a/source/source_estate/write_elecstat_pot.cpp +++ b/source/source_estate/write_elecstat_pot.cpp @@ -23,6 +23,7 @@ void write_elecstat_pot( const double* v_eff, const surchem& solvent, const std::vector* pcc_potential, + const std::vector* solvent_potential, const int precision, std::ofstream& ofs_running) { @@ -31,10 +32,15 @@ void write_elecstat_pot( std::vector v_elecstat(rho_basis->nrxx, 0.0); + if (solvent_potential != nullptr && solvent_potential->size() != static_cast(rho_basis->nrxx)) + { + ModuleBase::WARNING_QUIT("write_elecstat_pot", "Solvent potential does not match the output grid"); + } + const int nspin = PARAM.inp.nspin; const int efield = PARAM.inp.efield_flag; const int dip_corr = PARAM.inp.dip_cor_flag; - const bool imp_sol = PARAM.inp.imp_sol; + const int imp_sol = PARAM.inp.imp_sol; //========================================== // Hartree potential @@ -72,7 +78,11 @@ void write_elecstat_pot( { v_elecstat[ir] += (*pcc_potential)[ir]; } - if(imp_sol == true) + if (solvent_potential != nullptr) + { + v_elecstat[ir] += (*solvent_potential)[ir]; + } + if (imp_sol == 1) { v_elecstat[ir] += solvent.delta_phi[ir]; } diff --git a/source/source_estate/write_elecstat_pot.h b/source/source_estate/write_elecstat_pot.h index 79eb72cd698..15c7dc28662 100644 --- a/source/source_estate/write_elecstat_pot.h +++ b/source/source_estate/write_elecstat_pot.h @@ -35,6 +35,7 @@ void write_elecstat_pot( const double* v_eff_fixed, const surchem& solvent, const std::vector* pcc_potential, + const std::vector* solvent_potential, const int precision, std::ofstream& ofs_running); diff --git a/source/source_hamilt/CMakeLists.txt b/source/source_hamilt/CMakeLists.txt index 90fe2257c77..d0cdba7a1ef 100644 --- a/source/source_hamilt/CMakeLists.txt +++ b/source/source_hamilt/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(module_vdw) add_subdirectory(module_surchem) +add_subdirectory(module_sccs) add_subdirectory(module_xc) add_subdirectory(module_hcontainer) diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt new file mode 100644 index 00000000000..f76b46ec716 --- /dev/null +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -0,0 +1,16 @@ +add_library(sccs OBJECT + sccs_cavity.cpp + sccs_parameters.cpp + sccs_functional.cpp + sccs_ionic_charge.cpp + sccs_ionic_force.cpp + sccs_pw_coulomb.cpp + sccs_response.cpp) + +if(ENABLE_COVERAGE) + add_coverage(sccs) +endif() + +if(BUILD_TESTING) + add_subdirectory(test) +endif() diff --git a/source/source_hamilt/module_sccs/sccs_cavity.cpp b/source/source_hamilt/module_sccs/sccs_cavity.cpp new file mode 100644 index 00000000000..94f69ed0abc --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cavity.cpp @@ -0,0 +1,39 @@ +#include "sccs_cavity.h" +#include "source_base/constants.h" + +#include + +namespace ModuleSccs +{ +CavityPoint evaluate_cavity(double density, const CavityParameters& parameters) +{ + CavityPoint point; + if (density <= parameters.density_min) + { + point.epsilon = parameters.epsilon_bulk; + } + else if (density >= parameters.density_max) + { + point.solute = 1.0; + point.epsilon = 1.0; + } + else + { + const double density_ratio = parameters.density_max / parameters.density_min; + const double log_width = std::log(density_ratio); + const double relative_density = parameters.density_max / density; + const double log_relative_density = std::log(relative_density); + const double x = log_relative_density / log_width; + const double angle = ModuleBase::TWO_PI * x; + const double solvent = x - std::sin(angle) / ModuleBase::TWO_PI; + const double dsolvent_drho = -(1.0 - std::cos(angle)) / (log_width * density); + const double log_epsilon = std::log(parameters.epsilon_bulk); + const double log_epsilon_point = log_epsilon * solvent; + point.solute = 1.0 - solvent; + point.dsolute_drho = -dsolvent_drho; + point.epsilon = std::exp(log_epsilon_point); + point.depsilon_drho = point.epsilon * log_epsilon * dsolvent_drho; + } + return point; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_cavity.h b/source/source_hamilt/module_sccs/sccs_cavity.h new file mode 100644 index 00000000000..af3e4eeb9d9 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cavity.h @@ -0,0 +1,27 @@ +#ifndef SCCS_CAVITY_H +#define SCCS_CAVITY_H + +namespace ModuleSccs +{ +// Densities are in electrons/Bohr^3. The cavity uses the electronic density. +struct CavityParameters +{ + double density_min = 0.0; + double density_max = 0.0; + double epsilon_bulk = 0.0; +}; + +struct CavityPoint +{ + double solute = 0.0; + double dsolute_drho = 0.0; + double epsilon = 0.0; + double depsilon_drho = 0.0; +}; + +// Parameters are validated by the INPUT reader. Negative Fourier ringing +// uses the bulk limit. +CavityPoint evaluate_cavity(double density, const CavityParameters& parameters); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp new file mode 100644 index 00000000000..7547d7dae7b --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -0,0 +1,68 @@ +#include "sccs_functional.h" +#include "sccs_parameters.h" +#include "sccs_response.h" +#include "sccs_pw_coulomb.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include + +namespace ModuleSccs +{ +void evaluate_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result) +{ + const std::size_t size = charge.size(); + const double dv = basis.omega / basis.nxyz; + PeriodicCoulombOperator coulomb(basis, tpiba); + std::vector vacuum; + coulomb.apply_potential(charge, vacuum); + FunctionalResult candidate; + candidate.reaction_potential.resize(size); + candidate.electron_potential.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const double reaction = response.polarization.potential[i] - vacuum[i]; + candidate.reaction_potential[i] = reaction; + candidate.reaction_energy += 0.5 * charge[i] * reaction * dv; + candidate.electron_potential[i] = -reaction + response.cavity_potential[i]; + candidate.volume += response.solute[i] * dv; + } + Parallel_Reduce::reduce_pool(candidate.reaction_energy); + Parallel_Reduce::reduce_pool(candidate.volume); + candidate.volume_energy = config.pressure * candidate.volume; + + // Spectral surface derivative: dE/dn = (p - gamma div(grad s/|grad s|)) s'. + std::vector> solute_g(basis.npw); + basis.real2recip(response.solute.data(), solute_g.data()); + std::vector> unit_gradient(size); + XC_Functional::grad_rho(solute_g.data(), unit_gradient.data(), &basis, tpiba); + const double eta = config.surface_regularization; + for (std::size_t i = 0; i < size; ++i) + { + const double norm_squared = unit_gradient[i] * unit_gradient[i] + eta * eta; + const double norm = std::sqrt(norm_squared); + candidate.surface += (norm - eta) * dv; + unit_gradient[i] /= norm; + } + Parallel_Reduce::reduce_pool(candidate.surface); + candidate.surface_energy = config.surface_tension * candidate.surface; + std::vector divergence(size); + XC_Functional::grad_dot(unit_gradient.data(), divergence.data(), &basis, tpiba); + for (std::size_t i = 0; i < size; ++i) + { + const double nonel = (config.pressure - config.surface_tension * divergence[i]) + * response.dsolute_drho[i]; + candidate.electron_potential[i] += nonel; + } + result = std::move(candidate); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h new file mode 100644 index 00000000000..b0a021b9d47 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -0,0 +1,36 @@ +#ifndef SCCS_FUNCTIONAL_H +#define SCCS_FUNCTIONAL_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct SccsConfig; +struct SccsResponse; + +struct FunctionalResult +{ + double reaction_energy = 0.0; // Hartree + double surface_energy = 0.0; + double volume_energy = 0.0; + double surface = 0.0; // Bohr^2 + double volume = 0.0; // Bohr^3 + std::vector reaction_potential; // dielectric minus vacuum, Hartree + std::vector electron_potential; // derivative of all three energy terms +}; + +// Collective over the PW pool. +void evaluate_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp new file mode 100644 index 00000000000..a1015081797 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp @@ -0,0 +1,33 @@ +#include "sccs_ionic_charge.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/cell_tools.h" + +#include +#include + +namespace ModuleSccs +{ +void gaussian_ionic_density(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density) +{ + std::vector> charge(basis.npw, 0.0); + const double tpiba2 = tpiba * tpiba; + for (const unitcell::AtomData& atom : atoms) + { + for (int ig = 0; ig < basis.npw; ++ig) + { + const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; + const double phase = tpiba * (basis.gcar[ig] * atom.position); + const std::complex phase_factor = ModuleBase::NEG_IMAG_UNIT * phase; + charge[ig] += atom.valence_charge / basis.omega * std::exp(exponent) * std::exp(phase_factor); + } + } + density.resize(basis.nrxx); + basis.recip2real(charge.data(), density.data()); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.h b/source/source_hamilt/module_sccs/sccs_ionic_charge.h new file mode 100644 index 00000000000..bf29d2276fd --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_charge.h @@ -0,0 +1,29 @@ +#ifndef SCCS_IONIC_CHARGE_H +#define SCCS_IONIC_CHARGE_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} +namespace unitcell +{ +struct AtomData; +} + +namespace ModuleSccs +{ +// rho ~ exp(-r^2/spread^2), in Bohr; the Environ default ionic spread. +const double gaussian_ion_spread = 0.5; + +// AtomData is extracted by source_cell; positions are Cartesian Bohr. +// Collective over the PW pool. +void gaussian_ionic_density(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/sccs_ionic_force.cpp new file mode 100644 index 00000000000..0d6ee4f6726 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_force.cpp @@ -0,0 +1,45 @@ +#include "sccs_ionic_force.h" + +#include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/cell_tools.h" + +#include +#include + +namespace ModuleSccs +{ +void gaussian_ionic_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces) +{ + std::vector> potential_g(basis.npw); + basis.real2recip(reaction_potential.data(), potential_g.data()); + forces.assign(atoms.size(), ModuleBase::Vector3(0.0, 0.0, 0.0)); + const double tpiba2 = tpiba * tpiba; + // F = -integral(v_reaction d rho_ion/d R), using ABACUS's normalized FFT. + // The cell volume cancels the Z/omega in the Gaussian charge coefficients. + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + const unitcell::AtomData& atom = atoms[ia]; + for (int ig = 0; ig < basis.npw; ++ig) + { + const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; + const double gaussian = std::exp(exponent); + const double phase = tpiba * (basis.gcar[ig] * atom.position); + const std::complex phase_argument = ModuleBase::NEG_IMAG_UNIT * phase; + const std::complex phase_factor = std::exp(phase_argument); + const std::complex weighted = std::conj(potential_g[ig]) * phase_factor; + const double factor = -atom.valence_charge * tpiba * gaussian * weighted.imag(); + forces[ia] += basis.gcar[ig] * factor; + } + Parallel_Reduce::reduce_pool(forces[ia].x); + Parallel_Reduce::reduce_pool(forces[ia].y); + Parallel_Reduce::reduce_pool(forces[ia].z); + } +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.h b/source/source_hamilt/module_sccs/sccs_ionic_force.h new file mode 100644 index 00000000000..150781b03d8 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_force.h @@ -0,0 +1,30 @@ +#ifndef SCCS_IONIC_FORCE_H +#define SCCS_IONIC_FORCE_H + +#include + +namespace ModuleBase +{ +template class Vector3; +} +namespace ModulePW +{ +class PW_Basis; +} +namespace unitcell +{ +struct AtomData; +} +namespace ModuleSccs +{ +// Explicit Gaussian-ion derivative of the reaction energy at fixed electron +// density. reaction_potential is in Ha; force is in Ha/Bohr, atom-major. +// Collective over the PW pool. +void gaussian_ionic_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces); +} // namespace ModuleSccs +#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.cpp b/source/source_hamilt/module_sccs/sccs_parameters.cpp new file mode 100644 index 00000000000..740ce525d1b --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_parameters.cpp @@ -0,0 +1,79 @@ +#include "sccs_parameters.h" + +#include "source_base/tool_quit.h" + + +namespace ModuleSccs +{ +namespace +{ +// CODATA 2018 SI values. Environ 3.1.1 converts with 4.35974394e-18 J and +// 0.52917720859e-10 m, so the surface tension and pressure differ from its +// values by about 2e-7 (relative). ModuleBase's SI values differ from both. +const double hartree_joule = 4.3597447222071e-18; +const double bohr_metre = 5.29177210903e-11; +} + +double dyn_per_cm_to_hartree_per_bohr2(double value) +{ + const double joule_per_square_metre = value * 1e-3; + return joule_per_square_metre * bohr_metre * bohr_metre / hartree_joule; +} + +double gpa_to_hartree_per_bohr3(double value) +{ + const double pascal = value * 1e9; + return pascal * bohr_metre * bohr_metre * bohr_metre / hartree_joule; +} + +Preset parse_preset(const std::string& name) +{ + if (name == "custom") { return Preset::Custom; } + if (name == "vacuum") { return Preset::Vacuum; } + if (name == "water-neutral") { return Preset::WaterNeutral; } + if (name == "water-cation") { return Preset::WaterCation; } + if (name == "water-anion") { return Preset::WaterAnion; } + const std::string message = "Unknown SCCS preset: " + name; + ModuleBase::WARNING_QUIT("ModuleSccs::parse_preset", message); +} + +SccsConfig make_sccs_config(Preset preset) +{ + if (preset == Preset::Custom) + { + ModuleBase::WARNING_QUIT("ModuleSccs::make_sccs_config", + "SCCS custom configurations must be supplied explicitly"); + } + SccsConfig config; + config.cavity.epsilon_bulk = 78.3; + config.surface_regularization = 1e-8; + if (preset == Preset::Vacuum || preset == Preset::WaterNeutral) + { + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + if (preset == Preset::Vacuum) + { + config.cavity.epsilon_bulk = 1.0; + } + else + { + config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(47.9); + config.pressure = gpa_to_hartree_per_bohr3(-0.36); + } + } + else if (preset == Preset::WaterCation) + { + config.cavity.density_min = 2e-4; + config.cavity.density_max = 3.5e-3; + config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(5.0); + config.pressure = gpa_to_hartree_per_bohr3(0.125); + } + else + { + config.cavity.density_min = 2.4e-3; + config.cavity.density_max = 1.55e-2; + config.pressure = gpa_to_hartree_per_bohr3(0.45); + } + return config; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h new file mode 100644 index 00000000000..707c5c0c386 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -0,0 +1,36 @@ +#ifndef SCCS_PARAMETERS_H +#define SCCS_PARAMETERS_H + +#include "sccs_cavity.h" + +#include + +namespace ModuleSccs +{ +enum class Preset { Custom, Vacuum, WaterNeutral, WaterCation, WaterAnion }; + +struct SccsConfig +{ + CavityParameters cavity; + double surface_tension = 0.0; // Hartree/Bohr^2 + double pressure = 0.0; // Hartree/Bohr^3 + double surface_regularization = 0.0; +}; + +struct PolarizationSolverParameters +{ + int max_iterations = 200; + double tolerance_rms = 1e-10; + double tolerance_max = 1e-8; +}; + +// Exact, lower-case names as accepted by the INPUT reader; an unknown name +// stops the run. +Preset parse_preset(const std::string& name); +// Published presets only; custom configurations are supplied by callers. +SccsConfig make_sccs_config(Preset preset); +double dyn_per_cm_to_hartree_per_bohr2(double value); +double gpa_to_hartree_per_bohr3(double value); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp new file mode 100644 index 00000000000..62fbebabd20 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp @@ -0,0 +1,32 @@ +#include "sccs_pw_coulomb.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" + +namespace ModuleSccs +{ +PeriodicCoulombOperator::PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba) + : basis_(basis), tpiba_(tpiba) +{ +} + +void PeriodicCoulombOperator::apply_potential(const std::vector& charge, std::vector& potential) +{ + reciprocal_work_.resize(basis_.npw); + basis_.real2recip(charge.data(), reciprocal_work_.data()); + const double tpiba2 = tpiba_ * tpiba_; + for (int ig = 0; ig < basis_.npw; ++ig) + { + if (ig == basis_.ig_gge0 || basis_.gg[ig] == 0.0) + { + reciprocal_work_[ig] = 0.0; + } + else + { + reciprocal_work_[ig] = ModuleBase::FOUR_PI * reciprocal_work_[ig] / (tpiba2 * basis_.gg[ig]); + } + } + potential.resize(basis_.nrxx); + basis_.recip2real(reciprocal_work_.data(), potential.data()); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h new file mode 100644 index 00000000000..41245a45608 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h @@ -0,0 +1,31 @@ +#ifndef SCCS_PW_COULOMB_H +#define SCCS_PW_COULOMB_H + +#include +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +// Potential in Hartree for a signed charge density in e/Bohr^3. Calls are +// collective over the ABACUS pool used to initialize basis. G=0 is removed: +// this is the periodic Coulomb inverse with a uniform background. +class PeriodicCoulombOperator +{ +public: + PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba); + void apply_potential(const std::vector& charge, std::vector& potential); + +private: + const ModulePW::PW_Basis& basis_; + const double tpiba_; + // Scratch only, overwritten on each call; PW_Basis FFTs are also non-reentrant. + std::vector> reciprocal_work_; +}; +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp new file mode 100644 index 00000000000..58d247bce21 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -0,0 +1,316 @@ +#include "sccs_response.h" +#include "sccs_parameters.h" +#include "sccs_pw_coulomb.h" + +#include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_quit.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include +#include +#include + +namespace ModuleSccs +{ +namespace +{ +void prepare_cavity(const std::vector& density, + const CavityParameters& cavity, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& response, + std::vector& coefficient) +{ + const std::size_t size = density.size(); + response.solute.resize(size); + response.dsolute_drho.resize(size); + response.epsilon.resize(size); + response.depsilon_drho.resize(size); + response.grad_log_epsilon.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const CavityPoint point = evaluate_cavity(density[i], cavity); + response.solute[i] = point.solute; + response.dsolute_drho[i] = point.dsolute_drho; + response.epsilon[i] = point.epsilon; + response.depsilon_drho[i] = point.depsilon_drho; + } + + std::vector> density_g(basis.npw); + basis.real2recip(density.data(), density_g.data()); + std::vector> gradient(size); + std::vector laplacian(size); + // These XC helpers use G in units of tpiba and ABACUS's normalized FFT. + XC_Functional::grad_rho(density_g.data(), gradient.data(), &basis, tpiba); + XC_Functional::laplacian_rho(density_g.data(), laplacian.data(), &basis, tpiba); + const double density_ratio = cavity.density_max / cavity.density_min; + const double width = std::log(density_ratio); + const double log_bulk = std::log(cavity.epsilon_bulk); + coefficient.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + double second_log = 0.0; + if (density[i] > cavity.density_min && density[i] < cavity.density_max) + { + const double local_ratio = cavity.density_max / density[i]; + const double x = std::log(local_ratio) / width; + const double angle = ModuleBase::TWO_PI * x; + second_log = log_bulk * (1.0 - std::cos(angle) + ModuleBase::TWO_PI * std::sin(angle) / width) + / (width * density[i] * density[i]); + } + const double first_log = response.depsilon_drho[i] / response.epsilon[i]; + const ModuleBase::Vector3& local_gradient = gradient[i]; + const double gradient_square = local_gradient * local_gradient; + response.grad_log_epsilon[i] = local_gradient * first_log; + const double lap_log = first_log * laplacian[i] + second_log * gradient_square; + coefficient[i] = response.epsilon[i] * (0.5 * lap_log + 0.25 * first_log * first_log * gradient_square) + / ModuleBase::FOUR_PI; + } +} + +// Pool-reduced norms, so every rank takes the same convergence decision. A +// non-finite residual fails convergence and stops at the next CG step. +void residual_norms(const std::vector& values, + const ModulePW::PW_Basis& basis, + double& rms, + double& maximum) +{ + double square = 0.0; + maximum = 0.0; + for (double value : values) + { + square += value * value; + const double magnitude = std::abs(value); + maximum = std::max(maximum, magnitude); + } + Parallel_Reduce::reduce_pool(square); + Parallel_Reduce::reduce_max_pool(basis.poolnproc, maximum); + const double mean_square = square / basis.nxyz; + rms = std::sqrt(mean_square); +} + +bool converged(const PolarizationResult& result, const PolarizationSolverParameters& solver) +{ + return result.residual_rms <= solver.tolerance_rms && result.residual_max <= solver.tolerance_max; +} + +double grid_dot(const std::vector& left, + const std::vector& right, + const ModulePW::PW_Basis& basis) +{ + double value = 0.0; + for (std::size_t i = 0; i < left.size(); ++i) + { + value += left[i] * right[i]; + } + Parallel_Reduce::reduce_pool(value); + return value * basis.omega / basis.nxyz; +} + +// P r = eps^-1/2 G eps^-1/2 r; only FFT scratch survives an application. +void apply_preconditioner(const std::vector& rhs, + const std::vector& invsqrt, + PeriodicCoulombOperator& coulomb, + std::vector& weighted, + std::vector& value) +{ + const std::size_t size = rhs.size(); + weighted.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + weighted[i] = rhs[i] * invsqrt[i]; + } + coulomb.apply_potential(weighted, value); + for (std::size_t i = 0; i < size; ++i) + { + value[i] *= invsqrt[i]; + } +} + +// The preconditioned sqrt-CG system of one response. +struct SqrtCgOperator +{ + const std::vector& coefficient; + const std::vector& invsqrt; + const ModulePW::PW_Basis& basis; + PeriodicCoulombOperator& coulomb; +}; + +// Potential and charge residual of the sqrt-CG, with preconditioner scratch. +struct SqrtCgState +{ + std::vector potential; + std::vector residual; + std::vector weighted; + std::vector preconditioned; +}; + +// Replace the cold start by one fixed-point update of the previous solution, +// phi = P (charge - c phi_old), when that lowers the RMS charge residual. +void try_warm_start(const SqrtCgOperator& system, + const std::vector& charge, + const std::vector& initial_potential, + SqrtCgState& state, + PolarizationResult& polarization) +{ + const std::size_t size = charge.size(); + const std::vector& coefficient = system.coefficient; + std::vector& z = state.preconditioned; + std::vector guess_residual(size); + for (std::size_t i = 0; i < size; ++i) + { + guess_residual[i] = charge[i] - coefficient[i] * initial_potential[i]; + } + apply_preconditioner(guess_residual, system.invsqrt, system.coulomb, state.weighted, z); + for (std::size_t i = 0; i < size; ++i) + { + guess_residual[i] = coefficient[i] * (initial_potential[i] - z[i]); + } + double guess_rms = 0.0; + double guess_max = 0.0; + residual_norms(guess_residual, system.basis, guess_rms, guess_max); + if (guess_rms < polarization.residual_rms) + { + state.potential.swap(z); + state.residual.swap(guess_residual); + polarization.warm_started = true; + polarization.residual_rms = guess_rms; + polarization.residual_max = guess_max; + } +} + +// Preconditioned CG until both residual tolerances hold; a breakdown or a +// missed tolerance stops the run with WARNING_QUIT. +void iterate_sqrt_cg(const SqrtCgOperator& system, + const PolarizationSolverParameters& solver, + SqrtCgState& state, + PolarizationResult& polarization) +{ + const ModulePW::PW_Basis& basis = system.basis; + const std::vector& coefficient = system.coefficient; + std::vector& potential = state.potential; + std::vector& residual = state.residual; + std::vector& z = state.preconditioned; + const std::size_t size = residual.size(); + std::vector direction(size, 0.0); + std::vector image(size, 0.0); + double old_rz = 0.0; + for (int iteration = 1; !converged(polarization, solver) && iteration <= solver.max_iterations; ++iteration) + { + apply_preconditioner(residual, system.invsqrt, system.coulomb, state.weighted, z); + const double rz = grid_dot(residual, z, basis); + if (!std::isfinite(rz) || std::abs(rz) < 1e-30) + { + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", + "SCCS sqrt-CG has a null or nonfinite preconditioned residual"); + } + const double beta = std::abs(old_rz) > 1e-30 ? rz / old_rz : 0.0; + old_rz = rz; + for (std::size_t i = 0; i < size; ++i) + { + direction[i] = z[i] + beta * direction[i]; + image[i] = coefficient[i] * z[i] + residual[i] + beta * image[i]; + } + const double curvature = grid_dot(direction, image, basis); + if (!std::isfinite(curvature) || curvature == 0.0) + { + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", "SCCS sqrt-CG has invalid curvature"); + } + const double alpha = rz / curvature; + for (std::size_t i = 0; i < size; ++i) + { + potential[i] += alpha * direction[i]; + residual[i] -= alpha * image[i]; + } + polarization.iterations = iteration; + residual_norms(residual, basis, polarization.residual_rms, polarization.residual_max); + } + if (!converged(polarization, solver)) + { + std::ostringstream message; + message << "SCCS sqrt-CG did not reach both residual tolerances within sccs_maxiter = " + << solver.max_iterations << " iterations (RMS " << polarization.residual_rms << ", MAX " + << polarization.residual_max << ")"; + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", message.str()); + } +} + +// Keep the unshifted solution for warm starts, then store the zero-mean +// potential, its gradient and the cavity potential. +void finish_response(const ModulePW::PW_Basis& basis, + double tpiba, + std::vector& potential, + SccsResponse& response) +{ + PolarizationResult& polarization = response.polarization; + const std::size_t size = potential.size(); + response.restart_potential = potential; + double mean = 0.0; + for (double value : potential) + { + mean += value; + } + Parallel_Reduce::reduce_pool(mean); + mean /= basis.nxyz; + for (double& value : potential) + { + 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); + polarization.potential.swap(potential); + 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 + +void solve_sccs_response(const std::vector& density, + const std::vector& charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result) +{ + // Determine start mode collectively, including ranks with no real-space points. + double initial_count = initial_potential.size(); + Parallel_Reduce::reduce_pool(initial_count); + const bool warm_start = initial_count != 0.0; + + SccsResponse candidate; + std::vector coefficient; + prepare_cavity(density, cavity, basis, tpiba, candidate, coefficient); + const std::size_t size = density.size(); + std::vector invsqrt(size); + for (std::size_t i = 0; i < size; ++i) + { + invsqrt[i] = 1.0 / std::sqrt(candidate.epsilon[i]); + } + PeriodicCoulombOperator coulomb(basis, tpiba); + const SqrtCgOperator system = {coefficient, invsqrt, basis, coulomb}; + SqrtCgState state; + state.residual = charge; + state.potential.assign(size, 0.0); + PolarizationResult& polarization = candidate.polarization; + residual_norms(state.residual, basis, polarization.residual_rms, polarization.residual_max); + if (!converged(polarization, solver) && warm_start) + { + try_warm_start(system, charge, initial_potential, state, polarization); + } + iterate_sqrt_cg(system, solver, state, polarization); + finish_response(basis, tpiba, state.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 new file mode 100644 index 00000000000..0e9b5a89a98 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -0,0 +1,56 @@ +#ifndef SCCS_RESPONSE_H +#define SCCS_RESPONSE_H + +#include "sccs_cavity.h" +#include "source_base/vector3.h" + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct PolarizationSolverParameters; + +struct PolarizationResult +{ + int iterations = 0; + double residual_rms = 0.0; + double residual_max = 0.0; + bool warm_started = false; + std::vector potential; // Hartree, shifted to periodic zero mean + std::vector> gradient; // Hartree/Bohr +}; + +struct SccsResponse +{ + std::vector solute; + std::vector dsolute_drho; + std::vector epsilon; + std::vector depsilon_drho; + std::vector> grad_log_epsilon; + PolarizationResult polarization; + std::vector cavity_potential; // -epsilon' |grad v|^2/(8 pi), Hartree + std::vector restart_potential; // unshifted solution for explicit warm starts +}; + +// Periodic sqrt-CG response with chain-rule cavity derivatives, as Environ's +// CG solver with the sqrt preconditioner; charge is ions minus electrons. +// All ranks in the PW pool must call together with the same configuration. +// initial_potential is empty on all ranks for a cold start, otherwise a local grid. +// A solver breakdown or missed tolerance stops the run with WARNING_QUIT. +// No state or INPUT globals are read or retained. +void solve_sccs_response(const std::vector& cavity_density, + const std::vector& solute_charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt new file mode 100644 index 00000000000..722c01ebd01 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -0,0 +1,27 @@ +AddTest(TARGET MODULE_HAMILT_sccs_cavity + SOURCES test_sccs_cavity.cpp ../sccs_cavity.cpp) +AddTest(TARGET MODULE_HAMILT_sccs_parameters + LIBS base device + SOURCES test_sccs_parameters.cpp ../sccs_parameters.cpp ../sccs_cavity.cpp) + +AddTest(TARGET MODULE_HAMILT_sccs_poisson + LIBS parameter planewave device base + SOURCES test_sccs_pw_coulomb.cpp test_sccs_ionic_charge.cpp test_sccs_ionic_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + ../../../source_estate/test/test_pot_sccs.cpp + ../../../source_estate/module_pot/pot_sccs.cpp + ../../../source_cell/cell_tools.cpp ../../../source_cell/atom_spec.cpp + ../../../source_cell/atom_pseudo.cpp ../../../source_cell/pseudo.cpp + ../sccs_functional.cpp + ../sccs_parameters.cpp + ../sccs_ionic_charge.cpp ../sccs_ionic_force.cpp ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp + ../../module_xc/xc_grad_utils.cpp) +set_tests_properties(MODULE_HAMILT_sccs_poisson PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60) + +if(ENABLE_MPI) + add_test(NAME MODULE_HAMILT_sccs_poisson_mpi + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS}) + set_tests_properties(MODULE_HAMILT_sccs_poisson_mpi PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60 PROCESSORS 2) +endif() diff --git a/source/source_hamilt/module_sccs/test/sccs_test.h b/source/source_hamilt/module_sccs/test/sccs_test.h new file mode 100644 index 00000000000..6a1570cca4f --- /dev/null +++ b/source/source_hamilt/module_sccs/test/sccs_test.h @@ -0,0 +1,60 @@ +#ifndef SCCS_TEST_H +#define SCCS_TEST_H + +#include "source_base/constants.h" +#include "source_base/matrix3.h" +#include "source_base/parallel_comm.h" +#include "source_basis/module_pw/pw_basis.h" +#include + +#include + +namespace SccsTest +{ +extern int pool_size; +extern int pool_rank; + +class PwTest : public testing::Test +{ +protected: + PwTest() : basis("cpu", "double") {} + + void SetUp() override + { +#ifdef __MPI + basis.initmpi(pool_size, pool_rank, POOL_WORLD); +#endif + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + basis.initgrids(length, lattice, 20.0); + basis.initparameters(false, 20.0, 1, false); + basis.setuptransform(); + basis.collect_local_pw(); + tpiba = ModuleBase::TWO_PI / length; + } + + // A cosine along x or y, in the PW local real-grid ordering. + std::vector cosine_mode(int direction) const + { + std::vector values(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const int iy = ir / basis.nplane - ix * basis.ny; + const int coordinate = direction == 0 ? ix : iy; + const int count = direction == 0 ? basis.nx : basis.ny; + const double angle = ModuleBase::TWO_PI * coordinate / count; + values[ir] = std::cos(angle); + } + return values; + } + + const double length = 10.0; + double tpiba = 0.0; + ModulePW::PW_Basis basis; + std::string error; +}; +} // namespace SccsTest + +#endif diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp new file mode 100644 index 00000000000..c4d8fc15c4d --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp @@ -0,0 +1,110 @@ +#include "../sccs_cavity.h" +#include + +#include + +namespace +{ +class SccsCavityTest : public ::testing::Test +{ +protected: + ModuleSccs::CavityParameters parameters; + + void SetUp() override + { + parameters.density_min = 1e-4; + parameters.density_max = 5e-3; + parameters.epsilon_bulk = 78.3; + } +}; + +// Bulk and solute limits, continuity at both thresholds, and no dielectric +// response in vacuum. +TEST_F(SccsCavityTest, LimitsThresholdsAndVacuum) +{ + const double bulk_densities[] = {-1e-6, 0.0, 1e-4}; + for (double density : bulk_densities) + { + ModuleSccs::CavityPoint point; + point = ModuleSccs::evaluate_cavity(density, parameters); + EXPECT_DOUBLE_EQ(point.solute, 0.0); + EXPECT_DOUBLE_EQ(point.epsilon, 78.3); + EXPECT_DOUBLE_EQ(point.dsolute_drho, 0.0); + EXPECT_DOUBLE_EQ(point.depsilon_drho, 0.0); + } + const double solute_densities[] = {5e-3, 0.1}; + for (double density : solute_densities) + { + ModuleSccs::CavityPoint point; + point = ModuleSccs::evaluate_cavity(density, parameters); + EXPECT_DOUBLE_EQ(point.solute, 1.0); + EXPECT_DOUBLE_EQ(point.epsilon, 1.0); + EXPECT_DOUBLE_EQ(point.dsolute_drho, 0.0); + EXPECT_DOUBLE_EQ(point.depsilon_drho, 0.0); + } + const double near_min = parameters.density_min * (1.0 + 1e-8); + const double near_max = parameters.density_max * (1.0 - 1e-8); + ModuleSccs::CavityPoint lower; + ModuleSccs::CavityPoint upper; + lower = ModuleSccs::evaluate_cavity(near_min, parameters); + upper = ModuleSccs::evaluate_cavity(near_max, parameters); + EXPECT_NEAR(lower.solute, 0.0, 1e-14); + EXPECT_NEAR(upper.solute, 1.0, 1e-14); + EXPECT_NEAR(lower.epsilon, 78.3, 1e-12); + EXPECT_NEAR(upper.epsilon, 1.0, 1e-12); + EXPECT_NEAR(lower.dsolute_drho, 0.0, 1e-10); + EXPECT_NEAR(upper.dsolute_drho, 0.0, 1e-10); + EXPECT_NEAR(lower.depsilon_drho, 0.0, 1e-8); + EXPECT_NEAR(upper.depsilon_drho, 0.0, 1e-8); + parameters.epsilon_bulk = 1.0; + const ModuleSccs::CavityPoint vacuum = ModuleSccs::evaluate_cavity(1e-3, parameters); + EXPECT_GT(vacuum.solute, 0.0); + EXPECT_LT(vacuum.solute, 1.0); + EXPECT_DOUBLE_EQ(vacuum.epsilon, 1.0); + EXPECT_DOUBLE_EQ(vacuum.depsilon_drho, 0.0); +} + +// Analytic values at the logarithmic midpoint and finite differences across +// the transition. +TEST_F(SccsCavityTest, MidpointAndDerivativesAcrossTransition) +{ + const double product = parameters.density_min * parameters.density_max; + const double midpoint_density = std::sqrt(product); + const double ratio = parameters.density_max / parameters.density_min; + const double width = std::log(ratio); + const double epsilon = std::sqrt(parameters.epsilon_bulk); + const double log_epsilon = std::log(parameters.epsilon_bulk); + const double dsolute = 2.0 / (width * midpoint_density); + const double depsilon = -epsilon * log_epsilon * dsolute; + const ModuleSccs::CavityPoint midpoint = ModuleSccs::evaluate_cavity(midpoint_density, parameters); + EXPECT_NEAR(midpoint.solute, 0.5, 1e-14); + EXPECT_NEAR(midpoint.epsilon, epsilon, 1e-13); + EXPECT_NEAR(midpoint.dsolute_drho, dsolute, 1e-10); + EXPECT_NEAR(midpoint.depsilon_drho, depsilon, 1e-9); + + const double fractions[] = {0.2, 0.5, 0.8}; + for (double fraction : fractions) + { + const double scale = std::pow(ratio, fraction); + const double density = parameters.density_min * scale; + const double step = density * 1e-5; + const double lower_density = density - step; + const double upper_density = density + step; + ModuleSccs::CavityPoint point; + ModuleSccs::CavityPoint lower; + ModuleSccs::CavityPoint upper; + point = ModuleSccs::evaluate_cavity(density, parameters); + lower = ModuleSccs::evaluate_cavity(lower_density, parameters); + upper = ModuleSccs::evaluate_cavity(upper_density, parameters); + const double solute_fd = (upper.solute - lower.solute) / (2.0 * step); + const double epsilon_fd = (upper.epsilon - lower.epsilon) / (2.0 * step); + const double solute_ratio = solute_fd / point.dsolute_drho; + const double epsilon_ratio = epsilon_fd / point.depsilon_drho; + EXPECT_NEAR(solute_ratio, 1.0, 1e-8); + EXPECT_NEAR(epsilon_ratio, 1.0, 1e-8); + EXPECT_GT(point.dsolute_drho, 0.0); + EXPECT_LT(point.depsilon_drho, 0.0); + } +} + +} // namespace diff --git a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp new file mode 100644 index 00000000000..c25cde778d4 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -0,0 +1,84 @@ +#include "sccs_test.h" +#include "../sccs_functional.h" +#include "../sccs_parameters.h" +#include "../sccs_response.h" +#include "source_base/parallel_reduce.h" + +using SccsFunctionalTest = SccsTest::PwTest; + +TEST_F(SccsFunctionalTest, UniformDielectricEnergyAndElectronDerivative) +{ + const std::vector mode = cosine_mode(0); + std::vector charge = mode; + for (double& value : charge) + { + value *= 1e-3; + } + 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; + ModuleSccs::PolarizationSolverParameters solver; + const std::vector density(basis.nrxx, 0.0); + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold, + basis, tpiba, response); + ModuleSccs::FunctionalResult result; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, result); + const double kernel = ModuleBase::FOUR_PI / (tpiba * tpiba); + const double expected = 0.25 * basis.omega * 1e-6 * kernel * (0.2 - 1.0); + EXPECT_NEAR(result.reaction_energy, expected, 1e-13); + EXPECT_DOUBLE_EQ(result.surface_energy, 0.0); + EXPECT_DOUBLE_EQ(result.volume_energy, 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double potential = 0.8 * kernel * charge[ir]; + EXPECT_NEAR(result.electron_potential[ir], potential, 1e-13); + } +} + +TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) +{ + ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + config.surface_tension = 0.01; + config.pressure = -0.001; + config.surface_regularization = 0.02; + const std::vector mode = cosine_mode(0); + ModuleSccs::SccsResponse response; + response.solute.resize(basis.nrxx); + response.dsolute_drho.assign(basis.nrxx, 1.0); + response.polarization.potential.assign(basis.nrxx, 0.0); + response.cavity_potential.assign(basis.nrxx, 0.0); + const std::vector charge(basis.nrxx, 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + response.solute[ir] = 0.5 + 0.2 * mode[ir]; + } + ModuleSccs::FunctionalResult center; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, center); + // Avoid cancellation of grid-integrated energies at excessively small steps. + const double step = 1e-5; + double predicted = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + predicted += center.electron_potential[ir] * mode[ir]; + response.solute[ir] += step * mode[ir]; + } + Parallel_Reduce::reduce_pool(predicted); + predicted *= basis.omega / basis.nxyz; + ModuleSccs::FunctionalResult plus; + ModuleSccs::FunctionalResult minus; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, plus); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + response.solute[ir] -= 2.0 * step * mode[ir]; + } + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, minus); + const double fd = (plus.surface_energy + plus.volume_energy - minus.surface_energy - minus.volume_energy) + / (2.0 * step); + EXPECT_NEAR(fd, predicted, 1e-8); + const double expected_volume = 0.5 * basis.omega; + EXPECT_NEAR(center.volume, expected_volume, 1e-10); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp new file mode 100644 index 00000000000..27c2b795746 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp @@ -0,0 +1,53 @@ +#include "sccs_test.h" +#include "../sccs_ionic_charge.h" + +#include "source_base/parallel_reduce.h" +#include "source_cell/cell_tools.h" + +using SccsIonicChargeTest = SccsTest::PwTest; + +TEST_F(SccsIonicChargeTest, ValenceNormalizationAndFourierPhase) +{ + std::vector atoms(2); + atoms[0].position.x = length / 4.0; + atoms[0].valence_charge = 8.0; + atoms[1].valence_charge = 1.0; + const double spread = ModuleSccs::gaussian_ion_spread; + std::vector density; + ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, density); + double integral = 0.0; + for (double value : density) + { + integral += value; + } + Parallel_Reduce::reduce_pool(integral); + integral *= basis.omega / basis.nxyz; + EXPECT_NEAR(integral, 9.0, 1e-12); + + std::vector> charge_g(basis.npw); + basis.real2recip(density.data(), charge_g.data()); + const double exponent = -spread * spread * tpiba * tpiba / 4.0; + const double amplitude = 8.0 / basis.omega * std::exp(exponent); + const double second_amplitude = 1.0 / basis.omega * std::exp(exponent); + double selected = 0.0; + for (int ig = 0; ig < basis.npw; ++ig) + { + const ModuleBase::Vector3& g = basis.gdirect[ig]; + if (g.x == 1.0 && g.y == 0.0 && g.z == 0.0) + { + selected = 1.0; + EXPECT_NEAR(charge_g[ig].real(), second_amplitude, 1e-14); + const double expected_imaginary = -amplitude; + EXPECT_NEAR(charge_g[ig].imag(), expected_imaginary, 1e-14); + } + } + Parallel_Reduce::reduce_pool(selected); + EXPECT_DOUBLE_EQ(selected, 1.0); + atoms[0].position.x += length; + std::vector translated; + ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, translated); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(translated[ir], density[ir], 1e-13); + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp new file mode 100644 index 00000000000..9b482dec854 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp @@ -0,0 +1,97 @@ +#include "sccs_test.h" +#include "../sccs_ionic_force.h" +#include "../sccs_ionic_charge.h" +#include "../sccs_response.h" +#include "../sccs_parameters.h" +#include "../sccs_functional.h" + +#include "source_cell/cell_tools.h" + +using SccsIonicForceTest = SccsTest::PwTest; + +TEST_F(SccsIonicForceTest, AnalyticFourierForceAndGaugeInvariance) +{ + std::vector atoms(1); + atoms[0].valence_charge = 2.0; + atoms[0].position.x = length / 4.0; + std::vector potential = cosine_mode(0); + const double spread = ModuleSccs::gaussian_ion_spread; + std::vector> forces; + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + const double exponent = -0.25 * spread * spread * tpiba * tpiba; + const double expected = 2.0 * tpiba * std::exp(exponent); + EXPECT_NEAR(forces[0].x, expected, 1e-12); + EXPECT_NEAR(forces[0].y, 0.0, 1e-12); + EXPECT_NEAR(forces[0].z, 0.0, 1e-12); + atoms[0].position.x += length; + for (double& value : potential) { value += 4.0; } + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + EXPECT_NEAR(forces[0].x, expected, 1e-12); + potential.assign(basis.nrxx, 3.0); + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + EXPECT_NEAR(forces[0].x, 0.0, 1e-12); + EXPECT_NEAR(forces[0].y, 0.0, 1e-12); + EXPECT_NEAR(forces[0].z, 0.0, 1e-12); +} + +TEST_F(SccsIonicForceTest, ReactionEnergyFiniteDifferenceAtFixedNonuniformCavity) +{ + std::vector atoms(2); + atoms[0].valence_charge = 1.0; + atoms[0].position = ModuleBase::Vector3(2.1, 3.2, 4.3); + atoms[1].valence_charge = 1.0; + atoms[1].position = ModuleBase::Vector3(6.4, 5.1, 2.7); + std::vector density = cosine_mode(0); + for (double& value : density) { value = (2.0 / basis.omega) * (1.0 + 0.1 * value); } + ModuleSccs::SccsConfig config; + config.cavity.density_min = 1e-3; + config.cavity.density_max = 3e-3; + config.cavity.epsilon_bulk = 5.0; + config.surface_regularization = 1e-8; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-13; + solver.tolerance_max = 1e-12; + const double spread = ModuleSccs::gaussian_ion_spread; + const std::vector cold_start; + auto evaluate = [&](const std::vector& displaced, + ModuleSccs::FunctionalResult& functional) { + std::vector charge; + ModuleSccs::gaussian_ionic_density(displaced, basis, tpiba, spread, charge); + for (int ir = 0; ir < basis.nrxx; ++ir) { charge[ir] -= density[ir]; } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, functional); + }; + ModuleSccs::FunctionalResult baseline; + evaluate(atoms, baseline); + std::vector> forces; + ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); + 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 finite_difference = -(positive.reaction_energy - negative.reaction_energy) / (2.0 * step); + EXPECT_NEAR(forces[ia][axis], finite_difference, 1e-8); + } + } + config.cavity.epsilon_bulk = 1.0; + evaluate(atoms, baseline); + ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); + for (const ModuleBase::Vector3& force : forces) + { + EXPECT_NEAR(force.x, 0.0, 1e-12); + EXPECT_NEAR(force.y, 0.0, 1e-12); + EXPECT_NEAR(force.z, 0.0, 1e-12); + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp new file mode 100644 index 00000000000..aeedf8f4447 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -0,0 +1,43 @@ +#include "../sccs_parameters.h" +#include + +TEST(SccsParameters, ExactPresetNamesOnly) +{ + EXPECT_EQ(ModuleSccs::parse_preset("water-neutral"), ModuleSccs::Preset::WaterNeutral); + EXPECT_EQ(ModuleSccs::parse_preset("water-cation"), ModuleSccs::Preset::WaterCation); + EXPECT_EQ(ModuleSccs::parse_preset("water-anion"), ModuleSccs::Preset::WaterAnion); + EXPECT_EQ(ModuleSccs::parse_preset("vacuum"), ModuleSccs::Preset::Vacuum); + EXPECT_EQ(ModuleSccs::parse_preset("custom"), ModuleSccs::Preset::Custom); + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleSccs::parse_preset("unknown"), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(ModuleSccs::parse_preset("Water-Neutral"), ::testing::ExitedWithCode(1), ""); + testing::internal::GetCapturedStdout(); +} + +// The water presets in Hartree atomic units, which also fixes the dyn/cm and +// GPa conversions, and the vacuum preset. +TEST(SccsParameters, PresetValues) +{ + const ModuleSccs::Preset presets[] = {ModuleSccs::Preset::WaterNeutral, + ModuleSccs::Preset::WaterCation, + ModuleSccs::Preset::WaterAnion}; + const double minima[] = {1e-4, 2e-4, 2.4e-3}; + const double maxima[] = {5e-3, 3.5e-3, 1.55e-2}; + const double tensions[] = {3.076640259577260e-5, 3.211524279308205e-6, 0.0}; + const double pressures[] = {-1.223615131827540e-5, 4.248663652178958e-6, 1.529518914784425e-5}; + for (int index = 0; index < 3; ++index) + { + const ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(presets[index]); + EXPECT_DOUBLE_EQ(config.cavity.density_min, minima[index]); + EXPECT_DOUBLE_EQ(config.cavity.density_max, maxima[index]); + EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); + EXPECT_DOUBLE_EQ(config.surface_regularization, 1e-8); + EXPECT_NEAR(config.surface_tension, tensions[index], 1e-18); + EXPECT_NEAR(config.pressure, pressures[index], 1e-18); + } + const ModuleSccs::SccsConfig vacuum = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + EXPECT_DOUBLE_EQ(vacuum.cavity.epsilon_bulk, 1.0); + EXPECT_DOUBLE_EQ(vacuum.surface_tension, 0.0); + EXPECT_DOUBLE_EQ(vacuum.pressure, 0.0); +} + diff --git a/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp new file mode 100644 index 00000000000..a36b7247fef --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp @@ -0,0 +1,68 @@ +#include "sccs_test.h" +#include "../sccs_pw_coulomb.h" + +#include "source_base/parallel_global.h" + +namespace SccsTest +{ +int pool_size = 1; +int pool_rank = 0; +} + +using SccsPwCoulombTest = SccsTest::PwTest; + +TEST_F(SccsPwCoulombTest, FourierModeAndConstantBackground) +{ + const std::vector mode = cosine_mode(0); + std::vector charge = mode; + for (double& value : charge) + { + value += 0.3; + } + ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); + std::vector potential; + coulomb.apply_potential(charge, potential); + const double kernel = ModuleBase::FOUR_PI / (tpiba * tpiba); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = kernel * mode[ir]; + EXPECT_NEAR(potential[ir], expected, 1e-11); + } + charge.assign(basis.nrxx, 1.0); + coulomb.apply_potential(charge, potential); + for (double value : potential) + { + EXPECT_NEAR(value, 0.0, 1e-12); + } +} + +int main(int argc, char** argv) +{ + int threads = 1; + Parallel_Global::read_pal_param(argc, argv, SccsTest::pool_size, threads, SccsTest::pool_rank); + int band_size = 1; + int band_rank = 0; + int band_group = 0; + int pool = 0; +#ifdef __MPI + // These tests create only pool/band groups, not diagonalization/grid groups. + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; + Parallel_Global::init_pools(SccsTest::pool_size, + SccsTest::pool_rank, + 1, + 1, + band_size, + band_rank, + band_group, + SccsTest::pool_size, + SccsTest::pool_rank, + pool); +#endif + testing::InitGoogleTest(&argc, argv); + const int status = RUN_ALL_TESTS(); +#ifdef __MPI + Parallel_Global::finalize_mpi(); +#endif + return status; +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp new file mode 100644 index 00000000000..1ff9242bc88 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -0,0 +1,149 @@ +#include "sccs_test.h" +#include "../sccs_response.h" +#include "../sccs_parameters.h" + +#include "source_base/parallel_reduce.h" + +#include + +using SccsResponseTest = SccsTest::PwTest; + +TEST_F(SccsResponseTest, VacuumAndUniformDielectricAnalyticLimits) +{ + const double dielectrics[] = {1.0, 5.0}; + const std::vector charge = cosine_mode(0); + const std::vector density(basis.nrxx, 0.0); + const std::vector cold_start; + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + for (double epsilon : dielectrics) + { + cavity.epsilon_bulk = epsilon; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); + const double kernel = ModuleBase::FOUR_PI / (epsilon * tpiba * tpiba); + EXPECT_EQ(response.polarization.iterations, 1); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = kernel * charge[ir]; + EXPECT_NEAR(response.polarization.potential[ir], expected, 1e-11); + EXPECT_DOUBLE_EQ(response.cavity_potential[ir], 0.0); + } + } +} + +TEST_F(SccsResponseTest, ManufacturedNonuniformDielectricSolution) +{ + const std::vector mode_x = cosine_mode(0); + const std::vector mode_y = cosine_mode(1); + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 78.3; + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 + 2e-5 * mode_x[ir]; + const ModuleSccs::CavityPoint point = ModuleSccs::evaluate_cavity(density[ir], cavity); + // v=cos(ky), eps=eps(x): -div(eps grad v)/(4 pi)=eps k^2 cos(ky)/(4 pi). + charge[ir] = point.epsilon * tpiba * tpiba * mode_y[ir] / ModuleBase::FOUR_PI; + } + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + const std::vector cold_start; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); + double maximum_error = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double difference = response.polarization.potential[ir] - mode_y[ir]; + const double magnitude = std::abs(difference); + maximum_error = std::max(maximum_error, magnitude); + EXPECT_NEAR(response.polarization.gradient[ir].x, 0.0, 1e-9); + EXPECT_NEAR(response.polarization.gradient[ir].z, 0.0, 1e-9); + } + Parallel_Reduce::reduce_max_pool(basis.poolnproc, maximum_error); + EXPECT_LT(maximum_error, 1e-9); + EXPECT_GT(response.polarization.iterations, 1); +} + +TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) +{ + const std::vector mode_x = cosine_mode(0); + std::vector charge = cosine_mode(1); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 0.009 + 0.008 * mode_x[ir]; + charge[ir] *= 1e-3; + } + ModuleSccs::CavityParameters cavity; + cavity.density_min = 0.0024; + cavity.density_max = 0.0155; + cavity.epsilon_bulk = 78.3; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-11; + solver.tolerance_max = 1e-10; + const std::vector cold_start; + ModuleSccs::SccsResponse response; + // A missed tolerance stops the run; a forked death test cannot share the + // pool collectives of a multi-rank run. + if (SccsTest::pool_size == 1) + { + solver.max_iterations = 1; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, + response), + ::testing::ExitedWithCode(1), ""); + testing::internal::GetCapturedStdout(); + } + solver.max_iterations = 200; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, response); + EXPECT_LE(response.polarization.residual_rms, solver.tolerance_rms); + EXPECT_LE(response.polarization.residual_max, solver.tolerance_max); + for (double& value : density) + { + value *= 1.001; + } + ModuleSccs::SccsResponse cold; + ModuleSccs::SccsResponse warm; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, cold); + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, response.restart_potential, + basis, tpiba, warm); + EXPECT_TRUE(warm.polarization.warm_started); + EXPECT_LT(warm.polarization.iterations, cold.polarization.iterations); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(warm.polarization.potential[ir], cold.polarization.potential[ir], 1e-9); + } + + // Each of the two tolerances alone drives the iteration count. + solver.tolerance_rms = 1e-5; + solver.tolerance_max = 1e-5; + ModuleSccs::SccsResponse loose; + ModuleSccs::SccsResponse tight_maximum; + ModuleSccs::SccsResponse tight_rms; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, loose); + solver.tolerance_max = 1e-11; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_maximum); + solver.tolerance_rms = 1e-11; + solver.tolerance_max = 1.0; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_rms); + EXPECT_LE(tight_maximum.polarization.residual_max, 1e-11); + EXPECT_LE(tight_rms.polarization.residual_rms, 1e-11); + EXPECT_GT(tight_maximum.polarization.iterations, loose.polarization.iterations); + EXPECT_GT(tight_rms.polarization.iterations, loose.polarization.iterations); +} + diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index 76da6b23eba..e3c52fbd540 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -86,6 +86,21 @@ void Potential::cal_fixed_v(double* vl_pseudo) { } +void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = 0.0; + non_electrostatic = 0.0; +} + +void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ +} + +const std::vector* Potential::solvent_electrostatic_potential() const +{ + return nullptr; +} + } // namespace elecstate #include "source_cell/klist.h" diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 488235c6b1e..afa0b154b90 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -131,6 +131,7 @@ add_library( module_parameter/read_inp_dfpt.cpp module_parameter/read_inp_deepks.cpp module_parameter/read_inp_model.cpp + module_parameter/read_inp_sccs.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_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 22f78479a30..b070d9db8bc 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -167,6 +167,7 @@ void ctrl_output_fp(UnitCell& ucell, fn += geom_block + ".cube"; const std::vector* pcc_potential = pelec->pot->pcc_electron_potential(); + const std::vector* solvent_potential = pelec->pot->solvent_electrostatic_potential(); ModuleIO::write_elecstat_pot( #ifdef __MPI pw_big->bz, @@ -180,6 +181,7 @@ void ctrl_output_fp(UnitCell& ucell, pelec->pot->get_fixed_v(), solvent, pcc_potential, + solvent_potential, inp.out_pot[1], ofs_running); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 8f07b087a3f..036adc51fe3 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -576,12 +576,21 @@ struct Input_para double block_height = 0.1; ///< height of the block // implicit solvation model Menglin Sun added on 2022-04-04 - bool imp_sol = false; ///< true: implicit solvation correction; false: - ///< vacuum calculation(default) + int imp_sol = 0; ///< 0: vacuum, 1: legacy solvation, 2: SCCS double eb_k = 80; ///< the relative permittivity of the bulk solvent double tau = 1.0798e-05; ///< the effective surface tension parameter double sigma_k = 0.6; ///< the width of the diffuse cavity double nc_k = 0.00037; ///< the cut-off charge density + std::string sccs_preset = "custom"; + double sccs_epsilon = 1.0; // Environ env_static_permittivity default + double sccs_rho_min = 1e-4; + double sccs_rho_max = 5e-3; + double sccs_gamma = 0.0; // dyn/cm + double sccs_pressure = 0.0; // GPa + int sccs_maxiter = 200; + double sccs_tol_rms = 1e-10; + double sccs_tol_max = 1e-8; + double sccs_surface_eta = 1e-8; // Bohr^-1 // ============== #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 aa1d232feb5..892b235744a 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -2,6 +2,7 @@ #include "source_base/tool_quit.h" #include "read_input.h" #include "read_input_tool.h" +#include "read_inp_sccs.h" namespace ModuleIO { void ReadInput::item_model() @@ -181,16 +182,25 @@ void ReadInput::item_model() this->add_item(item); } - // imlicit_solvation + // Implicit solvation: imp_sol, the original model, then SCCS. { Input_Item item("imp_sol"); - item.annotation = "calculate implicit solvation correction or not"; + item.annotation = "implicit solvent model"; item.category = "Implicit solvation model"; - item.type = "Boolean"; - item.description = "Calculate implicit solvation correction"; - item.default_value = "False"; - item.unit = ""; - read_sync_bool(input.imp_sol); + item.type = "Integer"; + item.description = "Select 0 for vacuum, 1 for the original ABACUS implicit solvation model, or 2 for SCCS."; + item.default_value = "0"; + item.read_value = [](const Input_Item& item, Parameter& para) { + std::string error; + const bool valid = parse_solvation_model(item.str_values[0], para.input.imp_sol, error); + if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } + }; + sync_int(input.imp_sol); + item.check_value = [](const Input_Item&, const Parameter& para) { + std::string error; + const bool valid = validate_sccs_input(para.input, error); + if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } + }; this->add_item(item); } { @@ -201,7 +211,7 @@ void ReadInput::item_model() item.description = "The relative permittivity of the bulk solvent, 80 for water"; item.default_value = "80"; item.unit = ""; - item.set_availability("imp_sol==true"); + item.set_availability("imp_sol==1"); read_sync_double(input.eb_k); this->add_item(item); } @@ -238,6 +248,7 @@ void ReadInput::item_model() read_sync_double(input.nc_k); this->add_item(item); } + this->item_sccs(); // vdW Correction { diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp new file mode 100644 index 00000000000..5b317d9975d --- /dev/null +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -0,0 +1,220 @@ +#include "read_inp_sccs.h" +#include "read_input.h" +#include "read_input_tool.h" +#include "source_base/tool_quit.h" + +#include +#include + +namespace ModuleIO +{ +bool parse_solvation_model(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 sccs_context_error(const Input_para& input) +{ + std::string error; + if (input.assume_isolated != "none") { error = "SCCS currently requires assume_isolated none"; } + 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; +} + +// Cavity and surface model parameters, starting with the preset. +std::string sccs_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 sccs_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 validate_sccs_input(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 = sccs_context_error(input); + if (error.empty()) { error = sccs_model_error(input); } + if (error.empty()) { error = sccs_solver_error(input); } + return error.empty(); +} + +void ReadInput::item_sccs() +{ + { + Input_Item item("sccs_preset"); + item.annotation = "SCCS parameter preset"; + item.category = "Implicit solvation model"; + item.type = "String"; + item.description = "SCCS parameter preset. Allowed values: custom, vacuum, water-neutral, water-cation, water-anion. Non-custom presets override sccs_epsilon, sccs_rho_min, sccs_rho_max, sccs_gamma and sccs_pressure; solver controls and surface regularization remain user-controlled."; + item.default_value = "custom"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_string(input.sccs_preset); + this->add_item(item); + } + { + Input_Item item("sccs_epsilon"); + item.annotation = "Bulk dielectric constant >= 1 for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Bulk dielectric constant >= 1 for the custom preset. The default 1 (no dielectric) follows the Environ env_static_permittivity default for environ_type input. Water presets use 78.3 and vacuum uses 1."; + item.default_value = "1.0"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_epsilon); + this->add_item(item); + } + { + Input_Item item("sccs_rho_min"); + item.annotation = "Lower electronic-density cavity threshold for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Lower electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 1e-4; water-cation: 2e-4; water-anion: 2.4e-3."; + item.default_value = "1e-4"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_rho_min); + this->add_item(item); + } + { + Input_Item item("sccs_rho_max"); + item.annotation = "Upper electronic-density cavity threshold for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Upper electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 5e-3; water-cation: 3.5e-3; water-anion: 1.55e-2."; + item.default_value = "5e-3"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_rho_max); + this->add_item(item); + } + { + Input_Item item("sccs_gamma"); + item.annotation = "Surface coefficient for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Surface coefficient for the custom preset. Water-neutral: 47.9; water-cation: 5; water-anion/vacuum: 0."; + item.default_value = "0.0"; + item.unit = "dyn/cm"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_gamma); + this->add_item(item); + } + { + Input_Item item("sccs_pressure"); + item.annotation = "Volume coefficient for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Volume coefficient for the custom preset. Water-neutral: -0.36; water-cation: 0.125; water-anion: 0.45; vacuum: 0."; + item.default_value = "0.0"; + item.unit = "GPa"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_pressure); + this->add_item(item); + } + { + Input_Item item("sccs_maxiter"); + item.annotation = "Positive maximum number of SCCS sqrt-CG iterations"; + item.category = "Implicit solvation model"; + item.type = "Integer"; + item.description = "Positive maximum number of SCCS sqrt-CG iterations. Failure to satisfy both residual tolerances terminates the calculation."; + item.default_value = "200"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_int(input.sccs_maxiter); + this->add_item(item); + } + { + Input_Item item("sccs_tol_rms"); + item.annotation = "Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver."; + item.default_value = "1e-10"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_tol_rms); + this->add_item(item); + } + { + Input_Item item("sccs_tol_max"); + item.annotation = "Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver."; + item.default_value = "1e-8"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_tol_max); + this->add_item(item); + } + { + Input_Item item("sccs_surface_eta"); + item.annotation = "Positive regularization of the SCCS surface gradient norm"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive regularization of the SCCS surface gradient norm."; + item.default_value = "1e-8"; + item.unit = "bohr^-1"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_surface_eta); + 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 new file mode 100644 index 00000000000..8736ed5768e --- /dev/null +++ b/source/source_io/module_parameter/read_inp_sccs.h @@ -0,0 +1,10 @@ +#ifndef READ_INP_SCCS_H +#define READ_INP_SCCS_H +#include +struct Input_para; +namespace ModuleIO +{ +bool parse_solvation_model(const std::string& value, int& model, std::string& error); +bool validate_sccs_input(const Input_para& input, std::string& error); +} +#endif diff --git a/source/source_io/module_parameter/read_input.h b/source/source_io/module_parameter/read_input.h index ac180dcd4ea..941dde06734 100644 --- a/source/source_io/module_parameter/read_input.h +++ b/source/source_io/module_parameter/read_input.h @@ -146,6 +146,7 @@ class ReadInput void item_postprocess(); // items for some models void item_model(); + void item_sccs(); // items for exx void item_exx(); // items for dft+u diff --git a/source/source_io/test_serial/CMakeLists.txt b/source/source_io/test_serial/CMakeLists.txt index 3238fb743d4..0201fb14e01 100644 --- a/source/source_io/test_serial/CMakeLists.txt +++ b/source/source_io/test_serial/CMakeLists.txt @@ -17,6 +17,7 @@ add_library( ../module_parameter/read_inp_dfpt.cpp ../module_parameter/read_inp_deepks.cpp ../module_parameter/read_inp_model.cpp + ../module_parameter/read_inp_sccs.cpp ../module_parameter/read_inp_postproc.cpp ../module_parameter/read_inp_exx_dftu.cpp ../module_parameter/read_inp_other.cpp @@ -40,7 +41,7 @@ AddTest( AddTest( TARGET MODULE_IO_read_item_serial LIBS parameter base device io_input_serial - SOURCES read_input_item_test.cpp + SOURCES read_input_item_test.cpp test_read_inp_sccs.cpp ) AddTest( diff --git a/source/source_io/test_serial/test_read_inp_sccs.cpp b/source/source_io/test_serial/test_read_inp_sccs.cpp new file mode 100644 index 00000000000..94cb6bad718 --- /dev/null +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -0,0 +1,56 @@ +#include "source_io/module_parameter/read_inp_sccs.h" +#include "source_io/module_parameter/input_parameter.h" +#include + +TEST(ReadInpSccs, LegacyBooleansAndModelSelection) +{ + std::string error; + const std::string yes[] = {"true", "TRUE", "yes", "On", ".true.", "1", "t", "y"}; + const std::string no[] = {"false", "FALSE", "no", "Off", ".false.", "0", "f", "n"}; + int model = 0; + for (const std::string& value : yes) + { + ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + EXPECT_EQ(model, 1); + } + for (const std::string& value : no) + { + ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + EXPECT_EQ(model, 0); + } + ASSERT_TRUE(ModuleIO::parse_solvation_model("2", model, error)); + EXPECT_EQ(model, 2); + EXPECT_FALSE(ModuleIO::parse_solvation_model("3", model, error)); + EXPECT_EQ(model, 2); + EXPECT_FALSE(ModuleIO::parse_solvation_model("2.0", model, error)); +} + +TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) +{ + Input_para input; + input.imp_sol = 2; + input.device = "cpu"; + input.basis_type = "pw"; + std::string error; + ASSERT_TRUE(ModuleIO::validate_sccs_input(input, error)) << error; + input.cal_force = true; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + input.cal_stress = true; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.cal_stress = false; + input.cal_force = false; + input.assume_isolated = "pcc_0d"; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.assume_isolated = "none"; + input.sccs_rho_min = input.sccs_rho_max; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.sccs_rho_min = 1e-4; + input.sccs_preset = "unknown"; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.sccs_preset = "vacuum"; + input.basis_type = "lcao"; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + input.imp_sol = 1; + input.sccs_preset = "unused"; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); +} diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 4bd110ec348..4eb668ae8f0 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -190,6 +190,8 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, { parts.fpcc.create(ucell.nat, 3); pelec->pot->add_pcc_force(ucell, parts.fpcc); + if (parts.fsol.nr == 0) { parts.fsol.create(ucell.nat, 3); } + pelec->pot->add_solvation_force(ucell, parts.fsol); } // DFT+U force/stress diff --git a/source/source_lcao/force_stress_terms.cpp b/source/source_lcao/force_stress_terms.cpp index 6daa6f59d26..2dde8ff9ec7 100644 --- a/source/source_lcao/force_stress_terms.cpp +++ b/source/source_lcao/force_stress_terms.cpp @@ -93,7 +93,7 @@ void cal_external_field_forces(UnitCell& ucell, } //! atomic forces from implicit solvation model - if (PARAM.inp.imp_sol) + if (PARAM.inp.imp_sol == 1) { parts.fsol.create(ucell.nat, 3); solvent.cal_force_sol(ucell, rhopw, locpp.vloc, PARAM.inp.nspin, parts.fsol); diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index cda9d9ecacd..6c31e1f1f6a 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -88,7 +88,8 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, if (inp.vion_in_h) { pot_register_in.push_back("local"); } if (inp.vh_in_h) { pot_register_in.push_back("hartree"); } pot_register_in.push_back("xc"); - if (inp.imp_sol) { pot_register_in.push_back("surchem"); } + if (inp.imp_sol == 1) { pot_register_in.push_back("surchem"); } + if (inp.imp_sol == 2) { pot_register_in.push_back("sccs"); } if (inp.assume_isolated == "pcc_0d") { pot_register_in.push_back("pcc"); } if (inp.assume_isolated == "pcc_2d") { diff --git a/source/source_lcao/hamilt_lcao_factory.cpp b/source/source_lcao/hamilt_lcao_factory.cpp index 4f4f859901a..d1fd465e0f1 100644 --- a/source/source_lcao/hamilt_lcao_factory.cpp +++ b/source/source_lcao/hamilt_lcao_factory.cpp @@ -168,7 +168,7 @@ LcaoOpsBundle build_gamma_ops(const UnitCell& ucell, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, inp); // effective potential term Operator* veff = new Veff>(hsk, kv->kvec_d, pot_in, @@ -230,7 +230,7 @@ LcaoOpsBundle build_multik_ops(const UnitCell& ucell, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, inp); // Veff term ops = new Veff>(hsk, kv->kvec_d, pot_in, diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index f1d0d8eac95..d81997b89af 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -137,13 +137,21 @@ void Forces::cal_force(const int nspin, } } + const int solvation_model = PARAM.inp.imp_sol; ModuleBase::matrix forcesol; ModuleBase::matrix forcepcc(this->nat, 3); elec.pot->add_pcc_force(ucell, forcepcc); - if (PARAM.inp.imp_sol) + if (solvation_model != 0) { forcesol.create(this->nat, 3); - solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, nspin, forcesol); + if (solvation_model == 1) + { + solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, nspin, forcesol); + } + else + { + elec.pot->add_solvation_force(ucell, forcesol); + } if (PARAM.inp.test_force) { ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); @@ -177,7 +185,7 @@ void Forces::cal_force(const int nspin, force(iat, ipol) = force(iat, ipol) + force_gate(iat, ipol); } - if (PARAM.inp.imp_sol) + if (solvation_model != 0) { force(iat, ipol) = force(iat, ipol) + forcesol(iat, ipol); } @@ -261,7 +269,7 @@ void Forces::cal_force(const int nspin, force_gate, false); } - if (PARAM.inp.imp_sol) + if (solvation_model != 0) { ModuleIO::print_force(GlobalV::ofs_running, ucell, diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index 13977cd9f6e..d3d284643f5 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -83,10 +83,11 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, } // no variable can choose xc, maybe it is necessary pot_register_in.push_back("xc"); - if (input.imp_sol) + if (input.imp_sol == 1) { pot_register_in.push_back("surchem"); } + if (input.imp_sol == 2) { pot_register_in.push_back("sccs"); } if (input.assume_isolated == "pcc_0d") { pot_register_in.push_back("pcc"); @@ -122,7 +123,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, input); Operator* veff = new Veff>(isk, pot_in->get_veff_smooth_data(), pot_in->get_veff_smooth().nr, diff --git a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp index 7023d90719b..78dfc1a44aa 100644 --- a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp +++ b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp @@ -6,6 +6,9 @@ #include elecstate::Potential::~Potential(){} +void elecstate::Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void elecstate::Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* elecstate::Potential::solvent_electrostatic_potential() const { return nullptr; } void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase::matrix&){} void elecstate::Potential::cal_fixed_v(double*){} diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT new file mode 100644 index 00000000000..8eb0415ef14 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT @@ -0,0 +1,21 @@ +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 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README new file mode 100644 index 00000000000..44f4976649c --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README @@ -0,0 +1,2 @@ +PW neutral water with periodic SCCS (water-neutral preset), ecutwfc 20 and ecutrho 240. +Checks total energy, forces and electrostatic, surface and volume solvation terms. diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU new file mode 100644 index 00000000000..ed974dd43a5 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/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/220_PW_SCCS_PERIODIC_WATER/result.ref b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref new file mode 100644 index 00000000000..a1c90a54215 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref @@ -0,0 +1,6 @@ +etotref -441.6484359463214560 +etotperatomref -147.2161453154 +totalforceref 16.644169 +esolelref -0.2883000404 +esolcavref 0.0807425541 +totaltimeref 3.33 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index b0d6cade3e8..629c6782d34 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -125,3 +125,4 @@ scf_deltaspin2 scf_deltaspin4 221_PW_PCC0D_H3O 222_PW_PCC2D_H3O +220_PW_SCCS_PERIODIC_WATER diff --git a/tests/integrate/validation_tools/props_basic.sh b/tests/integrate/validation_tools/props_basic.sh index a73902d0978..992bf596816 100644 --- a/tests/integrate/validation_tools/props_basic.sh +++ b/tests/integrate/validation_tools/props_basic.sh @@ -94,7 +94,7 @@ run_basic_props_post_ml(){ #-------------------------------------------- # implicit solvation model #-------------------------------------------- -if ! test -z "$imp_sol" && [ $imp_sol == 1 ]; then +if ! test -z "$imp_sol" && { [ "$imp_sol" = 1 ] || [ "$imp_sol" = 2 ]; }; then esol_el=`grep E_sol_el $running_path | awk '{print $3}'` esol_cav=`grep E_sol_cav $running_path | awk '{print $3}'` echo "esolelref $esol_el" >>$props_result_file From 25f8f40911537163c31166643cbcb4bad7187a10 Mon Sep 17 00:00:00 2001 From: lyt Date: Sun, 18 Oct 2026 16:36:16 +0800 Subject: [PATCH 2/5] Address review comments for SCCS step 1 --- source/source_esolver/esolver_double_xc.cpp | 2 +- source/source_esolver/esolver_ks_lcao.cpp | 3 +- .../module_sccs/sccs_response.cpp | 5 + .../source_hamilt/module_sccs/sccs_response.h | 4 +- .../module_sccs/test/test_sccs_parameters.cpp | 1 - .../module_sccs/test/test_sccs_response.cpp | 106 ++++++++++++++++++ .../module_parameter/read_inp_model.cpp | 11 +- .../module_parameter/read_inp_sccs.cpp | 23 ++++ .../test_serial/test_read_inp_sccs.cpp | 38 +++++++ source/source_lcao/force_stress_lcao.cpp | 8 +- source/source_lcao/force_stress_lcao.h | 1 + tests/02_NAO_Gamma/CASES_CPU.txt | 1 + tests/02_NAO_Gamma/scf_sccs_water/INPUT | 27 +++++ tests/02_NAO_Gamma/scf_sccs_water/KPT | 4 + tests/02_NAO_Gamma/scf_sccs_water/README | 6 + tests/02_NAO_Gamma/scf_sccs_water/STRU | 28 +++++ tests/02_NAO_Gamma/scf_sccs_water/result.ref | 6 + 17 files changed, 267 insertions(+), 7 deletions(-) create mode 100644 tests/02_NAO_Gamma/scf_sccs_water/INPUT create mode 100644 tests/02_NAO_Gamma/scf_sccs_water/KPT create mode 100644 tests/02_NAO_Gamma/scf_sccs_water/README create mode 100644 tests/02_NAO_Gamma/scf_sccs_water/STRU create mode 100644 tests/02_NAO_Gamma/scf_sccs_water/result.ref diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 36b8009f764..5afe888d977 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -427,7 +427,7 @@ void ESolver_DoubleXC::cal_force(BaseCell& basecell, ModuleBase::matrix& const FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h, this->inp_->sc_mag_switch, this->inp_->device, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, - PARAM.globalv.gamma_only_pw}; + PARAM.globalv.gamma_only_pw, this->inp_->imp_sol}; fsl.getForceStress(ucell, this->get_vdw_result(), diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index e6c047a40c2..bb07dc88947 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -362,7 +362,8 @@ void ESolver_KS_LCAO::cal_force(BaseCell& basecell, ModuleBase::matrix& PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, - PARAM.globalv.gamma_only_pw}; + PARAM.globalv.gamma_only_pw, + this->inp_->imp_sol}; fsl.getForceStress(ucell, this->get_vdw_result(), diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index 58d247bce21..00fe82fc18b 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -265,6 +265,11 @@ void finish_response(const ModulePW::PW_Basis& basis, polarization.gradient.resize(size); XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); polarization.potential.swap(potential); + // Continuum derivative of the reaction energy, as in Environ. It is not the + // exact derivative of the discrete energy: at uniform density that is + // eps'/(16 pi) (2 v L v - L v^2) with the spectral Laplacian L, and v^2 + // aliases on the FFT grid. The two agree once v^2 is resolved by the grid + // (test CavityDerivativeConvergesWithGridRefinement). response.cavity_potential.resize(size); for (std::size_t i = 0; i < size; ++i) { diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index 0e9b5a89a98..4aea75c2428 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -33,7 +33,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 + // -epsilon' |grad v|^2/(8 pi), Hartree: the continuum derivative, exact only + // when v^2 is resolved by the grid (see finish_response). + std::vector cavity_potential; std::vector restart_potential; // unshifted solution for explicit warm starts }; diff --git a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp index aeedf8f4447..46a4c4bf975 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -40,4 +40,3 @@ TEST(SccsParameters, PresetValues) EXPECT_DOUBLE_EQ(vacuum.surface_tension, 0.0); EXPECT_DOUBLE_EQ(vacuum.pressure, 0.0); } - diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp index 1ff9242bc88..7a4b161b3d7 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -1,10 +1,13 @@ #include "sccs_test.h" +#include "../sccs_functional.h" #include "../sccs_response.h" #include "../sccs_parameters.h" #include "source_base/parallel_reduce.h" #include +#include +#include using SccsResponseTest = SccsTest::PwTest; @@ -147,3 +150,106 @@ TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) EXPECT_GT(tight_rms.polarization.iterations, loose.polarization.iterations); } +namespace +{ +// cos(2 pi harmonic x / L) on a basis' local real-space grid. +std::vector x_harmonic(const ModulePW::PW_Basis& basis, int harmonic) +{ + std::vector values(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const double angle = ModuleBase::TWO_PI * harmonic * ix / basis.nx; + values[ir] = std::cos(angle); + } + return values; +} + +// Relative difference between the analytic cavity potential and a central +// finite difference of the discrete reaction energy along a varying cavity +// density, on the grid set by gridecut. The base density is uniform inside the +// switching interval, so v = 0.1 cos(6 g x) solves the base response exactly. +double cavity_derivative_mismatch(double gridecut) +{ + ModulePW::PW_Basis basis("cpu", "double"); +#ifdef __MPI + basis.initmpi(SccsTest::pool_size, SccsTest::pool_rank, POOL_WORLD); +#endif + const double length = 10.0; + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + basis.initgrids(length, lattice, gridecut); + basis.initparameters(false, gridecut, 1, false); + basis.setuptransform(); + basis.collect_local_pw(); + const double tpiba = ModuleBase::TWO_PI / length; + + ModuleSccs::SccsConfig config; + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + config.cavity.epsilon_bulk = 78.3; + config.surface_regularization = 1e-8; + // On the perturbed cavity the CG stalls at a uniform residual (~2e-9 here) + // that the zero-mean Coulomb preconditioner cannot reduce; a uniform + // residual does not change the zero-mean potential or the energy, so 1e-8 + // is tight enough and the mismatches below are unchanged by it. + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-8; + solver.tolerance_max = 1e-8; + const double base_density = 7e-4; + const double base_epsilon = ModuleSccs::evaluate_cavity(base_density, config.cavity).epsilon; + const double wavevector = 6.0 * tpiba; + const std::vector mode6 = x_harmonic(basis, 6); + const std::vector mode3 = x_harmonic(basis, 3); + std::vector charge(basis.nrxx); + std::vector direction(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + charge[ir] = 0.1 * base_epsilon * wavevector * wavevector * mode6[ir] / ModuleBase::FOUR_PI; + direction[ir] = 1.0 + mode3[ir]; + } + const std::vector cold_start; + auto reaction_energy = [&](double step, ModuleSccs::SccsResponse& response) -> double { + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = base_density + step * direction[ir]; + } + ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, functional); + return functional.reaction_energy; + }; + ModuleSccs::SccsResponse center; + ModuleSccs::SccsResponse plus; + ModuleSccs::SccsResponse minus; + reaction_energy(0.0, center); + const double step = 1e-7; + const double backward_step = -step; + const double energy_plus = reaction_energy(step, plus); + const double energy_minus = reaction_energy(backward_step, minus); + const double finite_difference = (energy_plus - energy_minus) / (2.0 * step); + double analytic = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + analytic += center.cavity_potential[ir] * direction[ir]; + } + Parallel_Reduce::reduce_pool(analytic); + analytic *= basis.omega / basis.nxyz; + const double difference = finite_difference - analytic; + return std::abs(difference) / std::abs(finite_difference); +} +} // namespace + +// The cavity potential is the continuum derivative -eps'|grad v|^2/(8 pi). On +// the 15^3 grid (gridecut 20 Ry) v^2 = cos^2(6 g x) aliases, and the discrete +// reaction energy differs from it by about 65% along 1 + cos(3 g x); on the +// 30^3 grid (gridecut 80 Ry) v^2 is resolved and the two agree to O(step^2). +TEST(SccsResponseRefinement, CavityDerivativeConvergesWithGridRefinement) +{ + const double coarse = cavity_derivative_mismatch(20.0); + const double fine = cavity_derivative_mismatch(80.0); + EXPECT_GT(coarse, 0.1); + EXPECT_LT(fine, 1e-6); +} diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index 892b235744a..1a5a1a302ef 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -188,7 +188,16 @@ void ReadInput::item_model() item.annotation = "implicit solvent model"; item.category = "Implicit solvation model"; item.type = "Integer"; - item.description = "Select 0 for vacuum, 1 for the original ABACUS implicit solvation model, or 2 for SCCS."; + item.description = "Implicit solvation model:\n" + "* 0: no implicit solvent.\n" + "* 1: the original ABACUS implicit solvation model.\n" + "* 2: the self-consistent continuum solvation model (SCCS). Supported scope: CPU KS-DFT " + "(esolver_type ksdft, device cpu) with basis_type pw or lcao; calculation scf energies and " + "forces, no stress (cal_stress is rejected); nspin 1 or 2; periodic cells with " + "assume_isolated none; no electric/gate field, DFT-1/2, DeePKS or dm_to_rho. SCCS is meant " + "for neutral cells: 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" + "Legacy Boolean values (true/false) are read as 1/0."; item.default_value = "0"; item.read_value = [](const Input_Item& item, Parameter& para) { std::string error; diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 5b317d9975d..0a201083554 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -5,6 +5,8 @@ #include #include +#include +#include namespace ModuleIO { @@ -52,6 +54,26 @@ std::string sccs_context_error(const Input_para& input) 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 sccs_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 sccs_model_error(const Input_para& input) { @@ -89,6 +111,7 @@ bool validate_sccs_input(const Input_para& input, std::string& error) } if (input.imp_sol != 2) { return true; } error = sccs_context_error(input); + if (error.empty()) { error = sccs_nonfinite_error(input); } if (error.empty()) { error = sccs_model_error(input); } if (error.empty()) { error = sccs_solver_error(input); } return error.empty(); 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 94cb6bad718..c8411484db3 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,10 @@ #include "source_io/module_parameter/input_parameter.h" #include +#include +#include +#include + TEST(ReadInpSccs, LegacyBooleansAndModelSelection) { std::string error; @@ -54,3 +58,37 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) input.sccs_preset = "unused"; EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); } + +TEST(ReadInpSccs, RejectsNonfiniteReals) +{ + Input_para input; + input.imp_sol = 2; + input.device = "cpu"; + input.basis_type = "pw"; + std::string error; + // Signed surface and volume coefficients remain valid. + input.sccs_gamma = -1.0; + input.sccs_pressure = -0.36; + ASSERT_TRUE(ModuleIO::validate_sccs_input(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}; + const double nonfinite[] = {std::nan(""), std::numeric_limits::infinity(), + -std::numeric_limits::infinity()}; + for (double* value : reals) + { + const double saved = *value; + for (double bad : nonfinite) + { + *value = bad; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + const std::size_t position = error.find("must be finite"); + EXPECT_NE(position, std::string::npos) << error; + } + *value = saved; + } + // The parser path: std::stod accepts "nan". + input.sccs_surface_eta = std::stod("nan"); + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); +} diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 4eb668ae8f0..21f99dda951 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -190,8 +190,12 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, { parts.fpcc.create(ucell.nat, 3); pelec->pot->add_pcc_force(ucell, parts.fpcc); - if (parts.fsol.nr == 0) { parts.fsol.create(ucell.nat, 3); } - pelec->pot->add_solvation_force(ucell, parts.fsol); + // imp_sol 1 fills fsol in cal_vdw_fields_fs; only SCCS uses the component hook. + if (cfg.imp_sol == 2) + { + parts.fsol.create(ucell.nat, 3); + pelec->pot->add_solvation_force(ucell, parts.fsol); + } } // DFT+U force/stress diff --git a/source/source_lcao/force_stress_lcao.h b/source/source_lcao/force_stress_lcao.h index f40ddb27689..827dfc3db7a 100644 --- a/source/source_lcao/force_stress_lcao.h +++ b/source/source_lcao/force_stress_lcao.h @@ -45,6 +45,7 @@ struct FSCalcConfig bool domag_z; int gga_grad; bool gamma_only_pw; + int imp_sol; // 2 adds the SCCS component force }; // Force/stress component matrices assembled by getForceStress. Grouping them diff --git a/tests/02_NAO_Gamma/CASES_CPU.txt b/tests/02_NAO_Gamma/CASES_CPU.txt index 9ba3fa93de9..37a3e8404fc 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -40,3 +40,4 @@ scf_upf100 test.sum scf_pcc0d_water scf_pcc2d_water +scf_sccs_water diff --git a/tests/02_NAO_Gamma/scf_sccs_water/INPUT b/tests/02_NAO_Gamma/scf_sccs_water/INPUT new file mode 100644 index 00000000000..47c0b63503b --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water/INPUT @@ -0,0 +1,27 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +calculation scf +basis_type lcao +gamma_only 1 +ks_solver genelpa +nbands 8 +ecutwfc 60 +ecutrho 240 +scf_nmax 100 +scf_thr 1.0e-8 +smearing_method fixed +mixing_type pulay +mixing_beta 0.4 +symmetry 0 +init_wfc atomic +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 diff --git a/tests/02_NAO_Gamma/scf_sccs_water/KPT b/tests/02_NAO_Gamma/scf_sccs_water/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/02_NAO_Gamma/scf_sccs_water/README b/tests/02_NAO_Gamma/scf_sccs_water/README new file mode 100644 index 00000000000..a6df6c49e12 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water/README @@ -0,0 +1,6 @@ +LCAO (gamma-only) neutral water with periodic SCCS (water-neutral preset), same geometry as +01_PW/220_PW_SCCS_PERIODIC_WATER. Covers the LCAO path end to end: SCCS potential in the LCAO +Hamiltonian, local Pulay force terms and final force assembly of the SCCS ionic force. +Checks total energy, forces and the electrostatic and cavity solvation terms. +ecutwfc is 60 because LCAO requires ecutrho = 4 * ecutwfc; ecutrho 240 keeps the density grid of +the PW case. result.ref was generated on one MPI process and also passes on four. diff --git a/tests/02_NAO_Gamma/scf_sccs_water/STRU b/tests/02_NAO_Gamma/scf_sccs_water/STRU new file mode 100644 index 00000000000..63a4455cda1 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water/STRU @@ -0,0 +1,28 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +H_gga_8au_100Ry_2s1p.orb +O_gga_7au_100Ry_2s2p1d.orb + +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 1 1 1 +4.567 7.111 6.0 1 1 1 +O +0.0 +1 +6.0 6.0 6.0 1 1 1 diff --git a/tests/02_NAO_Gamma/scf_sccs_water/result.ref b/tests/02_NAO_Gamma/scf_sccs_water/result.ref new file mode 100644 index 00000000000..6c1895e23c2 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water/result.ref @@ -0,0 +1,6 @@ +etotref -466.8169430311206 +etotperatomref -155.6056476770 +totalforceref 3.691279 +esolelref -0.4602106631 +esolcavref 0.0746968006 +totaltimeref 11.42 From 5ba345f7eba4be5fa9059a1fc790dfc798b19754 Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Mon, 19 Oct 2026 21:09:05 +0800 Subject: [PATCH 3/5] Trim SCCS step 1 review tests and comments - Drop RejectsNonfiniteReals, as already done on the upper layers; the finiteness check itself stays. - Keep only the resolved-grid half of the cavity-derivative refinement test: the coarse-grid half asserted the known aliasing mismatch rather than a requirement. - Shorten the continuum-derivative comments. --- .../module_sccs/sccs_response.cpp | 7 ++-- .../source_hamilt/module_sccs/sccs_response.h | 4 +-- .../module_sccs/test/test_sccs_response.cpp | 35 ++++++------------ .../test_serial/test_read_inp_sccs.cpp | 36 ------------------- 4 files changed, 14 insertions(+), 68 deletions(-) diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index 00fe82fc18b..43519550a6b 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -265,11 +265,8 @@ void finish_response(const ModulePW::PW_Basis& basis, polarization.gradient.resize(size); XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); polarization.potential.swap(potential); - // Continuum derivative of the reaction energy, as in Environ. It is not the - // exact derivative of the discrete energy: at uniform density that is - // eps'/(16 pi) (2 v L v - L v^2) with the spectral Laplacian L, and v^2 - // aliases on the FFT grid. The two agree once v^2 is resolved by the grid - // (test CavityDerivativeConvergesWithGridRefinement). + // 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) { diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index 4aea75c2428..0e9b5a89a98 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -33,9 +33,7 @@ struct SccsResponse std::vector depsilon_drho; std::vector> grad_log_epsilon; PolarizationResult polarization; - // -epsilon' |grad v|^2/(8 pi), Hartree: the continuum derivative, exact only - // when v^2 is resolved by the grid (see finish_response). - std::vector cavity_potential; + std::vector cavity_potential; // -epsilon' |grad v|^2/(8 pi), Hartree std::vector restart_potential; // unshifted solution for explicit warm starts }; diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp index 7a4b161b3d7..b3294e9d571 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -164,18 +164,20 @@ std::vector x_harmonic(const ModulePW::PW_Basis& basis, int harmonic) } return values; } +} // namespace -// Relative difference between the analytic cavity potential and a central -// finite difference of the discrete reaction energy along a varying cavity -// density, on the grid set by gridecut. The base density is uniform inside the -// switching interval, so v = 0.1 cos(6 g x) solves the base response exactly. -double cavity_derivative_mismatch(double gridecut) +// The cavity potential -eps'|grad v|^2/(8 pi) against a central finite +// difference of the reaction energy along 1 + cos(3 g x). The base density is +// uniform inside the switching interval, so v = 0.1 cos(6 g x) solves the base +// response exactly; the 80 Ry grid resolves |grad v|^2. +TEST(SccsResponseRefinement, CavityDerivativeMatchesFiniteDifferenceOnAResolvedGrid) { ModulePW::PW_Basis basis("cpu", "double"); #ifdef __MPI basis.initmpi(SccsTest::pool_size, SccsTest::pool_rank, POOL_WORLD); #endif const double length = 10.0; + const double gridecut = 80.0; const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); @@ -190,10 +192,8 @@ double cavity_derivative_mismatch(double gridecut) config.cavity.density_max = 5e-3; config.cavity.epsilon_bulk = 78.3; config.surface_regularization = 1e-8; - // On the perturbed cavity the CG stalls at a uniform residual (~2e-9 here) - // that the zero-mean Coulomb preconditioner cannot reduce; a uniform - // residual does not change the zero-mean potential or the energy, so 1e-8 - // is tight enough and the mismatches below are unchanged by it. + // The CG stalls at a uniform residual (~2e-9) outside the range of the + // zero-mean preconditioner; it does not change the energy. ModuleSccs::PolarizationSolverParameters solver; solver.tolerance_rms = 1e-8; solver.tolerance_max = 1e-8; @@ -237,19 +237,6 @@ double cavity_derivative_mismatch(double gridecut) } Parallel_Reduce::reduce_pool(analytic); analytic *= basis.omega / basis.nxyz; - const double difference = finite_difference - analytic; - return std::abs(difference) / std::abs(finite_difference); -} -} // namespace - -// The cavity potential is the continuum derivative -eps'|grad v|^2/(8 pi). On -// the 15^3 grid (gridecut 20 Ry) v^2 = cos^2(6 g x) aliases, and the discrete -// reaction energy differs from it by about 65% along 1 + cos(3 g x); on the -// 30^3 grid (gridecut 80 Ry) v^2 is resolved and the two agree to O(step^2). -TEST(SccsResponseRefinement, CavityDerivativeConvergesWithGridRefinement) -{ - const double coarse = cavity_derivative_mismatch(20.0); - const double fine = cavity_derivative_mismatch(80.0); - EXPECT_GT(coarse, 0.1); - EXPECT_LT(fine, 1e-6); + const double tolerance = 1e-6 * std::abs(finite_difference); + EXPECT_NEAR(analytic, finite_difference, tolerance); } 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 c8411484db3..d61292c6ce5 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -2,8 +2,6 @@ #include "source_io/module_parameter/input_parameter.h" #include -#include -#include #include TEST(ReadInpSccs, LegacyBooleansAndModelSelection) @@ -58,37 +56,3 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) input.sccs_preset = "unused"; EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); } - -TEST(ReadInpSccs, RejectsNonfiniteReals) -{ - Input_para input; - input.imp_sol = 2; - input.device = "cpu"; - input.basis_type = "pw"; - std::string error; - // Signed surface and volume coefficients remain valid. - input.sccs_gamma = -1.0; - input.sccs_pressure = -0.36; - ASSERT_TRUE(ModuleIO::validate_sccs_input(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}; - const double nonfinite[] = {std::nan(""), std::numeric_limits::infinity(), - -std::numeric_limits::infinity()}; - for (double* value : reals) - { - const double saved = *value; - for (double bad : nonfinite) - { - *value = bad; - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); - const std::size_t position = error.find("must be finite"); - EXPECT_NE(position, std::string::npos) << error; - } - *value = saved; - } - // The parser path: std::stod accepts "nan". - input.sccs_surface_eta = std::stod("nan"); - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); -} From 7367e0ac1c960e679f557aab38249596ec247200 Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Tue, 20 Oct 2026 10:33:05 +0800 Subject: [PATCH 4/5] Address mohanchen's review of SCCS step 1 - Use ModuleBase::HARTREE_SI and ModuleBase::BOHR_RADIUS_SI from constants.h for the SCCS unit conversions instead of local constants. The preset surface tension and pressure change by about 6e-7 (relative); the SCCS integration energies move by about 2e-8 eV. - Keep the solvation interface out of Potential's vtable. Potential overrode the three PotBase solvation hooks, so every unit test that stubs Potential had to define them: the SDFT Hamiltonian test, the elecstate base/PW tests and the hsolver PW mock. Potential now aggregates its components through non-virtual members (component_solvation_energy, add_component_solvation_force, component_solvent_electrostatic_potential), and those tests are back to upstream. - Rename the PW test cases: 220_PW_SCCS_PERIODIC_WATER -> scf_sccs, 221_PW_PCC0D_H3O -> scf_pcc_0d, 222_PW_PCC2D_H3O -> scf_pcc_2d. --- source/source_estate/elecstate_energy.cpp | 2 +- source/source_estate/module_pot/potential_new.cpp | 6 +++--- source/source_estate/module_pot/potential_new.h | 8 +++++--- source/source_estate/test/elecstate_base_test.cpp | 4 ---- .../source_estate/test/elecstate_energy_test.cpp | 4 +--- source/source_estate/test/elecstate_pw_test.cpp | 4 ---- source/source_estate/test/potential_new_test.cpp | 4 ++-- .../source_hamilt/module_sccs/sccs_parameters.cpp | 15 ++++----------- .../module_sccs/test/test_sccs_parameters.cpp | 4 ++-- .../test/hsolver_supplementary_mock.h | 15 --------------- source/source_io/module_ctrl/ctrl_output_fp.cpp | 2 +- source/source_lcao/force_stress_lcao.cpp | 2 +- source/source_pw/module_pwdft/force_pw.cpp | 2 +- .../module_stodft/test/test_sto_hamilt_pw.cpp | 3 --- tests/01_PW/CASES_CPU.txt | 6 +++--- .../01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/INPUT | 0 .../KPT | 0 .../{221_PW_PCC0D_H3O => scf_pcc_0d}/README.md | 0 tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/STRU | 0 .../{221_PW_PCC0D_H3O => scf_pcc_0d}/result.ref | 0 .../01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/INPUT | 0 tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_2d}/KPT | 0 .../{222_PW_PCC2D_H3O => scf_pcc_2d}/README.md | 0 tests/01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/STRU | 0 .../{222_PW_PCC2D_H3O => scf_pcc_2d}/result.ref | 0 .../INPUT | 0 tests/01_PW/{222_PW_PCC2D_H3O => scf_sccs}/KPT | 0 .../README | 0 .../{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/STRU | 0 .../result.ref | 0 .../threshold | 0 tests/02_NAO_Gamma/scf_sccs_water/README | 2 +- 32 files changed, 25 insertions(+), 58 deletions(-) rename tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/INPUT (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_pcc_0d}/KPT (100%) rename tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/README.md (100%) rename tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/STRU (100%) rename tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_0d}/result.ref (100%) rename tests/01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/INPUT (100%) rename tests/01_PW/{221_PW_PCC0D_H3O => scf_pcc_2d}/KPT (100%) rename tests/01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/README.md (100%) rename tests/01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/STRU (100%) rename tests/01_PW/{222_PW_PCC2D_H3O => scf_pcc_2d}/result.ref (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/INPUT (100%) rename tests/01_PW/{222_PW_PCC2D_H3O => scf_sccs}/KPT (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/README (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/STRU (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/result.ref (100%) rename tests/01_PW/{220_PW_SCCS_PERIODIC_WATER => scf_sccs}/threshold (100%) diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp index af77832329c..0cc903c6968 100644 --- a/source/source_estate/elecstate_energy.cpp +++ b/source/source_estate/elecstate_energy.cpp @@ -338,7 +338,7 @@ void ElecState::cal_energies(const int type, } else if (imp_sol == 2) { - this->pot->get_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav); + this->pot->component_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav); } else { diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index 51dbced99f1..2a71bba6cc1 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -403,7 +403,7 @@ void Potential::add_pcc_force(const UnitCell& cell, ModuleBase::matrix& force) c } } -void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +void Potential::add_component_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const { for (const PotBase* component : this->components) { @@ -411,7 +411,7 @@ void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& fo } } -void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +void Potential::component_solvation_energy(double& electrostatic, double& non_electrostatic) const { electrostatic = 0.0; non_electrostatic = 0.0; @@ -425,7 +425,7 @@ void Potential::get_solvation_energy(double& electrostatic, double& non_electros } } -const std::vector* Potential::solvent_electrostatic_potential() const +const std::vector* Potential::component_solvent_electrostatic_potential() const { for (const PotBase* component : this->components) { diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index 9bf4683a139..032fdd59cc7 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -83,9 +83,11 @@ class Potential : public PotBase PotBase* get_pot_type(const std::string& pot_type, const Input_para& input); - void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; - void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; - const std::vector* solvent_electrostatic_potential() const override; + // Solvation terms of the registered components; not PotBase overrides, so + // modules that stub Potential need not know about solvation. + void component_solvation_energy(double& electrostatic, double& non_electrostatic) const; + void add_component_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const; + const std::vector* component_solvent_electrostatic_potential() const; /** * @brief Inject the shared RT-TDDFT field state before potential setup. diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index 8bcae6b32de..317ef2805a2 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -37,10 +37,6 @@ int main(int argc, char** argv) // Mock functions for testing elecstate.cpp namespace elecstate { -void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} -void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } -const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } - void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/elecstate_energy_test.cpp b/source/source_estate/test/elecstate_energy_test.cpp index 0ccadc08dbf..87fc2fce5b8 100644 --- a/source/source_estate/test/elecstate_energy_test.cpp +++ b/source/source_estate/test/elecstate_energy_test.cpp @@ -13,12 +13,10 @@ bool XC_Functional::ked_flag = false; namespace elecstate { Potential::~Potential() {} -void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} void Potential::cal_v_eff(const Charge*, const UnitCell*, ModuleBase::matrix&) {} void Potential::cal_fixed_v(double*) {} -const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } -void Potential::get_solvation_energy(double& el, double& cav) const +void Potential::component_solvation_energy(double& el, double& cav) const { el = 0.8; cav = 0.9; diff --git a/source/source_estate/test/elecstate_pw_test.cpp b/source/source_estate/test/elecstate_pw_test.cpp index fc74f035816..1c47c8356b8 100644 --- a/source/source_estate/test/elecstate_pw_test.cpp +++ b/source/source_estate/test/elecstate_pw_test.cpp @@ -26,10 +26,6 @@ class TestParameters int XC_Functional::func_type = 1; namespace elecstate { -void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} -void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } -const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } - void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/potential_new_test.cpp b/source/source_estate/test/potential_new_test.cpp index 026b2222d70..87d4f42f3c6 100644 --- a/source/source_estate/test/potential_new_test.cpp +++ b/source/source_estate/test/potential_new_test.cpp @@ -431,13 +431,13 @@ TEST_F(PotentialNewTest, SolvationForcesUseRegisteredComponentsOnce) { for (int axis = 0; axis < 3; ++axis) { force(ia, axis) = 7.0; } } - potential->add_solvation_force(*ucell, force); + potential->add_component_solvation_force(*ucell, force); for (int ia = 0; ia < 2; ++ia) { for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } } potential->pot_register({"fixed"}); - potential->add_solvation_force(*ucell, force); + potential->add_component_solvation_force(*ucell, force); for (int ia = 0; ia < 2; ++ia) { for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } diff --git a/source/source_hamilt/module_sccs/sccs_parameters.cpp b/source/source_hamilt/module_sccs/sccs_parameters.cpp index 740ce525d1b..8bcdced285e 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.cpp +++ b/source/source_hamilt/module_sccs/sccs_parameters.cpp @@ -1,29 +1,22 @@ #include "sccs_parameters.h" +#include "source_base/constants.h" #include "source_base/tool_quit.h" - namespace ModuleSccs { -namespace -{ -// CODATA 2018 SI values. Environ 3.1.1 converts with 4.35974394e-18 J and -// 0.52917720859e-10 m, so the surface tension and pressure differ from its -// values by about 2e-7 (relative). ModuleBase's SI values differ from both. -const double hartree_joule = 4.3597447222071e-18; -const double bohr_metre = 5.29177210903e-11; -} double dyn_per_cm_to_hartree_per_bohr2(double value) { const double joule_per_square_metre = value * 1e-3; - return joule_per_square_metre * bohr_metre * bohr_metre / hartree_joule; + return joule_per_square_metre * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI / ModuleBase::HARTREE_SI; } double gpa_to_hartree_per_bohr3(double value) { const double pascal = value * 1e9; - return pascal * bohr_metre * bohr_metre * bohr_metre / hartree_joule; + return pascal * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI + / ModuleBase::HARTREE_SI; } Preset parse_preset(const std::string& name) diff --git a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp index 46a4c4bf975..60bf39c97f2 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -23,8 +23,8 @@ TEST(SccsParameters, PresetValues) ModuleSccs::Preset::WaterAnion}; const double minima[] = {1e-4, 2e-4, 2.4e-3}; const double maxima[] = {5e-3, 3.5e-3, 1.55e-2}; - const double tensions[] = {3.076640259577260e-5, 3.211524279308205e-6, 0.0}; - const double pressures[] = {-1.223615131827540e-5, 4.248663652178958e-6, 1.529518914784425e-5}; + const double tensions[] = {3.076638359192053e-5, 3.211522295607572e-6, 0.0}; + const double pressures[] = {-1.223613888352364e-5, 4.248659334556821e-6, 1.529517360440455e-5}; for (int index = 0; index < 3; ++index) { const ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(presets[index]); diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index e3c52fbd540..76da6b23eba 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -86,21 +86,6 @@ void Potential::cal_fixed_v(double* vl_pseudo) { } -void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const -{ - electrostatic = 0.0; - non_electrostatic = 0.0; -} - -void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const -{ -} - -const std::vector* Potential::solvent_electrostatic_potential() const -{ - return nullptr; -} - } // namespace elecstate #include "source_cell/klist.h" diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 39a9daa287d..14ca304ab65 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -161,7 +161,7 @@ void ctrl_output_fp(UnitCell& ucell, fn += geom_block + ".cube"; const std::vector* pcc_potential = pelec->pot->pcc_electron_potential(); - const std::vector* solvent_potential = pelec->pot->solvent_electrostatic_potential(); + const std::vector* solvent_potential = pelec->pot->component_solvent_electrostatic_potential(); ModuleIO::write_elecstat_pot( #ifdef __MPI pw_big->bz, diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 50d59d9a711..cb391474872 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -194,7 +194,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, if (cfg.imp_sol == 2) { parts.fsol.create(ucell.nat, 3); - pelec->pot->add_solvation_force(ucell, parts.fsol); + pelec->pot->add_component_solvation_force(ucell, parts.fsol); } } diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index d81997b89af..9b92aaff1a5 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -150,7 +150,7 @@ void Forces::cal_force(const int nspin, } else { - elec.pot->add_solvation_force(ucell, forcesol); + elec.pot->add_component_solvation_force(ucell, forcesol); } if (PARAM.inp.test_force) { diff --git a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp index 78dfc1a44aa..7023d90719b 100644 --- a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp +++ b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp @@ -6,9 +6,6 @@ #include elecstate::Potential::~Potential(){} -void elecstate::Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} -void elecstate::Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } -const std::vector* elecstate::Potential::solvent_electrostatic_potential() const { return nullptr; } void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase::matrix&){} void elecstate::Potential::cal_fixed_v(double*){} diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 19d024103a9..e8721997956 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -125,6 +125,6 @@ scf_deltaspin2 scf_deltaspin4 nscf_hse_same nscf_hse_path -221_PW_PCC0D_H3O -222_PW_PCC2D_H3O -220_PW_SCCS_PERIODIC_WATER +scf_pcc_0d +scf_pcc_2d +scf_sccs diff --git a/tests/01_PW/221_PW_PCC0D_H3O/INPUT b/tests/01_PW/scf_pcc_0d/INPUT similarity index 100% rename from tests/01_PW/221_PW_PCC0D_H3O/INPUT rename to tests/01_PW/scf_pcc_0d/INPUT diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT b/tests/01_PW/scf_pcc_0d/KPT similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT rename to tests/01_PW/scf_pcc_0d/KPT diff --git a/tests/01_PW/221_PW_PCC0D_H3O/README.md b/tests/01_PW/scf_pcc_0d/README.md similarity index 100% rename from tests/01_PW/221_PW_PCC0D_H3O/README.md rename to tests/01_PW/scf_pcc_0d/README.md diff --git a/tests/01_PW/221_PW_PCC0D_H3O/STRU b/tests/01_PW/scf_pcc_0d/STRU similarity index 100% rename from tests/01_PW/221_PW_PCC0D_H3O/STRU rename to tests/01_PW/scf_pcc_0d/STRU diff --git a/tests/01_PW/221_PW_PCC0D_H3O/result.ref b/tests/01_PW/scf_pcc_0d/result.ref similarity index 100% rename from tests/01_PW/221_PW_PCC0D_H3O/result.ref rename to tests/01_PW/scf_pcc_0d/result.ref diff --git a/tests/01_PW/222_PW_PCC2D_H3O/INPUT b/tests/01_PW/scf_pcc_2d/INPUT similarity index 100% rename from tests/01_PW/222_PW_PCC2D_H3O/INPUT rename to tests/01_PW/scf_pcc_2d/INPUT diff --git a/tests/01_PW/221_PW_PCC0D_H3O/KPT b/tests/01_PW/scf_pcc_2d/KPT similarity index 100% rename from tests/01_PW/221_PW_PCC0D_H3O/KPT rename to tests/01_PW/scf_pcc_2d/KPT diff --git a/tests/01_PW/222_PW_PCC2D_H3O/README.md b/tests/01_PW/scf_pcc_2d/README.md similarity index 100% rename from tests/01_PW/222_PW_PCC2D_H3O/README.md rename to tests/01_PW/scf_pcc_2d/README.md diff --git a/tests/01_PW/222_PW_PCC2D_H3O/STRU b/tests/01_PW/scf_pcc_2d/STRU similarity index 100% rename from tests/01_PW/222_PW_PCC2D_H3O/STRU rename to tests/01_PW/scf_pcc_2d/STRU diff --git a/tests/01_PW/222_PW_PCC2D_H3O/result.ref b/tests/01_PW/scf_pcc_2d/result.ref similarity index 100% rename from tests/01_PW/222_PW_PCC2D_H3O/result.ref rename to tests/01_PW/scf_pcc_2d/result.ref diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT b/tests/01_PW/scf_sccs/INPUT similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT rename to tests/01_PW/scf_sccs/INPUT diff --git a/tests/01_PW/222_PW_PCC2D_H3O/KPT b/tests/01_PW/scf_sccs/KPT similarity index 100% rename from tests/01_PW/222_PW_PCC2D_H3O/KPT rename to tests/01_PW/scf_sccs/KPT diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README b/tests/01_PW/scf_sccs/README similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README rename to tests/01_PW/scf_sccs/README diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU b/tests/01_PW/scf_sccs/STRU similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU rename to tests/01_PW/scf_sccs/STRU diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref b/tests/01_PW/scf_sccs/result.ref similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref rename to tests/01_PW/scf_sccs/result.ref diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold b/tests/01_PW/scf_sccs/threshold similarity index 100% rename from tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold rename to tests/01_PW/scf_sccs/threshold diff --git a/tests/02_NAO_Gamma/scf_sccs_water/README b/tests/02_NAO_Gamma/scf_sccs_water/README index a6df6c49e12..ba2c529b6ff 100644 --- a/tests/02_NAO_Gamma/scf_sccs_water/README +++ b/tests/02_NAO_Gamma/scf_sccs_water/README @@ -1,5 +1,5 @@ LCAO (gamma-only) neutral water with periodic SCCS (water-neutral preset), same geometry as -01_PW/220_PW_SCCS_PERIODIC_WATER. Covers the LCAO path end to end: SCCS potential in the LCAO +01_PW/scf_sccs. Covers the LCAO path end to end: SCCS potential in the LCAO Hamiltonian, local Pulay force terms and final force assembly of the SCCS ionic force. Checks total energy, forces and the electrostatic and cavity solvation terms. ecutwfc is 60 because LCAO requires ecutrho = 4 * ecutwfc; ecutrho 240 keeps the density grid of From fe8b8783b1535e71ccbda4311d35127d34d1455a Mon Sep 17 00:00:00 2001 From: LKFEIYI Date: Tue, 20 Oct 2026 14:37:06 +0800 Subject: [PATCH 5/5] Shorten SCCS function and file names to at most 15 characters Follows mohanchen's review: rename the new SCCS functions and files whose names were longer than 15 characters. Continuation lines are re-aligned. - PotBase hooks: get_solvation_energy -> sol_energy, add_solvation_force -> add_sol_force, solvent_electrostatic_potential -> sol_elec_pot. - Potential aggregators: component_solvation_energy -> sum_sol_energy, add_component_solvation_force -> sum_sol_force, component_solvent_electrostatic_potential -> find_sol_pot, make_sccs_potential -> make_sccs_pot. - ModuleSccs: dyn_per_cm_to_hartree_per_bohr2 -> dyn_cm_to_au, gpa_to_hartree_per_bohr3 -> gpa_to_au, make_sccs_config -> preset_config, gaussian_ionic_density -> gauss_ion_rho, gaussian_ionic_force -> gauss_ion_force, solve_sccs_response -> solve_response, evaluate_functional -> eval_functional. - elecstate: make_sccs_config_from_input -> sccs_cfg_input, check_sccs_charge -> check_charge. - ModuleIO: parse_solvation_model -> parse_imp_sol, validate_sccs_input -> check_sccs_inp. - File-local helpers: apply_preconditioner -> apply_precond, sccs_context_error -> context_error, sccs_model_error -> model_error, sccs_nonfinite_error -> nonfinite_error, sccs_solver_error -> solver_error. - Files: sccs_ionic_charge -> sccs_ion_charge, sccs_ionic_force -> sccs_ion_force, with their tests. --- source/Makefile.Objects | 4 +- source/source_estate/elecstate_energy.cpp | 2 +- source/source_estate/module_pot/pot_base.h | 6 +-- source/source_estate/module_pot/pot_sccs.cpp | 40 +++++++++---------- source/source_estate/module_pot/pot_sccs.h | 14 +++---- .../module_pot/potential_new.cpp | 12 +++--- .../source_estate/module_pot/potential_new.h | 8 ++-- .../module_pot/potential_types.cpp | 10 ++--- .../test/elecstate_energy_test.cpp | 2 +- .../source_estate/test/potential_new_test.cpp | 6 +-- source/source_estate/test/test_pot_sccs.cpp | 22 +++++----- .../source_hamilt/module_sccs/CMakeLists.txt | 4 +- .../module_sccs/sccs_functional.cpp | 12 +++--- .../module_sccs/sccs_functional.h | 12 +++--- ...s_ionic_charge.cpp => sccs_ion_charge.cpp} | 12 +++--- ...{sccs_ionic_charge.h => sccs_ion_charge.h} | 14 +++---- ...ccs_ionic_force.cpp => sccs_ion_force.cpp} | 14 +++---- .../module_sccs/sccs_ion_force.h | 30 ++++++++++++++ .../module_sccs/sccs_ionic_force.h | 30 -------------- .../module_sccs/sccs_parameters.cpp | 18 ++++----- .../module_sccs/sccs_parameters.h | 6 +-- .../module_sccs/sccs_response.cpp | 36 ++++++++--------- .../source_hamilt/module_sccs/sccs_response.h | 16 ++++---- .../module_sccs/test/CMakeLists.txt | 4 +- .../module_sccs/test/test_sccs_functional.cpp | 14 +++---- ...ic_charge.cpp => test_sccs_ion_charge.cpp} | 6 +-- ...onic_force.cpp => test_sccs_ion_force.cpp} | 24 +++++------ .../module_sccs/test/test_sccs_parameters.cpp | 4 +- .../module_sccs/test/test_sccs_response.cpp | 38 +++++++++--------- .../source_io/module_ctrl/ctrl_output_fp.cpp | 2 +- .../module_parameter/read_inp_model.cpp | 4 +- .../module_parameter/read_inp_sccs.cpp | 20 +++++----- .../module_parameter/read_inp_sccs.h | 4 +- .../test_serial/test_read_inp_sccs.cpp | 26 ++++++------ source/source_lcao/force_stress_lcao.cpp | 2 +- source/source_pw/module_pwdft/force_pw.cpp | 2 +- 36 files changed, 240 insertions(+), 240 deletions(-) rename source/source_hamilt/module_sccs/{sccs_ionic_charge.cpp => sccs_ion_charge.cpp} (74%) rename source/source_hamilt/module_sccs/{sccs_ionic_charge.h => sccs_ion_charge.h} (52%) rename source/source_hamilt/module_sccs/{sccs_ionic_force.cpp => sccs_ion_force.cpp} (80%) create mode 100644 source/source_hamilt/module_sccs/sccs_ion_force.h delete mode 100644 source/source_hamilt/module_sccs/sccs_ionic_force.h rename source/source_hamilt/module_sccs/test/{test_sccs_ionic_charge.cpp => test_sccs_ion_charge.cpp} (89%) rename source/source_hamilt/module_sccs/test/{test_sccs_ionic_force.cpp => test_sccs_ion_force.cpp} (79%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 5f31dd8913c..f393a5eb316 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -633,8 +633,8 @@ OBJS_SURCHEM=surchem.o\ OBJS_SCCS=sccs_cavity.o\ sccs_functional.o\ sccs_parameters.o\ - sccs_ionic_charge.o\ - sccs_ionic_force.o\ + sccs_ion_charge.o\ + sccs_ion_force.o\ sccs_pw_coulomb.o\ sccs_response.o diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp index 0cc903c6968..61109aa912d 100644 --- a/source/source_estate/elecstate_energy.cpp +++ b/source/source_estate/elecstate_energy.cpp @@ -338,7 +338,7 @@ void ElecState::cal_energies(const int type, } else if (imp_sol == 2) { - this->pot->component_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav); + this->pot->sum_sol_energy(this->f_en.esol_el, this->f_en.esol_cav); } else { diff --git a/source/source_estate/module_pot/pot_base.h b/source/source_estate/module_pot/pot_base.h index 658962f276e..63b096ca2b8 100644 --- a/source/source_estate/module_pot/pot_base.h +++ b/source/source_estate/module_pot/pot_base.h @@ -31,15 +31,15 @@ class PotBase virtual double get_energy() const { return 0.0; } // Solvation contributions in Ry, added once after the band-energy potential subtraction. - virtual void get_solvation_energy(double& electrostatic, double& non_electrostatic) const + virtual void sol_energy(double& electrostatic, double& non_electrostatic) const { electrostatic = 0.0; non_electrostatic = 0.0; } // Explicit ionic solvation force in Ry/Bohr, added to existing contributions. - virtual void add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} + virtual void add_sol_force(const UnitCell&, ModuleBase::matrix&) const {} // Electronic electrostatic correction only; excludes cavity/surface derivatives. - virtual const std::vector* solvent_electrostatic_potential() const { return nullptr; } + virtual const std::vector* sol_elec_pot() const { return nullptr; } bool fixed_mode = 0; bool dynamic_mode = 0; diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp index 7d5dd323b19..6624e17624a 100644 --- a/source/source_estate/module_pot/pot_sccs.cpp +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -5,8 +5,8 @@ #include "source_base/tool_quit.h" #include "source_cell/cell_tools.h" #include "source_hamilt/module_sccs/sccs_functional.h" -#include "source_hamilt/module_sccs/sccs_ionic_charge.h" -#include "source_hamilt/module_sccs/sccs_ionic_force.h" +#include "source_hamilt/module_sccs/sccs_ion_charge.h" +#include "source_hamilt/module_sccs/sccs_ion_force.h" #include "source_hamilt/module_sccs/sccs_response.h" #include "source_io/module_parameter/input_parameter.h" @@ -16,9 +16,9 @@ namespace elecstate { -void make_sccs_config_from_input(const Input_para& input, - ModuleSccs::SccsConfig& config, - ModuleSccs::PolarizationSolverParameters& solver) +void sccs_cfg_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver) { const ModuleSccs::Preset preset = ModuleSccs::parse_preset(input.sccs_preset); if (preset == ModuleSccs::Preset::Custom) @@ -27,12 +27,12 @@ void make_sccs_config_from_input(const Input_para& input, config.cavity.density_min = input.sccs_rho_min; config.cavity.density_max = input.sccs_rho_max; config.cavity.epsilon_bulk = input.sccs_epsilon; - config.surface_tension = ModuleSccs::dyn_per_cm_to_hartree_per_bohr2(input.sccs_gamma); - config.pressure = ModuleSccs::gpa_to_hartree_per_bohr3(input.sccs_pressure); + config.surface_tension = ModuleSccs::dyn_cm_to_au(input.sccs_gamma); + config.pressure = ModuleSccs::gpa_to_au(input.sccs_pressure); } else { - config = ModuleSccs::make_sccs_config(preset); + config = ModuleSccs::preset_config(preset); } config.surface_regularization = input.sccs_surface_eta; solver.max_iterations = input.sccs_maxiter; @@ -40,7 +40,7 @@ void make_sccs_config_from_input(const Input_para& input, solver.tolerance_max = input.sccs_tol_max; } -std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count) +std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count) { double ionic_charge = 0.0; for (int it = 0; it < cell.ntype; ++it) @@ -69,7 +69,7 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: const ModulePW::PW_Basis& basis = *this->rho_basis_; const std::vector atoms = unitcell::get_atom_data(cell->atoms, cell->ntype, cell->lat0); std::vector ions; - ModuleSccs::gaussian_ionic_density(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions); + ModuleSccs::gauss_ion_rho(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions); std::vector density(basis.nrxx, 0.0); std::vector solute_charge(basis.nrxx); double grid_charges[2] = {0.0, 0.0}; // electrons, ions @@ -99,10 +99,10 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: ModuleBase::WARNING("PotSccs::cal_v_eff", message.str()); } ModuleSccs::SccsResponse response; - ModuleSccs::solve_sccs_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, - cell->tpiba, response); + ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, + cell->tpiba, response); ModuleSccs::FunctionalResult functional; - ModuleSccs::evaluate_functional(solute_charge, response, config_, basis, cell->tpiba, functional); + ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, functional); electrostatic_rydberg_ = 2.0 * functional.reaction_energy; non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); electrostatic_potential_.resize(basis.nrxx); @@ -116,9 +116,9 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: ModuleBase::timer::end("PotSccs", "cal_v_eff"); } -void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +void PotSccs::add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const { - ModuleBase::timer::start("PotSccs", "add_solvation_force"); + ModuleBase::timer::start("PotSccs", "add_sol_force"); const ModulePW::PW_Basis& basis = *this->rho_basis_; const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); std::vector reaction(basis.nrxx); @@ -127,8 +127,8 @@ void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& forc reaction[ir] = -0.5 * electrostatic_potential_[ir]; } std::vector> ionic_force; - ModuleSccs::gaussian_ionic_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, - ionic_force); + ModuleSccs::gauss_ion_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, + ionic_force); for (int ia = 0; ia < cell.nat; ++ia) { for (int axis = 0; axis < 3; ++axis) @@ -136,7 +136,7 @@ void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& forc force(ia, axis) += 2.0 * ionic_force[ia][axis]; } } - ModuleBase::timer::end("PotSccs", "add_solvation_force"); + ModuleBase::timer::end("PotSccs", "add_sol_force"); } double PotSccs::get_energy() const @@ -144,13 +144,13 @@ double PotSccs::get_energy() const return electrostatic_rydberg_ + non_electrostatic_rydberg_; } -void PotSccs::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +void PotSccs::sol_energy(double& electrostatic, double& non_electrostatic) const { electrostatic = electrostatic_rydberg_; non_electrostatic = non_electrostatic_rydberg_; } -const std::vector* PotSccs::solvent_electrostatic_potential() const +const std::vector* PotSccs::sol_elec_pot() const { return &electrostatic_potential_; } diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h index 00373447a83..ab68bb85010 100644 --- a/source/source_estate/module_pot/pot_sccs.h +++ b/source/source_estate/module_pot/pot_sccs.h @@ -10,15 +10,15 @@ struct Input_para; namespace elecstate { // Map INPUT values, already validated by the INPUT reader, onto the solver. -void make_sccs_config_from_input(const Input_para& input, - ModuleSccs::SccsConfig& config, - ModuleSccs::PolarizationSolverParameters& solver); +void sccs_cfg_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver); // Returns a warning (empty if none) for a charged cell in a dielectric: the // periodic Poisson solver drops the G = 0 component of the net charge, so the // energy depends on the cell size. electron_count is the resolved nelec (the // INPUT value, or the valence count when it is not set). -std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count); +std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count); class PotSccs : public PotBase { @@ -32,9 +32,9 @@ class PotSccs : public PotBase double electron_count); void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override; double get_energy() const override; - void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; - void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; - const std::vector* solvent_electrostatic_potential() const override; + void add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const override; + void sol_energy(double& electrostatic, double& non_electrostatic) const override; + const std::vector* sol_elec_pot() const override; private: const ModuleSccs::SccsConfig config_; diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index 2a71bba6cc1..aecf60006ed 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -403,15 +403,15 @@ void Potential::add_pcc_force(const UnitCell& cell, ModuleBase::matrix& force) c } } -void Potential::add_component_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +void Potential::sum_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const { for (const PotBase* component : this->components) { - component->add_solvation_force(cell, force); + component->add_sol_force(cell, force); } } -void Potential::component_solvation_energy(double& electrostatic, double& non_electrostatic) const +void Potential::sum_sol_energy(double& electrostatic, double& non_electrostatic) const { electrostatic = 0.0; non_electrostatic = 0.0; @@ -419,17 +419,17 @@ void Potential::component_solvation_energy(double& electrostatic, double& non_el { double el = 0.0; double nonel = 0.0; - component->get_solvation_energy(el, nonel); + component->sol_energy(el, nonel); electrostatic += el; non_electrostatic += nonel; } } -const std::vector* Potential::component_solvent_electrostatic_potential() const +const std::vector* Potential::find_sol_pot() const { for (const PotBase* component : this->components) { - const std::vector* value = component->solvent_electrostatic_potential(); + const std::vector* value = component->sol_elec_pot(); if (value != nullptr) { return value; } } return nullptr; diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index 032fdd59cc7..fdfaa190f08 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -85,9 +85,9 @@ class Potential : public PotBase // Solvation terms of the registered components; not PotBase overrides, so // modules that stub Potential need not know about solvation. - void component_solvation_energy(double& electrostatic, double& non_electrostatic) const; - void add_component_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const; - const std::vector* component_solvent_electrostatic_potential() const; + void sum_sol_energy(double& electrostatic, double& non_electrostatic) const; + void sum_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const; + const std::vector* find_sol_pot() const; /** * @brief Inject the shared RT-TDDFT field state before potential setup. @@ -237,7 +237,7 @@ class Potential : public PotBase // the registered PCC component, or nullptr when PCC is off const PotPcc* pcc_component() const; - PotBase* make_sccs_potential(const Input_para& input); + PotBase* make_sccs_pot(const Input_para& input); std::vector v_eff_fixed; ModuleBase::matrix v_eff; diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 6890d083bdc..0cd6048e667 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -22,15 +22,15 @@ namespace elecstate { -PotBase* Potential::make_sccs_potential(const Input_para& input) +PotBase* Potential::make_sccs_pot(const Input_para& input) { ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; - make_sccs_config_from_input(input, config, solver); + sccs_cfg_input(input, config, solver); if (this->ucell_ != nullptr) { - const std::string warning = check_sccs_charge(config, *this->ucell_, input.nelec); - if (!warning.empty()) { ModuleBase::WARNING("Potential::make_sccs_potential", warning); } + const std::string warning = check_charge(config, *this->ucell_, input.nelec); + if (!warning.empty()) { ModuleBase::WARNING("Potential::make_sccs_pot", warning); } } return new PotSccs(this->rho_basis_, config, solver, input.nelec); } @@ -62,7 +62,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& } else if (pot_type == "sccs") { - return this->make_sccs_potential(input); + return this->make_sccs_pot(input); } else if (pot_type == "surchem") { diff --git a/source/source_estate/test/elecstate_energy_test.cpp b/source/source_estate/test/elecstate_energy_test.cpp index 87fc2fce5b8..4145163b383 100644 --- a/source/source_estate/test/elecstate_energy_test.cpp +++ b/source/source_estate/test/elecstate_energy_test.cpp @@ -16,7 +16,7 @@ Potential::~Potential() {} void Potential::cal_v_eff(const Charge*, const UnitCell*, ModuleBase::matrix&) {} void Potential::cal_fixed_v(double*) {} -void Potential::component_solvation_energy(double& el, double& cav) const +void Potential::sum_sol_energy(double& el, double& cav) const { el = 0.8; cav = 0.9; diff --git a/source/source_estate/test/potential_new_test.cpp b/source/source_estate/test/potential_new_test.cpp index 87d4f42f3c6..e85130cadd4 100644 --- a/source/source_estate/test/potential_new_test.cpp +++ b/source/source_estate/test/potential_new_test.cpp @@ -130,7 +130,7 @@ class MockPotComponent : public PotBase } } - void add_solvation_force(const UnitCell&, ModuleBase::matrix& force) const override + void add_sol_force(const UnitCell&, ModuleBase::matrix& force) const override { if (type_ != "solvation") { return; } for (int ia = 0; ia < force.nr; ++ia) @@ -431,13 +431,13 @@ TEST_F(PotentialNewTest, SolvationForcesUseRegisteredComponentsOnce) { for (int axis = 0; axis < 3; ++axis) { force(ia, axis) = 7.0; } } - potential->add_component_solvation_force(*ucell, force); + potential->sum_sol_force(*ucell, force); for (int ia = 0; ia < 2; ++ia) { for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } } potential->pot_register({"fixed"}); - potential->add_component_solvation_force(*ucell, force); + potential->sum_sol_force(*ucell, force); for (int ia = 0; ia < 2; ++ia) { for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp index 4bf0aa26977..9ed4b9c89f3 100644 --- a/source/source_estate/test/test_pot_sccs.cpp +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -52,7 +52,7 @@ class PotSccsTest : public SccsTest::PwTest ModuleBase::matrix potential(1, basis.nrxx); component.cal_v_eff(&charge, &cell, potential); ModuleBase::matrix force(2, 3); - component.add_solvation_force(cell, force); + component.add_sol_force(cell, force); const double step = 1e-4; for (int ia = 0; ia < 2; ++ia) { @@ -79,7 +79,7 @@ class PotSccsTest : public SccsTest::PwTest { for (int axis = 0; axis < 3; ++axis) { accumulated(ia, axis) = 7.0; } } - component.add_solvation_force(cell, accumulated); + component.add_sol_force(cell, accumulated); for (int ia = 0; ia < 2; ++ia) { for (int axis = 0; axis < 3; ++axis) @@ -113,7 +113,7 @@ TEST_F(PotSccsTest, IndependentInstancesAndSpinChannels) Input_para input; ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; - elecstate::make_sccs_config_from_input(input, config, solver); + elecstate::sccs_cfg_input(input, config, solver); config.cavity.epsilon_bulk = 5.0; elecstate::PotSccs first(&basis, config, solver, 1.0); config.cavity.epsilon_bulk = 1.0; @@ -125,7 +125,7 @@ TEST_F(PotSccsTest, IndependentInstancesAndSpinChannels) EXPECT_LT(first_energy, 0.0); double el = 0.0; double cav = 0.0; - first.get_solvation_energy(el, cav); + first.sol_energy(el, cav); EXPECT_DOUBLE_EQ(el, first_energy); EXPECT_DOUBLE_EQ(cav, 0.0); ModuleBase::matrix zero(1, basis.nrxx); @@ -153,7 +153,7 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) Input_para input; ModuleSccs::SccsConfig defaults; ModuleSccs::PolarizationSolverParameters default_solver; - elecstate::make_sccs_config_from_input(input, defaults, default_solver); + elecstate::sccs_cfg_input(input, defaults, default_solver); EXPECT_EQ(input.sccs_preset, "custom"); EXPECT_DOUBLE_EQ(defaults.cavity.epsilon_bulk, 1.0); EXPECT_DOUBLE_EQ(defaults.cavity.density_min, 1e-4); @@ -167,7 +167,7 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) input.sccs_surface_eta = 1e-6; ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; - elecstate::make_sccs_config_from_input(input, config, solver); + elecstate::sccs_cfg_input(input, config, solver); EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); EXPECT_DOUBLE_EQ(config.cavity.density_min, 1e-4); EXPECT_LT(config.pressure, 0.0); @@ -183,7 +183,7 @@ TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) input.sccs_tol_max = 1e-12; ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; - elecstate::make_sccs_config_from_input(input, config, solver); + elecstate::sccs_cfg_input(input, config, solver); check_solvation_force(config, solver); } @@ -203,13 +203,13 @@ TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) Input_para input; ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; - elecstate::make_sccs_config_from_input(input, config, solver); + elecstate::sccs_cfg_input(input, config, solver); config.cavity.epsilon_bulk = 5.0; - EXPECT_TRUE(elecstate::check_sccs_charge(config, cell, 1.0).empty()); - EXPECT_NE(elecstate::check_sccs_charge(config, cell, 0.5).find("G = 0"), std::string::npos); + EXPECT_TRUE(elecstate::check_charge(config, cell, 1.0).empty()); + EXPECT_NE(elecstate::check_charge(config, cell, 0.5).find("G = 0"), std::string::npos); ModuleSccs::SccsConfig vacuum_config = config; vacuum_config.cavity.epsilon_bulk = 1.0; - EXPECT_TRUE(elecstate::check_sccs_charge(vacuum_config, cell, 0.5).empty()); + EXPECT_TRUE(elecstate::check_charge(vacuum_config, cell, 0.5).empty()); const double density_value = 0.5 / basis.omega; std::vector density(basis.nrxx, density_value); diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt index f76b46ec716..f0186d67228 100644 --- a/source/source_hamilt/module_sccs/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -2,8 +2,8 @@ add_library(sccs OBJECT sccs_cavity.cpp sccs_parameters.cpp sccs_functional.cpp - sccs_ionic_charge.cpp - sccs_ionic_force.cpp + sccs_ion_charge.cpp + sccs_ion_force.cpp sccs_pw_coulomb.cpp sccs_response.cpp) diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp index 7547d7dae7b..cf7ce44ba59 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -13,12 +13,12 @@ namespace ModuleSccs { -void evaluate_functional(const std::vector& charge, - const SccsResponse& response, - const SccsConfig& config, - const ModulePW::PW_Basis& basis, - double tpiba, - FunctionalResult& result) +void eval_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result) { const std::size_t size = charge.size(); const double dv = basis.omega / basis.nxyz; diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h index b0a021b9d47..5d6f7d04b99 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.h +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -25,12 +25,12 @@ struct FunctionalResult }; // Collective over the PW pool. -void evaluate_functional(const std::vector& charge, - const SccsResponse& response, - const SccsConfig& config, - const ModulePW::PW_Basis& basis, - double tpiba, - FunctionalResult& result); +void eval_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/sccs_ion_charge.cpp similarity index 74% rename from source/source_hamilt/module_sccs/sccs_ionic_charge.cpp rename to source/source_hamilt/module_sccs/sccs_ion_charge.cpp index a1015081797..0876ce5e0c2 100644 --- a/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp +++ b/source/source_hamilt/module_sccs/sccs_ion_charge.cpp @@ -1,4 +1,4 @@ -#include "sccs_ionic_charge.h" +#include "sccs_ion_charge.h" #include "source_base/constants.h" #include "source_basis/module_pw/pw_basis.h" @@ -9,11 +9,11 @@ namespace ModuleSccs { -void gaussian_ionic_density(const std::vector& atoms, - const ModulePW::PW_Basis& basis, - double tpiba, - double spread, - std::vector& density) +void gauss_ion_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density) { std::vector> charge(basis.npw, 0.0); const double tpiba2 = tpiba * tpiba; diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.h b/source/source_hamilt/module_sccs/sccs_ion_charge.h similarity index 52% rename from source/source_hamilt/module_sccs/sccs_ionic_charge.h rename to source/source_hamilt/module_sccs/sccs_ion_charge.h index bf29d2276fd..4cbebd72af8 100644 --- a/source/source_hamilt/module_sccs/sccs_ionic_charge.h +++ b/source/source_hamilt/module_sccs/sccs_ion_charge.h @@ -1,5 +1,5 @@ -#ifndef SCCS_IONIC_CHARGE_H -#define SCCS_IONIC_CHARGE_H +#ifndef SCCS_ION_CHARGE_H +#define SCCS_ION_CHARGE_H #include @@ -19,11 +19,11 @@ const double gaussian_ion_spread = 0.5; // AtomData is extracted by source_cell; positions are Cartesian Bohr. // Collective over the PW pool. -void gaussian_ionic_density(const std::vector& atoms, - const ModulePW::PW_Basis& basis, - double tpiba, - double spread, - std::vector& density); +void gauss_ion_rho(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/sccs_ion_force.cpp similarity index 80% rename from source/source_hamilt/module_sccs/sccs_ionic_force.cpp rename to source/source_hamilt/module_sccs/sccs_ion_force.cpp index 0d6ee4f6726..d76daa99d20 100644 --- a/source/source_hamilt/module_sccs/sccs_ionic_force.cpp +++ b/source/source_hamilt/module_sccs/sccs_ion_force.cpp @@ -1,4 +1,4 @@ -#include "sccs_ionic_force.h" +#include "sccs_ion_force.h" #include "source_base/constants.h" #include "source_base/parallel_reduce.h" @@ -10,12 +10,12 @@ namespace ModuleSccs { -void gaussian_ionic_force(const std::vector& atoms, - const std::vector& reaction_potential, - const ModulePW::PW_Basis& basis, - double tpiba, - double spread, - std::vector>& forces) +void gauss_ion_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces) { std::vector> potential_g(basis.npw); basis.real2recip(reaction_potential.data(), potential_g.data()); diff --git a/source/source_hamilt/module_sccs/sccs_ion_force.h b/source/source_hamilt/module_sccs/sccs_ion_force.h new file mode 100644 index 00000000000..a405ac4e551 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ion_force.h @@ -0,0 +1,30 @@ +#ifndef SCCS_ION_FORCE_H +#define SCCS_ION_FORCE_H + +#include + +namespace ModuleBase +{ +template class Vector3; +} +namespace ModulePW +{ +class PW_Basis; +} +namespace unitcell +{ +struct AtomData; +} +namespace ModuleSccs +{ +// Explicit Gaussian-ion derivative of the reaction energy at fixed electron +// density. reaction_potential is in Ha; force is in Ha/Bohr, atom-major. +// Collective over the PW pool. +void gauss_ion_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces); +} // namespace ModuleSccs +#endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.h b/source/source_hamilt/module_sccs/sccs_ionic_force.h deleted file mode 100644 index 150781b03d8..00000000000 --- a/source/source_hamilt/module_sccs/sccs_ionic_force.h +++ /dev/null @@ -1,30 +0,0 @@ -#ifndef SCCS_IONIC_FORCE_H -#define SCCS_IONIC_FORCE_H - -#include - -namespace ModuleBase -{ -template class Vector3; -} -namespace ModulePW -{ -class PW_Basis; -} -namespace unitcell -{ -struct AtomData; -} -namespace ModuleSccs -{ -// Explicit Gaussian-ion derivative of the reaction energy at fixed electron -// density. reaction_potential is in Ha; force is in Ha/Bohr, atom-major. -// Collective over the PW pool. -void gaussian_ionic_force(const std::vector& atoms, - const std::vector& reaction_potential, - const ModulePW::PW_Basis& basis, - double tpiba, - double spread, - std::vector>& forces); -} // namespace ModuleSccs -#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.cpp b/source/source_hamilt/module_sccs/sccs_parameters.cpp index 8bcdced285e..30be0dca433 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.cpp +++ b/source/source_hamilt/module_sccs/sccs_parameters.cpp @@ -6,13 +6,13 @@ namespace ModuleSccs { -double dyn_per_cm_to_hartree_per_bohr2(double value) +double dyn_cm_to_au(double value) { const double joule_per_square_metre = value * 1e-3; return joule_per_square_metre * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI / ModuleBase::HARTREE_SI; } -double gpa_to_hartree_per_bohr3(double value) +double gpa_to_au(double value) { const double pascal = value * 1e9; return pascal * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI @@ -30,11 +30,11 @@ Preset parse_preset(const std::string& name) ModuleBase::WARNING_QUIT("ModuleSccs::parse_preset", message); } -SccsConfig make_sccs_config(Preset preset) +SccsConfig preset_config(Preset preset) { if (preset == Preset::Custom) { - ModuleBase::WARNING_QUIT("ModuleSccs::make_sccs_config", + ModuleBase::WARNING_QUIT("ModuleSccs::preset_config", "SCCS custom configurations must be supplied explicitly"); } SccsConfig config; @@ -50,22 +50,22 @@ SccsConfig make_sccs_config(Preset preset) } else { - config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(47.9); - config.pressure = gpa_to_hartree_per_bohr3(-0.36); + config.surface_tension = dyn_cm_to_au(47.9); + config.pressure = gpa_to_au(-0.36); } } else if (preset == Preset::WaterCation) { config.cavity.density_min = 2e-4; config.cavity.density_max = 3.5e-3; - config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(5.0); - config.pressure = gpa_to_hartree_per_bohr3(0.125); + config.surface_tension = dyn_cm_to_au(5.0); + config.pressure = gpa_to_au(0.125); } else { config.cavity.density_min = 2.4e-3; config.cavity.density_max = 1.55e-2; - config.pressure = gpa_to_hartree_per_bohr3(0.45); + config.pressure = gpa_to_au(0.45); } return config; } diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h index 707c5c0c386..0c01362923c 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.h +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -28,9 +28,9 @@ struct PolarizationSolverParameters // stops the run. Preset parse_preset(const std::string& name); // Published presets only; custom configurations are supplied by callers. -SccsConfig make_sccs_config(Preset preset); -double dyn_per_cm_to_hartree_per_bohr2(double value); -double gpa_to_hartree_per_bohr3(double value); +SccsConfig preset_config(Preset preset); +double dyn_cm_to_au(double value); +double gpa_to_au(double value); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index 43519550a6b..b9d52725f61 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -112,11 +112,11 @@ double grid_dot(const std::vector& left, } // P r = eps^-1/2 G eps^-1/2 r; only FFT scratch survives an application. -void apply_preconditioner(const std::vector& rhs, - const std::vector& invsqrt, - PeriodicCoulombOperator& coulomb, - std::vector& weighted, - std::vector& value) +void apply_precond(const std::vector& rhs, + const std::vector& invsqrt, + PeriodicCoulombOperator& coulomb, + std::vector& weighted, + std::vector& value) { const std::size_t size = rhs.size(); weighted.resize(size); @@ -165,7 +165,7 @@ void try_warm_start(const SqrtCgOperator& system, { guess_residual[i] = charge[i] - coefficient[i] * initial_potential[i]; } - apply_preconditioner(guess_residual, system.invsqrt, system.coulomb, state.weighted, z); + apply_precond(guess_residual, system.invsqrt, system.coulomb, state.weighted, z); for (std::size_t i = 0; i < size; ++i) { guess_residual[i] = coefficient[i] * (initial_potential[i] - z[i]); @@ -201,11 +201,11 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, double old_rz = 0.0; for (int iteration = 1; !converged(polarization, solver) && iteration <= solver.max_iterations; ++iteration) { - apply_preconditioner(residual, system.invsqrt, system.coulomb, state.weighted, z); + apply_precond(residual, system.invsqrt, system.coulomb, state.weighted, z); const double rz = grid_dot(residual, z, basis); if (!std::isfinite(rz) || std::abs(rz) < 1e-30) { - ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", + ModuleBase::WARNING_QUIT("ModuleSccs::solve_response", "SCCS sqrt-CG has a null or nonfinite preconditioned residual"); } const double beta = std::abs(old_rz) > 1e-30 ? rz / old_rz : 0.0; @@ -218,7 +218,7 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, const double curvature = grid_dot(direction, image, basis); if (!std::isfinite(curvature) || curvature == 0.0) { - ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", "SCCS sqrt-CG has invalid curvature"); + ModuleBase::WARNING_QUIT("ModuleSccs::solve_response", "SCCS sqrt-CG has invalid curvature"); } const double alpha = rz / curvature; for (std::size_t i = 0; i < size; ++i) @@ -235,7 +235,7 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, message << "SCCS sqrt-CG did not reach both residual tolerances within sccs_maxiter = " << solver.max_iterations << " iterations (RMS " << polarization.residual_rms << ", MAX " << polarization.residual_max << ")"; - ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", message.str()); + ModuleBase::WARNING_QUIT("ModuleSccs::solve_response", message.str()); } } @@ -277,14 +277,14 @@ void finish_response(const ModulePW::PW_Basis& basis, } } // namespace -void solve_sccs_response(const std::vector& density, - const std::vector& charge, - const CavityParameters& cavity, - const PolarizationSolverParameters& solver, - const std::vector& initial_potential, - const ModulePW::PW_Basis& basis, - double tpiba, - SccsResponse& result) +void solve_response(const std::vector& density, + const std::vector& charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result) { // Determine start mode collectively, including ranks with no real-space points. double initial_count = initial_potential.size(); diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index 0e9b5a89a98..71be6fcf9c3 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -43,14 +43,14 @@ struct SccsResponse // initial_potential is empty on all ranks for a cold start, otherwise a local grid. // A solver breakdown or missed tolerance stops the run with WARNING_QUIT. // No state or INPUT globals are read or retained. -void solve_sccs_response(const std::vector& cavity_density, - const std::vector& solute_charge, - const CavityParameters& cavity, - const PolarizationSolverParameters& solver, - const std::vector& initial_potential, - const ModulePW::PW_Basis& basis, - double tpiba, - SccsResponse& result); +void solve_response(const std::vector& cavity_density, + const std::vector& solute_charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt index 722c01ebd01..6801ed87606 100644 --- a/source/source_hamilt/module_sccs/test/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -6,14 +6,14 @@ AddTest(TARGET MODULE_HAMILT_sccs_parameters AddTest(TARGET MODULE_HAMILT_sccs_poisson LIBS parameter planewave device base - SOURCES test_sccs_pw_coulomb.cpp test_sccs_ionic_charge.cpp test_sccs_ionic_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + SOURCES 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/module_pot/pot_sccs.cpp ../../../source_cell/cell_tools.cpp ../../../source_cell/atom_spec.cpp ../../../source_cell/atom_pseudo.cpp ../../../source_cell/pseudo.cpp ../sccs_functional.cpp ../sccs_parameters.cpp - ../sccs_ionic_charge.cpp ../sccs_ionic_force.cpp ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp + ../sccs_ion_charge.cpp ../sccs_ion_force.cpp ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp ../../module_xc/xc_grad_utils.cpp) set_tests_properties(MODULE_HAMILT_sccs_poisson PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60) 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 c25cde778d4..35bb7982070 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -23,10 +23,10 @@ TEST_F(SccsFunctionalTest, UniformDielectricEnergyAndElectronDerivative) const std::vector density(basis.nrxx, 0.0); const std::vector cold; ModuleSccs::SccsResponse response; - ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold, - basis, tpiba, response); + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + basis, tpiba, response); ModuleSccs::FunctionalResult result; - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, result); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, result); const double kernel = ModuleBase::FOUR_PI / (tpiba * tpiba); const double expected = 0.25 * basis.omega * 1e-6 * kernel * (0.2 - 1.0); EXPECT_NEAR(result.reaction_energy, expected, 1e-13); @@ -41,7 +41,7 @@ TEST_F(SccsFunctionalTest, UniformDielectricEnergyAndElectronDerivative) TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) { - ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::Vacuum); config.surface_tension = 0.01; config.pressure = -0.001; config.surface_regularization = 0.02; @@ -57,7 +57,7 @@ TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) response.solute[ir] = 0.5 + 0.2 * mode[ir]; } ModuleSccs::FunctionalResult center; - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, center); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, center); // Avoid cancellation of grid-integrated energies at excessively small steps. const double step = 1e-5; double predicted = 0.0; @@ -70,12 +70,12 @@ TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) predicted *= basis.omega / basis.nxyz; ModuleSccs::FunctionalResult plus; ModuleSccs::FunctionalResult minus; - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, plus); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, plus); for (int ir = 0; ir < basis.nrxx; ++ir) { response.solute[ir] -= 2.0 * step * mode[ir]; } - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, minus); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, minus); const double fd = (plus.surface_energy + plus.volume_energy - minus.surface_energy - minus.volume_energy) / (2.0 * step); EXPECT_NEAR(fd, predicted, 1e-8); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp similarity index 89% rename from source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp rename to source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp index 27c2b795746..19fedc4cc11 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_ion_charge.cpp @@ -1,5 +1,5 @@ #include "sccs_test.h" -#include "../sccs_ionic_charge.h" +#include "../sccs_ion_charge.h" #include "source_base/parallel_reduce.h" #include "source_cell/cell_tools.h" @@ -14,7 +14,7 @@ TEST_F(SccsIonicChargeTest, ValenceNormalizationAndFourierPhase) atoms[1].valence_charge = 1.0; const double spread = ModuleSccs::gaussian_ion_spread; std::vector density; - ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, density); + ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, spread, density); double integral = 0.0; for (double value : density) { @@ -45,7 +45,7 @@ TEST_F(SccsIonicChargeTest, ValenceNormalizationAndFourierPhase) EXPECT_DOUBLE_EQ(selected, 1.0); atoms[0].position.x += length; std::vector translated; - ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, translated); + ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, spread, translated); for (int ir = 0; ir < basis.nrxx; ++ir) { EXPECT_NEAR(translated[ir], density[ir], 1e-13); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp similarity index 79% rename from source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp rename to source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp index 9b482dec854..851a1ce74f4 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp @@ -1,6 +1,6 @@ #include "sccs_test.h" -#include "../sccs_ionic_force.h" -#include "../sccs_ionic_charge.h" +#include "../sccs_ion_force.h" +#include "../sccs_ion_charge.h" #include "../sccs_response.h" #include "../sccs_parameters.h" #include "../sccs_functional.h" @@ -17,7 +17,7 @@ TEST_F(SccsIonicForceTest, AnalyticFourierForceAndGaugeInvariance) std::vector potential = cosine_mode(0); const double spread = ModuleSccs::gaussian_ion_spread; std::vector> forces; - ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + ModuleSccs::gauss_ion_force(atoms, potential, basis, tpiba, spread, forces); const double exponent = -0.25 * spread * spread * tpiba * tpiba; const double expected = 2.0 * tpiba * std::exp(exponent); EXPECT_NEAR(forces[0].x, expected, 1e-12); @@ -25,10 +25,10 @@ TEST_F(SccsIonicForceTest, AnalyticFourierForceAndGaugeInvariance) EXPECT_NEAR(forces[0].z, 0.0, 1e-12); atoms[0].position.x += length; for (double& value : potential) { value += 4.0; } - ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + ModuleSccs::gauss_ion_force(atoms, potential, basis, tpiba, spread, forces); EXPECT_NEAR(forces[0].x, expected, 1e-12); potential.assign(basis.nrxx, 3.0); - ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + ModuleSccs::gauss_ion_force(atoms, potential, basis, tpiba, spread, forces); EXPECT_NEAR(forces[0].x, 0.0, 1e-12); EXPECT_NEAR(forces[0].y, 0.0, 1e-12); EXPECT_NEAR(forces[0].z, 0.0, 1e-12); @@ -56,17 +56,17 @@ TEST_F(SccsIonicForceTest, ReactionEnergyFiniteDifferenceAtFixedNonuniformCavity auto evaluate = [&](const std::vector& displaced, ModuleSccs::FunctionalResult& functional) { std::vector charge; - ModuleSccs::gaussian_ionic_density(displaced, basis, tpiba, spread, charge); + ModuleSccs::gauss_ion_rho(displaced, basis, tpiba, spread, charge); for (int ir = 0; ir < basis.nrxx; ++ir) { charge[ir] -= density[ir]; } ModuleSccs::SccsResponse response; - ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, functional); + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, functional); }; ModuleSccs::FunctionalResult baseline; evaluate(atoms, baseline); std::vector> forces; - ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, - tpiba, spread, forces); + ModuleSccs::gauss_ion_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); const double step = 1e-4; for (std::size_t ia = 0; ia < atoms.size(); ++ia) { @@ -86,8 +86,8 @@ TEST_F(SccsIonicForceTest, ReactionEnergyFiniteDifferenceAtFixedNonuniformCavity } config.cavity.epsilon_bulk = 1.0; evaluate(atoms, baseline); - ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, - tpiba, spread, forces); + ModuleSccs::gauss_ion_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); for (const ModuleBase::Vector3& force : forces) { EXPECT_NEAR(force.x, 0.0, 1e-12); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp index 60bf39c97f2..5c6e5a65261 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -27,7 +27,7 @@ TEST(SccsParameters, PresetValues) const double pressures[] = {-1.223613888352364e-5, 4.248659334556821e-6, 1.529517360440455e-5}; for (int index = 0; index < 3; ++index) { - const ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(presets[index]); + const ModuleSccs::SccsConfig config = ModuleSccs::preset_config(presets[index]); EXPECT_DOUBLE_EQ(config.cavity.density_min, minima[index]); EXPECT_DOUBLE_EQ(config.cavity.density_max, maxima[index]); EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); @@ -35,7 +35,7 @@ TEST(SccsParameters, PresetValues) EXPECT_NEAR(config.surface_tension, tensions[index], 1e-18); EXPECT_NEAR(config.pressure, pressures[index], 1e-18); } - const ModuleSccs::SccsConfig vacuum = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + const ModuleSccs::SccsConfig vacuum = ModuleSccs::preset_config(ModuleSccs::Preset::Vacuum); EXPECT_DOUBLE_EQ(vacuum.cavity.epsilon_bulk, 1.0); EXPECT_DOUBLE_EQ(vacuum.surface_tension, 0.0); EXPECT_DOUBLE_EQ(vacuum.pressure, 0.0); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp index b3294e9d571..91613008c91 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -27,8 +27,8 @@ TEST_F(SccsResponseTest, VacuumAndUniformDielectricAnalyticLimits) { cavity.epsilon_bulk = epsilon; ModuleSccs::SccsResponse response; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, response); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); const double kernel = ModuleBase::FOUR_PI / (epsilon * tpiba * tpiba); EXPECT_EQ(response.polarization.iterations, 1); for (int ir = 0; ir < basis.nrxx; ++ir) @@ -62,8 +62,8 @@ TEST_F(SccsResponseTest, ManufacturedNonuniformDielectricSolution) solver.tolerance_max = 1e-12; const std::vector cold_start; ModuleSccs::SccsResponse response; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, response); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); double maximum_error = 0.0; for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -103,13 +103,13 @@ TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) { solver.max_iterations = 1; testing::internal::CaptureStdout(); - EXPECT_EXIT(ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, - response), + EXPECT_EXIT(ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, basis, tpiba, + response), ::testing::ExitedWithCode(1), ""); testing::internal::GetCapturedStdout(); } solver.max_iterations = 200; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, basis, tpiba, response); EXPECT_LE(response.polarization.residual_rms, solver.tolerance_rms); EXPECT_LE(response.polarization.residual_max, solver.tolerance_max); for (double& value : density) @@ -118,10 +118,10 @@ TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) } ModuleSccs::SccsResponse cold; ModuleSccs::SccsResponse warm; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, cold); - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, response.restart_potential, - basis, tpiba, warm); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, cold); + ModuleSccs::solve_response(density, charge, cavity, solver, response.restart_potential, + basis, tpiba, warm); EXPECT_TRUE(warm.polarization.warm_started); EXPECT_LT(warm.polarization.iterations, cold.polarization.iterations); for (int ir = 0; ir < basis.nrxx; ++ir) @@ -135,15 +135,15 @@ TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) ModuleSccs::SccsResponse loose; ModuleSccs::SccsResponse tight_maximum; ModuleSccs::SccsResponse tight_rms; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, loose); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, loose); solver.tolerance_max = 1e-11; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, tight_maximum); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_maximum); solver.tolerance_rms = 1e-11; solver.tolerance_max = 1.0; - ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, - basis, tpiba, tight_rms); + ModuleSccs::solve_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_rms); EXPECT_LE(tight_maximum.polarization.residual_max, 1e-11); EXPECT_LE(tight_rms.polarization.residual_rms, 1e-11); EXPECT_GT(tight_maximum.polarization.iterations, loose.polarization.iterations); @@ -216,9 +216,9 @@ TEST(SccsResponseRefinement, CavityDerivativeMatchesFiniteDifferenceOnAResolvedG { density[ir] = base_density + step * direction[ir]; } - ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); ModuleSccs::FunctionalResult functional; - ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, functional); + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, functional); return functional.reaction_energy; }; ModuleSccs::SccsResponse center; diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 14ca304ab65..73e0d679e5b 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -161,7 +161,7 @@ void ctrl_output_fp(UnitCell& ucell, fn += geom_block + ".cube"; const std::vector* pcc_potential = pelec->pot->pcc_electron_potential(); - const std::vector* solvent_potential = pelec->pot->component_solvent_electrostatic_potential(); + const std::vector* solvent_potential = pelec->pot->find_sol_pot(); ModuleIO::write_elecstat_pot( #ifdef __MPI pw_big->bz, diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index d217ca01b1e..cdfd4d2ebed 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -213,13 +213,13 @@ void ReadInput::item_model() item.default_value = "0"; item.read_value = [](const Input_Item& item, Parameter& para) { std::string error; - const bool valid = parse_solvation_model(item.str_values[0], para.input.imp_sol, error); + const bool valid = parse_imp_sol(item.str_values[0], para.input.imp_sol, error); if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } }; sync_int(input.imp_sol); item.check_value = [](const Input_Item&, const Parameter& para) { std::string error; - const bool valid = validate_sccs_input(para.input, error); + const bool valid = check_sccs_inp(para.input, error); if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } }; this->add_item(item); diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 0a201083554..fa19b401f30 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -10,7 +10,7 @@ namespace ModuleIO { -bool parse_solvation_model(const std::string& value, int& model, std::string& error) +bool parse_imp_sol(const std::string& value, int& model, std::string& error) { error.clear(); std::string normalized = value; @@ -36,7 +36,7 @@ bool parse_solvation_model(const std::string& value, int& model, std::string& er namespace { // Calculation settings that SCCS supports; empty when all are supported. -std::string sccs_context_error(const Input_para& input) +std::string context_error(const Input_para& input) { std::string error; if (input.assume_isolated != "none") { error = "SCCS currently requires assume_isolated none"; } @@ -56,7 +56,7 @@ std::string sccs_context_error(const Input_para& input) // 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 sccs_nonfinite_error(const Input_para& input) +std::string nonfinite_error(const Input_para& input) { const std::vector> values = { {"sccs_epsilon", input.sccs_epsilon}, @@ -75,7 +75,7 @@ std::string sccs_nonfinite_error(const Input_para& input) } // Cavity and surface model parameters, starting with the preset. -std::string sccs_model_error(const Input_para& input) +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()) @@ -91,7 +91,7 @@ std::string sccs_model_error(const Input_para& input) } // sqrt-CG solver controls. -std::string sccs_solver_error(const Input_para& input) +std::string solver_error(const Input_para& input) { std::string error; if (input.sccs_maxiter <= 0) { error = "sccs_maxiter must be positive"; } @@ -101,7 +101,7 @@ std::string sccs_solver_error(const Input_para& input) } } // namespace -bool validate_sccs_input(const Input_para& input, std::string& error) +bool check_sccs_inp(const Input_para& input, std::string& error) { error.clear(); if (input.imp_sol < 0 || input.imp_sol > 2) @@ -110,10 +110,10 @@ bool validate_sccs_input(const Input_para& input, std::string& error) return false; } if (input.imp_sol != 2) { return true; } - error = sccs_context_error(input); - if (error.empty()) { error = sccs_nonfinite_error(input); } - if (error.empty()) { error = sccs_model_error(input); } - if (error.empty()) { error = sccs_solver_error(input); } + 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(); } diff --git a/source/source_io/module_parameter/read_inp_sccs.h b/source/source_io/module_parameter/read_inp_sccs.h index 8736ed5768e..cb7c926c52c 100644 --- a/source/source_io/module_parameter/read_inp_sccs.h +++ b/source/source_io/module_parameter/read_inp_sccs.h @@ -4,7 +4,7 @@ struct Input_para; namespace ModuleIO { -bool parse_solvation_model(const std::string& value, int& model, std::string& error); -bool validate_sccs_input(const Input_para& input, std::string& error); +bool parse_imp_sol(const std::string& value, int& model, std::string& error); +bool check_sccs_inp(const Input_para& input, std::string& error); } #endif 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 d61292c6ce5..01ea319b8be 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -12,19 +12,19 @@ TEST(ReadInpSccs, LegacyBooleansAndModelSelection) int model = 0; for (const std::string& value : yes) { - ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + ASSERT_TRUE(ModuleIO::parse_imp_sol(value, model, error)); EXPECT_EQ(model, 1); } for (const std::string& value : no) { - ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + ASSERT_TRUE(ModuleIO::parse_imp_sol(value, model, error)); EXPECT_EQ(model, 0); } - ASSERT_TRUE(ModuleIO::parse_solvation_model("2", model, error)); + ASSERT_TRUE(ModuleIO::parse_imp_sol("2", model, error)); EXPECT_EQ(model, 2); - EXPECT_FALSE(ModuleIO::parse_solvation_model("3", model, error)); + EXPECT_FALSE(ModuleIO::parse_imp_sol("3", model, error)); EXPECT_EQ(model, 2); - EXPECT_FALSE(ModuleIO::parse_solvation_model("2.0", model, error)); + EXPECT_FALSE(ModuleIO::parse_imp_sol("2.0", model, error)); } TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) @@ -34,25 +34,25 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) input.device = "cpu"; input.basis_type = "pw"; std::string error; - ASSERT_TRUE(ModuleIO::validate_sccs_input(input, error)) << error; + ASSERT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; input.cal_force = true; - EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); input.cal_stress = true; - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); input.cal_stress = false; input.cal_force = false; input.assume_isolated = "pcc_0d"; - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); input.assume_isolated = "none"; input.sccs_rho_min = input.sccs_rho_max; - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); input.sccs_rho_min = 1e-4; input.sccs_preset = "unknown"; - EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); input.sccs_preset = "vacuum"; input.basis_type = "lcao"; - EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); input.imp_sol = 1; input.sccs_preset = "unused"; - EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); } diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index cb391474872..1e2759a0ea5 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -194,7 +194,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, if (cfg.imp_sol == 2) { parts.fsol.create(ucell.nat, 3); - pelec->pot->add_component_solvation_force(ucell, parts.fsol); + pelec->pot->sum_sol_force(ucell, parts.fsol); } } diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index 9b92aaff1a5..5bf444dd0fb 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -150,7 +150,7 @@ void Forces::cal_force(const int nspin, } else { - elec.pot->add_component_solvation_force(ucell, forcesol); + elec.pot->sum_sol_force(ucell, forcesol); } if (PARAM.inp.test_force) {