diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index f613fbe3e01..6317fe3f421 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -751,6 +751,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 29d0be32766..abe149acbdd 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:\ @@ -105,6 +106,7 @@ ${OBJS_PSI_INITIALIZER}\ ${OBJS_PW}\ ${OBJS_RELAXATION}\ ${OBJS_SURCHEM}\ +${OBJS_SCCS}\ ${OBJS_SYMMETRY}\ ${OBJS_XC}\ ${OBJS_IO}\ @@ -280,6 +282,7 @@ OBJS_ELECSTAT=elecstate.o\ pcc_0d.o\ pcc_2d.o\ pot_pcc.o\ + pot_sccs.o\ elecstate_op.o\ efield.o\ gatefield.o\ @@ -623,6 +626,14 @@ OBJS_SURCHEM=surchem.o\ minimize_cg.o\ sol_force.o\ +OBJS_SCCS=sccs_cavity.o\ + sccs_functional.o\ + sccs_parameters.o\ + sccs_ionic_charge.o\ + sccs_ionic_force.o\ + sccs_pw_coulomb.o\ + sccs_response.o + OBJS_SYMMETRY=symm_other.o\ symm_rot_spin.o\ symm_basic.o\ @@ -750,6 +761,7 @@ OBJS_IO=module_parameter/input_conv.o\ module_parameter/read_inp_bse.o\ module_parameter/read_inp_deepks.o\ module_parameter/read_inp_model.o\ + module_parameter/read_inp_sccs.o\ module_parameter/read_inp_postproc.o\ module_parameter/read_inp_exx_dftu.o\ module_parameter/read_inp_dfpt.o\ 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_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..00fe82fc18b --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -0,0 +1,321 @@ +#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); + // 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) + { + 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..4aea75c2428 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -0,0 +1,58 @@ +#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; + // -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 +}; + +// 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..46a4c4bf975 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -0,0 +1,42 @@ +#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..7a4b161b3d7 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -0,0 +1,255 @@ +#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; + +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); +} + +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_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 cca93932acc..101c4af261f 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -132,6 +132,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 05bc5e6f6af..39a9daa287d 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -161,6 +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(); ModuleIO::write_elecstat_pot( #ifdef __MPI pw_big->bz, @@ -174,6 +175,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 d86993b63d9..d852df7fdc9 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..1a5a1a302ef 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,34 @@ 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 = "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; + 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 +220,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 +257,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..0a201083554 --- /dev/null +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -0,0 +1,243 @@ +#include "read_inp_sccs.h" +#include "read_input.h" +#include "read_input_tool.h" +#include "source_base/tool_quit.h" + +#include +#include +#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; +} + +// 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) +{ + 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_nonfinite_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..c8411484db3 --- /dev/null +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -0,0 +1,94 @@ +#include "source_io/module_parameter/read_inp_sccs.h" +#include "source_io/module_parameter/input_parameter.h" +#include + +#include +#include +#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)); +} + +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 4bd110ec348..21f99dda951 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -190,6 +190,12 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, { parts.fpcc.create(ucell.nat, 3); pelec->pot->add_pcc_force(ucell, parts.fpcc); + // 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/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/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 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