diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 77add8f43b..c61c7c3184 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -288,6 +288,10 @@ OBJS_ELECSTAT=elecstate.o\ pcc_2d.o\ pot_pcc.o\ pot_sccs.o\ + potential_sccs.o\ + sccs_fixed_src.o\ + sccs_output.o\ + sccs_collect.o\ sccs_pcc_0d.o\ sccs_pcc_2d.o\ elecstate_op.o\ @@ -651,6 +655,7 @@ OBJS_SURCHEM=surchem.o\ OBJS_SCCS=sccs_cavity.o\ sccs_cav_deriv.o\ + sccs_diag.o\ sccs_functional.o\ sccs_lowpass.o\ sccs_parameters.o\ diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index d81f3bf435..dcb2a3c61f 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -314,9 +314,12 @@ bool ESolver_KS::defer_iter_output(const int istep) const void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver) { const int iter_limit = this->scf_iter_limit(istep); + const bool two_fermi = PARAM.globalv.two_fermi; + const bool model_changed + = module_charge::update_pot_model(iter, this->drho, this->kv, *this->pelec->pot, this->p_chgmix); // 1.1) print out band gap - if (!PARAM.globalv.two_fermi) + if (!two_fermi) { this->pelec->cal_bandgap(); } @@ -336,7 +339,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool & // 2.1) compute magnetization, only for spin==2 ucell.magnet.compute_mag(ucell.omega, this->chr.nrxx, this->chr.nxyz, this->chr.rho, - this->inp_->nspin, PARAM.globalv.two_fermi, this->inp_->nelec, + this->inp_->nspin, two_fermi, this->inp_->nelec, this->pelec->nelec_spin.data()); // 2.2) charge mixing @@ -354,6 +357,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool & ctx.iter_limit = iter_limit; ctx.hsolver_error = hsolver_error; ctx.scf_thr = this->scf_thr; + ctx.model_changed = model_changed; ctx.scf_ene_thr = this->scf_ene_thr; ctx.converged_u = converged_u; ctx.ks_run = PARAM.globalv.ks_run; @@ -435,7 +439,9 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool & // print energies elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, iter_limit, drho, - dkin, duration, *this->inp_, PARAM.globalv.two_fermi, diag_ethr, 0, true, this->ds_rms_); + dkin, duration, *this->inp_, two_fermi, diag_ethr, 0, true, this->ds_rms_); + + module_charge::print_corr_iter(this->kv, this->inp_->sccs_debug, this->drho, *this->pelec->pot); if (!this->defer_iter_output(istep)) diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index e1e481a5ac..3dd1a4425b 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -22,6 +22,10 @@ list(APPEND objects module_pot/pot_local.cpp module_pot/pot_pcc.cpp module_pot/pot_sccs.cpp + module_pot/potential_sccs.cpp + module_pot/sccs_fixed_src.cpp + module_pot/sccs_output.cpp + module_pot/sccs_collect.cpp module_pot/sccs_pcc_0d.cpp module_pot/sccs_pcc_2d.cpp module_pot/potential_new.cpp diff --git a/source/source_estate/module_charge/chg_mix.cpp b/source/source_estate/module_charge/chg_mix.cpp index 45d93180a4..cda51fc2e2 100644 --- a/source/source_estate/module_charge/chg_mix.cpp +++ b/source/source_estate/module_charge/chg_mix.cpp @@ -194,6 +194,10 @@ void Charge_Mixing::mix_reset() { this->mixing->reset(); this->rho_mdata.reset(); + // mix_dmr and mix_uom reuse the Pulay coefficients of the rho history, so + // the density-matrix and DFT+U occupation histories must restart with it. + this->dmr_mdata.reset(); + this->uom_mdata.reset(); // initailize tau_mdata if (cfg_.mixing_tau) { @@ -209,6 +213,12 @@ void Charge_Mixing::reset_scf_history(int iter_limit) this->mixing_restart_count = 0; } +void Charge_Mixing::restart_for_new_model(const int iter) +{ + this->mix_reset(); + this->mixing_restart_last = iter; +} + bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold) { diff --git a/source/source_estate/module_charge/chg_mix.h b/source/source_estate/module_charge/chg_mix.h index db9d7a3249..3b4358f644 100644 --- a/source/source_estate/module_charge/chg_mix.h +++ b/source/source_estate/module_charge/chg_mix.h @@ -77,7 +77,15 @@ class Charge_Mixing /** @brief Reset restart counters and residual history for an independently limited SCF trial. */ void reset_scf_history(int iter_limit); - + + /** + * @brief restart mixing after the potential model changed during the SCF (SCCS delayed start) + * + * The drho of earlier iterations belongs to the old model, so the window of + * if_scf_oscillate also restarts at iteration iter, as after a mixing restart. + */ + void restart_for_new_model(const int iter); + // extracting parameters normally these parameters will not be used outside charge mixing // while Exx is using them as well as some other places const std::string& get_mixing_mode() const {return cfg_.mixing_mode;} diff --git a/source/source_estate/module_charge/chg_routine.cpp b/source/source_estate/module_charge/chg_routine.cpp index 020af1c3cd..98ef710d0e 100644 --- a/source/source_estate/module_charge/chg_routine.cpp +++ b/source/source_estate/module_charge/chg_routine.cpp @@ -2,8 +2,70 @@ #include "source_estate/module_charge/chg_dmr.h" #include "source_base/parallel_comm.h" +#include "source_base/parallel_common.h" #include "source_pw/module_pwdft/dftu_base.h" // Plus_U_Base members used below #include "source_estate/update_pot.h" +#include "source_cell/klist.h" + +#include + +namespace +{ +// SCCS reports and the start notice come from the first rank of the first pool. +bool reports_here(const K_Vectors& kv) +{ + return kv.para_k.my_pool == 0 && kv.para_k.rank_in_pool == 0; +} + +// sccs_debug reports go to the running log of that rank. +std::ostream* corr_log(const K_Vectors& kv, const int sccs_debug) +{ + if (sccs_debug > 0 && reports_here(kv)) + { + return &GlobalV::ofs_running; + } + return nullptr; +} +} // namespace + +bool module_charge::update_pot_model(const int iter, + const double drho, + const K_Vectors& kv, + elecstate::Potential& pot, + Charge_Mixing* p_chgmix) +{ + // Only while a potential component waits to change its model: ordinary + // runs skip the broadcast. The state is synchronized, so every rank takes + // the same branch. + if (!pot.model_pending()) + { + return false; + } + // Band-parallel groups may differ slightly in drho until mixing broadcasts + // it; broadcast first so every rank changes the model in the same iteration. + double density_residual = drho; + Parallel_Common::bcast_double(density_residual); + const bool reports = reports_here(kv); + std::ostream* const notice = reports ? &std::cout : nullptr; + const bool model_changed = pot.update_model(iter, density_residual, notice); + if (model_changed) + { + p_chgmix->restart_for_new_model(iter); + } + return model_changed; +} + +void module_charge::print_corr_iter(const K_Vectors& kv, + const int sccs_debug, + const double drho, + const elecstate::Potential& pot) +{ + std::ostream* const log = corr_log(kv, sccs_debug); + if (log != nullptr) + { + pot.write_corr_iter(*log, sccs_debug, drho); + } +} void module_charge::chgmixing_ks(const int iter, UnitCell& ucell, @@ -42,7 +104,7 @@ void module_charge::chgmixing_ks(const int iter, // not ground state bool not_restart_step = !(iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0); - conv_esolver = (drho < scf_thr && not_restart_step && converged_u); + conv_esolver = (drho < scf_thr && not_restart_step && converged_u && !ctx.model_changed); // add energy threshold for SCF convergence if (scf_ene_thr > 0.0) diff --git a/source/source_estate/module_charge/chg_routine.h b/source/source_estate/module_charge/chg_routine.h index 2bfe5c83c3..f70d9c8420 100644 --- a/source/source_estate/module_charge/chg_routine.h +++ b/source/source_estate/module_charge/chg_routine.h @@ -10,6 +10,7 @@ // Plus_U_Base forward declaration, full definition in source_pw/module_pwdft/dftu_base.h class Plus_U_Base; +class K_Vectors; namespace module_charge { @@ -23,11 +24,28 @@ struct ScfMixingCtx double scf_ene_thr; ///< energy convergence threshold bool converged_u; ///< whether DFT+U has converged bool ks_run; ///< whether the current run is a KS calculation (PARAM.globalv.ks_run) + bool model_changed = false; ///< a potential component changed its model in this iteration (SCCS delayed start); the SCF is not converged double drho; ///< charge density deviation (in/out) bool oscillate_esolver; ///< whether esolver oscillates (out) bool conv_esolver; ///< whether esolver converged (out) }; +/// Lets a potential component change its model before mixing (SCCS delayed +/// start); true only in the iteration in which it does, which then restarts the +/// mixing. The first rank of the first pool announces the change on the screen. +bool update_pot_model(const int iter, + const double drho, + const K_Vectors& kv, + elecstate::Potential& pot, + Charge_Mixing* p_chgmix); + +/// sccs_debug report of the SCCS correction in this SCF iteration, in the +/// running log of the first rank of the first pool. +void print_corr_iter(const K_Vectors& kv, + const int sccs_debug, + const double drho, + const elecstate::Potential& pot); + void chgmixing_ks(const int iter, UnitCell& ucell, elecstate::ElecState* pelec, diff --git a/source/source_estate/module_charge/unittests/test_chg_mix.cpp b/source/source_estate/module_charge/unittests/test_chg_mix.cpp index 597c447ef9..2d8a34bd37 100644 --- a/source/source_estate/module_charge/unittests/test_chg_mix.cpp +++ b/source/source_estate/module_charge/unittests/test_chg_mix.cpp @@ -1,6 +1,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #include "../chg_mix.h" +#include "../chg_dmr.h" #include "../chg_drho.h" #include "../chg_drho_detail.h" #include "../chg_precond.h" @@ -193,6 +194,97 @@ class ChargeMixingTest : public ::testing::Test } }; +// A model change (SCCS delayed start) raises drho; the oscillation check must +// not compare it with the drho of the old model. +TEST_F(ChargeMixingTest, ModelChangeRestartsTheOscillationWindow) +{ + const int scf_nmax = 20; + cfg.scf_nmax = scf_nmax; + Charge_Mixing plain; + Charge_Mixing restarted; + plain.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba); + restarted.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba); + plain.init_mixing(); + restarted.init_mixing(); + const int scf_os_ndim = 3; + const double scf_os_thr = -0.05; + const int model_change = 6; + // Steady decrease, a jump when the new model starts, then a steady decrease. + const double drho[9] = {1e-2, 3e-3, 9e-4, 2.7e-4, 8.1e-5, 2.43e-5, 1e-2, 5e-3, 2.5e-3}; + for (int iter = 1; iter <= 9; ++iter) + { + if (iter == model_change) + { + restarted.restart_for_new_model(iter); + } + const double residual = drho[iter - 1]; + const bool plain_oscillates = plain.if_scf_oscillate(iter, residual, scf_os_ndim, scf_os_thr); + const bool restarted_oscillates = restarted.if_scf_oscillate(iter, residual, scf_os_ndim, scf_os_thr); + EXPECT_FALSE(restarted_oscillates) << "iteration " << iter; + if (iter == model_change + 1) + { + EXPECT_TRUE(plain_oscillates); + } + } + EXPECT_EQ(restarted.mixing_restart_last, model_change); +} + +// mix_uom reuses the coefficients of the rho history, so after a model change +// (SCCS delayed start) the first occupation update must be a single simple +// mixing step, not a combination of the old occupations. +TEST_F(ChargeMixingTest, ModelChangeRestartsUomHistory) +{ + cfg.mixing_mode = "pulay"; + Charge_Mixing charge_mixing; + charge_mixing.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba); + charge_mixing.init_mixing(); + const int uom_size = 4; + charge_mixing.allocate_mixing_uom(uom_size); + for (int step = 0; step < 2; ++step) + { + const double uom_value = 1.0 + step; + const double uom_save_value = 0.5 + step; + std::vector uom(uom_size, uom_value); + std::vector uom_save(uom_size, uom_save_value); + charge_mixing.mix_uom(uom, uom_save); + } + + const int model_change = 3; + charge_mixing.restart_for_new_model(model_change); + std::vector uom = {2.0, 3.0, 4.0, 5.0}; + std::vector uom_save = {1.0, 1.0, 1.0, 1.0}; + charge_mixing.mix_uom(uom, uom_save); + + const double beta = cfg.mixing_beta; + for (int i = 0; i < uom_size; ++i) + { + const double expected = 1.0 + beta * (2.0 + i - 1.0); + EXPECT_NEAR(uom[i], expected, 1e-12) << "element " << i; + } +} + +TEST_F(ChargeMixingTest, MixResetClearsDmrHistory) +{ + cfg.mixing_dmr = true; + cfg.scf_thr_type = 2; + Charge_Mixing charge_mixing; + charge_mixing.set_mixing(make_cfg(), &pw_basis, &pw_basis, ucell.omega, ucell.tpiba); + charge_mixing.init_mixing(); + + Base_Mixing::Mixing_Data& dmr_mdata = charge_mixing.get_dmr_mdata(); + module_charge::init_mixing_dmr(charge_mixing.get_mixing(), dmr_mdata, 4, + charge_mixing.get_mixing_config()); + const double dmr[4] = {1.0, 2.0, 3.0, 4.0}; + dmr_mdata.push(dmr); + ASSERT_EQ(dmr_mdata.ndim_history, 1); + + charge_mixing.mix_reset(); + + EXPECT_EQ(dmr_mdata.ndim_use, 0); + EXPECT_EQ(dmr_mdata.ndim_history, 0); + EXPECT_EQ(dmr_mdata.start, -1); +} + TEST_F(ChargeMixingTest, SetMixingTest) { #ifdef _OPENMP diff --git a/source/source_estate/module_pot/pot_base.h b/source/source_estate/module_pot/pot_base.h index 63b096ca2b..1967369c2e 100644 --- a/source/source_estate/module_pot/pot_base.h +++ b/source/source_estate/module_pot/pot_base.h @@ -9,6 +9,8 @@ namespace elecstate { +struct SolventGridField; + /** This class is the base class of Potential module 1. Main class Potential is derived from it. 2. components of potentials on real space grids can derived from it and will be registered into Potential. @@ -40,6 +42,8 @@ class PotBase virtual void add_sol_force(const UnitCell&, ModuleBase::matrix&) const {} // Electronic electrostatic correction only; excludes cavity/surface derivatives. virtual const std::vector* sol_elec_pot() const { return nullptr; } + // Implicit-solvent fields of the last correction for out_sol; none by default. + virtual void add_sol_fields(std::vector&) const {} bool fixed_mode = 0; bool dynamic_mode = 0; diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp index c021964caf..d00f1ea916 100644 --- a/source/source_estate/module_pot/pot_sccs.cpp +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -1,4 +1,5 @@ #include "pot_sccs.h" +#include "sol_grid_field.h" #include "sccs_pcc_0d.h" #include "sccs_pcc_2d.h" @@ -44,6 +45,9 @@ void sccs_cfg_input(const Input_para& input, { config = ModuleSccs::preset_config(preset); } + config.start_drho = input.sccs_start_drho; + config.start_nmax = input.sccs_start_nmax; + config.collect_report = input.sccs_debug >= 2; config.core_electrons = input.sccs_solvent_mode == "full"; config.core_spreads = input.sccs_corespread; config.surface_regularization = input.sccs_surface_eta; @@ -120,20 +124,58 @@ std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& c PotSccs::PotSccs(const ModulePW::PW_Basis* basis, const ModuleSccs::SccsConfig& config, const ModuleSccs::PolarizationSolverParameters& solver, - double electron_count) + double electron_count, + const SccsResume& resume) : config_(config), solver_(solver), electron_count_(electron_count) { this->rho_basis_ = basis; this->dynamic_mode = true; + sccs_active_ = config_.start_drho <= 0.0 || resume.active; + restart_potential_ = resume.restart_potential; +} + +SccsResume PotSccs::resume_state() const +{ + SccsResume resume; + resume.active = sccs_active_; + resume.restart_potential = restart_potential_; + return resume; +} + +bool PotSccs::try_activate(int electronic_iteration, double density_residual) +{ + if (sccs_active_) + { + return false; + } + if (density_residual > config_.start_drho && electronic_iteration < config_.start_nmax) + { + return false; + } + sccs_active_ = true; + restart_potential_.clear(); + return 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_; + if (!sccs_active_) + { + electrostatic_rydberg_ = 0.0; + non_electrostatic_rydberg_ = 0.0; + electrostatic_potential_.assign(basis.nrxx, 0.0); + cavity_potential_.clear(); + ModuleBase::timer::end("PotSccs", "cal_v_eff"); + return; + } const std::vector atoms = unitcell::get_atom_data(cell->atoms, cell->ntype, cell->lat0); - std::vector ions; - ModuleSccs::gauss_ion_rho(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions); + // The ion and core densities are rebuilt only when the atoms changed. + const bool reused_fixed_sources = fixed_sources_.update(*cell, basis, atoms, config_); + output_.reused_fixed_sources = reused_fixed_sources; + const std::vector& ions = fixed_sources_.ionic_density(); + const std::vector>& grid_positions = fixed_sources_.positions(); std::vector density(basis.nrxx, 0.0); std::vector solute_charge(basis.nrxx); double grid_charges[2] = {0.0, 0.0}; // electrons, ions @@ -166,11 +208,11 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: std::unique_ptr coulomb; if (config_.boundary == ModuleSccs::Boundary::Pcc0d) { - pcc_0d_operator(*cell, basis, atoms, coulomb); + pcc_0d_operator(*cell, basis, atoms, grid_positions, coulomb); } else if (config_.boundary == ModuleSccs::Boundary::Pcc2d) { - pcc_2d_operator(*cell, basis, atoms, config_.pcc_2d_axis, coulomb); + pcc_2d_operator(*cell, basis, atoms, config_.pcc_2d_axis, grid_positions, coulomb); } else { @@ -178,8 +220,7 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: } if (config_.core_electrons) { - std::vector core_density; - ModuleSccs::gauss_core_rho(atoms, basis, cell->tpiba, config_.core_spreads, core_density); + const std::vector& core_density = fixed_sources_.core_density(); for (int ir = 0; ir < basis.nrxx; ++ir) { density[ir] += core_density[ir]; @@ -188,8 +229,24 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: ModuleSccs::SccsResponse response; ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, cell->tpiba, *coulomb, response); + // Only the Coulomb-operator FFTs of the sqrt-CG solve, its main cost; the + // functional below and the cavity derivatives transform outside this count. + output_.transforms = coulomb->fft_counts(); ModuleSccs::FunctionalResult functional; ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, *coulomb, functional); + output_.valid = true; + output_.iterations = response.polarization.iterations; + output_.warm_started = response.polarization.warm_started; + output_.residual_rms = response.polarization.residual_rms; + output_.residual_max = response.polarization.residual_max; + output_.reaction_energy = functional.reaction_energy; + output_.volume = functional.volume; + output_.surface = functional.surface; + output_.far_field_charge = response.far_field_polarization_charge; + if (config_.collect_report && config_.boundary != ModuleSccs::Boundary::Periodic) + { + collect_output(*cell, basis, atoms, grid_positions, ions, solute_charge, config_, response, output_); + } electrostatic_rydberg_ = 2.0 * functional.reaction_energy; non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); electrostatic_potential_.resize(basis.nrxx); @@ -200,6 +257,8 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: for (int spin = 0; spin < charge->nspin; ++spin) { potential(spin, ir) += value; } } cavity_potential_ = std::move(functional.cavity_potential); + epsilon_ = std::move(response.epsilon); + solute_ = std::move(response.solute); restart_potential_ = std::move(response.restart_potential); ModuleBase::timer::end("PotSccs", "cal_v_eff"); } @@ -207,6 +266,11 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: void PotSccs::add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const { ModuleBase::timer::start("PotSccs", "add_sol_force"); + if (!sccs_active_) + { + ModuleBase::timer::end("PotSccs", "add_sol_force"); + return; + } 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); @@ -237,6 +301,22 @@ void PotSccs::add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) con ModuleBase::timer::end("PotSccs", "add_sol_force"); } +void PotSccs::write_iter(std::ostream& output, int level, double residual, double pcc_energy) const +{ + if (level == 0) + { + return; + } + if (!sccs_active_) + { + write_sccs_wait(output, residual, config_.start_drho, config_.start_nmax); + return; + } + const double energy = get_energy(); + const bool pcc = config_.boundary != ModuleSccs::Boundary::Periodic; + write_sccs_out(output, output_, level, energy, pcc, pcc_energy); +} + double PotSccs::get_energy() const { return electrostatic_rydberg_ + non_electrostatic_rydberg_; @@ -252,4 +332,22 @@ const std::vector* PotSccs::sol_elec_pot() const { return &electrostatic_potential_; } + +void PotSccs::add_sol_fields(std::vector& fields) const +{ + // Fields are written collectively: gate on output_.valid, set on every rank, + // not on epsilon_, which is empty on a rank without grid points. + if (!sccs_active_ || !output_.valid) + { + return; + } + SolventGridField dielectric; + dielectric.name = "eps"; + dielectric.values = epsilon_; + fields.push_back(dielectric); + SolventGridField cavity; + cavity.name = "cavity"; + cavity.values = solute_; + fields.push_back(cavity); +} } // namespace elecstate diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h index 77962f5c54..fd975c808f 100644 --- a/source/source_estate/module_pot/pot_sccs.h +++ b/source/source_estate/module_pot/pot_sccs.h @@ -2,8 +2,11 @@ #define POT_SCCS_H #include "pot_base.h" +#include "sccs_output.h" +#include "sccs_fixed_src.h" #include "source_hamilt/module_sccs/sccs_parameters.h" +#include #include struct Input_para; @@ -26,31 +29,59 @@ std::string check_sccs_cell(const ModuleSccs::SccsConfig& config, const UnitCell // INPUT value, or the valence count when it is not set). std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count); +// State carried from the component of the previous ionic step, which is +// rebuilt every ionic step: an activation already reached, and the last +// unshifted sqrt-CG solution as a warm start. The solver keeps the warm start +// only when it lowers the initial charge residual, so moved atoms are safe. +struct SccsResume +{ + bool active = false; + std::vector restart_potential; +}; + 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. + // resume carries the activation and warm start of the previous ionic step. PotSccs(const ModulePW::PW_Basis* basis, const ModuleSccs::SccsConfig& config, const ModuleSccs::PolarizationSolverParameters& solver, - double electron_count); + double electron_count, + const SccsResume& resume); void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override; double get_energy() const override; void add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const override; void sol_energy(double& electrostatic, double& non_electrostatic) const override; const std::vector* sol_elec_pot() const override; + // eps and the cavity s of the last active evaluation; none before the start. + void add_sol_fields(std::vector& fields) const override; + + // Delayed start: called once per electronic iteration while SCCS waits, + // with the density residual broadcast to every rank. Returns true only in + // the iteration that activates SCCS; activation is irreversible. + bool try_activate(int electronic_iteration, double density_residual); + bool is_active() const { return sccs_active_; } + SccsResume resume_state() const; + bool has_output() const { return output_.valid; } + void write_iter(std::ostream& output, int level, double residual, double pcc_energy) const; private: + bool sccs_active_ = false; const ModuleSccs::SccsConfig config_; const ModuleSccs::PolarizationSolverParameters solver_; const double electron_count_; + SccsOutput output_; + SccsFixedSources fixed_sources_; double electrostatic_rydberg_ = 0.0; double non_electrostatic_rydberg_ = 0.0; std::vector electrostatic_potential_; std::vector restart_potential_; std::vector cavity_potential_; + std::vector epsilon_; + std::vector solute_; }; } // namespace elecstate #endif diff --git a/source/source_estate/module_pot/pot_surchem.cpp b/source/source_estate/module_pot/pot_surchem.cpp index 113544aa61..8af3e4a30f 100644 --- a/source/source_estate/module_pot/pot_surchem.cpp +++ b/source/source_estate/module_pot/pot_surchem.cpp @@ -1,4 +1,5 @@ #include "pot_surchem.h" +#include "sol_grid_field.h" namespace elecstate { @@ -25,16 +26,18 @@ PotSurChem::~PotSurChem() void PotSurChem::cal_v_eff(const Charge* const chg, const UnitCell* const ucell, ModuleBase::matrix& v_eff) { + // nspin from the charge: v_eff has no rows on a rank without real-space grid points. + const int nspin = chg->nspin; if (!this->allocated) { - this->surchem_->allocate(this->rho_basis_->nrxx, v_eff.nr); + this->surchem_->allocate(this->rho_basis_->nrxx, nspin); this->allocated = true; } - ModuleBase::matrix v_sol_correction(v_eff.nr, this->rho_basis_->nrxx); + ModuleBase::matrix v_sol_correction(nspin, this->rho_basis_->nrxx); this->surchem_->v_correction(*ucell, *chg->pgrid, const_cast(this->rho_basis_), - v_eff.nr, + nspin, chg->rho, this->vlocal, this->structure_factors_, @@ -42,4 +45,35 @@ void PotSurChem::cal_v_eff(const Charge* const chg, const UnitCell* const ucell, v_eff += v_sol_correction; } +void PotSurChem::add_sol_fields(std::vector& fields) const +{ + // Fields are written collectively; last_epsilon() is empty on a rank without + // grid points, has_last_eps() is set on every rank. + if (!this->surchem_->has_last_eps()) + { + return; + } + const std::vector& epsilon = this->surchem_->last_epsilon(); + SolventGridField dielectric; + dielectric.name = "eps"; + dielectric.values = epsilon; + fields.push_back(dielectric); + // epsilon = 1 + (eb_k - 1) S with the solvent shape function S, so the + // cavity s = 1 - S follows from epsilon; without a dielectric it is undefined. + const double contrast = this->surchem_->bulk_epsilon() - 1.0; + if (contrast <= 0.0) + { + return; + } + SolventGridField cavity; + cavity.name = "cavity"; + cavity.values.resize(epsilon.size()); + for (std::size_t ir = 0; ir < epsilon.size(); ++ir) + { + const double shape = (epsilon[ir] - 1.0) / contrast; + cavity.values[ir] = 1.0 - shape; + } + fields.push_back(cavity); +} + } // namespace elecstate diff --git a/source/source_estate/module_pot/pot_surchem.h b/source/source_estate/module_pot/pot_surchem.h index a93685b9e2..fa0c62ad66 100644 --- a/source/source_estate/module_pot/pot_surchem.h +++ b/source/source_estate/module_pot/pot_surchem.h @@ -20,6 +20,8 @@ class PotSurChem : public PotBase // Passing an explicit output matrix makes the lifetime and allocation explicit and avoids hidden allocations. void cal_v_eff(const Charge* const chg, const UnitCell* const ucell, ModuleBase::matrix& v_eff) override; + // eps, and the cavity s = 1 - S of the shape function when eb_k > 1. + void add_sol_fields(std::vector& fields) const override; private: surchem* surchem_ = nullptr; diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index aecf60006e..48b8fb9782 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -1,5 +1,7 @@ #include "potential_new.h" #include "pot_pcc.h" +#include "pot_sccs.h" +#include "sol_grid_field.h" #include "pot_ml_exx.h" #include "source_base/global_function.h" @@ -69,6 +71,9 @@ void Potential::pot_register(const std::vector& components_list) void Potential::pot_register(const std::vector& components_list, const Input_para& input) { ModuleBase::TITLE("Potential", "pot_register"); + // The SCCS activation and warm start of an earlier ionic step, taken + // before the old components are deleted. + const SccsResume resume = this->get_resume(); // delete old components first. if (this->components.size() > 0) { @@ -85,7 +90,7 @@ void Potential::pot_register(const std::vector& components_list, co //--------------------------- for (auto comp: components_list) { - PotBase* tmp = this->get_pot_type(comp, input); + PotBase* tmp = this->get_pot_type(comp, input, resume); this->components.push_back(tmp); } @@ -435,4 +440,14 @@ const std::vector* Potential::find_sol_pot() const return nullptr; } +std::vector Potential::sol_out_fields() const +{ + std::vector fields; + for (const PotBase* component : this->components) + { + component->add_sol_fields(fields); + } + return fields; +} + } // namespace elecstate diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index e459864f11..3f75c0b275 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -7,6 +7,7 @@ #include "source_hamilt/module_surchem/surchem.h" #include "source_pw/module_pwdft/vsep_pw.h" +#include #include #include @@ -15,6 +16,9 @@ struct Input_para; namespace elecstate { class PotPcc; +class PotSccs; +struct SolventGridField; +struct SccsResume; class TDFieldManager; /** @@ -81,13 +85,25 @@ class Potential : public PotBase // 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, const Input_para& input); - + // resume is the SCCS state of the previous ionic step, used by "sccs" only. + PotBase* get_pot_type(const std::string& pot_type, const Input_para& input, const SccsResume& resume); + + /// True while a component waits to change its model during the SCF (the + /// SCCS delayed start); the state is the same on every rank. + bool model_pending() const; + /// Called before SCF mixing with the pool-consistent density residual. + /// True only in the iteration in which a component changes its model; + /// log, when not null, receives a one-line notice. + bool update_model(int electronic_iteration, double density_residual, std::ostream* log); + /// sccs_debug output of the SCCS correction. + void write_corr_iter(std::ostream& output, int level, double residual) const; // Solvation terms of the registered components; not PotBase overrides, so // modules that stub Potential need not know about solvation. void sum_sol_energy(double& electrostatic, double& non_electrostatic) const; void sum_sol_force(const UnitCell& cell, ModuleBase::matrix& force) const; const std::vector* find_sol_pot() const; + /// Implicit-solvent fields of the registered components, for out_sol. + std::vector sol_out_fields() const; /** * @brief Inject the shared RT-TDDFT field state before potential setup. @@ -237,7 +253,12 @@ class Potential : public PotBase // the registered PCC component, or nullptr when PCC is off const PotPcc* pcc_component() const; - PotBase* make_sccs_pot(const Input_para& input); + const PotSccs* sccs_component() const; + PotSccs* sccs_component(); + // The SCCS activation and warm start of an earlier ionic step are passed + // to the rebuilt component. + PotBase* make_sccs_pot(const Input_para& input, const SccsResume& resume); + SccsResume get_resume() const; std::vector v_eff_fixed; ModuleBase::matrix v_eff; diff --git a/source/source_estate/module_pot/potential_sccs.cpp b/source/source_estate/module_pot/potential_sccs.cpp new file mode 100644 index 0000000000..833da342c1 --- /dev/null +++ b/source/source_estate/module_pot/potential_sccs.cpp @@ -0,0 +1,70 @@ +#include "potential_new.h" +#include "pot_sccs.h" + +#include + +namespace elecstate +{ +bool Potential::update_model(int electronic_iteration, double density_residual, std::ostream* log) +{ + PotSccs* sccs = this->sccs_component(); + if (sccs == nullptr) + { + return false; + } + const bool activated = sccs->try_activate(electronic_iteration, density_residual); + if (activated && log != nullptr) + { + *log << " SCCS activated at electronic iteration " << electronic_iteration + << ", DRHO = " << density_residual << std::endl; + } + return activated; +} + +bool Potential::model_pending() const +{ + const PotSccs* sccs = this->sccs_component(); + return sccs != nullptr && !sccs->is_active(); +} + +void Potential::write_corr_iter(std::ostream& output, int level, double residual) const +{ + const PotSccs* sccs = this->sccs_component(); + if (level == 0 || sccs == nullptr) + { + return; + } + // Active SCCS reports the coupled PCC energy too; a deferred one its start condition. + const double pcc_energy = this->pcc_energy_rydberg(); + sccs->write_iter(output, level, residual, pcc_energy); +} + +const PotSccs* Potential::sccs_component() const +{ + for (const PotBase* component : this->components) + { + const PotSccs* sccs = dynamic_cast(component); + if (sccs != nullptr) + { + return sccs; + } + } + return nullptr; +} + +PotSccs* Potential::sccs_component() +{ + const Potential& self = *this; + return const_cast(self.sccs_component()); +} + +SccsResume Potential::get_resume() const +{ + const PotSccs* sccs = this->sccs_component(); + if (sccs == nullptr) + { + return SccsResume(); + } + return sccs->resume_state(); +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 488bec38ea..b93c59345c 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -22,7 +22,7 @@ namespace elecstate { -PotBase* Potential::make_sccs_pot(const Input_para& input) +PotBase* Potential::make_sccs_pot(const Input_para& input, const SccsResume& resume) { ModuleSccs::SccsConfig config; ModuleSccs::PolarizationSolverParameters solver; @@ -40,10 +40,10 @@ PotBase* Potential::make_sccs_pot(const Input_para& input) ModuleBase::WARNING("Potential::make_sccs_pot", charge_warning); } } - return new PotSccs(this->rho_basis_, config, solver, input.nelec); + return new PotSccs(this->rho_basis_, config, solver, input.nelec, resume); } -PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& input) +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& input, const SccsResume& resume) { ModuleBase::TITLE("Potential", "get_pot_type"); if (pot_type == "local") @@ -70,7 +70,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& } else if (pot_type == "sccs") { - return this->make_sccs_pot(input); + return this->make_sccs_pot(input, resume); } else if (pot_type == "surchem") { diff --git a/source/source_estate/module_pot/sccs_collect.cpp b/source/source_estate/module_pot/sccs_collect.cpp new file mode 100644 index 0000000000..6e2fe6d56e --- /dev/null +++ b/source/source_estate/module_pot/sccs_collect.cpp @@ -0,0 +1,85 @@ +#include "sccs_output.h" +#include "sccs_pcc_0d.h" +#include "sccs_pcc_2d.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/cell_geometry.h" +#include "source_cell/cell_tools.h" +#include "source_cell/unitcell.h" +#include "source_estate/pcc_0d.h" +#include "source_estate/pcc_2d.h" +#include "source_hamilt/module_sccs/sccs_diag.h" +#include "source_hamilt/module_sccs/sccs_parameters.h" +#include "source_hamilt/module_sccs/sccs_moments.h" +#include "source_hamilt/module_sccs/sccs_response.h" + +namespace elecstate +{ +void collect_output(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + const std::vector>& grid_positions, + const std::vector& ions, + const std::vector& charge, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::SccsResponse& response, + SccsOutput& result) +{ + std::vector polarization; + ModuleSccs::cont_pol_charge(charge, response, basis, cell.tpiba, polarization); + PccOutput candidate; + candidate.slab = config.boundary == ModuleSccs::Boundary::Pcc2d; + candidate.axis = config.pcc_2d_axis; + unitcell::OrthogonalCell geometry; + unitcell::SlabCell slab; + Pcc0dParameters parameters; + Pcc2dParameters slab_parameters; + std::vector> positions; + std::vector ionic_charges; + for (const unitcell::AtomData& atom : atoms) + { + positions.push_back(atom.position); + ionic_charges.push_back(atom.valence_charge); + } + std::vector> grid = grid_positions; + // The same geometry as the SCCS operator of this evaluation. + if (candidate.slab) + { + pcc_2d_geometry(cell, atoms, candidate.axis, slab, slab_parameters); + candidate.coordinate = slab.origin; + candidate.normal = slab.normal; + pcc_2d_project(slab, grid); + pcc_2d_project(slab, positions); + } + else + { + pcc_0d_geometry(cell, atoms, geometry, parameters); + candidate.origin = geometry.origin; + pcc_0d_relative(geometry, grid); + pcc_0d_relative(geometry, positions); + } + std::vector electrons(charge.size()); + for (std::size_t i = 0; i < charge.size(); ++i) + { + electrons[i] = charge[i] - ions[i]; + } + candidate.moments[0] = ModuleSccs::pool_moments(charge, grid, basis); + candidate.moments[1] = ModuleSccs::pool_moments(electrons, grid, basis); + candidate.moments[2] = ModuleSccs::pool_moments(polarization, grid, basis); + const int atom_count = atoms.size(); + const double* ionic_charge_data = ionic_charges.data(); + const ModuleBase::Vector3* ionic_position_data = positions.data(); + const ChargeMoments point_ions = charge_moments(ionic_charge_data, ionic_position_data, atom_count, 1.0); + candidate.moments[1] = add_charge_moments(candidate.moments[1], point_ions); + candidate.moments[3] = add_charge_moments(candidate.moments[0], candidate.moments[2]); + if (candidate.slab) + { + candidate.smooth_energy = pcc_2d_energy(candidate.moments[0], slab_parameters); + } + else + { + candidate.smooth_energy = pcc_0d_energy(candidate.moments[0], parameters); + } + result.pcc = candidate; + result.expected_charge = -(1.0 - 1.0 / config.cavity.epsilon_bulk) * candidate.moments[0].charge; +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_fixed_src.cpp b/source/source_estate/module_pot/sccs_fixed_src.cpp new file mode 100644 index 0000000000..567b295896 --- /dev/null +++ b/source/source_estate/module_pot/sccs_fixed_src.cpp @@ -0,0 +1,72 @@ +#include "sccs_fixed_src.h" +#include "source_base/timer.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_tools.h" +#include "source_cell/unitcell.h" +#include "source_hamilt/module_sccs/sccs_ion_charge.h" +#include "source_hamilt/module_sccs/sccs_parameters.h" + +namespace +{ +bool same_atoms(const std::vector& left, const std::vector& right) +{ + if (left.size() != right.size()) + { + return false; + } + for (std::size_t i = 0; i < left.size(); ++i) + { + const ModuleBase::Vector3 displacement = left[i].position - right[i].position; + if (displacement.norm2() != 0.0 || left[i].valence_charge != right[i].valence_charge + || left[i].atomic_number != right[i].atomic_number) + { + return false; + } + } + return true; +} +} +namespace elecstate +{ +struct SccsFixedSources::Data +{ + bool valid = false; + std::vector atoms; + std::vector ionic_density; + std::vector core_density; + std::vector> positions; +}; +SccsFixedSources::SccsFixedSources() +{ + data_.reset(new Data); +} +SccsFixedSources::~SccsFixedSources() = default; + +bool SccsFixedSources::update(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + const ModuleSccs::SccsConfig& config) +{ + ModuleBase::timer::start("SccsFixedSources", "update"); + // Atoms are replicated, so every pool rank takes the same decision. + const bool reuse = data_->valid && same_atoms(data_->atoms, atoms); + if (!reuse) + { + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, data_->positions); + ModuleSccs::gauss_ion_rho(atoms, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, + data_->ionic_density); + if (config.core_electrons) + { + ModuleSccs::gauss_core_rho(atoms, basis, cell.tpiba, config.core_spreads, data_->core_density); + } + data_->atoms = atoms; + data_->valid = true; + } + ModuleBase::timer::end("SccsFixedSources", "update"); + return reuse; +} +const std::vector& SccsFixedSources::ionic_density() const { return data_->ionic_density; } +const std::vector& SccsFixedSources::core_density() const { return data_->core_density; } +const std::vector>& SccsFixedSources::positions() const { return data_->positions; } +} // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_fixed_src.h b/source/source_estate/module_pot/sccs_fixed_src.h new file mode 100644 index 0000000000..c4ef8631a7 --- /dev/null +++ b/source/source_estate/module_pot/sccs_fixed_src.h @@ -0,0 +1,36 @@ +#ifndef SCCS_FIXED_SRC_H +#define SCCS_FIXED_SRC_H + +#include +#include + +class UnitCell; +namespace ModuleBase { template class Vector3; } +namespace ModulePW { class PW_Basis; } +namespace unitcell { struct AtomData; } +namespace ModuleSccs { struct SccsConfig; } +namespace elecstate +{ +// Owns ion/core Gaussian densities and Cartesian grid coordinates. The cell, +// grid and SCCS configuration are fixed for one PotSccs, which is rebuilt every +// ionic step; update rebuilds only for different atoms and returns true when +// the previous sources were reused. +class SccsFixedSources +{ +public: + SccsFixedSources(); + ~SccsFixedSources(); + bool update(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + const ModuleSccs::SccsConfig& config); + // Empty until the first update; views live until the next update. + const std::vector& ionic_density() const; + const std::vector& core_density() const; + const std::vector>& positions() const; +private: + struct Data; + std::unique_ptr data_; +}; +} // namespace elecstate +#endif diff --git a/source/source_estate/module_pot/sccs_output.cpp b/source/source_estate/module_pot/sccs_output.cpp new file mode 100644 index 0000000000..ccd314c9fc --- /dev/null +++ b/source/source_estate/module_pot/sccs_output.cpp @@ -0,0 +1,163 @@ +#include "sccs_output.h" +#include "source_base/constants.h" +#include +#include +#include +#include +#include + +namespace +{ +// The report is a grid: every cell (tag, key or value) has one width and is +// left-aligned, so every cell starts on a column edge. Energies keep 9 +// significant digits, other reals 8. +const int cell_width = 16; +const int energy_digits = 9; +const int value_digits = 8; + +std::string real_text(double value, int digits) +{ + std::ostringstream text; + const int mantissa_digits = digits - 1; + text << std::scientific << std::setprecision(mantissa_digits) << value; + return text.str(); +} +std::string int_text(int value) +{ + return std::to_string(value); +} +// One line of cells, without trailing blanks. A cell wider than the column +// keeps one blank after it. +void write_row(std::ostream& output, const std::vector& cells) +{ + const std::size_t width = cell_width; + std::string line = " "; + for (std::size_t i = 0; i + 1 < cells.size(); ++i) + { + line += cells[i]; + const std::size_t padding = cells[i].size() < width ? width - cells[i].size() : 1; + line.append(padding, ' '); + } + line += cells.back(); + output << line << '\n'; +} +void write_origin(std::ostream& output, const elecstate::PccOutput& result) +{ + if (result.slab) + { + const std::vector cells = {"PCC2D_ORIGIN", "AXIS", int_text(result.axis), + "COORDINATE/Bohr", real_text(result.coordinate, value_digits)}; + write_row(output, cells); + return; + } + const std::vector cells = {"PCC0D_ORIGIN", + "X/Bohr", real_text(result.origin.x, value_digits), + "Y/Bohr", real_text(result.origin.y, value_digits), + "Z/Bohr", real_text(result.origin.z, value_digits)}; + write_row(output, cells); +} +// Smooth and point solute, polarization and screened charge: one header row, +// then one row per kind. A slab reports the dipole and second moment along the normal. +void write_moments(std::ostream& output, const elecstate::PccOutput& result) +{ + const std::string tag = result.slab ? "PCC2D_MOMENTS" : "PCC0D_MOMENTS"; + std::vector header = {tag, "KIND", "Q/e"}; + if (result.slab) + { + header.push_back("PN/eBohr"); + header.push_back("QNN/eBohr2"); + } + else + { + header.push_back("PX/eBohr"); + header.push_back("PY/eBohr"); + header.push_back("PZ/eBohr"); + header.push_back("QRR/eBohr2"); + } + write_row(output, header); + const char* labels[] = {"SMOOTH", "POINT", "POLARIZATION", "SCREENED"}; + for (int i = 0; i < 4; ++i) + { + const elecstate::ChargeMoments& moments = result.moments[i]; + std::vector row = {tag, labels[i], real_text(moments.charge, value_digits)}; + if (result.slab) + { + const double dipole = moments.dipole * result.normal; + row.push_back(real_text(dipole, value_digits)); + } + else + { + row.push_back(real_text(moments.dipole.x, value_digits)); + row.push_back(real_text(moments.dipole.y, value_digits)); + row.push_back(real_text(moments.dipole.z, value_digits)); + } + row.push_back(real_text(moments.second_moment, value_digits)); + write_row(output, row); + } +} +} +namespace elecstate +{ +void write_sccs_wait(std::ostream& output, double residual, double start_drho, int start_nmax) +{ + const std::vector cells = {"SCCS_DEFERRED", "DRHO", real_text(residual, value_digits), + "START_DRHO", real_text(start_drho, value_digits), + "START_NMAX", int_text(start_nmax)}; + write_row(output, cells); +} + +void write_sccs_out(std::ostream& output, const SccsOutput& result, int level, + double solvation_energy, bool pcc, double pcc_energy) +{ + if (level == 0 || !result.valid) + { + return; + } + std::vector iteration = {"SCCS_ITER", "CG_STEPS", int_text(result.iterations), + "E_SOL/Ry", real_text(solvation_energy, energy_digits)}; + if (pcc) + { + iteration.push_back("E_PCC/Ry"); + iteration.push_back(real_text(pcc_energy, energy_digits)); + } + write_row(output, iteration); + if (level < 2) + { + return; + } + const std::vector residual = {"SCCS_RESIDUAL", "RMS", real_text(result.residual_rms, value_digits), + "MAX", real_text(result.residual_max, value_digits), + "WARM_START", int_text(result.warm_started)}; + write_row(output, residual); + if (pcc) + { + const std::vector gauss = {"SCCS_GAUSS", + "Q_FAR_FIELD/e", real_text(result.far_field_charge, value_digits), + "Q_POL/e", real_text(result.pcc.moments[2].charge, value_digits), + "Q_EXPECTED/e", real_text(result.expected_charge, value_digits)}; + write_row(output, gauss); + } + const std::vector cavity = {"SCCS_CAVITY", "VOLUME/Bohr3", real_text(result.volume, value_digits), + "SURFACE/Bohr2", real_text(result.surface, value_digits)}; + write_row(output, cavity); + const std::vector transforms = {"SCCS_FFT", + "R2G_CALLS", int_text(result.transforms.forward_calls), + "G2R_CALLS", int_text(result.transforms.inverse_calls), + "CACHED_SOURCES", int_text(result.reused_fixed_sources)}; + write_row(output, transforms); + if (!pcc) + { + return; + } + const PccOutput& data = result.pcc; + write_origin(output, data); + write_moments(output, data); + // PCC energies are in Ry; ModuleBase::e2 converts Hartree to Rydberg. + const double point_energy = pcc_energy / ModuleBase::e2; + const std::string tag = data.slab ? "PCC2D_ENERGY" : "PCC0D_ENERGY"; + const std::vector energies = {tag, "REACTION/Ha", real_text(result.reaction_energy, energy_digits), + "PCC_SMOOTH/Ha", real_text(data.smooth_energy, energy_digits), + "PCC_POINT/Ha", real_text(point_energy, energy_digits)}; + write_row(output, energies); +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_output.h b/source/source_estate/module_pot/sccs_output.h new file mode 100644 index 0000000000..15253c2c3b --- /dev/null +++ b/source/source_estate/module_pot/sccs_output.h @@ -0,0 +1,63 @@ +#ifndef SCCS_OUTPUT_H +#define SCCS_OUTPUT_H + +#include "source_estate/charge_moments.h" +#include "source_hamilt/module_sccs/sccs_coulomb.h" +#include +#include + +class UnitCell; +namespace ModulePW { class PW_Basis; } +namespace unitcell { struct AtomData; } +namespace ModuleSccs +{ +struct SccsConfig; +struct SccsResponse; +struct FunctionalResult; +} +namespace elecstate +{ +// Scalar reporting snapshot only: no retained grid buffers or model controls. +struct PccOutput +{ + bool slab = false; + int axis = 2; + double coordinate = 0.0; + ModuleBase::Vector3 origin; + ModuleBase::Vector3 normal; + ChargeMoments moments[4]; // smooth, point, polarization, screened + double smooth_energy = 0.0; // Hartree +}; +struct SccsOutput +{ + bool valid = false; + bool reused_fixed_sources = false; + int iterations = 0; + bool warm_started = false; + double residual_rms = 0.0; + double residual_max = 0.0; + double reaction_energy = 0.0; + double volume = 0.0; + double surface = 0.0; + double far_field_charge = 0.0; + double expected_charge = 0.0; + ModuleSccs::CoulombTransformCounts transforms; + PccOutput pcc; +}; +// All ranks in the PW pool call together. +void collect_output(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + const std::vector>& grid_positions, + const std::vector& ions, + const std::vector& charge, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::SccsResponse& response, + SccsOutput& result); +// sccs_debug lines of one SCF iteration, tagged and column-aligned. +void write_sccs_out(std::ostream& output, const SccsOutput& result, int level, + double solvation_energy, bool pcc, double pcc_energy); +// The line of an SCF iteration while a delayed SCCS waits. +void write_sccs_wait(std::ostream& output, double residual, double start_drho, int start_nmax); +} // namespace elecstate +#endif diff --git a/source/source_estate/module_pot/sccs_pcc_0d.cpp b/source/source_estate/module_pot/sccs_pcc_0d.cpp index f2d1de6f3e..8014517936 100644 --- a/source/source_estate/module_pot/sccs_pcc_0d.cpp +++ b/source/source_estate/module_pot/sccs_pcc_0d.cpp @@ -2,7 +2,6 @@ #include "source_base/tool_quit.h" #include "source_basis/module_pw/pw_basis.h" -#include "source_basis/module_pw/pw_grid_geometry.h" #include "source_cell/cell_geometry.h" #include "source_cell/cell_tools.h" #include "source_cell/unitcell.h" @@ -45,13 +44,13 @@ void pcc_0d_relative(const unitcell::OrthogonalCell& geometry, void pcc_0d_operator(const UnitCell& cell, const ModulePW::PW_Basis& basis, const std::vector& atoms, + const std::vector>& grid_positions, std::unique_ptr& coulomb) { unitcell::OrthogonalCell geometry; Pcc0dParameters parameters; pcc_0d_geometry(cell, atoms, geometry, parameters); - std::vector> relative_positions; - ModulePW::grid_positions(basis, cell.latvec, cell.lat0, relative_positions); + std::vector> relative_positions = grid_positions; pcc_0d_relative(geometry, relative_positions); coulomb.reset(new ModuleSccs::Pcc0dCoulombOperator(basis, cell.tpiba, relative_positions, parameters)); } diff --git a/source/source_estate/module_pot/sccs_pcc_0d.h b/source/source_estate/module_pot/sccs_pcc_0d.h index 97a6b8cd01..6f2ab530ad 100644 --- a/source/source_estate/module_pot/sccs_pcc_0d.h +++ b/source/source_estate/module_pot/sccs_pcc_0d.h @@ -38,10 +38,12 @@ void pcc_0d_relative(const unitcell::OrthogonalCell& geometry, std::vector>& positions); // Assemble an SCCS grid-charge operator using the existing cell geometry and -// mass-weighted ionic origin. An unsupported cell stops the run. +// mass-weighted ionic origin, on the Cartesian grid coordinates grid_positions +// (e.g. from SccsFixedSources). An unsupported cell stops the run. void pcc_0d_operator(const UnitCell& cell, const ModulePW::PW_Basis& basis, const std::vector& atoms, + const std::vector>& grid_positions, std::unique_ptr& coulomb); } // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_pcc_2d.cpp b/source/source_estate/module_pot/sccs_pcc_2d.cpp index 44a403fefa..b0a38b135a 100644 --- a/source/source_estate/module_pot/sccs_pcc_2d.cpp +++ b/source/source_estate/module_pot/sccs_pcc_2d.cpp @@ -2,7 +2,6 @@ #include "source_base/tool_quit.h" #include "source_basis/module_pw/pw_basis.h" -#include "source_basis/module_pw/pw_grid_geometry.h" #include "source_cell/cell_geometry.h" #include "source_cell/cell_tools.h" #include "source_cell/unitcell.h" @@ -51,13 +50,13 @@ void pcc_2d_operator(const UnitCell& cell, const ModulePW::PW_Basis& basis, const std::vector& atoms, int open_axis, + const std::vector>& grid_positions, std::unique_ptr& coulomb) { unitcell::SlabCell geometry; Pcc2dParameters parameters; pcc_2d_geometry(cell, atoms, open_axis, geometry, parameters); - std::vector> projected_positions; - ModulePW::grid_positions(basis, cell.latvec, cell.lat0, projected_positions); + std::vector> projected_positions = grid_positions; pcc_2d_project(geometry, projected_positions); coulomb.reset(new ModuleSccs::Pcc2dCoulombOperator(basis, cell.tpiba, projected_positions, geometry.normal, parameters)); diff --git a/source/source_estate/module_pot/sccs_pcc_2d.h b/source/source_estate/module_pot/sccs_pcc_2d.h index dfb2a4f6f8..dc88e62d59 100644 --- a/source/source_estate/module_pot/sccs_pcc_2d.h +++ b/source/source_estate/module_pot/sccs_pcc_2d.h @@ -40,11 +40,13 @@ void pcc_2d_project(const unitcell::SlabCell& geometry, std::vector>& positions); // Assemble an SCCS grid-charge operator using the existing cell geometry and -// mass-weighted ionic origin. An unsupported cell stops the run. +// mass-weighted ionic origin, on the Cartesian grid coordinates grid_positions +// (e.g. from SccsFixedSources). An unsupported cell stops the run. void pcc_2d_operator(const UnitCell& cell, const ModulePW::PW_Basis& basis, const std::vector& atoms, int open_axis, + const std::vector>& grid_positions, std::unique_ptr& coulomb); } // namespace elecstate diff --git a/source/source_estate/module_pot/sol_grid_field.h b/source/source_estate/module_pot/sol_grid_field.h new file mode 100644 index 0000000000..0665824c38 --- /dev/null +++ b/source/source_estate/module_pot/sol_grid_field.h @@ -0,0 +1,18 @@ +#ifndef SOL_GRID_FIELD_H +#define SOL_GRID_FIELD_H + +#include +#include + +namespace elecstate +{ +/// A real-space field of an implicit-solvent model on this rank's density +/// grid, written by out_sol as sol_{name}.cube. +struct SolventGridField +{ + std::string name; + std::vector values; +}; +} // namespace elecstate + +#endif diff --git a/source/source_estate/unittests/CMakeLists.txt b/source/source_estate/unittests/CMakeLists.txt index 9cb03fbcc3..3a565c9abc 100644 --- a/source/source_estate/unittests/CMakeLists.txt +++ b/source/source_estate/unittests/CMakeLists.txt @@ -28,6 +28,7 @@ AddTest( if (ENABLE_MPI) set(pcc_test_sources + ../module_pot/sccs_output.cpp ../module_pot/pot_pcc.cpp ../pcc_0d.cpp ../pcc_2d.cpp ../charge_moments.cpp ../../source_cell/cell_geometry.cpp ../../source_cell/cell_tools.cpp ../../source_cell/atom_spec.cpp ../../source_cell/atom_pseudo.cpp ../../source_cell/pseudo.cpp diff --git a/source/source_estate/unittests/test_pot_sccs.cpp b/source/source_estate/unittests/test_pot_sccs.cpp index 2e3ab0ac00..1184ecab3d 100644 --- a/source/source_estate/unittests/test_pot_sccs.cpp +++ b/source/source_estate/unittests/test_pot_sccs.cpp @@ -9,6 +9,11 @@ #include #include +namespace +{ +const elecstate::SccsResume no_resume = elecstate::SccsResume(); +} + // Public storage fixtures: no INPUT initialization or privately owned atom maps. UnitCell::UnitCell() {} UnitCell::~UnitCell() {} @@ -48,7 +53,7 @@ class PotSccsTest : public SccsTest::PwTest Charge charge; charge.nspin = 1; charge.rho = channels; - elecstate::PotSccs component(&basis, config, solver, 2.0); + elecstate::PotSccs component(&basis, config, solver, 2.0, no_resume); ModuleBase::matrix potential(1, basis.nrxx); component.cal_v_eff(&charge, &cell, potential); ModuleBase::matrix force(2, 3); @@ -115,9 +120,9 @@ TEST_F(PotSccsTest, IndependentInstancesAndSpinChannels) ModuleSccs::PolarizationSolverParameters solver; elecstate::sccs_cfg_input(input, config, solver); config.cavity.epsilon_bulk = 5.0; - elecstate::PotSccs first(&basis, config, solver, 1.0); + elecstate::PotSccs first(&basis, config, solver, 1.0, no_resume); config.cavity.epsilon_bulk = 1.0; - elecstate::PotSccs vacuum(&basis, config, solver, 1.0); + elecstate::PotSccs vacuum(&basis, config, solver, 1.0, no_resume); EXPECT_DOUBLE_EQ(first.get_energy(), 0.0); ModuleBase::matrix potential(1, basis.nrxx); first.cal_v_eff(&charge, &cell, potential); @@ -209,6 +214,18 @@ TEST_F(PotSccsTest, InputMapsCavityModelKeys) EXPECT_EQ(config.core_spreads, input.sccs_corespread); // The cavity model keys leave the preset dielectric unchanged. EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); + EXPECT_DOUBLE_EQ(config.start_drho, 0.0); + EXPECT_FALSE(config.collect_report); + input.sccs_start_drho = 1e-3; + input.sccs_start_nmax = 3; + input.sccs_debug = 2; + elecstate::sccs_cfg_input(input, config, solver); + EXPECT_DOUBLE_EQ(config.start_drho, 1e-3); + EXPECT_EQ(config.start_nmax, 3); + EXPECT_TRUE(config.collect_report); + input.sccs_debug = 1; + elecstate::sccs_cfg_input(input, config, solver); + EXPECT_FALSE(config.collect_report); } TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) @@ -240,6 +257,79 @@ TEST_F(PotSccsTest, FullCavityForceMatchesTheEnergyDerivative) check_solvation_force(config, solver); } +TEST_F(PotSccsTest, DelayedActivationHasOneIrreversibleTransitionAndResumes) +{ + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + config.start_drho = 1e-3; + config.start_nmax = 3; + elecstate::PotSccs threshold(&basis, config, solver, 1.0, no_resume); + EXPECT_FALSE(threshold.is_active()); + EXPECT_FALSE(threshold.try_activate(1, 0.1)); + EXPECT_TRUE(threshold.try_activate(2, 1e-3)); + EXPECT_TRUE(threshold.is_active()); + EXPECT_FALSE(threshold.try_activate(3, 0.2)); + EXPECT_FALSE(threshold.try_activate(1, 0.2)); + elecstate::PotSccs forced(&basis, config, solver, 1.0, no_resume); + EXPECT_FALSE(forced.try_activate(2, 0.1)); + EXPECT_TRUE(forced.try_activate(3, 0.1)); + config.start_drho = 0.0; + elecstate::PotSccs immediate(&basis, config, solver, 1.0, no_resume); + EXPECT_TRUE(immediate.is_active()); + EXPECT_FALSE(immediate.try_activate(1, 0.0)); + // A later ionic step resumes an activation already reached. + config.start_drho = 1e-3; + elecstate::PotSccs first_ionic_step(&basis, config, solver, 1.0, no_resume); + EXPECT_FALSE(first_ionic_step.is_active()); + elecstate::SccsResume resumed; + resumed.active = true; + elecstate::PotSccs next_ionic_step(&basis, config, solver, 1.0, resumed); + EXPECT_TRUE(next_ionic_step.is_active()); + EXPECT_FALSE(next_ionic_step.try_activate(1, 0.1)); + EXPECT_DOUBLE_EQ(next_ionic_step.get_energy(), 0.0); + const std::vector* potential = next_ionic_step.sol_elec_pot(); + EXPECT_TRUE(potential->empty()); +} + +TEST_F(PotSccsTest, WaitingSccsContributesNoPotentialEnergyOrForce) +{ + 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 value = 1.0 / basis.omega; + std::vector density(basis.nrxx, value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + config.start_drho = 1e-3; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::PotSccs delayed(&basis, config, solver, 1.0, no_resume); + ModuleBase::matrix potential(1, basis.nrxx); + delayed.cal_v_eff(&charge, &cell, potential); + EXPECT_DOUBLE_EQ(delayed.get_energy(), 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_DOUBLE_EQ(potential(0, ir), 0.0); + } + ModuleBase::matrix force(1, 3); + force(0, 0) = 42.0; + delayed.add_sol_force(cell, force); + EXPECT_DOUBLE_EQ(force(0, 0), 42.0); + EXPECT_TRUE(delayed.try_activate(2, 1e-4)); + delayed.cal_v_eff(&charge, &cell, potential); + EXPECT_LT(delayed.get_energy(), 0.0); +} + TEST_F(PotSccsTest, FullCavityStructureChecksWarnOnceAndRejectWidthCount) { UnitCell cell; @@ -311,7 +401,7 @@ TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) 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); + elecstate::PotSccs matched(&basis, config, solver, 0.5, no_resume); ModuleBase::matrix matched_potential(1, basis.nrxx); matched.cal_v_eff(&charge, &cell, matched_potential); warning_log.flush(); @@ -321,7 +411,7 @@ TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) 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); + elecstate::PotSccs mismatched(&basis, config, solver, 1.0, no_resume); ModuleBase::matrix mismatched_potential(1, basis.nrxx); mismatched.cal_v_eff(&charge, &cell, mismatched_potential); warning_log.close(); diff --git a/source/source_estate/unittests/test_potential_new.cpp b/source/source_estate/unittests/test_potential_new.cpp index e85130cadd..83930ce75c 100644 --- a/source/source_estate/unittests/test_potential_new.cpp +++ b/source/source_estate/unittests/test_potential_new.cpp @@ -1,4 +1,5 @@ #include "source_estate/module_pot/potential_new.h" +#include "source_estate/module_pot/pot_sccs.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" @@ -162,12 +163,18 @@ int MockPotComponent::destroyed = 0; int MockPotComponent::fixed_calls = 0; int MockPotComponent::dynamic_calls = 0; -PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para&) +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para&, const SccsResume&) { const int grid_size = this->get_rho_basis() == nullptr ? 0 : this->get_rho_basis()->nrxx; return new MockPotComponent(pot_type, grid_size); } +// The mock components hold no SCCS state. +SccsResume Potential::get_resume() const +{ + return SccsResume(); +} + } // namespace elecstate class PotentialNewTest : public ::testing::Test diff --git a/source/source_estate/unittests/test_sccs_fixed_src.cpp b/source/source_estate/unittests/test_sccs_fixed_src.cpp new file mode 100644 index 0000000000..981404eb6f --- /dev/null +++ b/source/source_estate/unittests/test_sccs_fixed_src.cpp @@ -0,0 +1,196 @@ +#include "source_estate/module_pot/sccs_fixed_src.h" +#include "source_estate/module_pot/pot_sccs.h" +#include "source_cell/cell_tools.h" +#include "source_hamilt/module_sccs/sccs_ion_charge.h" +#include "source_hamilt/module_sccs/test/sccs_test.h" +#include +#include + +namespace +{ +const elecstate::SccsResume no_resume = elecstate::SccsResume(); + +// The value after a key in an sccs_debug report. +std::string report_value(const std::string& report, const std::string& key) +{ + std::istringstream tokens(report); + std::string token; + while (tokens >> token) + { + if (token == key && tokens >> token) + { + return token; + } + } + return std::string(); +} +} + +class SccsFixedSourcesTest : public SccsTest::PwTest +{ +protected: + void SetUp() override + { + SccsTest::PwTest::SetUp(); + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + unitcell::AtomData atom; + atom.position = {2.5, 2.5, 2.5}; + atom.mass = 1.0; + atom.valence_charge = 1.0; + atom.atomic_number = 1; + atoms.push_back(atom); + atom.position.x = 4.0; + atom.mass = 16.0; + atom.valence_charge = 6.0; + atom.atomic_number = 8; + atoms.push_back(atom); + config.core_electrons = true; + config.core_spreads = {0.5, 0.7}; + } + void rebuild_and_compare() + { + EXPECT_FALSE(cache.update(cell, basis, atoms, config)); + std::vector expected; + ModuleSccs::gauss_ion_rho(atoms, basis, cell.tpiba, + ModuleSccs::gaussian_ion_spread, expected); + EXPECT_EQ(cache.ionic_density(), expected); + ModuleSccs::gauss_core_rho(atoms, basis, cell.tpiba, config.core_spreads, expected); + EXPECT_EQ(cache.core_density(), expected); + EXPECT_TRUE(cache.update(cell, basis, atoms, config)); + } + UnitCell cell; + std::vector atoms; + ModuleSccs::SccsConfig config; + elecstate::SccsFixedSources cache; +}; + +TEST_F(SccsFixedSourcesTest, ReusesSourcesUntilTheAtomsChange) +{ + EXPECT_TRUE(cache.ionic_density().empty()); + rebuild_and_compare(); + const double* stored_density = cache.ionic_density().data(); + EXPECT_TRUE(cache.update(cell, basis, atoms, config)); + EXPECT_EQ(cache.ionic_density().data(), stored_density); + EXPECT_EQ(cache.positions().size(), static_cast(basis.nrxx)); + atoms[1].mass = 17.0; // Not an input to the Gaussians. + EXPECT_TRUE(cache.update(cell, basis, atoms, config)); + atoms[0].position.x += 0.1; + rebuild_and_compare(); + atoms[1].valence_charge = 5.0; + rebuild_and_compare(); + atoms[0].atomic_number = 2; // Changes which atoms a single sccs_corespread skips (Z == zv). + rebuild_and_compare(); +} + +TEST_F(SccsFixedSourcesTest, CachedSourcesMatchFreshFullCavityEnergyPotentialAndForce) +{ + Atom atom; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.mass = 1.0; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 0.6 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channel = density.data(); + Charge charge; + charge.nspin = 1; + charge.rho = &channel; + config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + config.boundary = ModuleSccs::Boundary::Pcc0d; + config.cavity.epsilon_bulk = 5.0; + config.core_electrons = true; + config.core_spreads = {0.5}; + config.collect_report = true; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + elecstate::PotSccs cached(&basis, config, solver, 0.6, no_resume); + for (int step = 0; step < 3; ++step) + { + if (step == 1) { for (double& value : density) { value *= 0.95; } } + if (step == 2) + { + atom.tau[0].x += 0.01; + } + elecstate::PotSccs fresh(&basis, config, solver, 0.6, no_resume); + ModuleBase::matrix first(1, basis.nrxx); + ModuleBase::matrix second(1, basis.nrxx); + cached.cal_v_eff(&charge, &cell, first); + fresh.cal_v_eff(&charge, &cell, second); + EXPECT_NEAR(cached.get_energy(), fresh.get_energy(), 1e-11); + for (int i = 0; i < basis.nrxx; ++i) + { + EXPECT_NEAR(first(0, i), second(0, i), 1e-9); + } + ModuleBase::matrix force1(1, 3); + ModuleBase::matrix force2(1, 3); + cached.add_sol_force(cell, force1); + fresh.add_sol_force(cell, force2); + for (int i = 0; i < 3; ++i) + { + EXPECT_NEAR(force1(0, i), force2(0, i), 1e-9); + } + std::ostringstream output; + cached.write_iter(output, 2, 1e-5, 0.0); + const std::string reused = step == 1 ? "1" : "0"; + EXPECT_EQ(report_value(output.str(), "CACHED_SOURCES"), reused); + } +} + +TEST_F(SccsFixedSourcesTest, NextIonicStepResumesTheWarmStart) +{ + Atom atom; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.mass = 1.0; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 0.6 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channel = density.data(); + Charge charge; + charge.nspin = 1; + charge.rho = &channel; + config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + config.boundary = ModuleSccs::Boundary::Pcc0d; + config.cavity.epsilon_bulk = 5.0; + config.core_electrons = true; + config.core_spreads = {0.5}; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + elecstate::PotSccs first_step(&basis, config, solver, 0.6, no_resume); + ModuleBase::matrix first_potential(1, basis.nrxx); + first_step.cal_v_eff(&charge, &cell, first_potential); + const elecstate::SccsResume resume = first_step.resume_state(); + EXPECT_TRUE(resume.active); + EXPECT_EQ(resume.restart_potential.size(), static_cast(basis.nrxx)); + + // The rebuilt component of the next ionic step starts from the previous + // solution; moved atoms do not discard it. + atom.tau[0].x += 0.002; + elecstate::PotSccs resumed(&basis, config, solver, 0.6, resume); + elecstate::PotSccs cold(&basis, config, solver, 0.6, no_resume); + ModuleBase::matrix resumed_potential(1, basis.nrxx); + ModuleBase::matrix cold_potential(1, basis.nrxx); + resumed.cal_v_eff(&charge, &cell, resumed_potential); + cold.cal_v_eff(&charge, &cell, cold_potential); + EXPECT_NEAR(resumed.get_energy(), cold.get_energy(), 1e-10); + for (int i = 0; i < basis.nrxx; ++i) + { + EXPECT_NEAR(resumed_potential(0, i), cold_potential(0, i), 1e-8); + } + std::ostringstream resumed_output; + std::ostringstream cold_output; + resumed.write_iter(resumed_output, 2, 1e-5, 0.0); + cold.write_iter(cold_output, 2, 1e-5, 0.0); + EXPECT_EQ(report_value(resumed_output.str(), "WARM_START"), "1"); + EXPECT_EQ(report_value(cold_output.str(), "WARM_START"), "0"); +} diff --git a/source/source_estate/unittests/test_sccs_output.cpp b/source/source_estate/unittests/test_sccs_output.cpp new file mode 100644 index 0000000000..bd2635c10c --- /dev/null +++ b/source/source_estate/unittests/test_sccs_output.cpp @@ -0,0 +1,322 @@ +#include "source_estate/module_pot/sccs_output.h" +#include "source_estate/module_pot/pot_sccs.h" +#include "source_estate/module_pot/pot_pcc.h" +#include "source_hamilt/module_sccs/test/sccs_test.h" +#include "source_estate/module_pot/sol_grid_field.h" +#include "source_io/module_parameter/input_parameter.h" +#include "source_base/parallel_reduce.h" + +#include +#include +#include +#include + +namespace +{ +const elecstate::SccsResume no_resume = elecstate::SccsResume(); +} + +using SccsOutputTest = SccsTest::PwTest; + +TEST_F(SccsOutputTest, OutputLevelsSelectTheReportsAndPreserveStreamState) +{ + elecstate::SccsOutput result; + result.valid = true; + result.iterations = 4; + result.reaction_energy = -0.1; + std::ostringstream output; + output.setf(std::ios::scientific, std::ios::floatfield); + output.precision(3); + const std::ios::fmtflags flags = output.flags(); + elecstate::write_sccs_out(output, result, 0, -0.2, true, 0.05); + EXPECT_TRUE(output.str().empty()); + elecstate::write_sccs_out(output, result, 1, -0.2, true, 0.05); + // Cells of 16 characters, each starting on a column edge. + const std::string iteration_line = " SCCS_ITER CG_STEPS 4 E_SOL/Ry " + "-2.00000000e-01 E_PCC/Ry 5.00000000e-02\n"; + EXPECT_EQ(output.str(), iteration_line); + EXPECT_EQ(output.str().find("SCCS_RESIDUAL"), std::string::npos); + output.str(""); + elecstate::write_sccs_out(output, result, 2, -0.2, true, 0.05); + // Every cell starts on the grid of 16 columns: no cell runs into the next, + // every cell begins on its column edge, and so do the table headers. + const std::size_t cell = 16; + std::istringstream lines(output.str()); + std::string line; + int table_rows = 0; + while (std::getline(lines, line)) + { + for (std::size_t start = 1 + cell; start < line.size(); start += cell) + { + EXPECT_EQ(line[start - 1], ' ') << line; + EXPECT_NE(line[start], ' ') << line; + } + if (line.find("PCC0D_MOMENTS") == std::string::npos) + { + continue; + } + ++table_rows; + if (line.find("KIND") != std::string::npos) + { + EXPECT_EQ(line.find("Q/e"), 1 + 2 * cell); + EXPECT_EQ(line.find("QRR/eBohr2"), 1 + 6 * cell); + } + } + EXPECT_EQ(table_rows, 5); + EXPECT_NE(output.str().find(" PCC0D_MOMENTS SCREENED"), std::string::npos); + EXPECT_NE(output.str().find("CACHED_SOURCES 0"), std::string::npos); + EXPECT_EQ(output.str().find("V_FILLED"), std::string::npos); + EXPECT_EQ(output.precision(), 3); + EXPECT_EQ(output.flags(), flags); + output.str(""); + result.pcc.slab = true; + result.pcc.normal.z = 1.0; + elecstate::write_sccs_out(output, result, 2, -0.2, true, 0.05); + EXPECT_NE(output.str().find("PCC2D_ORIGIN"), std::string::npos); + EXPECT_NE(output.str().find(" PCC2D_MOMENTS SCREENED"), std::string::npos); + EXPECT_EQ(output.precision(), 3); + EXPECT_EQ(output.flags(), flags); +} + +TEST_F(SccsOutputTest, CellAsWideAsTheColumnKeepsOneBlank) +{ + elecstate::SccsOutput result; + result.valid = true; + result.iterations = 4; + std::ostringstream output; + // A three-digit exponent fills the whole 16-character cell. + elecstate::write_sccs_out(output, result, 1, -1.5e-100, true, 0.05); + EXPECT_NE(output.str().find("-1.50000000e-100 E_PCC/Ry"), std::string::npos) << output.str(); +} + +TEST_F(SccsOutputTest, AdvancedDiagnosticsPreservePccPotentialEnergyAndForce) +{ + 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.mass = 1.0; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + std::vector density(basis.nrxx, 1e-6); + const double electron_count = 1e-6 * basis.omega; + double* channel = density.data(); + Charge charge; + charge.nspin = 1; + charge.rho = &channel; + const char* boundaries[] = {"pcc_0d", "pcc_2d"}; + for (const char* boundary : boundaries) + { + Input_para input; + input.assume_isolated = boundary; + input.sccs_epsilon = 5.0; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::sccs_cfg_input(input, config, solver); + elecstate::PotSccs plain(&basis, config, solver, electron_count, no_resume); + config.collect_report = true; + elecstate::PotSccs checked(&basis, config, solver, electron_count, no_resume); + elecstate::PotPcc::Dimension dimension = elecstate::PotPcc::Dimension::molecule; + if (config.boundary == ModuleSccs::Boundary::Pcc2d) + { + dimension = elecstate::PotPcc::Dimension::slab; + } + elecstate::PotPcc pcc(&basis, dimension, 2, electron_count); + ModuleBase::matrix vacuum(1, basis.nrxx); + pcc.cal_v_eff(&charge, &cell, vacuum); + ModuleBase::matrix first(1, basis.nrxx); + ModuleBase::matrix second(1, basis.nrxx); + plain.cal_v_eff(&charge, &cell, first); + checked.cal_v_eff(&charge, &cell, second); + EXPECT_DOUBLE_EQ(plain.get_energy(), checked.get_energy()); + for (int i = 0; i < basis.nrxx; ++i) + { + EXPECT_DOUBLE_EQ(first(0, i), second(0, i)); + } + ModuleBase::matrix force1(1, 3); + ModuleBase::matrix force2(1, 3); + plain.add_sol_force(cell, force1); + checked.add_sol_force(cell, force2); + for (int i = 0; i < 3; ++i) + { + EXPECT_DOUBLE_EQ(force1(0, i), force2(0, i)); + } + const double pcc_energy = pcc.get_energy(); + std::ostringstream summary; + checked.write_iter(summary, 2, 1e-6, pcc_energy); + EXPECT_NE(summary.str().find("SCCS_GAUSS"), std::string::npos); + EXPECT_TRUE(checked.has_output()); + config.start_drho = 1e-4; + elecstate::PotSccs deferred(&basis, config, solver, electron_count, no_resume); + EXPECT_FALSE(deferred.is_active()); + EXPECT_FALSE(deferred.has_output()); + std::ostringstream waiting; + deferred.write_iter(waiting, 1, 1e-3, pcc_energy); + EXPECT_NE(waiting.str().find("SCCS_DEFERRED"), std::string::npos); + } +} + +// out_sol fields of PotSccs, written by ctrl_output_fp with collective calls. +using SccsSolventFieldTest = SccsTest::PwTest; + +TEST_F(SccsSolventFieldTest, SolventFieldsFollowTheLastActiveEvaluation) +{ + 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)}; + // A density across the switching range gives a nontrivial cavity. + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double step = ir % 7; + density[ir] = 1e-5 * (1.0 + step); + } + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + Input_para input; + elecstate::sccs_cfg_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + config.cavity.density_min = 2e-5; + config.cavity.density_max = 5e-5; + config.start_drho = 1e-3; + elecstate::PotSccs delayed(&basis, config, solver, 1.0, no_resume); + std::vector fields; + delayed.add_sol_fields(fields); + EXPECT_TRUE(fields.empty()); + config.start_drho = 0.0; + elecstate::PotSccs active(&basis, config, solver, 1.0, no_resume); + active.add_sol_fields(fields); + EXPECT_TRUE(fields.empty()); + ModuleBase::matrix potential(1, basis.nrxx); + active.cal_v_eff(&charge, &cell, potential); + active.add_sol_fields(fields); + ASSERT_EQ(fields.size(), 2u); + EXPECT_EQ(fields[0].name, "eps"); + EXPECT_EQ(fields[1].name, "cavity"); + ASSERT_EQ(fields[0].values.size(), static_cast(basis.nrxx)); + const double log_bulk = std::log(config.cavity.epsilon_bulk); + bool transition = false; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double solute = fields[1].values[ir]; + EXPECT_GE(solute, 0.0); + EXPECT_LE(solute, 1.0); + const double expected = std::exp(log_bulk * (1.0 - solute)); + EXPECT_NEAR(fields[0].values[ir], expected, 1e-12); + transition = transition || (solute > 0.0 && solute < 1.0); + } + EXPECT_TRUE(transition); +} + +// With more pool ranks than z planes a rank owns no real-space grid points. +class PotSccsEmptyPartitionTest : public SccsTest::EmptyPartitionTest +{ +protected: + // One active evaluation, then the out_sol field list. The pool sum over the + // field count is collective, so a rank with a different list hangs. + void evaluate_and_collect(const ModuleSccs::SccsConfig& config, + double& energy, + std::vector& names, + std::vector& sums) + { + double empty_ranks = basis.nrxx == 0 ? 1.0 : 0.0; + Parallel_Reduce::reduce_pool(empty_ranks); + if (SccsTest::pool_size >= 3) + { + EXPECT_GE(empty_ranks, 1.0); + } + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.latvec = lattice; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.label = "H"; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.1)}; + const std::vector mode = cosine_mode(0); + const double mean_density = 3.5e-5; + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = mean_density * (1.0 + 0.8 * mode[ir]); + } + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + ModuleSccs::PolarizationSolverParameters solver; + const double electron_count = mean_density * basis.omega; + elecstate::PotSccs component(&basis, config, solver, electron_count, no_resume); + ModuleBase::matrix potential(1, basis.nrxx); + component.cal_v_eff(&charge, &cell, potential); + energy = component.get_energy(); + + std::vector fields; + component.add_sol_fields(fields); + const std::size_t field_count = fields.size(); + names.clear(); + sums.assign(field_count, 0.0); + for (std::size_t field = 0; field < field_count; ++field) + { + names.push_back(fields[field].name); + ASSERT_EQ(fields[field].values.size(), static_cast(basis.nrxx)); + for (double value : fields[field].values) + { + sums[field] += value; + } + } + // Collective: a rank with a different field count would hang here. + Parallel_Reduce::reduce_pool(sums.data(), static_cast(field_count)); + } + + ModuleSccs::SccsConfig active_config() const + { + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + Input_para input; + elecstate::sccs_cfg_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + config.cavity.density_min = 2e-5; + config.cavity.density_max = 5e-5; + config.start_drho = 0.0; + return config; + } +}; + +// Every rank lists eps and cavity, including one without grid points, and the +// energy and fields match the serial run. +TEST_F(PotSccsEmptyPartitionTest, SolventOutputIsCollectiveOnEmptyGridPartitions) +{ + double energy = 0.0; + std::vector names; + std::vector sums; + evaluate_and_collect(active_config(), energy, names, sums); + const std::vector expected_names = {"eps", "cavity"}; + ASSERT_EQ(names, expected_names); + // Serial reference values. + EXPECT_NEAR(energy, -0.748979082383938, 1e-10); + EXPECT_NEAR(sums[0], 795.834560078087, 1e-8); + EXPECT_NEAR(sums[1], 154.256465345260, 1e-8); +} diff --git a/source/source_estate/unittests/test_sccs_pcc_0d.cpp b/source/source_estate/unittests/test_sccs_pcc_0d.cpp index d9ab205bc7..e57a9c25ad 100644 --- a/source/source_estate/unittests/test_sccs_pcc_0d.cpp +++ b/source/source_estate/unittests/test_sccs_pcc_0d.cpp @@ -1,6 +1,7 @@ #include "source_estate/module_pot/pot_sccs.h" #include "source_estate/module_pot/pot_pcc.h" #include "source_estate/module_pot/sccs_pcc_0d.h" +#include "source_basis/module_pw/pw_grid_geometry.h" #include "source_cell/cell_tools.h" #include "source_hamilt/module_sccs/sccs_coulomb.h" #include "source_hamilt/module_sccs/sccs_functional.h" @@ -10,6 +11,11 @@ #include "sccs_cell_test.h" #include "source_io/module_parameter/input_parameter.h" +namespace +{ +const elecstate::SccsResume no_resume = elecstate::SccsResume(); +} + class SccsPcc0dTest : public SccsTest::CellTest { protected: @@ -34,8 +40,10 @@ TEST_F(SccsPcc0dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) ModuleSccs::PolarizationSolverParameters solver; elecstate::sccs_cfg_input(input, config, solver); const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::vector> grid; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, grid); std::unique_ptr coulomb; - elecstate::pcc_0d_operator(cell, basis, atoms, coulomb); + elecstate::pcc_0d_operator(cell, basis, atoms, grid, coulomb); std::vector ions; ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, ModuleSccs::gaussian_ion_spread, ions); @@ -51,7 +59,7 @@ TEST_F(SccsPcc0dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) ModuleSccs::FunctionalResult functional; ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, *coulomb, functional); - elecstate::PotSccs solvent(&basis, config, solver, 0.6); + elecstate::PotSccs solvent(&basis, config, solver, 0.6, no_resume); elecstate::PotPcc correction(&basis, 0.6); ModuleBase::matrix potential(1, basis.nrxx); correction.cal_v_eff(&charge, &cell, potential); diff --git a/source/source_estate/unittests/test_sccs_pcc_2d.cpp b/source/source_estate/unittests/test_sccs_pcc_2d.cpp index 6c0f86ef3e..0960210f82 100644 --- a/source/source_estate/unittests/test_sccs_pcc_2d.cpp +++ b/source/source_estate/unittests/test_sccs_pcc_2d.cpp @@ -1,6 +1,7 @@ #include "source_estate/module_pot/pot_sccs.h" #include "source_estate/module_pot/pot_pcc.h" #include "source_estate/module_pot/sccs_pcc_2d.h" +#include "source_basis/module_pw/pw_grid_geometry.h" #include "source_cell/cell_tools.h" #include "source_hamilt/module_sccs/sccs_coulomb.h" #include "source_hamilt/module_sccs/sccs_functional.h" @@ -10,6 +11,11 @@ #include "sccs_cell_test.h" #include "source_io/module_parameter/input_parameter.h" +namespace +{ +const elecstate::SccsResume no_resume = elecstate::SccsResume(); +} + class SccsPcc2dTest : public SccsTest::CellTest { @@ -37,8 +43,10 @@ TEST_F(SccsPcc2dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) ModuleSccs::PolarizationSolverParameters solver; elecstate::sccs_cfg_input(input, config, solver); const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::vector> grid; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, grid); std::unique_ptr coulomb; - elecstate::pcc_2d_operator(cell, basis, atoms, input.pcc_2d_axis, coulomb); + elecstate::pcc_2d_operator(cell, basis, atoms, input.pcc_2d_axis, grid, coulomb); std::vector ions; ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, ModuleSccs::gaussian_ion_spread, ions); @@ -54,7 +62,7 @@ TEST_F(SccsPcc2dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) ModuleSccs::FunctionalResult functional; ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, *coulomb, functional); - elecstate::PotSccs solvent(&basis, config, solver, 0.6); + elecstate::PotSccs solvent(&basis, config, solver, 0.6, no_resume); elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, input.pcc_2d_axis, 0.6); ModuleBase::matrix potential(1, basis.nrxx); correction.cal_v_eff(&charge, &cell, potential); @@ -96,15 +104,18 @@ TEST_F(SccsPcc2dTest, AcceptsThePccSlabAlignmentTolerance) // digits, passes both PCC and the SCCS operator; a real tilt stops both. cell.latvec.e31 = 1.5e-8; const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::vector> grid; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, grid); std::unique_ptr coulomb; - elecstate::pcc_2d_operator(cell, basis, atoms, 2, coulomb); + elecstate::pcc_2d_operator(cell, basis, atoms, 2, grid, coulomb); EXPECT_NE(coulomb.get(), nullptr); elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, 2, 0.6); ModuleBase::matrix potential(1, basis.nrxx); correction.cal_v_eff(&charge, &cell, potential); cell.latvec.e31 = 1.5e-4; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, grid); std::unique_ptr tilted; - EXPECT_EXIT(elecstate::pcc_2d_operator(cell, basis, atoms, 2, tilted), + EXPECT_EXIT(elecstate::pcc_2d_operator(cell, basis, atoms, 2, grid, tilted), testing::ExitedWithCode(1), ""); } diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt index 5760b9f5b5..8cb7628ed9 100644 --- a/source/source_hamilt/module_sccs/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -1,4 +1,5 @@ add_library(sccs OBJECT + sccs_diag.cpp sccs_cavity.cpp sccs_cav_deriv.cpp sccs_parameters.cpp diff --git a/source/source_hamilt/module_sccs/sccs_coul_0d.h b/source/source_hamilt/module_sccs/sccs_coul_0d.h index 1c90d7afd3..79d8c5b7d9 100644 --- a/source/source_hamilt/module_sccs/sccs_coul_0d.h +++ b/source/source_hamilt/module_sccs/sccs_coul_0d.h @@ -17,6 +17,7 @@ class Pcc0dCoulombOperator : public CoulombOperator const std::vector>& relative_positions, const elecstate::Pcc0dParameters& parameters); bool has_correction() const override; + CoulombTransformCounts fft_counts() const override { return periodic_.fft_counts(); } void apply_potential(const std::vector& charge, std::vector& potential) override; private: diff --git a/source/source_hamilt/module_sccs/sccs_coul_2d.h b/source/source_hamilt/module_sccs/sccs_coul_2d.h index 5b1f3a333c..c0e55199a8 100644 --- a/source/source_hamilt/module_sccs/sccs_coul_2d.h +++ b/source/source_hamilt/module_sccs/sccs_coul_2d.h @@ -18,6 +18,7 @@ class Pcc2dCoulombOperator : public CoulombOperator const ModuleBase::Vector3& normal, const elecstate::Pcc2dParameters& parameters); bool has_correction() const override; + CoulombTransformCounts fft_counts() const override { return periodic_.fft_counts(); } void apply_potential(const std::vector& charge, std::vector& potential) override; private: diff --git a/source/source_hamilt/module_sccs/sccs_coulomb.h b/source/source_hamilt/module_sccs/sccs_coulomb.h index 0878021579..fa365d59ae 100644 --- a/source/source_hamilt/module_sccs/sccs_coulomb.h +++ b/source/source_hamilt/module_sccs/sccs_coulomb.h @@ -5,6 +5,12 @@ namespace ModuleSccs { +struct CoulombTransformCounts +{ + int forward_calls = 0; + int inverse_calls = 0; +}; + // Potential in Hartree for signed charge density in e/Bohr^3. Applications // are collective over the PW pool. class CoulombOperator @@ -12,6 +18,7 @@ class CoulombOperator public: virtual ~CoulombOperator() {} virtual bool has_correction() const = 0; + virtual CoulombTransformCounts fft_counts() const { return {}; } virtual void apply_potential(const std::vector& charge, std::vector& potential) = 0; }; } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_diag.cpp b/source/source_hamilt/module_sccs/sccs_diag.cpp new file mode 100644 index 0000000000..07ae5bc730 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_diag.cpp @@ -0,0 +1,29 @@ +#include "sccs_diag.h" +#include "sccs_response.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include + +namespace ModuleSccs +{ +void cont_pol_charge(const std::vector& charge, + const SccsResponse& response, + const ModulePW::PW_Basis& basis, + double tpiba, + std::vector& density) +{ + std::vector> potential_g(basis.npw); + basis.real2recip(response.polarization.potential.data(), potential_g.data()); + std::vector> gradient(charge.size()); + XC_Functional::grad_rho(potential_g.data(), gradient.data(), &basis, tpiba); + density.resize(charge.size()); + for (std::size_t i = 0; i < charge.size(); ++i) + { + const double contraction = response.grad_log_epsilon[i] * gradient[i]; + density[i] = contraction / ModuleBase::FOUR_PI + charge[i] * (1.0 / response.epsilon[i] - 1.0); + } +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_diag.h b/source/source_hamilt/module_sccs/sccs_diag.h new file mode 100644 index 0000000000..568f77b2ee --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_diag.h @@ -0,0 +1,18 @@ +#ifndef SCCS_DIAG_H +#define SCCS_DIAG_H + +#include + +namespace ModulePW { class PW_Basis; } +namespace ModuleSccs +{ +struct SccsResponse; + +// Environ dielectric_of_potential; diagnostic density, not a new solve source. +void cont_pol_charge(const std::vector& charge, + const SccsResponse& response, + const ModulePW::PW_Basis& basis, + double tpiba, + std::vector& density); +} // namespace ModuleSccs +#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h index 145546b199..9277f26e3c 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.h +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -17,6 +17,9 @@ struct SccsConfig CavityParameters cavity; Boundary boundary = Boundary::Periodic; int pcc_2d_axis = 2; + double start_drho = 0.0; // SCF DRHO at which a delayed SCCS starts; zero starts immediately + int start_nmax = 30; // Latest electronic iteration of a delayed start + bool collect_report = false; // sccs_debug 2: PCC moments and Gauss-law check for the report bool core_electrons = false; std::vector core_spreads = {0.5}; // Bohr: singleton skips Z == zv; per-atom list overrides; <= 0 disables double surface_tension = 0.0; // Hartree/Bohr^2 diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp index 36f7e678a5..db2c4bdc35 100644 --- a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp @@ -19,6 +19,7 @@ void PeriodicCoulombOperator::apply_potential(const std::vector& charge, { reciprocal_work_.resize(basis_.npw); basis_.real2recip(charge.data(), reciprocal_work_.data()); + ++counts_.forward_calls; const double tpiba2 = tpiba_ * tpiba_; for (int ig = 0; ig < basis_.npw; ++ig) { @@ -33,5 +34,6 @@ void PeriodicCoulombOperator::apply_potential(const std::vector& charge, } potential.resize(basis_.nrxx); basis_.recip2real(reciprocal_work_.data(), potential.data()); + ++counts_.inverse_calls; } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h index adcce673e3..70c0756d5d 100644 --- a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h @@ -21,9 +21,11 @@ class PeriodicCoulombOperator : public CoulombOperator public: PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba); bool has_correction() const override; + CoulombTransformCounts fft_counts() const override { return counts_; } void apply_potential(const std::vector& charge, std::vector& potential) override; private: + CoulombTransformCounts counts_; const ModulePW::PW_Basis& basis_; const double tpiba_; // Scratch only, overwritten on each call; PW_Basis FFTs are also non-reentrant. diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt index ad98a0503c..bb6c64dccf 100644 --- a/source/source_hamilt/module_sccs/test/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -10,6 +10,11 @@ AddTest(TARGET MODULE_HAMILT_sccs_poisson ../../../source_estate/unittests/test_pot_sccs.cpp ../../../source_estate/unittests/test_sccs_pcc_0d.cpp ../../../source_estate/unittests/test_sccs_pcc_2d.cpp + ../../../source_estate/unittests/test_sccs_fixed_src.cpp + ../../../source_estate/module_pot/sccs_fixed_src.cpp + ../../../source_estate/unittests/test_sccs_output.cpp + ../../../source_estate/module_pot/sccs_output.cpp + ../../../source_estate/module_pot/sccs_collect.cpp ../../../source_estate/module_pot/pot_sccs.cpp ../../../source_estate/module_pot/sccs_pcc_0d.cpp ../../../source_estate/module_pot/sccs_pcc_2d.cpp @@ -20,6 +25,7 @@ AddTest(TARGET MODULE_HAMILT_sccs_poisson ../../../source_estate/charge_moments.cpp ../../../source_estate/pcc_0d.cpp ../sccs_coul_0d.cpp ../sccs_coul_2d.cpp + ../sccs_diag.cpp ../sccs_moments.cpp ../sccs_lowpass.cpp ../sccs_cav_deriv.cpp @@ -36,4 +42,11 @@ if(ENABLE_MPI) $ ${MPIEXEC_POSTFLAGS}) set_tests_properties(MODULE_HAMILT_sccs_poisson_mpi PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60 PROCESSORS 2) + # Three ranks on a grid with two z planes leave one rank with nrxx == 0; + # a rank that skips a collective hangs, and the timeout turns that into a failure. + add_test(NAME MODULE_HAMILT_sccs_empty_partition_mpi + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 3 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS} --gtest_filter=*EmptyPartition*) + set_tests_properties(MODULE_HAMILT_sccs_empty_partition_mpi PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60 PROCESSORS 3) endif() diff --git a/source/source_hamilt/module_sccs/test/sccs_test.h b/source/source_hamilt/module_sccs/test/sccs_test.h index c79017964e..75a535776a 100644 --- a/source/source_hamilt/module_sccs/test/sccs_test.h +++ b/source/source_hamilt/module_sccs/test/sccs_test.h @@ -49,6 +49,22 @@ class PwTest : public testing::Test basis.initmpi(pool_size, pool_rank, POOL_WORLD); #endif basis.initgrids(length, lattice, 20.0); + finish_basis(); + } + + // Explicit FFT grid. A pool with more ranks than nz leaves ranks without + // real-space grid points (nrxx == 0); they still own G vectors. + void set_up_grid(const ModuleBase::Matrix3& lattice, int nx, int ny, int nz) + { +#ifdef __MPI + basis.initmpi(pool_size, pool_rank, POOL_WORLD); +#endif + basis.initgrids(length, lattice, nx, ny, nz); + finish_basis(); + } + + void finish_basis() + { basis.initparameters(false, 20.0, 1, false); basis.setuptransform(); basis.collect_local_pw(); @@ -76,6 +92,20 @@ class PwTest : public testing::Test ModulePW::PW_Basis basis; std::string error; }; +// A cell short along c on a 12 x 12 x 2 grid: with three pool ranks the +// third owns no z plane, while the G sphere keeps sticks for every rank. +class EmptyPartitionTest : public PwTest +{ +protected: + void SetUp() override + { + set_up_grid(lattice, 12, 12, 2); + } + + const ModuleBase::Matrix3 lattice = ModuleBase::Matrix3(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 0.25); +}; } // namespace SccsTest #endif 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 index a36b7247fe..e8a5d52289 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp @@ -36,6 +36,21 @@ TEST_F(SccsPwCoulombTest, FourierModeAndConstantBackground) } } +// Each application is one forward and one inverse FFT, reported by sccs_debug 2. +TEST_F(SccsPwCoulombTest, FftCountsFollowTheApplications) +{ + ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); + const std::vector charge = cosine_mode(0); + std::vector potential; + EXPECT_EQ(coulomb.fft_counts().forward_calls, 0); + for (int application = 1; application <= 2; ++application) + { + coulomb.apply_potential(charge, potential); + EXPECT_EQ(coulomb.fft_counts().forward_calls, application); + EXPECT_EQ(coulomb.fft_counts().inverse_calls, application); + } +} + int main(int argc, char** argv) { int threads = 1; diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index 4a3077a2d6..efc315d256 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -95,6 +95,9 @@ void surchem::cal_vel(const UnitCell& cell, double *epsilon = new double[rho_basis->nrxx]; double *epsilon0 = new double[rho_basis->nrxx]; cal_epsilon(rho_basis, PS_TOTN_real, epsilon, epsilon0); + double* epsilon_end = epsilon + rho_basis->nrxx; + this->last_epsilon_.assign(epsilon, epsilon_end); + this->last_epsilon_set_ = true; std::complex *Sol_phi = new std::complex[rho_basis->npw]; std::complex *Sol_phi0 = new std::complex[rho_basis->npw]; diff --git a/source/source_hamilt/module_surchem/h_corr_pw.cpp b/source/source_hamilt/module_surchem/h_corr_pw.cpp index 9cdb503fa4..cd9014500e 100644 --- a/source/source_hamilt/module_surchem/h_corr_pw.cpp +++ b/source/source_hamilt/module_surchem/h_corr_pw.cpp @@ -18,8 +18,6 @@ void surchem::v_correction(const UnitCell& cell, ModuleBase::TITLE("surchem", "v_cor"); ModuleBase::timer::start("surchem", "v_cor"); - assert(rho_basis->nrxx>0); - double* porter = new double[rho_basis->nrxx]; for (int i = 0; i < rho_basis->nrxx; i++) { diff --git a/source/source_hamilt/module_surchem/surchem.h b/source/source_hamilt/module_surchem/surchem.h index 2c8f73157f..4208e06ab2 100644 --- a/source/source_hamilt/module_surchem/surchem.h +++ b/source/source_hamilt/module_surchem/surchem.h @@ -7,6 +7,8 @@ #include "source_basis/module_pw/pw_basis.h" #include "source_cell/unitcell.h" +#include + // forward-declared: used below only as pointer/reference class Parallel_Grid; class Structure_Factor; @@ -156,9 +158,18 @@ class surchem void induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, double* induced_rho) const; + // Dielectric function of the last cal_vel on the local grid; empty before it + // and on a rank without real-space grid points. + const std::vector& last_epsilon() const { return last_epsilon_; } + // True on every rank once cal_vel has run, whatever the local grid size. + bool has_last_eps() const { return last_epsilon_set_; } + double bulk_epsilon() const { return parameters_.eb_k; } + private: SurchemParameters parameters_; bool parameters_set_ = false; + std::vector last_epsilon_; + bool last_epsilon_set_ = false; }; #endif \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/test/cal_vel_test.cpp b/source/source_hamilt/module_surchem/test/cal_vel_test.cpp index 88678e487e..a71a73714f 100644 --- a/source/source_hamilt/module_surchem/test/cal_vel_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_vel_test.cpp @@ -235,6 +235,20 @@ TEST_F(cal_vel_test, cal_vel) EXPECT_NEAR(v_res(0, 0), 0.0532168705, 1e-10); EXPECT_NEAR(v_res(0, 1), 0.0447818244, 1e-10); + // out_sol: cal_vel keeps the dielectric function of this correction. + std::vector ps_totn_real(nrxx); + pwtest.recip2real(PS_TOTN, ps_totn_real.data()); + std::vector epsilon(nrxx); + std::vector epsilon0(nrxx); + solvent_model.cal_epsilon(&pwtest, ps_totn_real.data(), epsilon.data(), epsilon0.data()); + const std::vector& stored = solvent_model.last_epsilon(); + ASSERT_EQ(stored.size(), static_cast(nrxx)); + for (int ir = 0; ir < nrxx; ++ir) + { + EXPECT_EQ(stored[ir], epsilon[ir]); + } + EXPECT_DOUBLE_EQ(solvent_model.bulk_epsilon(), 80.0); + delete[] PS_TOTN; delete[] TOTN; } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index a57166a82b..57bed1d3a1 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -4,6 +4,7 @@ #include "source_estate/elecstate.h" // use elecstate::ElecState #include "source_estate/module_charge/charge.h" // use Charge #include "source_estate/module_charge/chg_symm.h" // use module_charge::cal_rhog_symm +#include "source_estate/module_pot/sol_grid_field.h" // use SolventGridField #include "source_estate/write_elecstat_pot.h" // use write_elecstat_pot #include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional #include "source_io/module_dipole/dipole_io.h" // use write_dipole @@ -180,6 +181,31 @@ void ctrl_output_fp(UnitCell& ucell, ofs_running); } + // 5b) write the implicit-solvent dielectric function and cavity + if (inp.out_sol[0] == 1 && should_output) + { + const std::vector fields = pelec->pot->sol_out_fields(); + for (const elecstate::SolventGridField& field : fields) + { + const std::string fn = global_out_dir + "sol_" + field.name + geom_block + ".cube"; + const std::string data_desc = "implicit solvent " + field.name; + ModuleIO::write_vdata_palgrid(para_grid, + field.values.data(), + -1, // spin-independent field + nspin, + istep, + fn, + 0.0, // efermi + &(ucell), + inp.out_sol[1], // precision + 0, // out_fermi + false, // two_fermi, unused without out_fermi + false, // reduce_all_pool + ofs_running, + data_desc); + } + } + // 6) write ELF if (inp.out_elf[0] > 0 && should_output) { diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1caf2ec975..677e21920f 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -465,6 +465,7 @@ struct Input_para ///< frequency controlled by out_freq_ion/out_freq_elec std::vector out_xc_r = {-1, 3}; ///< output xc(r). -1: no; >=0: output the order of xc(r) std::vector out_pot = {0, 8}; ///< output potential + std::vector out_sol = {0, 8}; ///< output the implicit-solvent dielectric function and cavity int out_wfc_pw = 0; ///< 0: no; 1: txt; 2: dat std::vector out_band = {0, 8}; ///< band calculation pengfei 2014-10-13 int out_dos = 0; ///< dos calculation. mohan add 20090909 @@ -612,6 +613,9 @@ struct Input_para double sccs_surface_eta = 1e-8; // Bohr^-1 double sccs_lowpass_p1 = -1.0; // PCC switching-derivative filter slope double sccs_lowpass_p2 = -1.0; // PCC switching-derivative filter offset + double sccs_start_drho = 0.0; + int sccs_start_nmax = 30; + int sccs_debug = 0; // 0: silent, 1: iteration summary, 2: full diagnostics std::string sccs_solvent_mode = "electronic"; std::vector sccs_corespread = {0.5}; // Bohr: one value or one per atom diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index 37aa6d915a..9899698050 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -192,13 +192,13 @@ void ReadInput::item_model() "* 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; assume_isolated none, pcc_0d or " - "pcc_2d; no electric/gate field, DFT-1/2, DeePKS or dm_to_rho. Use assume_isolated pcc_0d " - "or pcc_2d for charged systems: with assume_isolated none a charged cell only prints a " - "warning and runs, and its energy then depends on the cell size because the periodic " - "Poisson solver drops the G = 0 net charge. With pcc_0d or pcc_2d, see sccs_lowpass_p1 " - "for the accuracy of forces.\n\n" + "(esolver_type ksdft, device cpu) with basis_type pw or lcao; calculation scf (energies and " + "forces) or fixed-cell relax, no stress (cal_stress and cell-relax are rejected); nspin 1 " + "or 2; assume_isolated none, pcc_0d or pcc_2d; no electric/gate field, DFT-1/2, DeePKS or " + "dm_to_rho. Use assume_isolated pcc_0d or pcc_2d for charged systems: with assume_isolated " + "none a charged cell only prints a warning and runs, and its energy then depends on the " + "cell size because the periodic Poisson solver drops the G = 0 net charge. With pcc_0d or " + "pcc_2d, see sccs_lowpass_p1 for the accuracy of forces.\n\n" "Legacy Boolean values (true/false) are read as 1/0."; item.default_value = "0"; item.read_value = [](const Input_Item& item, Parameter& para) { diff --git a/source/source_io/module_parameter/read_inp_out.cpp b/source/source_io/module_parameter/read_inp_out.cpp index 1755de440c..9830989fa2 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -180,6 +180,68 @@ In molecular dynamics calculations, the output frequency is controlled by out_fr sync_intvec(input.out_pot, 2, 0); this->add_item(item); } + { + Input_Item item("out_sol"); + item.annotation = "output the solvent dielectric function and cavity (with precision 8)"; + item.category = "Output information"; + item.type = R"(Integer \[Integer\](optional))"; + item.description = R"(* 0: No output. +* 1: Output the implicit-solvent fields of the last solvent correction on real space grids into files in the folder OUT.{suffix}, for imp_sol 1 (the original implicit solvation model) and imp_sol 2 (SCCS): + * `sol_eps.cube`: the dielectric function epsilon(r), 1 in the solute and the bulk permittivity in the solvent. + * `sol_cavity.cube`: the solute boundary s(r), 1 in the solute and 0 in the solvent; for imp_sol 1 it is 1 minus the solvent shape function (written only when eb_k > 1). + With out_freq_ion > 0 the files are written every out_freq_ion ionic steps with the geometry index in the name, e.g. `sol_epsg1.cube`; with out_freq_ion = 0 they are overwritten at every step. Nothing is written while SCCS waits for its delayed start (sccs_start_drho). + +The optional second integer controls the output precision. If not provided, the default precision is 8.)"; + item.default_value = "0"; + item.unit = ""; + item.read_value = [](const Input_Item& item, Parameter& para) { + const size_t count = item.get_size(); + if (count < 1 || count > 2) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_sol takes 1 or 2 integers"); + } + // The whole token must be an integer: std::stoi alone accepts "6junk". + auto parse_int = [](const std::string& text) { + size_t parsed = 0; + int value = 0; + try + { + value = std::stoi(text, &parsed); + } + catch (const std::exception&) + { + parsed = 0; + } + if (parsed == 0 || parsed != text.size()) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_sol takes integers"); + } + return value; + }; + para.input.out_sol[0] = parse_int(item.str_values[0]); + para.input.out_sol[1] = 8; + if (count == 2) + { + para.input.out_sol[1] = parse_int(item.str_values[1]); + } + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.out_sol[0] != 0 && para.input.out_sol[0] != 1) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_sol must be 0 or 1"); + } + if (para.input.out_sol[0] == 1 && para.input.imp_sol == 0) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_sol 1 requires imp_sol 1 or 2"); + } + if (para.input.out_sol[1] <= 0) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_sol precision must be positive"); + } + }; + sync_intvec(input.out_sol, 2, 0); + this->add_item(item); + } { Input_Item item("out_dmk"); item.annotation = ">0 output density matrix DM(k) for each k-point"; diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 2c8761c9a5..3d48b6c504 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -183,5 +183,46 @@ void ReadInput::item_sccs() sync_doublevec(input.sccs_corespread, para.input.sccs_corespread.size(), 0.0); this->add_item(item); } + { + Input_Item item("sccs_start_drho"); + item.annotation = "SCCS delayed-start density threshold"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Delays the start of SCCS until DRHO is at or below this value; until then the SCF runs without the solvent, while PCC, if set, stays on. 0 (default) starts SCCS immediately. A positive value must be larger than scf_thr, so that the SCF cannot converge before SCCS starts. When SCCS starts, the charge mixing restarts, that iteration does not count as converged, and the iteration is printed on the screen. SCCS then stays on for the rest of the run, including later ionic steps. Not set by sccs_preset."; + item.default_value = "0.0"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_start_drho); + this->add_item(item); + } + { + Input_Item item("sccs_start_nmax"); + item.annotation = "SCCS delayed-start iteration limit"; + item.category = "Implicit solvation model"; + item.type = "Integer"; + item.description = "With sccs_start_drho > 0, starts SCCS at this electronic iteration at the latest, if DRHO has not yet reached sccs_start_drho. Must be positive, and smaller than scf_nmax when sccs_start_drho > 0. Default 30. Not set by sccs_preset."; + item.default_value = "30"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_int(input.sccs_start_nmax); + this->add_item(item); + } + { + Input_Item item("sccs_debug"); + item.annotation = "SCCS report level"; + item.category = "Implicit solvation model"; + item.type = "Integer"; + item.description = "Level of the SCCS reports, written to the running log (running_*.log in OUT.{suffix}): 0 (default) writes none; 1 writes the sqrt-CG iteration count and the correction energies of every SCF iteration; 2 also writes the residuals, warm start, cavity volume and surface, the FFT counts of the Coulomb operator in the sqrt-CG solve and the count of reused fixed sources, and with PCC the Gauss-law check (far-field polarization charge and its expected value), the charge moments and the energy components."; + item.default_value = "0"; + item.unit = ""; + read_sync_int(input.sccs_debug); + item.check_value = [](const Input_Item&, const Parameter& para) { + if (para.input.sccs_debug < 0 || para.input.sccs_debug > 2) + { + ModuleBase::WARNING_QUIT("ReadInput", "sccs_debug must be 0, 1, or 2"); + } + }; + this->add_item(item); + } } } // namespace ModuleIO diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 43adbfc816..3977bfef04 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -563,7 +563,7 @@ Available options are: * pcc_2d: self-consistent PCC correction for a slab, open along the lattice vector selected by pcc_2d_axis (default: the third) and periodic along the other two. The open vector must be perpendicular to the periodic plane (rewrite a tilted cell as the equivalent perpendicular one), and the k-point sampling along it must be Gamma only. The monopole term uses the open planar kernel that vanishes on the plane of the charge, -pi*Q*L/(3*A) with L the cell length along the open vector and A the periodic area, so the energy of a charged slab converges with the vacuum size; it is referenced to zero potential on the plane of the charge. ENVIRON uses -pi*Q/(3*L), so charged-slab energies agree with ENVIRON only for A = L^2. * makov-payne, m-p, mp: compute the Makov-Payne correction to the total energy and estimate a corrected vacuum level for eigenvalue alignment. This option is available only for cubic lattices (latname = sc, fcc, or bcc). -pcc_0d and pcc_2d contribute to the energy, the potential and the fixed-cell forces; the correction energy is printed as E_pcc. They require CPU KS-DFT (esolver_type ksdft) with basis_type pw or lcao, calculation scf or relax, and nspin 1 or 2, without efield_flag, gate_flag, cal_stress, DFT-1/2, deepks output or dm_to_rho. With imp_sol=2, pcc_0d and pcc_2d support SCF energies and forces; imp_sol=1 is not supported. +pcc_0d and pcc_2d contribute to the energy, the potential and the fixed-cell forces; the correction energy is printed as E_pcc. They require CPU KS-DFT (esolver_type ksdft) with basis_type pw or lcao, calculation scf or relax, and nspin 1 or 2, without efield_flag, gate_flag, cal_stress, DFT-1/2, deepks output or dm_to_rho. With imp_sol=2, pcc_0d and pcc_2d support SCF energies, forces and fixed-cell relax; imp_sol=1 is not supported. Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995); PCC: O. Andreussi and N. Marzari, Phys. Rev. B 90, 245101 (2014).)"; item.default_value = "none"; diff --git a/source/source_io/module_parameter/sccs_inp_check.cpp b/source/source_io/module_parameter/sccs_inp_check.cpp index d056775321..670fa3634e 100644 --- a/source/source_io/module_parameter/sccs_inp_check.cpp +++ b/source/source_io/module_parameter/sccs_inp_check.cpp @@ -83,9 +83,9 @@ std::string context_error(const Input_para& input) { error = "SCCS requires basis_type pw or lcao"; } - else if (input.calculation != "scf" || input.cal_stress) + else if ((input.calculation != "scf" && input.calculation != "relax") || input.cal_stress) { - error = "SCCS currently supports SCF energies and forces without stress"; + error = "SCCS supports calculation scf or fixed-cell relax without stress"; } else if (input.nspin != 1 && input.nspin != 2) { @@ -120,7 +120,8 @@ std::string nonfinite_error(const Input_para& input) {"sccs_tol_max", input.sccs_tol_max}, {"sccs_surface_eta", input.sccs_surface_eta}, {"sccs_lowpass_p1", input.sccs_lowpass_p1}, - {"sccs_lowpass_p2", input.sccs_lowpass_p2}}; + {"sccs_lowpass_p2", input.sccs_lowpass_p2}, + {"sccs_start_drho", input.sccs_start_drho}}; for (const std::pair& value : values) { if (!std::isfinite(value.second)) @@ -174,9 +175,10 @@ std::string model_error(const Input_para& input) return error; } -// sqrt-CG solver controls. +// sqrt-CG solver controls and the delayed start. std::string solver_error(const Input_para& input) { + const bool delayed_start = input.sccs_start_drho > 0.0; std::string error; if (input.sccs_maxiter <= 0) { @@ -190,6 +192,18 @@ std::string solver_error(const Input_para& input) { error = "sccs_tol_max must be positive"; } + else if (input.sccs_start_drho < 0.0) + { + error = "sccs_start_drho must be non-negative"; + } + else if (delayed_start && input.sccs_start_drho <= input.scf_thr) + { + error = "sccs_start_drho must be zero or larger than scf_thr"; + } + else if (input.sccs_start_nmax <= 0 || (delayed_start && input.sccs_start_nmax >= input.scf_nmax)) + { + error = "sccs_start_nmax must be positive and smaller than scf_nmax when delayed start is enabled"; + } return error; } } // namespace diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index b1221ff0b7..354f0925fc 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -723,6 +723,58 @@ TEST_F(InputTest, Item_test) EXPECT_EQ(TestParameters::input(param).out_pot[0], 1); EXPECT_EQ(TestParameters::input(param).out_pot[1], 2); } + { // out_sol + auto it = find_label("out_sol", readinput.input_lists); + it->second.str_values = {"1"}; + it->second.read_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).out_sol[0], 1); + EXPECT_EQ(TestParameters::input(param).out_sol[1], 8); + it->second.str_values = {"1", "5"}; + it->second.read_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).out_sol[1], 5); + + TestParameters::input(param).imp_sol = 2; + check_item(it->second, param); + TestParameters::input(param).imp_sol = 0; + testing::internal::CaptureStdout(); + EXPECT_EXIT(check_item(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + TestParameters::input(param).imp_sol = 1; + TestParameters::input(param).out_sol[0] = 2; + testing::internal::CaptureStdout(); + EXPECT_EXIT(check_item(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + TestParameters::input(param).out_sol[0] = 1; + TestParameters::input(param).out_sol[1] = 0; + testing::internal::CaptureStdout(); + EXPECT_EXIT(check_item(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + it->second.str_values = {"x"}; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.read_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("out_sol takes integers")); + it->second.str_values = {"6junk"}; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.read_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("out_sol takes integers")); + it->second.str_values = {"1", "8x"}; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.read_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("out_sol takes integers")); + it->second.str_values = {"1", "8", "3"}; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.read_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("out_sol takes 1 or 2 integers")); + TestParameters::input(param).out_sol = {0, 8}; + TestParameters::input(param).imp_sol = 0; + } { // out_dos auto it = find_label("out_dos", readinput.input_lists); TestParameters::input(param).calculation = "get_wf"; diff --git a/source/source_io/test_serial/test_read_inp_sccs.cpp b/source/source_io/test_serial/test_read_inp_sccs.cpp index 37dcec861b..18cbabb607 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -54,7 +54,7 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); input.cal_force = false; input.calculation = "relax"; - EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); input.calculation = "scf"; input.assume_isolated = "none"; input.sccs_rho_min = input.sccs_rho_max; @@ -110,6 +110,63 @@ TEST(ReadInpSccs, FullCavityWidthsAndModeValidation) EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); } +TEST(ReadInpSccs, DelayedStartRequiresRoomBeforeScfConvergenceAndIterationLimit) +{ + Input_para input; + input.imp_sol = 2; + input.scf_thr = 1e-9; + input.scf_nmax = 10; + std::string error; + ASSERT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_drho = 1e-3; + input.sccs_start_nmax = 3; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_drho = input.scf_thr; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_drho = -1.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_drho = 1e-3; + input.sccs_start_nmax = 10; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_drho = 0.0; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.sccs_start_nmax = 0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); +} + +TEST(ReadInpSccs, FixedCellRelaxSupportsAllBoundariesAndKeepsStressOutOfScope) +{ + Input_para input; + input.imp_sol = 2; + input.calculation = "relax"; + input.cal_force = true; + input.sccs_solvent_mode = "full"; + input.sccs_corespread = {0.0, 0.5}; + input.sccs_start_drho = 1e-3; + input.sccs_start_nmax = 3; + std::string error; + const std::string boundaries[] = {"none", "pcc_0d", "pcc_2d"}; + const std::string bases[] = {"pw", "lcao"}; + for (const std::string& basis : bases) + { + input.basis_type = basis; + for (const std::string& boundary : boundaries) + { + input.assume_isolated = boundary; + ASSERT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; + input.cal_stress = true; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.cal_stress = false; + } + } + const std::string unsupported[] = {"cell-relax", "md", "nscf"}; + for (const std::string& calculation : unsupported) + { + input.calculation = calculation; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + } +} + TEST(ReadInpSccs, CoreSpreadWordsMustBeFiniteNumbers) { std::vector spreads = {0.5}; @@ -145,7 +202,7 @@ TEST(ReadInpSccs, RejectsNonfiniteReals) double* const reals[] = {&input.sccs_epsilon, &input.sccs_rho_min, &input.sccs_rho_max, &input.sccs_gamma, &input.sccs_pressure, &input.sccs_tol_rms, &input.sccs_tol_max, &input.sccs_surface_eta, &input.sccs_lowpass_p1, - &input.sccs_lowpass_p2}; + &input.sccs_lowpass_p2, &input.sccs_start_drho}; const double infinity = std::numeric_limits::infinity(); const double nonfinite[] = {std::nan(""), infinity, -infinity}; for (double* value : reals) diff --git a/source/source_io/unittests/support/INPUT b/source/source_io/unittests/support/INPUT index 1e44965f10..ccd94790cd 100644 --- a/source/source_io/unittests/support/INPUT +++ b/source/source_io/unittests/support/INPUT @@ -59,6 +59,7 @@ init_chg atomic #start charge is from 'atomic' or file chg_extrap atomic #atomic; first-order; second-order; dm:coefficients of SIA out_chg 0 #>0 output charge density for selected electron steps out_pot 2 #output realspace potential +out_sol 0 6 #output solvent dielectric function and cavity out_wfc_pw 0 #output wave functions out_dos 0 #output energy and dos out_ldos 1 #output local density of states, second parameter controls the precision diff --git a/source/source_io/unittests/test_read_input_parallel.cpp b/source/source_io/unittests/test_read_input_parallel.cpp index 3c8ce4c54b..7476f0082a 100644 --- a/source/source_io/unittests/test_read_input_parallel.cpp +++ b/source/source_io/unittests/test_read_input_parallel.cpp @@ -200,6 +200,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.deepks_bandgap, 0); EXPECT_EQ(param.inp.deepks_out_unittest, 0); EXPECT_EQ(param.inp.out_pot[0], 2); + EXPECT_EQ(param.inp.out_sol[0], 0); + EXPECT_EQ(param.inp.out_sol[1], 6); EXPECT_EQ(param.inp.out_wfc_pw, 0); EXPECT_EQ(param.inp.out_dos, 0); EXPECT_EQ(param.inp.out_ldos[0], 1); diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 0b97a36781..1b8261b7a8 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -68,6 +68,7 @@ scf_in_chgdat 082_PW_gatefield scf_sol_H2 scf_sol_H2O +scf_sol_H2O_out_sol 085_PW_get_pchg 086_PW_get_wf 087_PW_get_pchg_kpar_bndpar @@ -132,3 +133,4 @@ scf_sccs_pcc_2d scf_sccs_pcc_0d scf_sccs_pcc_0d_lowpass scf_sccs_full +relax_sccs diff --git a/tests/01_PW/relax_sccs/INPUT b/tests/01_PW/relax_sccs/INPUT new file mode 100644 index 0000000000..1ad788d83b --- /dev/null +++ b/tests/01_PW/relax_sccs/INPUT @@ -0,0 +1,25 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation relax +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 +relax_nmax 2 +sccs_start_drho 1.0e-3 +sccs_debug 1 +out_sol 1 diff --git a/tests/01_PW/relax_sccs/KPT b/tests/01_PW/relax_sccs/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/01_PW/relax_sccs/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/relax_sccs/README b/tests/01_PW/relax_sccs/README new file mode 100644 index 0000000000..f034c87904 --- /dev/null +++ b/tests/01_PW/relax_sccs/README @@ -0,0 +1 @@ +PW neutral water, fixed-cell relax with periodic SCCS, two ionic steps: SCCS starts with a delay in the first step, and the warm start and started state carry over to the rebuilt potential of the second step. out_sol 1 also checks the cavity volume integrated from sol_cavity.cube. ecutwfc 20 and ecutrho 240. diff --git a/tests/01_PW/relax_sccs/STRU b/tests/01_PW/relax_sccs/STRU new file mode 100644 index 0000000000..48d9dbb11e --- /dev/null +++ b/tests/01_PW/relax_sccs/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 1 1 1 +4.567 7.111 6.0 1 1 1 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/relax_sccs/result.ref b/tests/01_PW/relax_sccs/result.ref new file mode 100644 index 0000000000..0dc864fd2b --- /dev/null +++ b/tests/01_PW/relax_sccs/result.ref @@ -0,0 +1,7 @@ +etotref -442.0103622303894895 +etotperatomref -147.3367874101 +totalforceref 4.865466 +esolelref -0.3032039609 +esolcavref 0.0827906974 +sol_cavity_integralref 271.0889616 +totaltimeref 1.95 diff --git a/tests/01_PW/relax_sccs/threshold b/tests/01_PW/relax_sccs/threshold new file mode 100644 index 0000000000..a125b952d6 --- /dev/null +++ b/tests/01_PW/relax_sccs/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/scf_sol_H2O_out_sol/INPUT b/tests/01_PW/scf_sol_H2O_out_sol/INPUT new file mode 100644 index 0000000000..e4f93e95d1 --- /dev/null +++ b/tests/01_PW/scf_sol_H2O_out_sol/INPUT @@ -0,0 +1,24 @@ +INPUT_PARAMETERS + +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +cal_force 1 +#cal_stress 1 # implicit solvation stress is not implemented yet +nbands 20 +ecutwfc 20 +scf_nmax 50 + +#Parameters (Solvation Model) +imp_sol 1 +eb_k 80 +tau 0.000010798 +sigma_k 0.6 +nc_k 0.00037 + + +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +out_sol 1 diff --git a/tests/01_PW/scf_sol_H2O_out_sol/KPT b/tests/01_PW/scf_sol_H2O_out_sol/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/01_PW/scf_sol_H2O_out_sol/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sol_H2O_out_sol/README b/tests/01_PW/scf_sol_H2O_out_sol/README new file mode 100644 index 0000000000..b6e7d9e6a5 --- /dev/null +++ b/tests/01_PW/scf_sol_H2O_out_sol/README @@ -0,0 +1,2 @@ +test implicit solvation model correction and force for H2O with out_sol 1 +Same input as scf_sol_H2O plus out_sol; also checks the cavity volume integrated from sol_cavity.cube. diff --git a/tests/01_PW/scf_sol_H2O_out_sol/STRU b/tests/01_PW/scf_sol_H2O_out_sol/STRU new file mode 100644 index 0000000000..9aa513c379 --- /dev/null +++ b/tests/01_PW/scf_sol_H2O_out_sol/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 + +LATTICE_VECTORS +10 0 0 +0 10 0 +0 0 10 + +ATOMIC_POSITIONS +Cartesian # Cartesian(Unit is LATTICE_CONSTANT) + +H +0.0 +2 +-23.196290411877 12.932693996280 6.850457923460 0 0 0 +-24.006711594306 11.850811692381 9.324876689886 0 0 0 +O +0.0 +1 +-22.515669775626 12.546754901616 8.620909072323 0 0 0 diff --git a/tests/01_PW/scf_sol_H2O_out_sol/result.ref b/tests/01_PW/scf_sol_H2O_out_sol/result.ref new file mode 100644 index 0000000000..2201e0c477 --- /dev/null +++ b/tests/01_PW/scf_sol_H2O_out_sol/result.ref @@ -0,0 +1,7 @@ +etotref -442.8368758671132923 +etotperatomref -147.6122919557 +totalforceref 16.736445 +esolelref -0.9245824197 +esolcavref 0.0325548859 +sol_cavity_integralref 299.7598532 +totaltimeref 0.75 diff --git a/tests/02_NAO_Gamma/CASES_CPU.txt b/tests/02_NAO_Gamma/CASES_CPU.txt index 2bdb478b1c..ce8e970bdb 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -42,3 +42,5 @@ scf_pcc0d_water scf_pcc2d_water scf_sccs_pcc2d scf_sccs_water +scf_sccs_water_delay +relax_sccs_pcc0d_h3o diff --git a/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/INPUT b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/INPUT new file mode 100644 index 0000000000..6aef8f6af6 --- /dev/null +++ b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/INPUT @@ -0,0 +1,31 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +calculation relax +basis_type lcao +gamma_only 1 +ks_solver genelpa +nbands 8 +ecutwfc 60 +ecutrho 240 +scf_nmax 100 +scf_thr 1.0e-7 +smearing_method fixed +mixing_type pulay +mixing_beta 0.6 +symmetry 0 +init_wfc atomic +cal_force 1 +nelec_delta -1 + +imp_sol 2 +assume_isolated pcc_0d +sccs_preset water-cation +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +sccs_lowpass_p1 10 +sccs_lowpass_p2 5 +relax_nmax 2 diff --git a/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/KPT b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/README b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/README new file mode 100644 index 0000000000..91bf3965f0 --- /dev/null +++ b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/README @@ -0,0 +1,4 @@ +gamma-only LCAO H3O+ (nelec_delta -1) fixed-cell relax with SCCS (water-cation preset), 0D PCC and the cavity +low-pass filter (p1 10, p2 5), two ionic steps: the warm start of a charged solute carries over to the rebuilt +potential of the second step. ecutwfc 60 and ecutrho 240 as in scf_sccs_water; scf_thr 1e-7 and mixing_beta 0.6 +keep the run short. diff --git a/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/STRU b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/STRU new file mode 100644 index 0000000000..9bea9bb547 --- /dev/null +++ b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/STRU @@ -0,0 +1,29 @@ +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 +3 +7.80775 5.6969 6.29183 1 1 1 +4.94391 4.55719 6.49787 1 1 1 +5.64292 6.58583 4.27607 1 1 1 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/result.ref b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/result.ref new file mode 100644 index 0000000000..e15e090d95 --- /dev/null +++ b/tests/02_NAO_Gamma/relax_sccs_pcc0d_h3o/result.ref @@ -0,0 +1,7 @@ +etotref -478.6011205443832 +etotperatomref -119.6502801361 +totalforceref 1.982263 +esolelref -4.8068734558 +esolcavref 0.0335325829 +epccref 3.4923138796 +totaltimeref 3.93 diff --git a/tests/02_NAO_Gamma/scf_sccs_water_delay/INPUT b/tests/02_NAO_Gamma/scf_sccs_water_delay/INPUT new file mode 100644 index 0000000000..29e2f2a1aa --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water_delay/INPUT @@ -0,0 +1,30 @@ +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 +mixing_restart 1.0e-2 +mixing_dmr 1 +sccs_start_drho 1.0e-3 diff --git a/tests/02_NAO_Gamma/scf_sccs_water_delay/KPT b/tests/02_NAO_Gamma/scf_sccs_water_delay/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water_delay/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_delay/README b/tests/02_NAO_Gamma/scf_sccs_water_delay/README new file mode 100644 index 0000000000..f2588ea58d --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water_delay/README @@ -0,0 +1,3 @@ +LCAO (gamma-only) neutral water with periodic SCCS (water-neutral preset) and a delayed start, same input as +scf_sccs_water plus mixing_restart 1e-2, mixing_dmr 1 and sccs_start_drho 1e-3: the density-matrix mixing +starts at the mixing restart; when SCCS starts later, the mixing restarts, including the density-matrix history. diff --git a/tests/02_NAO_Gamma/scf_sccs_water_delay/STRU b/tests/02_NAO_Gamma/scf_sccs_water_delay/STRU new file mode 100644 index 0000000000..63a4455cda --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water_delay/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_delay/result.ref b/tests/02_NAO_Gamma/scf_sccs_water_delay/result.ref new file mode 100644 index 0000000000..6ed7b8f3c1 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_water_delay/result.ref @@ -0,0 +1,6 @@ +etotref -466.8169430490147 +etotperatomref -155.6056476830 +totalforceref 3.691279 +esolelref -0.4602106579 +esolcavref 0.0746967829 +totaltimeref 1.67 diff --git a/tests/integrate/validation_tools/props_basic.sh b/tests/integrate/validation_tools/props_basic.sh index ea6ce37248..0431108328 100644 --- a/tests/integrate/validation_tools/props_basic.sh +++ b/tests/integrate/validation_tools/props_basic.sh @@ -101,6 +101,21 @@ if ! test -z "$imp_sol" && { [ "$imp_sol" = 1 ] || [ "$imp_sol" = 2 ]; }; then echo "esolcavref $esol_cav" >>$props_result_file fi +# out_sol 1: both cubes must exist (the cases set imp_sol 2, or imp_sol 1 with +# eb_k > 1, so sol_cavity.cube is written). The cavity volume int s dV of +# sol_cavity.cube is compared; sol_eps.cube is only checked for presence, since +# its integral is dominated by the bulk permittivity times the volume. +if [[ "$out_sol" == "1" ]]; then + for cube in sol_eps.cube sol_cavity.cube; do + if [ ! -f "OUT.autotest/$cube" ]; then + echo "Error: out_sol 1 but OUT.autotest/$cube is missing" + exit 1 + fi + done + integral=$(python3 "$CUBE_TOOL" integrate "OUT.autotest/sol_cavity.cube") + echo "sol_cavity_integralref $integral" >>$props_result_file +fi + #-------------------------------------------- # vacuum point-countercharge correction #-------------------------------------------- diff --git a/tests/integrate/validation_tools/props_common.sh b/tests/integrate/validation_tools/props_common.sh index bca1717e6e..983d2f1b6d 100644 --- a/tests/integrate/validation_tools/props_common.sh +++ b/tests/integrate/validation_tools/props_common.sh @@ -127,6 +127,7 @@ props_init(){ out_mul=$(get_input_key_value "out_mul" "INPUT") gamma_only=$(get_input_key_value "gamma_only" "INPUT") imp_sol=$(get_input_key_value "imp_sol" "INPUT") + out_sol=$(get_input_key_value "out_sol" "INPUT") assume_isolated=$(get_input_key_value "assume_isolated" "INPUT") run_rpa=$(get_input_key_value "rpa" "INPUT") out_pot=$(get_input_key_value "out_pot" "INPUT")