diff --git a/source/Makefile.Objects b/source/Makefile.Objects index eb0b99bd9b3..3aa07144a0c 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -677,7 +677,8 @@ OBJS_SCCS=sccs_cavity.o\ sccs_coul_2d.o\ sccs_moments.o\ sccs_pw_coulomb.o\ - sccs_response.o + sccs_response.o\ + sccs_sa.o OBJS_SYMMETRY=symm_other.o\ symm_rot_spin.o\ diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp index d00f1ea9162..3ba439fd850 100644 --- a/source/source_estate/module_pot/pot_sccs.cpp +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -11,6 +11,7 @@ #include "source_hamilt/module_sccs/sccs_ion_charge.h" #include "source_hamilt/module_sccs/sccs_ion_force.h" #include "source_hamilt/module_sccs/sccs_response.h" +#include "source_hamilt/module_sccs/sccs_thread_sum.h" #include "source_hamilt/module_sccs/sccs_pw_coulomb.h" #include "source_io/module_parameter/input_parameter.h" @@ -53,6 +54,11 @@ void sccs_cfg_input(const Input_para& input, config.surface_regularization = input.sccs_surface_eta; config.cavity.lowpass_p1 = input.sccs_lowpass_p1; config.cavity.lowpass_p2 = input.sccs_lowpass_p2; + config.cavity.solvent_aware.solvent_radius = input.sccs_solvent_radius; + config.cavity.solvent_aware.radial_scale = input.sccs_radial_scale; + config.cavity.solvent_aware.radial_spread = input.sccs_radial_spread; + config.cavity.solvent_aware.filling_thr = input.sccs_filling_thr; + config.cavity.solvent_aware.filling_spread = input.sccs_filling_spread; if (input.assume_isolated == "pcc_0d") { config.boundary = ModuleSccs::Boundary::Pcc0d; @@ -132,6 +138,10 @@ PotSccs::PotSccs(const ModulePW::PW_Basis* basis, this->dynamic_mode = true; sccs_active_ = config_.start_drho <= 0.0 || resume.active; restart_potential_ = resume.restart_potential; + if (resume.probe_ready) + { + fixed_sources_.reuse_probe(resume.probe_kernel); + } } SccsResume PotSccs::resume_state() const @@ -139,6 +149,8 @@ SccsResume PotSccs::resume_state() const SccsResume resume; resume.active = sccs_active_; resume.restart_potential = restart_potential_; + resume.probe_ready = fixed_sources_.probe_ready(); + resume.probe_kernel = fixed_sources_.probe_kernel(); return resume; } @@ -178,14 +190,14 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: 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 - for (int ir = 0; ir < basis.nrxx; ++ir) - { + const auto add_point = [&](std::size_t ir, std::array& sums) { for (int spin = 0; spin < charge->nspin; ++spin) { density[ir] += charge->rho[spin][ir]; } solute_charge[ir] = ions[ir] - density[ir]; - grid_charges[0] += density[ir]; - grid_charges[1] += ions[ir]; - } + sums[0] += density[ir]; + sums[1] += ions[ir]; + }; + const std::array sums = ModuleSccs::thread_sums<2>(basis.nrxx, add_point); + double grid_charges[2] = {sums[0], sums[1]}; // electrons, ions // Pool-reduced, so every rank takes the same decision. Parallel_Reduce::reduce_pool(grid_charges, 2); const double dv = basis.omega / basis.nxyz; @@ -221,19 +233,24 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: if (config_.core_electrons) { const std::vector& core_density = fixed_sources_.core_density(); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ir = 0; ir < basis.nrxx; ++ir) { density[ir] += core_density[ir]; } } + const std::vector& probe_kernel = fixed_sources_.probe_kernel(); ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, - cell->tpiba, *coulomb, response); + ModuleSccs::solve_response(density, solute_charge, config_.cavity, probe_kernel, 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); + ModuleSccs::eval_functional(solute_charge, response, config_, probe_kernel, basis, cell->tpiba, *coulomb, + functional); output_.valid = true; output_.iterations = response.polarization.iterations; output_.warm_started = response.polarization.warm_started; @@ -242,6 +259,8 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: output_.reaction_energy = functional.reaction_energy; output_.volume = functional.volume; output_.surface = functional.surface; + output_.filled = response.solvent_aware.enabled(); + output_.filled_volume = response.solvent_aware.volume; output_.far_field_charge = response.far_field_polarization_charge; if (config_.collect_report && config_.boundary != ModuleSccs::Boundary::Periodic) { @@ -250,6 +269,9 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: electrostatic_rydberg_ = 2.0 * functional.reaction_energy; non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); electrostatic_potential_.resize(basis.nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ir = 0; ir < basis.nrxx; ++ir) { electrostatic_potential_[ir] = -2.0 * functional.reaction_potential[ir]; @@ -259,6 +281,8 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: cavity_potential_ = std::move(functional.cavity_potential); epsilon_ = std::move(response.epsilon); solute_ = std::move(response.solute); + local_solute_ = std::move(response.solvent_aware.local); + filled_fraction_ = std::move(response.solvent_aware.filling.fraction); restart_potential_ = std::move(response.restart_potential); ModuleBase::timer::end("PotSccs", "cal_v_eff"); } @@ -274,6 +298,9 @@ void PotSccs::add_sol_force(const UnitCell& cell, ModuleBase::matrix& force) con 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); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ir = 0; ir < basis.nrxx; ++ir) { reaction[ir] = -0.5 * electrostatic_potential_[ir]; @@ -349,5 +376,18 @@ void PotSccs::add_sol_fields(std::vector& fields) const cavity.name = "cavity"; cavity.values = solute_; fields.push_back(cavity); + // output_.filled comes from FilledCavity::enabled(), set on every rank. + if (!output_.filled) + { + return; + } + SolventGridField local; + local.name = "cavity_local"; + local.values = local_solute_; + fields.push_back(local); + SolventGridField fraction; + fraction.name = "filled_fraction"; + fraction.values = filled_fraction_; + fields.push_back(fraction); } } // namespace elecstate diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h index fd975c808f1..14912581f00 100644 --- a/source/source_estate/module_pot/pot_sccs.h +++ b/source/source_estate/module_pot/pot_sccs.h @@ -30,13 +30,16 @@ std::string check_sccs_cell(const ModuleSccs::SccsConfig& config, const UnitCell 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. +// rebuilt every ionic step: an activation already reached, the last +// unshifted sqrt-CG solution as a warm start, and the solvent-aware probe +// kernel of the fixed cell and grid. 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; + bool probe_ready = false; // same on every pool rank + std::vector probe_kernel; // local G vectors }; class PotSccs : public PotBase @@ -56,7 +59,9 @@ class PotSccs : public PotBase 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. + // eps and the cavity s of the last active evaluation, plus the local s + // (cavity_local) and the probe solute fraction c (filled_fraction) of a + // solvent-aware cavity; none before the start. void add_sol_fields(std::vector& fields) const override; // Delayed start: called once per electronic iteration while SCCS waits, @@ -82,6 +87,8 @@ class PotSccs : public PotBase std::vector cavity_potential_; std::vector epsilon_; std::vector solute_; + std::vector local_solute_; + std::vector filled_fraction_; }; } // namespace elecstate #endif diff --git a/source/source_estate/module_pot/sccs_fixed_src.cpp b/source/source_estate/module_pot/sccs_fixed_src.cpp index 567b295896f..f724f9c5311 100644 --- a/source/source_estate/module_pot/sccs_fixed_src.cpp +++ b/source/source_estate/module_pot/sccs_fixed_src.cpp @@ -6,6 +6,7 @@ #include "source_cell/unitcell.h" #include "source_hamilt/module_sccs/sccs_ion_charge.h" #include "source_hamilt/module_sccs/sccs_parameters.h" +#include "source_hamilt/module_sccs/sccs_sa.h" namespace { @@ -36,6 +37,9 @@ struct SccsFixedSources::Data std::vector ionic_density; std::vector core_density; std::vector> positions; + std::vector probe_kernel; + // Set on every rank; probe_kernel is collective to compute. + bool probe_kernel_ready = false; }; SccsFixedSources::SccsFixedSources() { @@ -51,6 +55,14 @@ bool SccsFixedSources::update(const UnitCell& cell, 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); + // The probe depends only on the cell and grid, which are fixed here. + const bool filled = ModuleSccs::uses_sa(config.cavity.solvent_aware); + if (filled && !data_->probe_kernel_ready) + { + data_->probe_kernel = ModuleSccs::sa_probe_kernel(basis, cell.latvec, cell.lat0, + config.cavity.solvent_aware); + data_->probe_kernel_ready = true; + } if (!reuse) { ModulePW::grid_positions(basis, cell.latvec, cell.lat0, data_->positions); @@ -69,4 +81,11 @@ bool SccsFixedSources::update(const UnitCell& cell, 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; } +const std::vector& SccsFixedSources::probe_kernel() const { return data_->probe_kernel; } +bool SccsFixedSources::probe_ready() const { return data_->probe_kernel_ready; } +void SccsFixedSources::reuse_probe(const std::vector& kernel) +{ + data_->probe_kernel = kernel; + data_->probe_kernel_ready = true; +} } // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_fixed_src.h b/source/source_estate/module_pot/sccs_fixed_src.h index c4ef8631a73..3d5e1d090df 100644 --- a/source/source_estate/module_pot/sccs_fixed_src.h +++ b/source/source_estate/module_pot/sccs_fixed_src.h @@ -11,10 +11,12 @@ namespace unitcell { struct AtomData; } namespace ModuleSccs { struct SccsConfig; } namespace elecstate { -// Owns ion/core Gaussian densities and Cartesian grid coordinates. The cell, +// Owns ion/core Gaussian densities, Cartesian grid coordinates and the +// solvent-aware probe kernel (empty without the filling). 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. +// the previous sources were reused. The probe depends only on the fixed cell +// and grid, so a kernel of the previous ionic step can be handed over. class SccsFixedSources { public: @@ -28,6 +30,12 @@ class SccsFixedSources const std::vector& ionic_density() const; const std::vector& core_density() const; const std::vector>& positions() const; + const std::vector& probe_kernel() const; + // Same value on every pool rank; the kernel itself is empty on a rank + // without local G vectors. + bool probe_ready() const; + // Called with the same ready flag on every pool rank, before update. + void reuse_probe(const std::vector& kernel); private: struct Data; std::unique_ptr data_; diff --git a/source/source_estate/module_pot/sccs_output.cpp b/source/source_estate/module_pot/sccs_output.cpp index ccd314c9fcd..acc38e02100 100644 --- a/source/source_estate/module_pot/sccs_output.cpp +++ b/source/source_estate/module_pot/sccs_output.cpp @@ -137,8 +137,14 @@ void write_sccs_out(std::ostream& output, const SccsOutput& result, int level, "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)}; + std::vector cavity = {"SCCS_CAVITY", "VOLUME/Bohr3", real_text(result.volume, value_digits), + "SURFACE/Bohr2", real_text(result.surface, value_digits)}; + if (result.filled) + { + const std::string filled_text = real_text(result.filled_volume, value_digits); + cavity.push_back("V_FILLED/Bohr3"); + cavity.push_back(filled_text); + } write_row(output, cavity); const std::vector transforms = {"SCCS_FFT", "R2G_CALLS", int_text(result.transforms.forward_calls), diff --git a/source/source_estate/module_pot/sccs_output.h b/source/source_estate/module_pot/sccs_output.h index 15253c2c3ba..5666b98ccfe 100644 --- a/source/source_estate/module_pot/sccs_output.h +++ b/source/source_estate/module_pot/sccs_output.h @@ -39,6 +39,8 @@ struct SccsOutput double reaction_energy = 0.0; double volume = 0.0; double surface = 0.0; + bool filled = false; // solvent-aware cavity + double filled_volume = 0.0; double far_field_charge = 0.0; double expected_charge = 0.0; ModuleSccs::CoulombTransformCounts transforms; diff --git a/source/source_estate/unittests/test_pot_sccs.cpp b/source/source_estate/unittests/test_pot_sccs.cpp index 1184ecab3df..7bb16f61984 100644 --- a/source/source_estate/unittests/test_pot_sccs.cpp +++ b/source/source_estate/unittests/test_pot_sccs.cpp @@ -214,6 +214,19 @@ 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_FALSE(ModuleSccs::uses_sa(config.cavity.solvent_aware)); + input.sccs_solvent_radius = 3.0; + input.sccs_radial_scale = 1.5; + input.sccs_radial_spread = 0.4; + input.sccs_filling_thr = 0.6; + input.sccs_filling_spread = 0.05; + elecstate::sccs_cfg_input(input, config, solver); + const ModuleSccs::SolventAwareParameters& filling = config.cavity.solvent_aware; + EXPECT_DOUBLE_EQ(filling.solvent_radius, 3.0); + EXPECT_DOUBLE_EQ(filling.radial_scale, 1.5); + EXPECT_DOUBLE_EQ(filling.radial_spread, 0.4); + EXPECT_DOUBLE_EQ(filling.filling_thr, 0.6); + EXPECT_DOUBLE_EQ(filling.filling_spread, 0.05); EXPECT_DOUBLE_EQ(config.start_drho, 0.0); EXPECT_FALSE(config.collect_report); input.sccs_start_drho = 1e-3; @@ -291,6 +304,56 @@ TEST_F(PotSccsTest, DelayedActivationHasOneIrreversibleTransitionAndResumes) EXPECT_TRUE(potential->empty()); } +// The probe kernel of the fixed cell and grid is handed to the next ionic +// step instead of being rebuilt: a seeded identity kernel (c = s) survives the +// evaluation unchanged. +TEST_F(PotSccsTest, ResumeCarriesTheSolventAwareProbeKernel) +{ + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 1.0 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + Input_para input; + input.sccs_epsilon = 5.0; + input.sccs_tol_rms = 1e-13; + input.sccs_tol_max = 1e-12; + input.sccs_solvent_radius = 0.5; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::sccs_cfg_input(input, config, solver); + elecstate::PotSccs first(&basis, config, solver, 1.0, no_resume); + EXPECT_FALSE(first.resume_state().probe_ready); + ModuleBase::matrix potential(1, basis.nrxx); + first.cal_v_eff(&charge, &cell, potential); + const elecstate::SccsResume resume = first.resume_state(); + EXPECT_TRUE(resume.probe_ready); + EXPECT_EQ(resume.probe_kernel.size(), static_cast(basis.npw)); + elecstate::PotSccs next(&basis, config, solver, 1.0, resume); + potential.zero_out(); + next.cal_v_eff(&charge, &cell, potential); + EXPECT_NEAR(next.get_energy(), first.get_energy(), 1e-10); + elecstate::SccsResume seeded = resume; + seeded.probe_kernel.assign(basis.npw, 1.0); + elecstate::PotSccs identity(&basis, config, solver, 1.0, seeded); + potential.zero_out(); + identity.cal_v_eff(&charge, &cell, potential); + const std::vector kept = identity.resume_state().probe_kernel; + EXPECT_EQ(kept, seeded.probe_kernel); +} + TEST_F(PotSccsTest, WaitingSccsContributesNoPotentialEnergyOrForce) { UnitCell cell; diff --git a/source/source_estate/unittests/test_sccs_output.cpp b/source/source_estate/unittests/test_sccs_output.cpp index bd2635c10c2..d46d9b450f2 100644 --- a/source/source_estate/unittests/test_sccs_output.cpp +++ b/source/source_estate/unittests/test_sccs_output.cpp @@ -229,12 +229,10 @@ TEST_F(SccsSolventFieldTest, SolventFieldsFollowTheLastActiveEvaluation) 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) + // One active evaluation, then the out_sol field list of this rank. + void evaluate_fields(const ModuleSccs::SccsConfig& config, + double& energy, + std::vector& fields) { double empty_ranks = basis.nrxx == 0 ? 1.0 : 0.0; Parallel_Reduce::reduce_pool(empty_ranks); @@ -272,9 +270,19 @@ class PotSccsEmptyPartitionTest : public SccsTest::EmptyPartitionTest ModuleBase::matrix potential(1, basis.nrxx); component.cal_v_eff(&charge, &cell, potential); energy = component.get_energy(); + fields.clear(); + component.add_sol_fields(fields); + } + // 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) + { std::vector fields; - component.add_sol_fields(fields); + evaluate_fields(config, energy, fields); const std::size_t field_count = fields.size(); names.clear(); sums.assign(field_count, 0.0); @@ -320,3 +328,44 @@ TEST_F(PotSccsEmptyPartitionTest, SolventOutputIsCollectiveOnEmptyGridPartitions EXPECT_NEAR(sums[0], 795.834560078087, 1e-8); EXPECT_NEAR(sums[1], 154.256465345260, 1e-8); } + +// The solvent-aware cavity adds cavity_local and filled_fraction on every rank, +// including one without grid points. The probe has unit integral, so the +// probe solute fraction has the same sum as the local boundary; the sums of +// squares tell the two fields apart. Serial reference values. +TEST_F(PotSccsEmptyPartitionTest, FilledSolventOutputIsCollectiveOnEmptyGridPartitions) +{ + ModuleSccs::SccsConfig config = active_config(); + config.cavity.solvent_aware.solvent_radius = 0.5; + config.cavity.solvent_aware.radial_spread = 0.3; + config.cavity.solvent_aware.filling_thr = 0.6; + config.cavity.solvent_aware.filling_spread = 0.1; + double energy = 0.0; + std::vector fields; + evaluate_fields(config, energy, fields); + const std::size_t field_count = fields.size(); + std::vector names; + std::vector sums(field_count, 0.0); + std::vector squares(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; + squares[field] += value * value; + } + } + Parallel_Reduce::reduce_pool(sums.data(), static_cast(field_count)); + Parallel_Reduce::reduce_pool(squares.data(), static_cast(field_count)); + const std::vector expected_names = {"eps", "cavity", "cavity_local", "filled_fraction"}; + ASSERT_EQ(names, expected_names); + EXPECT_NEAR(energy, -0.931004969822491, 1e-10); + EXPECT_NEAR(sums[0], 781.635153985345, 1e-8); + EXPECT_NEAR(sums[1], 160.406546243018, 1e-8); + EXPECT_NEAR(sums[2], 154.256465345260, 1e-8); + EXPECT_NEAR(sums[3], 154.256465345260, 1e-8); + EXPECT_NEAR(squares[2], 144.378502438537, 1e-8); + EXPECT_NEAR(squares[3], 133.256718776487, 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 e57a9c25ad1..f25c02b912d 100644 --- a/source/source_estate/unittests/test_sccs_pcc_0d.cpp +++ b/source/source_estate/unittests/test_sccs_pcc_0d.cpp @@ -54,10 +54,10 @@ TEST_F(SccsPcc0dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) } const std::vector cold; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, solute_charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, solute_charge, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, *coulomb, response); ModuleSccs::FunctionalResult functional; - ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, + ModuleSccs::eval_functional(solute_charge, response, config, SccsTest::no_probe, basis, tpiba, *coulomb, functional); elecstate::PotSccs solvent(&basis, config, solver, 0.6, no_resume); elecstate::PotPcc correction(&basis, 0.6); diff --git a/source/source_estate/unittests/test_sccs_pcc_2d.cpp b/source/source_estate/unittests/test_sccs_pcc_2d.cpp index 0960210f822..c7ceebcbbbc 100644 --- a/source/source_estate/unittests/test_sccs_pcc_2d.cpp +++ b/source/source_estate/unittests/test_sccs_pcc_2d.cpp @@ -57,10 +57,10 @@ TEST_F(SccsPcc2dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) } const std::vector cold; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, solute_charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, solute_charge, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, *coulomb, response); ModuleSccs::FunctionalResult functional; - ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, + ModuleSccs::eval_functional(solute_charge, response, config, SccsTest::no_probe, basis, tpiba, *coulomb, functional); elecstate::PotSccs solvent(&basis, config, solver, 0.6, no_resume); elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, input.pcc_2d_axis, 0.6); diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt index 1c476064308..d43b9b5c5ef 100644 --- a/source/source_hamilt/module_sccs/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -12,7 +12,8 @@ add_library( sccs_coul_2d.cpp sccs_moments.cpp sccs_pw_coulomb.cpp - sccs_response.cpp) + sccs_response.cpp + sccs_sa.cpp) if(ENABLE_COVERAGE) add_coverage(sccs) diff --git a/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp index 54c920bf8ce..d78631e63ff 100644 --- a/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp +++ b/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp @@ -2,8 +2,12 @@ #include "sccs_cavity.h" #include "sccs_lowpass.h" #include "sccs_response.h" +#include "sccs_sa.h" +#include "sccs_thread_sum.h" #include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_quit.h" #include "source_basis/module_pw/pw_basis.h" #include "source_hamilt/module_xc/xc_functional.h" @@ -14,9 +18,71 @@ namespace ModuleSccs { namespace { +// Gradient and Laplacian of the field with coefficients values_g. +void spec_derivs(const std::vector>& values_g, + const ModulePW::PW_Basis& basis, + double tpiba, + std::vector>& gradient, + std::vector& laplacian) +{ + gradient.resize(basis.nrxx); + laplacian.resize(basis.nrxx); + XC_Functional::grad_rho(values_g.data(), gradient.data(), &basis, tpiba); + XC_Functional::laplacian_rho(values_g.data(), laplacian.data(), &basis, tpiba); +} + +// Gradient and Laplacian on the FFT grid, multiplied by filter unless it is empty. +void fft_derivatives(const std::vector& values, + const std::vector& filter, + const ModulePW::PW_Basis& basis, + double tpiba, + std::vector>& gradient, + std::vector& laplacian) +{ + std::vector> values_g(basis.npw); + basis.real2recip(values.data(), values_g.data()); + if (!filter.empty()) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < basis.npw; ++ig) + { + values_g[ig] *= filter[ig]; + } + } + spec_derivs(values_g, basis, tpiba, gradient, laplacian); +} + +// FFT derivatives of the local s and of the probe solute fraction c = p * s from +// one transform of s: the probe is diagonal in G, so p * grad s = grad c. +void fft_fill_derivs(const std::vector& probe_kernel, + const ModulePW::PW_Basis& basis, + double tpiba, + const SccsResponse& response, + std::vector>& gradient, + std::vector& laplacian) +{ + std::vector> local_g(basis.npw); + const FilledCavity& filled = response.solvent_aware; + basis.real2recip(filled.local.data(), local_g.data()); + spec_derivs(local_g, basis, tpiba, gradient, laplacian); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < basis.npw; ++ig) + { + local_g[ig] *= probe_kernel[ig]; + } + std::vector> fraction_gradient; + std::vector fraction_laplacian; + spec_derivs(local_g, basis, tpiba, fraction_gradient, fraction_laplacian); + sa_chain_derivs(filled.local, filled.filling, fraction_gradient, fraction_laplacian, gradient, + laplacian); +} + // Environ eps_of_boundary: eps = exp(L (1 - s)), L = ln(eps_bulk), and -// the sqrt-CG coefficient eps (lapl ln eps / 2 + |grad ln eps|^2 / 4) / (4 pi) -// from grad s and lapl s. +// the sqrt-CG coefficient eps (lapl ln eps / 2 + |grad ln eps|^2 / 4) / (4 pi). void eps_of_boundary(const CavityParameters& cavity, const std::vector>& gradient, const std::vector& laplacian, @@ -26,6 +92,9 @@ void eps_of_boundary(const CavityParameters& cavity, const std::size_t size = response.solute.size(); const double log_bulk = std::log(cavity.epsilon_bulk); coefficient.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { const double log_epsilon = log_bulk * (1.0 - response.solute[i]); @@ -38,10 +107,82 @@ void eps_of_boundary(const CavityParameters& cavity, / ModuleBase::FOUR_PI; } } + +// Environ solvent_aware_boundary: keep s as the local boundary and make solute +// the filled s_sa. The chain grad n, grad c = p * s' grad n and s'' serve the +// chain surface with every boundary condition; density_g is the transform of +// the cavity density. +void fill_cavity(const std::vector& density, + const CavityParameters& cavity, + const std::vector& probe_kernel, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& response, + std::vector>& density_g) +{ + const std::size_t size = density.size(); + FilledCavity& filled = response.solvent_aware; + filled.active = true; + // Move s out and s_sa in without copies; filling.boundary is left empty. + filled.local.swap(response.solute); + filled.filling = sa_boundary(filled.local, probe_kernel, cavity.solvent_aware, basis); + response.solute.swap(filled.filling.boundary); + const std::vector& boundary = response.solute; + const std::vector& local = filled.local; + const auto add_filling = [&boundary, &local](std::size_t i, std::array& sum) { + sum[0] += boundary[i] - local[i]; + }; + double filled_volume = thread_sums<1>(size, add_filling)[0]; + Parallel_Reduce::reduce_pool(filled_volume); + filled.volume = filled_volume * basis.omega / basis.nxyz; + + density_g.resize(basis.npw); + basis.real2recip(density.data(), density_g.data()); + filled.density_gradient.resize(size); + XC_Functional::grad_rho(density_g.data(), filled.density_gradient.data(), &basis, tpiba); + std::vector> local_gradient(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + local_gradient[i] = filled.density_gradient[i] * response.dsolute_drho[i]; + } + filled.fraction_gradient = convolve_grad(probe_kernel, basis, local_gradient); +} + +// Environ deriv_method 'chain': grad s = s' grad n and +// lapl s = s' lapl n + s'' |grad n|^2 from the FFT derivatives of n. +void chain_derivs(const std::vector>& density_g, + const std::vector>& density_gradient, + const std::vector& dsolute_drho, + const std::vector& d2solute_drho2, + const ModulePW::PW_Basis& basis, + double tpiba, + std::vector>& gradient, + std::vector& laplacian) +{ + const std::size_t size = density_gradient.size(); + std::vector density_laplacian(size); + // These XC helpers use G in units of tpiba and ABACUS's normalized FFT. + XC_Functional::laplacian_rho(density_g.data(), density_laplacian.data(), &basis, tpiba); + gradient.resize(size); + laplacian.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + const double gradient_square = density_gradient[i] * density_gradient[i]; + gradient[i] = density_gradient[i] * dsolute_drho[i]; + laplacian[i] = dsolute_drho[i] * density_laplacian[i] + d2solute_drho2[i] * gradient_square; + } +} } // namespace void make_cav_derivs(const std::vector& density, const CavityParameters& cavity, + const std::vector& probe_kernel, const ModulePW::PW_Basis& basis, double tpiba, bool open_boundary, @@ -50,66 +191,94 @@ void make_cav_derivs(const std::vector& density, { derivatives.filter = switch_filter(cavity, basis); const bool lowpass = uses_lowpass(cavity); + // The configuration, not the local array sizes, decides the filling. + const bool filled = uses_sa(cavity.solvent_aware); + const std::size_t kernel_size = basis.npw; + if (filled && probe_kernel.size() != kernel_size) + { + ModuleBase::WARNING_QUIT("ModuleSccs::make_cav_derivs", + "solvent-aware SCCS needs the probe kernel on the local G vectors"); + } derivatives.gradient.clear(); std::vector& coefficient = derivatives.coefficient; const std::size_t size = density.size(); response.solute.resize(size); response.dsolute_drho.resize(size); response.epsilon.resize(size); - response.depsilon_drho.resize(size); response.grad_log_epsilon.resize(size); + response.solvent_aware = FilledCavity(); std::vector d2solute_drho2(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { const CavityPoint point = evaluate_cavity(density[i], cavity); response.solute[i] = point.solute; response.dsolute_drho[i] = point.dsolute_drho; response.epsilon[i] = point.epsilon; - response.depsilon_drho[i] = point.depsilon_drho; d2solute_drho2[i] = point.d2solute_drho2; } + std::vector> density_g; + if (filled) + { + response.solvent_aware.d2solute_drho2.swap(d2solute_drho2); + fill_cavity(density, cavity, probe_kernel, basis, tpiba, response, density_g); + } - std::vector> gradient(size); - std::vector laplacian(size); if (open_boundary) { - // Original PCC path: differentiate the switching boundary on the FFT - // grid, since its chain coefficient is not converged at the cavity edge. - std::vector> boundary_g(basis.npw); - basis.real2recip(response.solute.data(), boundary_g.data()); + // Open boundaries differentiate the switching boundary on the FFT grid, + // since its chain coefficient is not converged at the cavity edge. + // Unfiltered FFT derivatives of s_sa ring where the filling switches + // within a grid spacing and stall the SCF, so without lowpass the + // filling is applied analytically to those of the local s. With + // lowpass, s_sa is differentiated on the filtered grid (Environ fft + // with deriv_lowpass): the analytic filling would put unfiltered + // pointwise f'' and f''' into the exact cavity derivative. + std::vector> gradient; + std::vector laplacian; + if (filled && !lowpass) + { + fft_fill_derivs(probe_kernel, basis, tpiba, response, gradient, laplacian); + } + else + { + fft_derivatives(response.solute, derivatives.filter, basis, tpiba, gradient, laplacian); + } + eps_of_boundary(cavity, gradient, laplacian, response, coefficient); if (lowpass) { - for (int ig = 0; ig < basis.npw; ++ig) - { - boundary_g[ig] *= derivatives.filter[ig]; - } + derivatives.gradient.swap(gradient); } - XC_Functional::grad_rho(boundary_g.data(), gradient.data(), &basis, tpiba); - XC_Functional::laplacian_rho(boundary_g.data(), laplacian.data(), &basis, tpiba); + return; } - else + + // The filled cavity already holds the transform and gradient of n and s''. + std::vector> plain_gradient; + if (!filled) { - // Environ deriv_method 'chain': grad s = s' grad n and - // lapl s = s' lapl n + s'' |grad n|^2 from the FFT derivatives of n. - std::vector> density_g(basis.npw); + density_g.resize(basis.npw); basis.real2recip(density.data(), density_g.data()); - std::vector> density_gradient(size); - std::vector density_laplacian(size); - // These XC helpers use G in units of tpiba and ABACUS's normalized FFT. - XC_Functional::grad_rho(density_g.data(), density_gradient.data(), &basis, tpiba); - XC_Functional::laplacian_rho(density_g.data(), density_laplacian.data(), &basis, tpiba); - for (std::size_t i = 0; i < size; ++i) - { - const double gradient_square = density_gradient[i] * density_gradient[i]; - gradient[i] = density_gradient[i] * response.dsolute_drho[i]; - laplacian[i] = response.dsolute_drho[i] * density_laplacian[i] + d2solute_drho2[i] * gradient_square; - } + plain_gradient.resize(size); + XC_Functional::grad_rho(density_g.data(), plain_gradient.data(), &basis, tpiba); } - eps_of_boundary(cavity, gradient, laplacian, response, coefficient); - if (open_boundary && lowpass) + const FilledCavity& filling = response.solvent_aware; + const std::vector>& density_gradient + = filled ? filling.density_gradient : plain_gradient; + const std::vector& second_derivative = filled ? filling.d2solute_drho2 : d2solute_drho2; + std::vector> gradient; + std::vector laplacian; + chain_derivs(density_g, density_gradient, response.dsolute_drho, second_derivative, basis, tpiba, gradient, + laplacian); + if (filled) { - derivatives.gradient.swap(gradient); + // Environ deriv_method 'chain' through the filling. + const std::vector fraction_laplacian = convolve_probe(probe_kernel, basis, laplacian); + sa_chain_derivs(filling.local, filling.filling, filling.fraction_gradient, + fraction_laplacian, gradient, laplacian); } + eps_of_boundary(cavity, gradient, laplacian, response, coefficient); } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_cav_deriv.h b/source/source_hamilt/module_sccs/sccs_cav_deriv.h index d06fd8a8c90..e83dca05696 100644 --- a/source/source_hamilt/module_sccs/sccs_cav_deriv.h +++ b/source/source_hamilt/module_sccs/sccs_cav_deriv.h @@ -23,9 +23,14 @@ struct CavityDerivatives }; // Collective over the PW pool. Periodic cells use density chain derivatives; -// open boundaries differentiate solute on the FFT grid. +// open boundaries differentiate solute on the FFT grid. With the solvent-aware +// filling (positive cavity.solvent_aware.solvent_radius, with probe_kernel the +// sa_probe_kernel on the local G vectors), periodic cells carry the chain through the +// filling; open boundaries differentiate s_sa on the filtered grid with lowpass +// and otherwise apply the analytic filling to the FFT derivatives of the local s. void make_cav_derivs(const std::vector& density, const CavityParameters& cavity, + const std::vector& probe_kernel, const ModulePW::PW_Basis& basis, double tpiba, bool open_boundary, diff --git a/source/source_hamilt/module_sccs/sccs_cavity.cpp b/source/source_hamilt/module_sccs/sccs_cavity.cpp index 182ee672182..4af0fd0f579 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.cpp +++ b/source/source_hamilt/module_sccs/sccs_cavity.cpp @@ -10,6 +10,11 @@ bool uses_lowpass(const CavityParameters& parameters) return parameters.lowpass_p1 > 0.0 && parameters.lowpass_p2 > 0.0; } +bool uses_sa(const SolventAwareParameters& parameters) +{ + return parameters.solvent_radius > 0.0; +} + CavityPoint evaluate_cavity(double density, const CavityParameters& parameters) { CavityPoint point; diff --git a/source/source_hamilt/module_sccs/sccs_cavity.h b/source/source_hamilt/module_sccs/sccs_cavity.h index 9c36bd2a50e..1b50f31acbf 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.h +++ b/source/source_hamilt/module_sccs/sccs_cavity.h @@ -3,6 +3,18 @@ namespace ModuleSccs { +// Environ solvent-aware interface (Andreussi et al., JCTC 15, 1996 (2019)): +// cavity voids too small for a solvent molecule are filled with solute. +// Lengths are in Bohr; validated by the INPUT reader. +struct SolventAwareParameters +{ + double solvent_radius = 0.0; // zero disables the filling + double radial_scale = 2.0; // probe radius over solvent_radius + double radial_spread = 0.5; // erfc spread of the probe sphere + double filling_thr = 0.825; // probe solute fraction that fills a point + double filling_spread = 0.02; // erfc width of the filling step +}; + // Densities are in electrons/Bohr^3. The cavity uses the electronic density. struct CavityParameters { @@ -12,6 +24,7 @@ struct CavityParameters // Both positive enable the PCC switching-function derivative filter. double lowpass_p1 = -1.0; double lowpass_p2 = -1.0; + SolventAwareParameters solvent_aware; }; struct CavityPoint @@ -25,6 +38,8 @@ struct CavityPoint bool uses_lowpass(const CavityParameters& parameters); +bool uses_sa(const SolventAwareParameters& parameters); + // Parameters are validated by the INPUT reader. Negative Fourier ringing // uses the bulk limit. CavityPoint evaluate_cavity(double density, const CavityParameters& parameters); diff --git a/source/source_hamilt/module_sccs/sccs_coul_0d.cpp b/source/source_hamilt/module_sccs/sccs_coul_0d.cpp index 615c69c285a..a5095ffbbf3 100644 --- a/source/source_hamilt/module_sccs/sccs_coul_0d.cpp +++ b/source/source_hamilt/module_sccs/sccs_coul_0d.cpp @@ -24,6 +24,9 @@ void Pcc0dCoulombOperator::apply_potential(const std::vector& charge, st { const elecstate::ChargeMoments moments = pool_moments(charge, relative_positions_, basis_); periodic_.apply_potential(charge, potential); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ir = 0; ir < basis_.nrxx; ++ir) { potential[ir] += elecstate::pcc_0d_potential(moments, relative_positions_[ir], parameters_); diff --git a/source/source_hamilt/module_sccs/sccs_coul_2d.cpp b/source/source_hamilt/module_sccs/sccs_coul_2d.cpp index 64bd937523a..618ead7a6b1 100644 --- a/source/source_hamilt/module_sccs/sccs_coul_2d.cpp +++ b/source/source_hamilt/module_sccs/sccs_coul_2d.cpp @@ -26,6 +26,9 @@ void Pcc2dCoulombOperator::apply_potential(const std::vector& charge, st { const elecstate::ChargeMoments moments = pool_moments(charge, projected_positions_, basis_); periodic_.apply_potential(charge, potential); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ir = 0; ir < basis_.nrxx; ++ir) { const double coordinate = projected_positions_[ir] * normal_; diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp index a2b28a516b1..3e9ab0d3f32 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -1,6 +1,8 @@ #include "sccs_functional.h" #include "sccs_parameters.h" #include "sccs_response.h" +#include "sccs_sa.h" +#include "sccs_thread_sum.h" #include "sccs_pw_coulomb.h" #include "source_base/parallel_reduce.h" @@ -13,6 +15,74 @@ namespace ModuleSccs { +namespace +{ +// Spectral surface int(|grad s|_r - r), r = surface_regularization, and +// dE/ds = pressure - gamma div(grad s/|grad s|_r). +void spec_surface(const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + double& surface, + std::vector& boundary_potential) +{ + const std::size_t size = response.solute.size(); + const double dv = basis.omega / basis.nxyz; + std::vector> solute_g(basis.npw); + basis.real2recip(response.solute.data(), solute_g.data()); + std::vector> unit_gradient(size); + XC_Functional::grad_rho(solute_g.data(), unit_gradient.data(), &basis, tpiba); + const double eta = config.surface_regularization; + const auto add_point = [&](std::size_t i, std::array& sum) { + const double norm_squared = unit_gradient[i] * unit_gradient[i] + eta * eta; + const double norm = std::sqrt(norm_squared); + sum[0] += (norm - eta) * dv; + unit_gradient[i] /= norm; + }; + surface = thread_sums<1>(size, add_point)[0]; + Parallel_Reduce::reduce_pool(surface); + std::vector divergence(size); + XC_Functional::grad_dot(unit_gradient.data(), divergence.data(), &basis, tpiba); + boundary_potential.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + boundary_potential[i] = config.pressure - config.surface_tension * divergence[i]; + } +} + +// Surface of the filled boundary from the chain gradient of s_sa (Environ +// 'chain'); the spectral gradient of s_sa rings where the filling switches +// within a grid spacing, so every boundary condition uses it. density_potential +// is the exact derivative of the discrete surface and volume energies (see +// sa_nonelec), bounded where the gradient vanishes. +void filled_nonelec(const SccsResponse& response, + const SccsConfig& config, + const std::vector& probe_kernel, + const ModulePW::PW_Basis& basis, + double tpiba, + double& surface, + std::vector& density_potential) +{ + const double eta = config.surface_regularization; + SolventAwareNonelectrostatic filled = sa_nonelec( + response.solvent_aware, response.dsolute_drho, probe_kernel, basis, tpiba, eta, config.surface_tension, + config.pressure); + const double dv = basis.omega / basis.nxyz; + const std::vector>& gradient = filled.gradient; + const auto add_point = [&](std::size_t i, std::array& sum) { + const double norm_squared = gradient[i].norm2() + eta * eta; + const double norm = std::sqrt(norm_squared); + sum[0] += (norm - eta) * dv; + }; + surface = thread_sums<1>(gradient.size(), add_point)[0]; + Parallel_Reduce::reduce_pool(surface); + density_potential.swap(filled.density_potential); +} +} // namespace + void eval_functional(const std::vector& charge, const SccsResponse& response, const SccsConfig& config, @@ -21,12 +91,14 @@ void eval_functional(const std::vector& charge, FunctionalResult& result) { PeriodicCoulombOperator coulomb(basis, tpiba); - eval_functional(charge, response, config, basis, tpiba, coulomb, result); + const std::vector no_probe; + eval_functional(charge, response, config, no_probe, basis, tpiba, coulomb, result); } void eval_functional(const std::vector& charge, const SccsResponse& response, const SccsConfig& config, + const std::vector& probe_kernel, const ModulePW::PW_Basis& basis, double tpiba, CoulombOperator& coulomb, @@ -46,43 +118,41 @@ void eval_functional(const std::vector& charge, { candidate.cavity_potential = response.cavity_potential; } - for (std::size_t i = 0; i < size; ++i) - { + const auto add_point = [&](std::size_t i, std::array& sums) { const double reaction = response.polarization.potential[i] - vacuum[i]; candidate.reaction_potential[i] = reaction; - candidate.reaction_energy += 0.5 * charge[i] * reaction * dv; candidate.electron_potential[i] = -reaction + response.cavity_potential[i]; - candidate.volume += response.solute[i] * dv; - } + sums[0] += 0.5 * charge[i] * reaction * dv; + sums[1] += response.solute[i] * dv; + }; + const std::array sums = thread_sums<2>(size, add_point); + candidate.reaction_energy = sums[0]; + candidate.volume = sums[1]; Parallel_Reduce::reduce_pool(candidate.reaction_energy); Parallel_Reduce::reduce_pool(candidate.volume); candidate.volume_energy = config.pressure * candidate.volume; - // Spectral surface derivative: dE/dn = (p - gamma div(grad s/|grad s|)) s'. - std::vector> solute_g(basis.npw); - basis.real2recip(response.solute.data(), solute_g.data()); - std::vector> unit_gradient(size); - XC_Functional::grad_rho(solute_g.data(), unit_gradient.data(), &basis, tpiba); - const double eta = config.surface_regularization; - for (std::size_t i = 0; i < size; ++i) + std::vector nonel; + if (response.solvent_aware.enabled()) { - const double norm_squared = unit_gradient[i] * unit_gradient[i] + eta * eta; - const double norm = std::sqrt(norm_squared); - candidate.surface += (norm - eta) * dv; - unit_gradient[i] /= norm; + filled_nonelec(response, config, probe_kernel, basis, tpiba, candidate.surface, nonel); + } + else + { + std::vector boundary_potential; + spec_surface(response, config, basis, tpiba, candidate.surface, boundary_potential); + boundary_to_rho(response, probe_kernel, basis, boundary_potential, nonel); } - Parallel_Reduce::reduce_pool(candidate.surface); candidate.surface_energy = config.surface_tension * candidate.surface; - std::vector divergence(size); - XC_Functional::grad_dot(unit_gradient.data(), divergence.data(), &basis, tpiba); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { - const double nonel = (config.pressure - config.surface_tension * divergence[i]) - * response.dsolute_drho[i]; - candidate.electron_potential[i] += nonel; + candidate.electron_potential[i] += nonel[i]; if (core_electrons) { - candidate.cavity_potential[i] += nonel; + candidate.cavity_potential[i] += nonel[i]; } } result = std::move(candidate); diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h index 8912ec57f68..2464c80f2a2 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.h +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -27,7 +27,7 @@ struct FunctionalResult std::vector cavity_potential; // derivative w.r.t. cavity density, Hartree; full cavity only }; -// Collective over the PW pool. +// Collective over the PW pool; without the solvent-aware filling. void eval_functional(const std::vector& charge, const SccsResponse& response, const SccsConfig& config, @@ -36,9 +36,11 @@ void eval_functional(const std::vector& charge, FunctionalResult& result); // Use the same operator for the dielectric response and vacuum subtraction. +// probe_kernel is the one passed to solve_response for response. void eval_functional(const std::vector& charge, const SccsResponse& response, const SccsConfig& config, + const std::vector& probe_kernel, const ModulePW::PW_Basis& basis, double tpiba, CoulombOperator& coulomb, diff --git a/source/source_hamilt/module_sccs/sccs_lowpass.cpp b/source/source_hamilt/module_sccs/sccs_lowpass.cpp index 8f0758e318e..7b75953dd2b 100644 --- a/source/source_hamilt/module_sccs/sccs_lowpass.cpp +++ b/source/source_hamilt/module_sccs/sccs_lowpass.cpp @@ -19,6 +19,9 @@ std::vector switch_filter(const CavityParameters& cavity, const ModulePW if (lowpass) { filter.resize(basis.npw); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ig = 0; ig < basis.npw; ++ig) { // Environ deriv_lowpass: Gcut is the density cutoff, not ecutwfc. @@ -29,7 +32,7 @@ std::vector switch_filter(const CavityParameters& cavity, const ModulePW return filter; } -void lowpass_cav_pot(const std::vector& charge, +void lowpass_bnd_pot(const std::vector& charge, const CavityParameters& cavity, const CavityDerivatives& derivatives, const ModulePW::PW_Basis& basis, @@ -40,6 +43,9 @@ void lowpass_cav_pot(const std::vector& charge, const double log_bulk = std::log(cavity.epsilon_bulk); const std::vector& potential = response.polarization.potential; std::vector weight(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { weight[i] = response.epsilon[i] * potential[i] * potential[i] / (8.0 * ModuleBase::PI); @@ -49,6 +55,9 @@ void lowpass_cav_pot(const std::vector& charge, // The filtered derivative transposes share one inverse FFT. std::vector> transpose_g(basis.npw); basis.real2recip(weight.data(), transpose_g.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ig = 0; ig < basis.npw; ++ig) { transpose_g[ig] *= -tpiba * tpiba * basis.gg[ig] * derivatives.filter[ig]; @@ -57,11 +66,17 @@ void lowpass_cav_pot(const std::vector& charge, std::vector> component_g(basis.npw); for (int d = 0; d < 3; ++d) { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { component[i] = weight[i] * derivatives.gradient[i][d]; } basis.real2recip(component.data(), component_g.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ig = 0; ig < basis.npw; ++ig) { transpose_g[ig] += log_bulk * ModuleBase::IMAG_UNIT * tpiba * basis.gcar[ig][d] @@ -71,11 +86,13 @@ void lowpass_cav_pot(const std::vector& charge, std::vector transpose(size); basis.recip2real(transpose_g.data(), transpose.data()); std::vector candidate(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { - const double boundary_potential = 0.5 * log_bulk * (charge[i] * potential[i] + transpose[i]); - candidate[i] = boundary_potential * response.dsolute_drho[i]; + candidate[i] = 0.5 * log_bulk * (charge[i] * potential[i] + transpose[i]); } - response.cavity_potential = std::move(candidate); + response.boundary_potential = std::move(candidate); } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_lowpass.h b/source/source_hamilt/module_sccs/sccs_lowpass.h index 99cc5938ebc..821f66b6a32 100644 --- a/source/source_hamilt/module_sccs/sccs_lowpass.h +++ b/source/source_hamilt/module_sccs/sccs_lowpass.h @@ -19,9 +19,9 @@ struct SccsResponse; // erfc(p1 G^2 / Gcut^2 - p2) / 2, applied before every FFT derivative. std::vector switch_filter(const CavityParameters& cavity, const ModulePW::PW_Basis& basis); -// Exact derivative of the filtered, discrete sqrt-CG reaction energy through -// the electronic-density cavity; collective over the PW pool. -void lowpass_cav_pot(const std::vector& charge, +// Exact derivative of the filtered, discrete sqrt-CG reaction energy with +// respect to the boundary s; collective over the PW pool. +void lowpass_bnd_pot(const std::vector& charge, const CavityParameters& cavity, const CavityDerivatives& derivatives, const ModulePW::PW_Basis& basis, diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp index db2c4bdc359..914dc98214a 100644 --- a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp @@ -21,6 +21,9 @@ void PeriodicCoulombOperator::apply_potential(const std::vector& charge, basis_.real2recip(charge.data(), reciprocal_work_.data()); ++counts_.forward_calls; const double tpiba2 = tpiba_ * tpiba_; +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (int ig = 0; ig < basis_.npw; ++ig) { if (ig == basis_.ig_gge0 || basis_.gg[ig] == 0.0) diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index 6a94f5593bf..0339a7ccadf 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -3,6 +3,7 @@ #include "sccs_lowpass.h" #include "sccs_parameters.h" #include "sccs_pw_coulomb.h" +#include "sccs_thread_sum.h" #include "source_base/constants.h" #include "source_base/parallel_reduce.h" @@ -31,13 +32,13 @@ void finish_open_bc(const std::vector& charge, { const std::size_t size = charge.size(); const std::vector& potential = result.polarization.potential; - double source_sum = 0.0; - double solute_sum = 0.0; - for (std::size_t i = 0; i < size; ++i) - { - source_sum += (charge[i] - coefficient[i] * potential[i]) * invsqrt[i]; - solute_sum += charge[i]; - } + const auto add_point = [&](std::size_t i, std::array& sums) { + sums[0] += (charge[i] - coefficient[i] * potential[i]) * invsqrt[i]; + sums[1] += charge[i]; + }; + const std::array sums = thread_sums<2>(size, add_point); + double source_sum = sums[0]; + double solute_sum = sums[1]; Parallel_Reduce::reduce_pool(source_sum); Parallel_Reduce::reduce_pool(solute_sum); const double volume_element = basis.omega / basis.nxyz; @@ -52,14 +53,21 @@ void residual_norms(const std::vector& values, double& rms, double& maximum) { - double square = 0.0; - maximum = 0.0; - for (double value : values) + const auto add_square = [&values](std::size_t i, std::array& sum) { + sum[0] += values[i] * values[i]; + }; + double square = thread_sums<1>(values.size(), add_square)[0]; + double largest = 0.0; + const long size = values.size(); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) reduction(max : largest) +#endif + for (long i = 0; i < size; ++i) { - square += value * value; - const double magnitude = std::abs(value); - maximum = std::max(maximum, magnitude); + const double magnitude = std::abs(values[i]); + largest = std::max(largest, magnitude); } + maximum = largest; Parallel_Reduce::reduce_pool(square); Parallel_Reduce::reduce_max_pool(basis.poolnproc, maximum); const double mean_square = square / basis.nxyz; @@ -75,11 +83,10 @@ double grid_dot(const std::vector& left, const std::vector& right, const ModulePW::PW_Basis& basis) { - double value = 0.0; - for (std::size_t i = 0; i < left.size(); ++i) - { - value += left[i] * right[i]; - } + const auto add_product = [&left, &right](std::size_t i, std::array& sum) { + sum[0] += left[i] * right[i]; + }; + double value = thread_sums<1>(left.size(), add_product)[0]; Parallel_Reduce::reduce_pool(value); return value * basis.omega / basis.nxyz; } @@ -93,11 +100,17 @@ void apply_precond(const std::vector& rhs, { const std::size_t size = rhs.size(); weighted.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { weighted[i] = rhs[i] * invsqrt[i]; } coulomb.apply_potential(weighted, value); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { value[i] *= invsqrt[i]; @@ -134,11 +147,17 @@ void try_warm_start(const SqrtCgOperator& system, const std::vector& coefficient = system.coefficient; std::vector& z = state.preconditioned; std::vector guess_residual(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { guess_residual[i] = charge[i] - coefficient[i] * initial_potential[i]; } apply_precond(guess_residual, system.invsqrt, system.coulomb, state.weighted, z); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { guess_residual[i] = coefficient[i] * (initial_potential[i] - z[i]); @@ -183,6 +202,9 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, } const double beta = std::abs(old_rz) > 1e-30 ? rz / old_rz : 0.0; old_rz = rz; +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { direction[i] = z[i] + beta * direction[i]; @@ -194,6 +216,9 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, ModuleBase::WARNING_QUIT("ModuleSccs::solve_response", "SCCS sqrt-CG has invalid curvature"); } const double alpha = rz / curvature; +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { potential[i] += alpha * direction[i]; @@ -232,16 +257,18 @@ void finish_response(const SqrtCgOperator& system, // Preserve the physical gauge fixed by a boundary correction. if (!open_boundary) { - double mean = 0.0; - for (double value : potential) - { - mean += value; - } + const auto add_value = [&potential](std::size_t i, std::array& sum) { + sum[0] += potential[i]; + }; + double mean = thread_sums<1>(size, add_value)[0]; Parallel_Reduce::reduce_pool(mean); mean /= basis.nxyz; - for (double& value : potential) +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) { - value -= mean; + potential[i] -= mean; } } // The same condition under which make_cav_derivs keeps the filtered gradient. @@ -260,23 +287,52 @@ void finish_response(const SqrtCgOperator& system, } if (lowpass) { - lowpass_cav_pot(charge, cavity, derivatives, basis, tpiba, response); + lowpass_bnd_pot(charge, cavity, derivatives, basis, tpiba, response); } else { - // Continuum derivative, as in Environ: exact for the discrete energy only - // when |grad v|^2 is resolved by the FFT grid. - response.cavity_potential.resize(size); + // Continuum derivative -(d eps/ds) |grad v|^2/(8 pi), as in Environ: exact + // for the discrete energy only when |grad v|^2 is resolved by the FFT grid. + const double log_bulk = std::log(cavity.epsilon_bulk); + response.boundary_potential.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { const ModuleBase::Vector3& gradient = polarization.gradient[i]; const double gradient_square = gradient * gradient; - response.cavity_potential[i] = -response.depsilon_drho[i] * gradient_square / (8.0 * ModuleBase::PI); + response.boundary_potential[i] + = log_bulk * response.epsilon[i] * gradient_square / (8.0 * ModuleBase::PI); } } } } // namespace +void boundary_to_rho(const SccsResponse& response, + const std::vector& probe_kernel, + const ModulePW::PW_Basis& basis, + const std::vector& boundary_potential, + std::vector& density_potential) +{ + // Environ calculator: sa_de_dboundary, then de_dboundary * dscaled. + if (response.solvent_aware.enabled()) + { + density_potential = sa_adjoint(response.solvent_aware, probe_kernel, basis, boundary_potential); + } + else + { + density_potential = boundary_potential; + } +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < density_potential.size(); ++i) + { + density_potential[i] *= response.dsolute_drho[i]; + } +} + void solve_response(const std::vector& density, const std::vector& charge, const CavityParameters& cavity, @@ -287,12 +343,14 @@ void solve_response(const std::vector& density, SccsResponse& result) { PeriodicCoulombOperator coulomb(basis, tpiba); - solve_response(density, charge, cavity, solver, initial_potential, basis, tpiba, coulomb, result); + const std::vector no_probe; + solve_response(density, charge, cavity, no_probe, solver, initial_potential, basis, tpiba, coulomb, result); } void solve_response(const std::vector& density, const std::vector& charge, const CavityParameters& cavity, + const std::vector& probe_kernel, const PolarizationSolverParameters& solver, const std::vector& initial_potential, const ModulePW::PW_Basis& basis, @@ -308,10 +366,13 @@ void solve_response(const std::vector& density, const bool open_boundary = coulomb.has_correction(); SccsResponse candidate; CavityDerivatives derivatives; - make_cav_derivs(density, cavity, basis, tpiba, open_boundary, candidate, derivatives); + make_cav_derivs(density, cavity, probe_kernel, basis, tpiba, open_boundary, candidate, derivatives); const std::vector& coefficient = derivatives.coefficient; const std::size_t size = density.size(); std::vector invsqrt(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif for (std::size_t i = 0; i < size; ++i) { invsqrt[i] = 1.0 / std::sqrt(candidate.epsilon[i]); @@ -329,6 +390,8 @@ void solve_response(const std::vector& density, iterate_sqrt_cg(system, solver, state, polarization); std::vector& potential = state.potential; finish_response(system, charge, cavity, derivatives, tpiba, potential, candidate); + boundary_to_rho(candidate, probe_kernel, basis, candidate.boundary_potential, + candidate.cavity_potential); result = std::move(candidate); } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index eda88069a08..98c783a826f 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -2,6 +2,7 @@ #define SCCS_RESPONSE_H #include "sccs_cavity.h" +#include "sccs_sa.h" #include "source_base/vector3.h" #include @@ -29,21 +30,37 @@ struct PolarizationResult struct SccsResponse { + // Dielectric boundary; with the solvent-aware filling it is s_sa, while + // dsolute_drho stays ds/dn of the local boundary s. std::vector solute; std::vector dsolute_drho; + FilledCavity solvent_aware; // empty without the filling std::vector epsilon; - std::vector depsilon_drho; std::vector> grad_log_epsilon; PolarizationResult polarization; - // Hartree: discrete cavity derivative with lowpass, otherwise the continuum - // derivative -epsilon' |grad v|^2/(8 pi). + // Hartree: reaction-energy derivative with respect to the dielectric + // boundary (s_sa when filled); the discrete derivative with lowpass on an + // open boundary, otherwise the continuum derivative L eps |grad v|^2/(8 pi) + // with L = ln(eps_bulk). + std::vector boundary_potential; + // Hartree: the same derivative with respect to the cavity density. std::vector cavity_potential; std::vector restart_potential; // unshifted solution for explicit warm starts double far_field_polarization_charge = 0.0; // open-boundary screened charge minus solute charge }; +// Chain a derivative with respect to the boundary to the cavity density, +// through the solvent-aware adjoint when the boundary is filled; probe_kernel +// is the one the response was solved with. +void boundary_to_rho(const SccsResponse& response, + const std::vector& probe_kernel, + const ModulePW::PW_Basis& basis, + const std::vector& boundary_potential, + std::vector& density_potential); + // Periodic sqrt-CG response with chain-rule cavity derivatives, as Environ's -// CG solver with the sqrt preconditioner; charge is ions minus electrons. +// CG solver with the sqrt preconditioner, without the solvent-aware filling; +// charge is ions minus electrons. // All ranks in the PW pool must call together with the same configuration. // initial_potential is empty on all ranks for a cold start, otherwise a local grid. // A solver breakdown or missed tolerance stops the run with WARNING_QUIT. @@ -58,9 +75,13 @@ void solve_response(const std::vector& cavity_density, SccsResponse& result); // Use the same operator for the dielectric response and vacuum subtraction. +// cavity.solvent_aware decides on every pool rank whether the cavity is +// filled; probe_kernel is then its sa_probe_kernel, and is unused (may be +// empty) without the filling. void solve_response(const std::vector& cavity_density, const std::vector& solute_charge, const CavityParameters& cavity, + const std::vector& probe_kernel, const PolarizationSolverParameters& solver, const std::vector& initial_potential, const ModulePW::PW_Basis& basis, diff --git a/source/source_hamilt/module_sccs/sccs_sa.cpp b/source/source_hamilt/module_sccs/sccs_sa.cpp new file mode 100644 index 00000000000..bba3f396e30 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_sa.cpp @@ -0,0 +1,353 @@ +#include "sccs_sa.h" +#include "sccs_cavity.h" + +#include "source_base/constants.h" +#include "source_base/matrix3.h" +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_basis/module_pw/pw_basis.h" + +#include +#include + +namespace ModuleSccs +{ +namespace +{ +// Environ tools_math erfc threshold: dsfunct2 and d2sfunct2 vanish beyond it. +const double filling_argument_cutoff = 6.0; +// Environ function_erfc cuts the probe at width + 5 spread. +const double probe_cutoff_spreads = 5.0; + +// Fractional grid coordinate folded into [-1/2, 1/2). +double folded_fraction(int index, int count) +{ + double fraction = static_cast(index) / static_cast(count); + if (fraction >= 0.5) + { + fraction -= 1.0; + } + return fraction; +} + +// Periodic images n with |r0 + n a_k| <= cutoff for r0 in the folded cell: +// the cell height along a_k is the volume over the area of the other two. +int image_count(double cutoff, + const ModuleBase::Vector3& first, + const ModuleBase::Vector3& second, + double volume) +{ + const double reach = cutoff * (first ^ second).norm() / volume + 0.5; + return static_cast(std::ceil(reach)); +} +} // namespace + +std::vector sa_probe_kernel(const ModulePW::PW_Basis& basis, + const ModuleBase::Matrix3& lattice, + double lattice_scale, + const SolventAwareParameters& parameters) +{ + ModuleBase::timer::start("ModuleSccs", "sa_probe_kernel"); + const ModuleBase::Vector3 row1(lattice.e11, lattice.e12, lattice.e13); + const ModuleBase::Vector3 row2(lattice.e21, lattice.e22, lattice.e23); + const ModuleBase::Vector3 row3(lattice.e31, lattice.e32, lattice.e33); + const ModuleBase::Vector3 a1 = row1 * lattice_scale; + const ModuleBase::Vector3 a2 = row2 * lattice_scale; + const ModuleBase::Vector3 a3 = row3 * lattice_scale; + const double triple_product = a1 * (a2 ^ a3); + const double volume = std::abs(triple_product); + const double width = parameters.solvent_radius * parameters.radial_scale; + const double spread = parameters.radial_spread; + const double cutoff = width + probe_cutoff_spreads * spread; + const int count_1 = image_count(cutoff, a2, a3, volume); + const int count_2 = image_count(cutoff, a3, a1, volume); + const int count_3 = image_count(cutoff, a1, a2, volume); + + std::vector probe(basis.nrxx, 0.0); +#ifdef _OPENMP +#pragma omp parallel for schedule(dynamic, 64) +#endif + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const int iy = ir / basis.nplane - ix * basis.ny; + const int iz = ir % basis.nplane + basis.startz_current; + const ModuleBase::Vector3 folded = a1 * folded_fraction(ix, basis.nx) + + a2 * folded_fraction(iy, basis.ny) + + a3 * folded_fraction(iz, basis.nz); + double value = 0.0; + for (int n1 = -count_1; n1 <= count_1; ++n1) + { + for (int n2 = -count_2; n2 <= count_2; ++n2) + { + for (int n3 = -count_3; n3 <= count_3; ++n3) + { + const ModuleBase::Vector3 image = folded + a1 * static_cast(n1) + + a2 * static_cast(n2) + + a3 * static_cast(n3); + const double distance = image.norm(); + if (distance <= cutoff) + { + const double argument = (distance - width) / spread; + value += std::erfc(argument); + } + } + } + } + probe[ir] = value; + } + // The image sums dominate; a serial sum keeps the normalization independent + // of the thread count. + double integral = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + integral += probe[ir]; + } + Parallel_Reduce::reduce_pool(integral); + integral *= volume / static_cast(basis.nxyz); + // Environ convolution: fwfft(f1) fwfft(f2) omega with the 1/N forward + // transform, so the kernel is omega times the transform of the normalized + // probe. The probe is even on the grid, so its transform is real. + std::vector> probe_g(basis.npw); + basis.real2recip(probe.data(), probe_g.data()); + std::vector kernel(basis.npw); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < basis.npw; ++ig) + { + kernel[ig] = probe_g[ig].real() * volume / integral; + } + ModuleBase::timer::end("ModuleSccs", "sa_probe_kernel"); + return kernel; +} + +std::vector convolve_probe(const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector& values) +{ + std::vector> values_g(basis.npw); + basis.real2recip(values.data(), values_g.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < basis.npw; ++ig) + { + values_g[ig] *= kernel[ig]; + } + std::vector result(values.size()); + basis.recip2real(values_g.data(), result.data()); + return result; +} + +std::vector> convolve_grad( + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector>& gradient) +{ + const std::size_t size = gradient.size(); + std::vector component(size); + std::vector> result(size); + for (int d = 0; d < 3; ++d) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + component[i] = gradient[i][d]; + } + const std::vector convolved = convolve_probe(kernel, basis, component); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + result[i][d] = convolved[i]; + } + } + return result; +} + +SolventAwareBoundary sa_boundary(const std::vector& local, + const std::vector& kernel, + const SolventAwareParameters& parameters, + const ModulePW::PW_Basis& basis) +{ + SolventAwareBoundary filled; + filled.fraction = convolve_probe(kernel, basis, local); + const std::size_t size = local.size(); + const double threshold = parameters.filling_thr; + const double spread = parameters.filling_spread; + const double sqrt_pi = std::sqrt(ModuleBase::PI); + filled.filling.resize(size); + filled.dfilling.assign(size, 0.0); + filled.d2filling.assign(size, 0.0); + filled.boundary.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + // Environ: fill = 1 - sfunct2, dfill = -dsfunct2, d2fill = -d2sfunct2. + const double argument = (filled.fraction[i] - threshold) / spread; + filled.filling[i] = 1.0 - 0.5 * std::erfc(argument); + if (std::abs(argument) <= filling_argument_cutoff) + { + const double exponent = -argument * argument; + const double gaussian = std::exp(exponent) / sqrt_pi; + filled.dfilling[i] = gaussian / spread; + filled.d2filling[i] = -2.0 * argument * gaussian / (spread * spread); + } + filled.boundary[i] = local[i] + (1.0 - local[i]) * filled.filling[i]; + } + return filled; +} + +void sa_chain_derivs(const std::vector& local, + const SolventAwareBoundary& filled, + const std::vector>& fraction_gradient, + const std::vector& fraction_laplacian, + std::vector>& gradient, + std::vector& laplacian) +{ +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < local.size(); ++i) + { + const double empty = 1.0 - filled.filling[i]; + const double solvent = 1.0 - local[i]; + const double projection = gradient[i] * fraction_gradient[i]; + laplacian[i] = laplacian[i] * empty - 2.0 * projection * filled.dfilling[i] + + solvent * (filled.d2filling[i] * fraction_gradient[i].norm2() + + filled.dfilling[i] * fraction_laplacian[i]); + gradient[i] = gradient[i] * empty + fraction_gradient[i] * (solvent * filled.dfilling[i]); + } +} + +SolventAwareNonelectrostatic sa_nonelec(const FilledCavity& cavity, + const std::vector& dsolute_drho, + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + double tpiba, + double regularization, + double surface_tension, + double pressure) +{ + ModuleBase::timer::start("ModuleSccs", "sa_nonelec"); + const std::vector& local = cavity.local; + const SolventAwareBoundary& filled = cavity.filling; + const std::vector>& density_gradient = cavity.density_gradient; + const std::vector& d2solute_drho2 = cavity.d2solute_drho2; + const std::vector>& fraction_gradient = cavity.fraction_gradient; + const std::size_t size = local.size(); + const double regularization_square = regularization * regularization; + SolventAwareNonelectrostatic result; + result.gradient.resize(size); + std::vector> filled_part(size); + std::vector> fraction_part(size); + std::vector convolved(size); + std::vector pointwise(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + const ModuleBase::Vector3 local_gradient = density_gradient[i] * dsolute_drho[i]; + const double empty = 1.0 - filled.filling[i]; + const double solvent = 1.0 - local[i]; + const double fraction_weight = solvent * filled.dfilling[i]; + const ModuleBase::Vector3 g = local_gradient * empty + fraction_gradient[i] * fraction_weight; + result.gradient[i] = g; + const double norm_square = g.norm2() + regularization_square; + const ModuleBase::Vector3 unit = g / std::sqrt(norm_square); + const double local_projection = unit * local_gradient; + const double fraction_projection = unit * fraction_gradient[i]; + filled_part[i] = unit * empty; + fraction_part[i] = unit * fraction_weight; + // gamma alpha + p (1 - s) f' and gamma beta + p (1 - f): the convolved + // and pointwise coefficients of the variation of s. + const double alpha = -filled.dfilling[i] * local_projection + + solvent * filled.d2filling[i] * fraction_projection; + const double beta = -filled.dfilling[i] * fraction_projection; + convolved[i] = surface_tension * alpha + pressure * fraction_weight; + pointwise[i] = surface_tension * beta + pressure * empty; + } + // A = (1 - f) u + p * ((1 - s) f' u): the coefficient of the variation of a. + const std::vector> fraction_adjoint = convolve_grad(kernel, basis, + fraction_part); + const std::vector convolved_adjoint = convolve_probe(kernel, basis, convolved); + std::vector component(size); + std::vector> component_g(basis.npw); + std::vector> divergence_g(basis.npw, std::complex(0.0, 0.0)); + result.density_potential.resize(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + filled_part[i] += fraction_adjoint[i]; + const double curvature = d2solute_drho2[i] * (filled_part[i] * density_gradient[i]); + result.density_potential[i] = dsolute_drho[i] * (pointwise[i] + convolved_adjoint[i]) + + surface_tension * curvature; + } + // a = s' D n with the spectral gradient D, whose transpose is -div. + for (int d = 0; d < 3; ++d) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + component[i] = dsolute_drho[i] * filled_part[i][d]; + } + basis.real2recip(component.data(), component_g.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < basis.npw; ++ig) + { + divergence_g[ig] += ModuleBase::IMAG_UNIT * tpiba * basis.gcar[ig][d] * component_g[ig]; + } + } + std::vector divergence(size); + basis.recip2real(divergence_g.data(), divergence.data()); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + result.density_potential[i] -= surface_tension * divergence[i]; + } + ModuleBase::timer::end("ModuleSccs", "sa_nonelec"); + return result; +} + +std::vector sa_adjoint(const FilledCavity& cavity, + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector& boundary_potential) +{ + const std::vector& local = cavity.local; + const SolventAwareBoundary& filled = cavity.filling; + const std::size_t size = local.size(); + std::vector weighted(size); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + weighted[i] = (1.0 - local[i]) * filled.dfilling[i] * boundary_potential[i]; + } + std::vector result = convolve_probe(kernel, basis, weighted); +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + result[i] += (1.0 - filled.filling[i]) * boundary_potential[i]; + } + return result; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_sa.h b/source/source_hamilt/module_sccs/sccs_sa.h new file mode 100644 index 00000000000..dc9d0444435 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_sa.h @@ -0,0 +1,121 @@ +#ifndef SCCS_SA_H +#define SCCS_SA_H + +#include "source_base/vector3.h" + +#include + +namespace ModuleBase +{ +class Matrix3; +} +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct SolventAwareParameters; + +// Probe erfc((|r| - R)/radial_spread), R = solvent_radius * radial_scale, cut +// at R + 5 radial_spread, unit grid integral, on the local G vectors (Environ +// solvent_probe, but summing all periodic images within the cutoff, so narrow +// cells are not truncated). lattice in lattice_scale Bohr; collective. +std::vector sa_probe_kernel(const ModulePW::PW_Basis& basis, + const ModuleBase::Matrix3& lattice, + double lattice_scale, + const SolventAwareParameters& parameters); + +// Periodic convolution with the probe. +std::vector convolve_probe(const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector& values); + +// Component-wise convolution; of grad s it is grad c = p * grad s. +std::vector> convolve_grad( + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector>& gradient); + +// Filled boundary s_sa = s + (1 - s) f(c) of the probe solute fraction c = p * s, +// f(c) = 1 - erfc((c - filling_thr)/filling_spread)/2, with f' and f'' +// (Environ solvent_aware_boundary). +struct SolventAwareBoundary +{ + std::vector fraction; + std::vector filling; + std::vector dfilling; + std::vector d2filling; + std::vector boundary; +}; + +SolventAwareBoundary sa_boundary(const std::vector& local, + const std::vector& kernel, + const SolventAwareParameters& parameters, + const ModulePW::PW_Basis& basis); + +// Solvent-aware state of an SCCS cavity, empty without the filling: the +// local boundary s, its filling, the filled volume int(s_sa - s) in Bohr^3 +// and the chain data of the filled surface, s'' = d2s/dn2, grad n and the +// gradient of the probe solute fraction, grad c = p * (s' grad n). +// filling.boundary is empty: s_sa itself is SccsResponse::solute. +// active is set on every pool rank; the grid arrays are empty on a rank +// without real-space grid points. +struct FilledCavity +{ + bool enabled() const { return active; } + + bool active = false; + std::vector local; + SolventAwareBoundary filling; + double volume = 0.0; + std::vector d2solute_drho2; + std::vector> density_gradient; + std::vector> fraction_gradient; +}; + +// Environ solvent_aware_boundary, deriv_method 'chain': turn grad s and +// lapl s of the local boundary into those of s_sa, given the derivatives of +// the probe solute fraction, fraction_gradient = grad c = p * grad s and +// fraction_laplacian = lapl c = p * lapl s. +void sa_chain_derivs(const std::vector& local, + const SolventAwareBoundary& filled, + const std::vector>& fraction_gradient, + const std::vector& fraction_laplacian, + std::vector>& gradient, + std::vector& laplacian); + +// Filled surface and volume, a = s' grad n, b = p * a, g = (1 - f) a + (1 - s) f' b: +// S = sum_i (|g_i|_r - r), |g|_r = sqrt(|g|^2 + r^2); V = sum_i s_sa,i. +// With u = g/|g|_r, A = (1 - f) u + p * ((1 - s) f' u), alpha = -f' u.a + +// (1 - s) f'' u.b and beta = -f' u.b, the exact derivatives are +// dS/dn = s' (beta + p * alpha) + s'' A.grad n - div(s' A), +// dV/dn = s' ((1 - f) + p * ((1 - s) f')), +// bounded where g vanishes (|u| <= 1), unlike Environ's continuum curvature. +// density_potential = surface_tension dS/dn + pressure dV/dn; collective. +struct SolventAwareNonelectrostatic +{ + std::vector> gradient; + std::vector density_potential; +}; + +SolventAwareNonelectrostatic sa_nonelec(const FilledCavity& cavity, + const std::vector& dsolute_drho, + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + double tpiba, + double regularization, + double surface_tension, + double pressure); + +// Environ calc_solvent_aware_de_dboundary, from dE/ds_sa to dE/ds: +// (1 - f) dE/ds_sa + p * ((1 - s) f' dE/ds_sa). The probe is even, so the +// convolution is its own transpose and this is the exact adjoint. +std::vector sa_adjoint(const FilledCavity& cavity, + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector& boundary_potential); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_thread_sum.h b/source/source_hamilt/module_sccs/sccs_thread_sum.h new file mode 100644 index 00000000000..bcd6866d9e5 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_thread_sum.h @@ -0,0 +1,62 @@ +#ifndef SCCS_THREAD_SUM_H +#define SCCS_THREAD_SUM_H + +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +namespace ModuleSccs +{ +// Sums over the grid points [0, size): each thread runs accumulate(i, partial) +// over a fixed static share of the points in index order, starting from zero, +// and the partials are added in thread order. Unlike an OpenMP reduction, whose +// partials meet in run-dependent order, a given thread count gives the same +// sums on every run, and one thread gives the serial loop. The sqrt-CG +// iteration count depends on the last bits of its inner products. +template +std::array thread_sums(std::size_t size, const Accumulate& accumulate) +{ +#ifdef _OPENMP + const int thread_count = omp_get_max_threads(); +#else + const int thread_count = 1; +#endif + std::array zero; + zero.fill(0.0); + std::vector> partials(thread_count, zero); +#ifdef _OPENMP +#pragma omp parallel num_threads(thread_count) +#endif + { +#ifdef _OPENMP + const int thread = omp_get_thread_num(); +#else + const int thread = 0; +#endif + std::array partial = zero; +#ifdef _OPENMP +#pragma omp for schedule(static) +#endif + for (std::size_t i = 0; i < size; ++i) + { + accumulate(i, partial); + } + partials[thread] = partial; + } + std::array sums = zero; + for (int thread = 0; thread < thread_count; ++thread) + { + for (std::size_t k = 0; k < Count; ++k) + { + sums[k] += partials[thread][k]; + } + } + return sums; +} +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt index d3d56acc033..6bfb9313bca 100644 --- a/source/source_hamilt/module_sccs/test/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -17,6 +17,10 @@ addtest( test_sccs_ion_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + test_sccs_sa.cpp + test_sccs_sa_probe.cpp + sccs_sa_test.cpp + test_sccs_thread_sum.cpp ../../../source_estate/unittests/test_pot_sccs.cpp ../../../source_estate/unittests/test_sccs_pcc_0d.cpp ../../../source_estate/unittests/test_sccs_pcc_2d.cpp @@ -50,9 +54,33 @@ addtest( ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp + ../sccs_sa.cpp ../../module_xc/xc_grad_utils.cpp) -set_tests_properties(MODULE_HAMILT_sccs_poisson - PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60) +# All entries below run this binary in the same directory, where +# ChargedDielectricCellWarnsAndRuns writes and removes its warning log, so they +# share a lock and never run at the same time. +set_tests_properties( + MODULE_HAMILT_sccs_poisson PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT + 60 RESOURCE_LOCK sccs_poisson_workdir) +if(ENABLE_OPENMP) + # Four threads: the grid sums combine per-thread partial sums, and a fixed + # team must repeat a solve bit for bit. OMP_DYNAMIC=FALSE keeps the runtime + # from shrinking the team. + add_test( + NAME MODULE_HAMILT_sccs_poisson_omp4 + COMMAND $ + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties( + MODULE_HAMILT_sccs_poisson_omp4 + PROPERTIES ENVIRONMENT + "OMP_NUM_THREADS=4;OMP_DYNAMIC=FALSE" + TIMEOUT + 60 + PROCESSORS + 4 + RESOURCE_LOCK + sccs_poisson_workdir) +endif() if(ENABLE_MPI) add_test( @@ -60,8 +88,15 @@ if(ENABLE_MPI) COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} $ ${MPIEXEC_POSTFLAGS}) set_tests_properties( - MODULE_HAMILT_sccs_poisson_mpi PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" - TIMEOUT 60 PROCESSORS 2) + MODULE_HAMILT_sccs_poisson_mpi + PROPERTIES ENVIRONMENT + "OMP_NUM_THREADS=1" + TIMEOUT + 60 + PROCESSORS + 2 + RESOURCE_LOCK + sccs_poisson_workdir) # 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( @@ -72,5 +107,12 @@ if(ENABLE_MPI) --gtest_filter=*EmptyPartition*) set_tests_properties( MODULE_HAMILT_sccs_empty_partition_mpi - PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60 PROCESSORS 3) + PROPERTIES ENVIRONMENT + "OMP_NUM_THREADS=1" + TIMEOUT + 60 + PROCESSORS + 3 + RESOURCE_LOCK + sccs_poisson_workdir) endif() diff --git a/source/source_hamilt/module_sccs/test/sccs_sa_test.cpp b/source/source_hamilt/module_sccs/test/sccs_sa_test.cpp new file mode 100644 index 00000000000..c9b4e1dc679 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/sccs_sa_test.cpp @@ -0,0 +1,200 @@ +#include "sccs_sa_test.h" + +#include "source_base/parallel_reduce.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include + +namespace SccsTest +{ +void setup_basis(ModulePW::PW_Basis& basis, double scale, const ModuleBase::Matrix3& lattice, double ecut) +{ +#ifdef __MPI + basis.initmpi(pool_size, pool_rank, POOL_WORLD); +#endif + basis.initgrids(scale, lattice, ecut); + basis.initparameters(false, ecut, 1, false); + basis.setuptransform(); + basis.collect_local_pw(); +} + +void fold_into_cube(double length, ModuleBase::Vector3& r) +{ + for (int d = 0; d < 3; ++d) + { + const double fraction = r[d] / length; + r[d] -= length * std::round(fraction); + } +} + +ModuleBase::Vector3 cubic_displacement(const ModulePW::PW_Basis& basis, double length, int ir) +{ + const int ix = ir / (basis.ny * basis.nplane); + const int iy = ir / basis.nplane - ix * basis.ny; + const int iz = ir % basis.nplane + basis.startz_current; + const double x = length * ix / basis.nx; + const double y = length * iy / basis.ny; + const double z = length * iz / basis.nz; + ModuleBase::Vector3 r(x, y, z); + fold_into_cube(length, r); + return r; +} + +std::vector> spectral_gradient(const std::vector& values, + const ModulePW::PW_Basis& basis, + double tpiba) +{ + std::vector> values_g(basis.npw); + basis.real2recip(values.data(), values_g.data()); + std::vector> gradient(values.size()); + XC_Functional::grad_rho(values_g.data(), gradient.data(), &basis, tpiba); + return gradient; +} + +std::vector spectral_laplacian(const std::vector& values, + const ModulePW::PW_Basis& basis, + double tpiba) +{ + std::vector> values_g(basis.npw); + basis.real2recip(values.data(), values_g.data()); + std::vector laplacian(values.size()); + XC_Functional::laplacian_rho(values_g.data(), laplacian.data(), &basis, tpiba); + return laplacian; +} + +double pool_dot(const std::vector& left, const std::vector& right) +{ + double sum = 0.0; + for (std::size_t i = 0; i < left.size(); ++i) + { + sum += left[i] * right[i]; + } + Parallel_Reduce::reduce_pool(sum); + return sum; +} + +double pool_max(const ModulePW::PW_Basis& basis, double value) +{ + Parallel_Reduce::reduce_max_pool(basis.poolnproc, value); + return value; +} + +std::vector shifted(const std::vector& values, + const std::vector& direction, + double step) +{ + std::vector result(values.size()); + for (std::size_t i = 0; i < values.size(); ++i) + { + result[i] = values[i] + step * direction[i]; + } + return result; +} + +ModuleSccs::SolventAwareParameters probe_parameters(double solvent_radius) +{ + ModuleSccs::SolventAwareParameters parameters; + parameters.solvent_radius = solvent_radius; + return parameters; +} + +void SolventAwareTest::set_up_cube(double edge, double ecut, const ModuleSccs::SolventAwareParameters& probe) +{ + length = edge; + tpiba = ModuleBase::TWO_PI / edge; + parameters = probe; + setup_basis(basis, edge, cubic_lattice, ecut); + kernel = ModuleSccs::sa_probe_kernel(basis, cubic_lattice, edge, probe); +} + +std::vector SolventAwareTest::sphere(const ModuleBase::Vector3& center, double radius, double spread) const +{ + std::vector values(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + ModuleBase::Vector3 r = displacement(ir) - center; + fold_into_cube(length, r); + const double argument = (r.norm() - radius) / spread; + values[ir] = 0.5 * std::erfc(argument); + } + return values; +} + +std::vector SolventAwareTest::gaussian_pair(double offset) const +{ + std::vector density(basis.nrxx); + const ModuleBase::Vector3 first(offset, 0.0, 0.0); + const ModuleBase::Vector3 second(-offset, 0.0, 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const ModuleBase::Vector3 r = displacement(ir); + const double exponent_first = -(r - first).norm2(); + const double exponent_second = -(r - second).norm2(); + const double first_gaussian = std::exp(exponent_first); + const double second_gaussian = std::exp(exponent_second); + density[ir] = 4.8e-3 * (first_gaussian + second_gaussian); + } + return density; +} + +DensitySurface SolventAwareTest::surface_of_density(const std::vector& density, + const ModuleSccs::CavityParameters& cavity, + double regularization) const +{ + const std::size_t size = density.size(); + ModuleSccs::FilledCavity filled; + filled.active = true; + filled.local.resize(size); + filled.d2solute_drho2.resize(size); + std::vector dsolute(size); + for (std::size_t i = 0; i < size; ++i) + { + const ModuleSccs::CavityPoint point = ModuleSccs::evaluate_cavity(density[i], cavity); + filled.local[i] = point.solute; + dsolute[i] = point.dsolute_drho; + filled.d2solute_drho2[i] = point.d2solute_drho2; + } + filled.filling = fill(filled.local); + filled.density_gradient = spectral_gradient(density, basis, tpiba); + std::vector> local_gradient(size); + for (std::size_t i = 0; i < size; ++i) + { + local_gradient[i] = filled.density_gradient[i] * dsolute[i]; + } + filled.fraction_gradient = ModuleSccs::convolve_grad(kernel, basis, local_gradient); + // Unit surface tension and no pressure: density_potential is dS/dn. + const double surface_tension = 1.0; + const double pressure = 0.0; + DensitySurface result; + result.chain = ModuleSccs::sa_nonelec(filled, dsolute, kernel, basis, tpiba, regularization, + surface_tension, pressure); + for (std::size_t i = 0; i < size; ++i) + { + const double norm_square = result.chain.gradient[i].norm2() + regularization * regularization; + result.surface += std::sqrt(norm_square) - regularization; + } + Parallel_Reduce::reduce_pool(result.surface); + return result; +} + +double SolventAwareTest::surface_derivative_error(const std::vector& density, + const std::vector& direction, + const ModuleSccs::CavityParameters& cavity, + double regularization, + double step) const +{ + const DensitySurface center = surface_of_density(density, cavity, regularization); + const double backward = -step; + const std::vector plus = shifted(density, direction, step); + const std::vector minus = shifted(density, direction, backward); + const double surface_plus = surface_of_density(plus, cavity, regularization).surface; + const double surface_minus = surface_of_density(minus, cavity, regularization).surface; + const double finite_difference = (surface_plus - surface_minus) / (2.0 * step); + const double analytic = pool_dot(center.chain.density_potential, direction); + EXPECT_GT(std::abs(finite_difference), 1.0e-3); + const double miss = std::abs(analytic - finite_difference); + return miss / std::abs(finite_difference); +} + +} // namespace SccsTest diff --git a/source/source_hamilt/module_sccs/test/sccs_sa_test.h b/source/source_hamilt/module_sccs/test/sccs_sa_test.h new file mode 100644 index 00000000000..90a6366ae06 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/sccs_sa_test.h @@ -0,0 +1,93 @@ +#ifndef SCCS_SA_TEST_H +#define SCCS_SA_TEST_H + +#include "sccs_test.h" +#include "../sccs_cavity.h" +#include "../sccs_sa.h" + +#include + +// Grid helpers and the cubic-cell fixture of the solvent-aware tests +// (definitions in sccs_sa_test.cpp). +namespace SccsTest +{ +const ModuleBase::Matrix3 cubic_lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); +const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + +void setup_basis(ModulePW::PW_Basis& basis, double scale, const ModuleBase::Matrix3& lattice, double ecut); + +// Fold a displacement into the cube [-length/2, length/2). +void fold_into_cube(double length, ModuleBase::Vector3& r); + +// Minimum-image displacement of grid point ir from the origin in a cube. +ModuleBase::Vector3 cubic_displacement(const ModulePW::PW_Basis& basis, double length, int ir); + +std::vector> spectral_gradient(const std::vector& values, + const ModulePW::PW_Basis& basis, + double tpiba); + +std::vector spectral_laplacian(const std::vector& values, + const ModulePW::PW_Basis& basis, + double tpiba); + +double pool_dot(const std::vector& left, const std::vector& right); + +double pool_max(const ModulePW::PW_Basis& basis, double value); + +// values + step * direction, for central differences. +std::vector shifted(const std::vector& values, + const std::vector& direction, + double step); + +ModuleSccs::SolventAwareParameters probe_parameters(double solvent_radius); + +struct DensitySurface +{ + double surface = 0.0; + ModuleSccs::SolventAwareNonelectrostatic chain; +}; + +// A cubic cell and the probe kernel of the solvent-aware parameters. +class SolventAwareTest : public testing::Test +{ +protected: + SolventAwareTest() : basis("cpu", "double") {} + + void set_up_cube(double edge, double ecut, const ModuleSccs::SolventAwareParameters& probe); + + ModuleBase::Vector3 displacement(int ir) const { return cubic_displacement(basis, length, ir); } + + // Smooth sphere erfc((|r - center| - radius)/spread)/2. + std::vector sphere(const ModuleBase::Vector3& center, double radius, double spread) const; + + ModuleSccs::SolventAwareBoundary fill(const std::vector& local) const + { + return ModuleSccs::sa_boundary(local, kernel, parameters, basis); + } + + // Two Gaussians of density along x, centred at +-offset; their midpoint is a density saddle. + std::vector gaussian_pair(double offset) const; + + DensitySurface surface_of_density(const std::vector& density, + const ModuleSccs::CavityParameters& cavity, + double regularization) const; + + // Central difference of the discrete surface along direction against the + // projected analytic derivative; returns the relative error. + double surface_derivative_error(const std::vector& density, + const std::vector& direction, + const ModuleSccs::CavityParameters& cavity, + double regularization, + double step) const; + + double length = 0.0; + double tpiba = 0.0; + ModulePW::PW_Basis basis; + ModuleSccs::SolventAwareParameters parameters; + std::vector kernel; +}; +} // namespace SccsTest + +#endif diff --git a/source/source_hamilt/module_sccs/test/sccs_test.h b/source/source_hamilt/module_sccs/test/sccs_test.h index 75a535776ac..4d34938aa54 100644 --- a/source/source_hamilt/module_sccs/test/sccs_test.h +++ b/source/source_hamilt/module_sccs/test/sccs_test.h @@ -8,6 +8,7 @@ #include #include +#include namespace SccsTest { @@ -28,6 +29,8 @@ const double cavity_amplitude = 2e-5; extern int pool_size; extern int pool_rank; +// Probe kernel of a cavity without the solvent-aware filling. +const std::vector no_probe; class PwTest : public testing::Test { diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp index d23ef2d0374..abce0ffee2d 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp @@ -2,11 +2,51 @@ #include "../sccs_cav_deriv.h" #include "../sccs_cavity.h" #include "../sccs_response.h" +#include "../sccs_sa.h" + +#include "source_hamilt/module_xc/xc_functional.h" #include "source_base/parallel_reduce.h" using SccsCavityDerivativesTest = SccsTest::PwTest; +namespace +{ +std::vector> fft_gradient(const std::vector& values, + const std::vector& filter, + const ModulePW::PW_Basis& basis, + double tpiba) +{ + std::vector> values_g(basis.npw); + basis.real2recip(values.data(), values_g.data()); + for (int ig = 0; ig < basis.npw && !filter.empty(); ++ig) + { + values_g[ig] *= filter[ig]; + } + std::vector> gradient(values.size()); + XC_Functional::grad_rho(values_g.data(), gradient.data(), &basis, tpiba); + return gradient; +} + +// grad s_sa = (1 - f) grad s + (1 - s) f' p * grad s for a given grad s. +std::vector> filled_gradient(const ModuleSccs::SccsResponse& response, + const std::vector& kernel, + const ModulePW::PW_Basis& basis, + const std::vector>& local) +{ + const std::vector> fraction = ModuleSccs::convolve_grad(kernel, basis, + local); + std::vector> gradient(local.size()); + for (std::size_t i = 0; i < local.size(); ++i) + { + const double empty = 1.0 - response.solvent_aware.filling.filling[i]; + const double solvent = 1.0 - response.solvent_aware.local[i]; + gradient[i] = local[i] * empty + fraction[i] * (solvent * response.solvent_aware.filling.dfilling[i]); + } + return gradient; +} +} // namespace + TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFft) { ModuleSccs::CavityParameters cavity; @@ -42,7 +82,7 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf ModuleSccs::SccsResponse response; ModuleSccs::CavityDerivatives derivatives; const std::vector& coefficient = derivatives.coefficient; - ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, open_boundary, response, derivatives); const double log_bulk = std::log(cavity.epsilon_bulk); for (int ir = 0; ir < basis.nrxx; ++ir) @@ -68,7 +108,7 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFf ModuleSccs::SccsResponse periodic; ModuleSccs::CavityDerivatives periodic_derivatives; const std::vector& periodic_coefficient = periodic_derivatives.coefficient; - ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, periodic_boundary, + ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, periodic_boundary, periodic, periodic_derivatives); double difference = 0.0; for (int ir = 0; ir < basis.nrxx; ++ir) @@ -134,7 +174,7 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryLowpassFiltersSwitchingFunctionDer const bool open_boundary = true; ModuleSccs::SccsResponse response; ModuleSccs::CavityDerivatives derivatives; - ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, open_boundary, response, derivatives); const double log_bulk = std::log(cavity.epsilon_bulk); const double weight = 0.5 * std::erfc(0.5); @@ -151,3 +191,147 @@ TEST_F(SccsCavityDerivativesTest, OpenBoundaryLowpassFiltersSwitchingFunctionDer EXPECT_NEAR(derivatives.gradient[ir].x, gradient, 1e-12); } } + +// The filled dielectric: periodic cells chain through the density; open +// boundaries differentiate s_sa on the filtered grid with lowpass, otherwise +// fill the FFT gradient of the local s analytically. +TEST_F(SccsCavityDerivativesTest, FilledDielectricPathsFollowBoundaryAndLowpass) +{ + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + cavity.solvent_aware.solvent_radius = 1.5; + cavity.solvent_aware.filling_thr = 0.67; + cavity.solvent_aware.filling_spread = 0.1; + const ModuleBase::Matrix3 lattice; + const std::vector kernel = ModuleSccs::sa_probe_kernel(basis, lattice, length, cavity.solvent_aware); + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 * (1.0 + 0.3 * mode[ir]); + } + const double log_bulk = std::log(cavity.epsilon_bulk); + const std::vector no_filter; + for (int path = 0; path < 3; ++path) + { + const bool open_boundary = path > 0; + const bool lowpass = path == 2; + cavity.lowpass_p1 = lowpass ? basis.ggecut : -1.0; + cavity.lowpass_p2 = lowpass ? 0.5 : -1.0; + ModuleSccs::SccsResponse response; + ModuleSccs::CavityDerivatives derivatives; + ModuleSccs::make_cav_derivs(density, cavity, kernel, basis, tpiba, open_boundary, response, + derivatives); + double filled_volume = response.solvent_aware.volume; + EXPECT_GT(filled_volume, 50.0); + std::vector> expected; + if (lowpass) + { + expected = fft_gradient(response.solute, derivatives.filter, basis, tpiba); + } + else if (open_boundary) + { + const std::vector> local = fft_gradient(response.solvent_aware.local, no_filter, + basis, tpiba); + expected = filled_gradient(response, kernel, basis, local); + } + else + { + std::vector> local = fft_gradient(density, no_filter, basis, tpiba); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + local[ir] *= response.dsolute_drho[ir]; + } + expected = filled_gradient(response, kernel, basis, local); + } + double error = 0.0; + double scale = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const ModuleBase::Vector3 target = expected[ir] * (-log_bulk); + const ModuleBase::Vector3 difference = response.grad_log_epsilon[ir] - target; + const double target_norm = target.norm(); + error = std::max(error, difference.norm()); + scale = std::max(scale, target_norm); + const double epsilon_exponent = log_bulk * (1.0 - response.solute[ir]); + EXPECT_NEAR(response.epsilon[ir], std::exp(epsilon_exponent), 1e-12); + } + Parallel_Reduce::reduce_max_pool(basis.poolnproc, error); + Parallel_Reduce::reduce_max_pool(basis.poolnproc, scale); + EXPECT_GT(scale, 1e-2); + const double tolerance = 1e-12 * scale; + EXPECT_LT(error, tolerance); + EXPECT_EQ(derivatives.gradient.empty(), !lowpass); + } +} + +// The solvent radius alone switches the filling on, the same on every pool +// rank: the other solvent-aware parameters without a radius give the plain +// cavity on every path. +TEST_F(SccsCavityDerivativesTest, SolventAwareParametersWithoutARadiusLeaveTheCavityUnfilled) +{ + ModuleSccs::CavityParameters plain; + plain.density_min = 1e-4; + plain.density_max = 5e-3; + plain.epsilon_bulk = 5.0; + ModuleSccs::CavityParameters parameters_only = plain; + parameters_only.solvent_aware.radial_scale = 1.5; + parameters_only.solvent_aware.radial_spread = 0.3; + parameters_only.solvent_aware.filling_thr = 0.6; + parameters_only.solvent_aware.filling_spread = 0.05; + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 * (1.0 + 0.3 * mode[ir]); + } + for (int path = 0; path < 3; ++path) + { + const bool open_boundary = path > 0; + const bool lowpass = path == 2; + const double p1 = lowpass ? basis.ggecut : -1.0; + const double p2 = lowpass ? 0.5 : -1.0; + plain.lowpass_p1 = p1; + plain.lowpass_p2 = p2; + parameters_only.lowpass_p1 = p1; + parameters_only.lowpass_p2 = p2; + ModuleSccs::SccsResponse expected; + ModuleSccs::CavityDerivatives expected_derivatives; + ModuleSccs::make_cav_derivs(density, plain, SccsTest::no_probe, basis, tpiba, open_boundary, + expected, expected_derivatives); + ModuleSccs::SccsResponse response; + ModuleSccs::CavityDerivatives derivatives; + ModuleSccs::make_cav_derivs(density, parameters_only, SccsTest::no_probe, basis, tpiba, + open_boundary, response, derivatives); + EXPECT_FALSE(response.solvent_aware.enabled()); + EXPECT_EQ(response.solute, expected.solute); + EXPECT_EQ(derivatives.coefficient, expected_derivatives.coefficient); + } +} + +// A filled cavity without its probe kernel stops the run instead of reading +// past the kernel; a forked death test cannot share the pool collectives of a +// multi-rank run. +TEST_F(SccsCavityDerivativesTest, SolventAwareCavityWithoutAProbeKernelStops) +{ + if (SccsTest::pool_size != 1) + { + GTEST_SKIP(); + } + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + cavity.solvent_aware.solvent_radius = 1.0; + const std::vector density(basis.nrxx, 1e-3); + const bool open_boundary = false; + ModuleSccs::SccsResponse response; + ModuleSccs::CavityDerivatives derivatives; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, open_boundary, + response, derivatives), + ::testing::ExitedWithCode(1), ""); + testing::internal::GetCapturedStdout(); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp b/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp index 4acaf4e3a81..03cab6a1b77 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp @@ -113,14 +113,14 @@ TEST_F(SccsPcc0dCoulombTest, ChargedUniformDielectricPreservesGaugeAndVacuumSubt ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::Vacuum); config.cavity.epsilon_bulk = epsilon; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, charge, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, response); EXPECT_EQ(response.polarization.iterations, 1); const double solute_sum = integrate_product(charge, ones); const double screening = (1.0 / epsilon - 1.0) * solute_sum; EXPECT_NEAR(response.far_field_polarization_charge, screening, 1e-13); ModuleSccs::FunctionalResult result; - ModuleSccs::eval_functional(charge, response, config, basis, tpiba, + ModuleSccs::eval_functional(charge, response, config, SccsTest::no_probe, basis, tpiba, coulomb, result); for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -149,9 +149,9 @@ TEST_F(SccsPcc0dCoulombTest, ChargedUniformDielectricPreservesGaugeAndVacuumSubt } ModuleSccs::SccsResponse shifted_response; ModuleSccs::FunctionalResult shifted_result; - ModuleSccs::solve_response(density, shifted, config.cavity, solver, cold, + ModuleSccs::solve_response(density, shifted, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, shifted_response); - ModuleSccs::eval_functional(shifted, shifted_response, config, basis, tpiba, + ModuleSccs::eval_functional(shifted, shifted_response, config, SccsTest::no_probe, basis, tpiba, coulomb, shifted_result); energy[index++] = shifted_result.reaction_energy; } @@ -181,14 +181,14 @@ TEST_F(SccsPcc0dCoulombTest, NonuniformResponseSatisfiesFixedPointAndFarFieldCha solver.tolerance_max = 1e-12; const std::vector cold; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, charge, cavity, solver, cold, + ModuleSccs::solve_response(density, charge, cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, response); EXPECT_GT(response.polarization.iterations, 1); const bool open_boundary = true; ModuleSccs::SccsResponse derivatives; ModuleSccs::CavityDerivatives cavity_derivatives; const std::vector& coefficient = cavity_derivatives.coefficient; - ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, open_boundary, derivatives, cavity_derivatives); std::vector source(basis.nrxx); for (int ir = 0; ir < basis.nrxx; ++ir) @@ -209,7 +209,7 @@ TEST_F(SccsPcc0dCoulombTest, NonuniformResponseSatisfiesFixedPointAndFarFieldCha const double expected_screening = source_sum / std::sqrt(cavity.epsilon_bulk) - solute_sum; EXPECT_NEAR(response.far_field_polarization_charge, expected_screening, 1e-13); ModuleSccs::SccsResponse restarted; - ModuleSccs::solve_response(density, charge, cavity, solver, response.restart_potential, + ModuleSccs::solve_response(density, charge, cavity, SccsTest::no_probe, solver, response.restart_potential, basis, tpiba, coulomb, restarted); EXPECT_TRUE(restarted.polarization.warm_started); EXPECT_LT(restarted.polarization.iterations, response.polarization.iterations); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp b/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp index c4ced8cb89d..2ebcfddc97f 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp @@ -126,11 +126,11 @@ TEST_F(SccsPcc2dCoulombTest, UniformDielectricRetainsChargedGaugeAndReactionEner ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::Vacuum); config.cavity.epsilon_bulk = epsilon; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, charge, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, response); EXPECT_EQ(response.polarization.iterations, 1); ModuleSccs::FunctionalResult result; - ModuleSccs::eval_functional(charge, response, config, basis, tpiba, + ModuleSccs::eval_functional(charge, response, config, SccsTest::no_probe, basis, tpiba, coulomb, result); for (int ir = 0; ir < basis.nrxx; ++ir) { @@ -169,13 +169,13 @@ TEST_F(SccsPcc2dCoulombTest, NonuniformResponseSatisfiesFixedPointAndWarmStart) ModuleSccs::Pcc2dCoulombOperator coulomb(basis, tpiba, positions, normal, parameters); const std::vector cold; ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(density, charge, cavity, solver, cold, + ModuleSccs::solve_response(density, charge, cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, response); const bool open_boundary = true; ModuleSccs::SccsResponse derivatives; ModuleSccs::CavityDerivatives cavity_derivatives; const std::vector& coefficient = cavity_derivatives.coefficient; - ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + ModuleSccs::make_cav_derivs(density, cavity, SccsTest::no_probe, basis, tpiba, open_boundary, derivatives, cavity_derivatives); std::vector source(basis.nrxx); for (int ir = 0; ir < basis.nrxx; ++ir) @@ -191,7 +191,7 @@ TEST_F(SccsPcc2dCoulombTest, NonuniformResponseSatisfiesFixedPointAndWarmStart) EXPECT_NEAR(response.polarization.potential[ir], expected, 1e-10); } ModuleSccs::SccsResponse restarted; - ModuleSccs::solve_response(density, charge, cavity, solver, response.restart_potential, + ModuleSccs::solve_response(density, charge, cavity, SccsTest::no_probe, solver, response.restart_potential, basis, tpiba, coulomb, restarted); EXPECT_TRUE(restarted.polarization.warm_started); EXPECT_LT(restarted.polarization.iterations, response.polarization.iterations); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp index 1d335c5d0ec..8036b10c884 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -3,9 +3,82 @@ #include "../sccs_parameters.h" #include "../sccs_response.h" #include "../sccs_pw_coulomb.h" +#include "../sccs_sa.h" #include "source_base/parallel_reduce.h" -using SccsFunctionalTest = SccsTest::PwTest; +class SccsFunctionalTest : public SccsTest::PwTest +{ +protected: + // cos(2 pi x/L + 0.3) along x, whose flat planes miss the grid. + std::vector shifted_mode() const + { + std::vector values(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const double angle = ModuleBase::TWO_PI * ix / basis.nx + 0.3; + values[ir] = std::cos(angle); + } + return values; + } + + // Directional derivative of the non-electrostatic energy of a filled + // cavity in vacuum along a density cosine: electron-potential projection + // and central difference of the evaluated energy. + void filled_non_electrostatic_derivative(double surface_tension, + double pressure, + double& predicted, + double& finite_difference, + double& filled_volume) + { + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::Vacuum); + config.surface_tension = surface_tension; + config.pressure = pressure; + config.surface_regularization = 1e-8; + // s is about 0.67 +- 0.13; the probe solute fraction stays inside the filling step. + config.cavity.solvent_aware.solvent_radius = 1.5; + config.cavity.solvent_aware.filling_thr = 0.67; + config.cavity.solvent_aware.filling_spread = 0.1; + const ModuleBase::Matrix3 lattice; + const std::vector kernel = ModuleSccs::sa_probe_kernel(basis, lattice, length, + config.cavity.solvent_aware); + const std::vector mode = shifted_mode(); + const std::vector charge(basis.nrxx, 0.0); + const std::vector cold; + ModuleSccs::PolarizationSolverParameters solver; + ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); + auto evaluate = [&](double step, ModuleSccs::SccsResponse& response, ModuleSccs::FunctionalResult& result) { + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 * (1.0 + 0.3 * mode[ir]) + step * mode[ir]; + } + ModuleSccs::solve_response(density, charge, config.cavity, kernel, solver, cold, basis, tpiba, coulomb, + response); + ModuleSccs::eval_functional(charge, response, config, kernel, basis, tpiba, coulomb, result); + }; + ModuleSccs::SccsResponse center_response; + ModuleSccs::FunctionalResult center; + evaluate(0.0, center_response, center); + filled_volume = center_response.solvent_aware.volume; + predicted = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + predicted += center.electron_potential[ir] * mode[ir]; + } + Parallel_Reduce::reduce_pool(predicted); + predicted *= basis.omega / basis.nxyz; + const double step = 1e-7; + ModuleSccs::SccsResponse response; + ModuleSccs::FunctionalResult plus; + ModuleSccs::FunctionalResult minus; + evaluate(step, response, plus); + evaluate(-step, response, minus); + const double energy_plus = plus.surface_energy + plus.volume_energy; + const double energy_minus = minus.surface_energy + minus.volume_energy; + finite_difference = (energy_plus - energy_minus) / (2.0 * step); + } +}; TEST_F(SccsFunctionalTest, UniformDielectricEnergyAndElectronDerivative) { @@ -103,7 +176,7 @@ TEST_F(SccsFunctionalTest, ExplicitPeriodicOperatorPreservesResponseAndFunctiona basis, tpiba, legacy); ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); EXPECT_FALSE(coulomb.has_correction()); - ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, charge, config.cavity, SccsTest::no_probe, solver, cold, basis, tpiba, coulomb, explicit_response); EXPECT_EQ(legacy.polarization.iterations, explicit_response.polarization.iterations); EXPECT_EQ(legacy.polarization.potential, explicit_response.polarization.potential); @@ -113,7 +186,7 @@ TEST_F(SccsFunctionalTest, ExplicitPeriodicOperatorPreservesResponseAndFunctiona ModuleSccs::FunctionalResult explicit_result; ModuleSccs::eval_functional(charge, legacy, config, basis, tpiba, legacy_result); - ModuleSccs::eval_functional(charge, explicit_response, config, basis, tpiba, + ModuleSccs::eval_functional(charge, explicit_response, config, SccsTest::no_probe, basis, tpiba, coulomb, explicit_result); EXPECT_DOUBLE_EQ(legacy_result.reaction_energy, explicit_result.reaction_energy); EXPECT_DOUBLE_EQ(legacy_result.surface_energy, explicit_result.surface_energy); @@ -160,3 +233,22 @@ TEST_F(SccsFunctionalTest, CavityDensityDerivativeOnlyForTheFullCavity) EXPECT_NEAR(full.cavity_potential[ir], expected, 1e-14); } } + +// The filled surface and volume potential is the exact derivative of the +// discrete energies, with the pressure folded into the surface convolution. +TEST_F(SccsFunctionalTest, FilledNonElectrostaticPotentialDifferentiatesTheDiscreteEnergy) +{ + // Hartree/Bohr^2 and Hartree/Bohr^3, so that both terms enter the derivative. + const double surface_tension = 1e-2; + const double pressure = 1e-3; + double predicted = 0.0; + double finite_difference = 0.0; + double filled_volume = 0.0; + filled_non_electrostatic_derivative(surface_tension, pressure, predicted, finite_difference, filled_volume); + // The filling is switched on, and the derivative is far from zero, so the + // relative tolerance below tests it. + EXPECT_GT(filled_volume, 50.0); + EXPECT_GT(std::abs(finite_difference), 1e2); + const double tolerance = 1e-6 * std::abs(finite_difference); + EXPECT_NEAR(predicted, finite_difference, tolerance); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp index 0fb5ff96d3e..da16a73dd3d 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_ion_force.cpp @@ -4,6 +4,7 @@ #include "../sccs_response.h" #include "../sccs_parameters.h" #include "../sccs_functional.h" +#include "../sccs_sa.h" #include "source_cell/cell_tools.h" #include "source_basis/module_pw/pw_grid_geometry.h" @@ -195,9 +196,9 @@ TEST_F(SccsIonicForceTest, FullCavityPccLowpassTotalCorrectionForceFiniteDiffere ions[ir] -= electron_density; } ModuleSccs::SccsResponse response; - ModuleSccs::solve_response(cavity_density, ions, config.cavity, solver, cold, basis, tpiba, coulomb, - response); - ModuleSccs::eval_functional(ions, response, config, basis, tpiba, coulomb, functional); + ModuleSccs::solve_response(cavity_density, ions, config.cavity, SccsTest::no_probe, solver, cold, basis, + tpiba, coulomb, response); + ModuleSccs::eval_functional(ions, response, config, SccsTest::no_probe, basis, tpiba, coulomb, functional); }; ModuleSccs::FunctionalResult baseline; evaluate(atoms, baseline); @@ -228,3 +229,97 @@ TEST_F(SccsIonicForceTest, FullCavityPccLowpassTotalCorrectionForceFiniteDiffere } } } + +// The same force through the solvent-aware adjoint of a filled cavity. +TEST_F(SccsIonicForceTest, FilledFullCavityPccLowpassTotalCorrectionForceFiniteDifference) +{ + std::vector atoms(2); + atoms[0].valence_charge = 1.0; + atoms[0].position = ModuleBase::Vector3(3.1, 3.2, 4.3); + atoms[1].valence_charge = 1.0; + atoms[1].position = ModuleBase::Vector3(6.4, 5.1, 5.7); + const ModuleBase::Matrix3 lattice; + std::vector> positions; + ModulePW::grid_positions(basis, lattice, length, positions); + const double half = 0.5 * length; + const ModuleBase::Vector3 center(half, half, half); + for (ModuleBase::Vector3& position : positions) + { + position -= center; + } + elecstate::Pcc0dParameters parameters; + parameters.length = length; + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, parameters); + ModuleSccs::SccsConfig config; + config.boundary = ModuleSccs::Boundary::Pcc0d; + config.core_electrons = true; + config.core_spreads = {2.0, 0.0}; + config.cavity.epsilon_bulk = 5.0; + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + config.cavity.lowpass_p1 = basis.ggecut; + config.cavity.lowpass_p2 = 0.5; + config.surface_regularization = 1e-8; + config.surface_tension = 1e-5; + config.pressure = 1e-6; + // The uniform electrons alone give s = 0.92; the core Gaussian moves the + // probe solute fraction through the filling step. + config.cavity.solvent_aware.solvent_radius = 1.0; + config.cavity.solvent_aware.filling_thr = 0.93; + config.cavity.solvent_aware.filling_spread = 0.05; + const std::vector probe_kernel = ModuleSccs::sa_probe_kernel(basis, lattice, length, + config.cavity.solvent_aware); + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-13; + solver.tolerance_max = 1e-13; + const std::vector cold; + double filled_volume = 0.0; + auto evaluate = [&](const std::vector& displaced, + ModuleSccs::FunctionalResult& functional) { + std::vector ions; + std::vector core; + ModuleSccs::gauss_ion_rho(displaced, basis, tpiba, ModuleSccs::gaussian_ion_spread, ions); + ModuleSccs::gauss_core_rho(displaced, basis, tpiba, config.core_spreads, core); + std::vector cavity_density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double electron_density = 2.0 / basis.omega; + cavity_density[ir] = electron_density + core[ir]; + ions[ir] -= electron_density; + } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(cavity_density, ions, config.cavity, probe_kernel, solver, cold, basis, tpiba, + coulomb, response); + ModuleSccs::eval_functional(ions, response, config, probe_kernel, basis, tpiba, coulomb, functional); + filled_volume = response.solvent_aware.volume; + }; + ModuleSccs::FunctionalResult baseline; + evaluate(atoms, baseline); + EXPECT_GT(filled_volume, 1e-3); + std::vector> ionic; + std::vector> core; + ModuleSccs::gauss_ion_force(atoms, baseline.reaction_potential, basis, + tpiba, ModuleSccs::gaussian_ion_spread, ionic); + ModuleSccs::gauss_core_frc(atoms, baseline.cavity_potential, basis, + tpiba, config.core_spreads, core); + const double step = 1e-4; + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + std::vector plus = atoms; + std::vector minus = atoms; + plus[ia].position[axis] += step; + minus[ia].position[axis] -= step; + ModuleSccs::FunctionalResult positive; + ModuleSccs::FunctionalResult negative; + evaluate(plus, positive); + evaluate(minus, negative); + const double energy_plus = positive.reaction_energy + positive.surface_energy + positive.volume_energy; + const double energy_minus = negative.reaction_energy + negative.surface_energy + negative.volume_energy; + const double finite_difference = -(energy_plus - energy_minus) / (2.0 * step); + const double predicted = ionic[ia][axis] + core[ia][axis]; + EXPECT_NEAR(predicted, finite_difference, 1e-8); + } + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp index cc69a79a966..0b8c4cb57ed 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_lowpass.cpp @@ -6,6 +6,7 @@ #include "../sccs_coul_0d.h" #include "../sccs_coul_2d.h" #include "../sccs_response.h" +#include "../sccs_sa.h" #include "source_base/parallel_reduce.h" #include "source_basis/module_pw/pw_grid_geometry.h" @@ -15,6 +16,13 @@ class SccsLowpassTest : public SccsTest::PwTest { protected: void check_electron_derivative(ModuleSccs::CoulombOperator& coulomb, ModuleSccs::Boundary boundary) + { + const bool filled = false; + check_derivative(coulomb, boundary, filled); + } + + // A filled cavity puts the solvent-aware adjoint between dE/ds_sa and dE/dn. + void check_derivative(ModuleSccs::CoulombOperator& coulomb, ModuleSccs::Boundary boundary, bool filled) { ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); config.boundary = boundary; @@ -25,6 +33,16 @@ class SccsLowpassTest : public SccsTest::PwTest config.cavity.lowpass_p2 = 0.5; config.surface_tension = 0.0; config.pressure = 0.0; + std::vector probe_kernel; + if (filled) + { + // The probe solute fraction is about 0.67 everywhere, inside the filling step. + config.cavity.solvent_aware.solvent_radius = 1.5; + config.cavity.solvent_aware.filling_thr = 0.67; + config.cavity.solvent_aware.filling_spread = 0.05; + const ModuleBase::Matrix3 lattice; + probe_kernel = ModuleSccs::sa_probe_kernel(basis, lattice, length, config.cavity.solvent_aware); + } ModuleSccs::PolarizationSolverParameters solver; solver.tolerance_rms = 1e-13; solver.tolerance_max = 1e-13; @@ -40,11 +58,15 @@ class SccsLowpassTest : public SccsTest::PwTest } const std::vector cold; ModuleSccs::SccsResponse center_response; - ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + ModuleSccs::solve_response(density, charge, config.cavity, probe_kernel, solver, cold, basis, tpiba, coulomb, center_response); EXPECT_TRUE(center_response.polarization.gradient.empty()); + if (filled) + { + EXPECT_GT(center_response.solvent_aware.volume, 50.0); + } ModuleSccs::FunctionalResult center; - ModuleSccs::eval_functional(charge, center_response, config, basis, + ModuleSccs::eval_functional(charge, center_response, config, probe_kernel, basis, tpiba, coulomb, center); double predicted = 0.0; for (int ir = 0; ir < basis.nrxx; ++ir) @@ -70,9 +92,10 @@ class SccsLowpassTest : public SccsTest::PwTest } ModuleSccs::SccsResponse response; ModuleSccs::solve_response(displaced_density, displaced_charge, - config.cavity, solver, cold, basis, tpiba, coulomb, response); + config.cavity, probe_kernel, solver, cold, basis, tpiba, coulomb, + response); ModuleSccs::FunctionalResult functional; - ModuleSccs::eval_functional(displaced_charge, response, config, + ModuleSccs::eval_functional(displaced_charge, response, config, probe_kernel, basis, tpiba, coulomb, functional); energies[side] = functional.reaction_energy; } @@ -118,6 +141,17 @@ TEST_F(SccsLowpassTest, PccElectronPotentialDifferentiatesDiscreteEnergy) check_electron_derivative(slab_coulomb, ModuleSccs::Boundary::Pcc2d); } +// PCC 0D with the solvent-aware filling. +TEST_F(SccsLowpassTest, FilledPccElectronPotentialDifferentiatesDiscreteEnergy) +{ + elecstate::Pcc0dParameters molecule; + molecule.length = length; + const std::vector> positions = centered_positions(); + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, molecule); + const bool filled = true; + check_derivative(coulomb, ModuleSccs::Boundary::Pcc0d, filled); +} + TEST_F(SccsLowpassTest, FilterUsesDensityCutoff) { ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); 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 e8a5d52289e..d99f9cec74d 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 @@ -75,6 +75,14 @@ int main(int argc, char** argv) pool); #endif testing::InitGoogleTest(&argc, argv); + // A forked death-test child would inherit the OpenMP thread pool and + // deadlock in its first parallel region, so a single rank re-executes the + // binary instead. Under mpirun a re-executed child would join the MPI + // job; there the death tests that run OpenMP code are skipped. + if (SccsTest::pool_size == 1) + { + testing::GTEST_FLAG(death_test_style) = "threadsafe"; + } const int status = RUN_ALL_TESTS(); #ifdef __MPI Parallel_Global::finalize_mpi(); diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp index 91613008c91..21fbbe40d56 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -150,6 +150,55 @@ TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) EXPECT_GT(tight_rms.polarization.iterations, loose.polarization.iterations); } +// Two solves of the same problem with the same OpenMP team give the same bits: +// thread_sums adds the per-thread partial sums in a fixed order, and the +// sqrt-CG iteration count depends on the last bits of its inner products. +// MODULE_HAMILT_sccs_poisson_omp4 runs this with four threads. +TEST_F(SccsResponseTest, RepeatedSolveWithTheSameTeamIsBitIdentical) +{ + const std::vector mode_x = cosine_mode(0); + std::vector charge = cosine_mode(1); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 0.009 + 0.008 * mode_x[ir]; + charge[ir] *= 1e-3; + } + ModuleSccs::SccsConfig config; + config.cavity.density_min = 0.0024; + config.cavity.density_max = 0.0155; + config.cavity.epsilon_bulk = 78.3; + config.surface_tension = 1e-2; + config.pressure = 1e-3; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-11; + solver.tolerance_max = 1e-10; + const std::vector cold_start; + ModuleSccs::SccsResponse first; + ModuleSccs::SccsResponse second; + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, first); + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, second); + ModuleSccs::FunctionalResult first_energy; + ModuleSccs::FunctionalResult second_energy; + ModuleSccs::eval_functional(charge, first, config, basis, tpiba, first_energy); + ModuleSccs::eval_functional(charge, second, config, basis, tpiba, second_energy); + EXPECT_GT(first.polarization.iterations, 1); + EXPECT_EQ(first.polarization.iterations, second.polarization.iterations); + EXPECT_EQ(first.polarization.residual_rms, second.polarization.residual_rms); + double difference = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double change = first.polarization.potential[ir] - second.polarization.potential[ir]; + const double deviation = std::abs(change); + difference = std::max(difference, deviation); + } + Parallel_Reduce::reduce_max_pool(basis.poolnproc, difference); + EXPECT_EQ(difference, 0.0); + EXPECT_EQ(first_energy.reaction_energy, second_energy.reaction_energy); + EXPECT_EQ(first_energy.surface, second_energy.surface); + EXPECT_EQ(first_energy.volume, second_energy.volume); +} + namespace { // cos(2 pi harmonic x / L) on a basis' local real-space grid. diff --git a/source/source_hamilt/module_sccs/test/test_sccs_sa.cpp b/source/source_hamilt/module_sccs/test/test_sccs_sa.cpp new file mode 100644 index 00000000000..0d3ce78494b --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_sa.cpp @@ -0,0 +1,148 @@ +#include "sccs_sa_test.h" +#include "source_base/parallel_reduce.h" + +using SccsSolventAware = SccsTest::SolventAwareTest; +using SccsTest::DensitySurface; +using SccsTest::origin; +using SccsTest::pool_dot; +using SccsTest::pool_max; +using SccsTest::probe_parameters; +using SccsTest::shifted; +using SccsTest::spectral_gradient; +using SccsTest::spectral_laplacian; + +// E(s) = int g s_sa(s): the adjoint applied to g is dE/ds, including the +// nonlocal probe term that the pointwise (1 - f) g misses (by 27 percent). +TEST_F(SccsSolventAware, AdjointIsTheDerivativeOfTheFilledBoundary) +{ + ModuleSccs::SolventAwareParameters probe = probe_parameters(1.5); + probe.filling_thr = 0.45; + probe.filling_spread = 0.1; + set_up_cube(16.0, 20.0, probe); + const std::vector local = sphere(origin, 3.5, 0.6); + const ModuleBase::Vector3 center(2.5, 1.0, 0.0); + const std::vector direction = sphere(center, 2.0, 0.8); + std::vector weight(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const ModuleBase::Vector3 r = displacement(ir); + const double angle_x = ModuleBase::TWO_PI * r.x / length; + const double angle_y = ModuleBase::TWO_PI * r.y / length; + weight[ir] = 1.0 + 0.5 * std::cos(angle_x) + 0.3 * std::sin(angle_y); + } + ModuleSccs::FilledCavity cavity; + cavity.active = true; + cavity.local = local; + cavity.filling = fill(local); + const std::vector adjoint = ModuleSccs::sa_adjoint(cavity, kernel, basis, weight); + std::vector pointwise(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + pointwise[ir] = (1.0 - cavity.filling.filling[ir]) * weight[ir]; + } + const double step = 1.0e-5; + const double backward = -step; + const std::vector plus = shifted(local, direction, step); + const std::vector minus = shifted(local, direction, backward); + const ModuleSccs::SolventAwareBoundary filled_plus = fill(plus); + const ModuleSccs::SolventAwareBoundary filled_minus = fill(minus); + const double energy_plus = pool_dot(weight, filled_plus.boundary); + const double energy_minus = pool_dot(weight, filled_minus.boundary); + const double finite_difference = (energy_plus - energy_minus) / (2.0 * step); + const double analytic = pool_dot(adjoint, direction); + const double local_only = pool_dot(pointwise, direction); + const double scale = std::abs(finite_difference); + // Measured 9e-10 relative. + const double tolerance = 1.0e-8 * scale; + EXPECT_NEAR(analytic, finite_difference, tolerance); + const double local_miss = std::abs(local_only - finite_difference); + const double local_floor = 1.0e-2 * scale; + EXPECT_GT(local_miss, local_floor); +} + +// For a well-resolved boundary the Environ chain formulas agree with spectral +// derivatives of s_sa up to the grid error of the composition. The filling +// switches over filling_spread / |grad c|, which must span several grid +// points; with Environ's filling_spread 0.02 it usually does not. +TEST_F(SccsSolventAware, ChainDerivativesMatchSpectralDerivativesOfTheFilledBoundary) +{ + ModuleSccs::SolventAwareParameters probe = probe_parameters(1.5); + probe.filling_thr = 0.45; + probe.filling_spread = 0.3; + set_up_cube(16.0, 80.0, probe); + const std::vector local = sphere(origin, 3.5, 1.5); + const ModuleSccs::SolventAwareBoundary filled = fill(local); + std::vector> gradient = spectral_gradient(local, basis, tpiba); + std::vector laplacian = spectral_laplacian(local, basis, tpiba); + const std::vector> fraction_gradient + = ModuleSccs::convolve_grad(kernel, basis, gradient); + const std::vector fraction_laplacian = ModuleSccs::convolve_probe(kernel, basis, laplacian); + ModuleSccs::sa_chain_derivs(local, filled, fraction_gradient, fraction_laplacian, gradient, + laplacian); + + const std::vector> expected_gradient = spectral_gradient(filled.boundary, basis, + tpiba); + const std::vector expected_laplacian = spectral_laplacian(filled.boundary, basis, tpiba); + double gradient_error = 0.0; + double gradient_scale = 0.0; + double laplacian_error = 0.0; + double laplacian_scale = 0.0; + double filling_change = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const ModuleBase::Vector3 gradient_difference = gradient[ir] - expected_gradient[ir]; + const double gradient_norm = expected_gradient[ir].norm(); + const double laplacian_difference = std::abs(laplacian[ir] - expected_laplacian[ir]); + const double laplacian_norm = std::abs(expected_laplacian[ir]); + const double change = filled.boundary[ir] - local[ir]; + gradient_error = std::max(gradient_error, gradient_difference.norm()); + gradient_scale = std::max(gradient_scale, gradient_norm); + laplacian_error = std::max(laplacian_error, laplacian_difference); + laplacian_scale = std::max(laplacian_scale, laplacian_norm); + filling_change = std::max(filling_change, change); + } + EXPECT_GT(pool_max(basis, filling_change), 0.05); + // Measured 8e-6 and 1e-4 of the maxima. + const double gradient_tolerance = 5.0e-5 * pool_max(basis, gradient_scale); + const double laplacian_tolerance = 5.0e-4 * pool_max(basis, laplacian_scale); + EXPECT_LT(pool_max(basis, gradient_error), gradient_tolerance); + EXPECT_LT(pool_max(basis, laplacian_error), laplacian_tolerance); +} + +// A density saddle on a grid node inside the switching range: the continuum +// curvature gives -tr(H)/r there (about 2e8 at r = 1e-8); the discrete +// derivative stays of the size of its neighbours and remains exact. +TEST_F(SccsSolventAware, ChainSurfaceDerivativeStaysBoundedAtADensitySaddle) +{ + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 1.0; + cavity.solvent_aware = probe_parameters(0.5); + cavity.solvent_aware.filling_thr = 0.99; + set_up_cube(12.0, 60.0, cavity.solvent_aware); + const std::vector density = gaussian_pair(1.5); + const double regularization = 1.0e-8; + const DensitySurface surface = surface_of_density(density, cavity, regularization); + double saddle = 0.0; + double largest = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double value = std::abs(surface.chain.density_potential[ir]); + if (displacement(ir).norm() < 1.0e-12) + { + saddle = value; + } + else + { + largest = std::max(largest, value); + } + } + Parallel_Reduce::reduce_pool(saddle); + largest = pool_max(basis, largest); + EXPECT_LT(saddle, largest); + const ModuleBase::Vector3 center(0.3, 0.2, 0.1); + const std::vector direction = sphere(center, 1.0, 0.5); + const double error = surface_derivative_error(density, direction, cavity, regularization, 1.0e-7); + EXPECT_LT(error, 5.0e-6); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_sa_probe.cpp b/source/source_hamilt/module_sccs/test/test_sccs_sa_probe.cpp new file mode 100644 index 00000000000..9f4aa57d77a --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_sa_probe.cpp @@ -0,0 +1,213 @@ +#include "sccs_sa_test.h" +#include "source_base/parallel_reduce.h" + +#include + +using SccsSolventAware = SccsTest::SolventAwareTest; +using SccsTest::origin; +using SccsTest::pool_max; +using SccsTest::probe_parameters; +using SccsTest::setup_basis; + +namespace +{ +// Normalized transform of the isolated probe, 4 pi int r^2 p(r) j0(G r) dr +// over the same cut at width + 5 spread, divided by its G = 0 value. +double isolated_probe_transform(double g, double width, double spread) +{ + const double cutoff = width + 5.0 * spread; + const int intervals = 40000; + const double step = cutoff / intervals; + double value = 0.0; + double norm = 0.0; + for (int i = 0; i <= intervals; ++i) + { + const double r = i * step; + const double weight = (i == 0 || i == intervals) ? 1.0 : (i % 2 == 1 ? 4.0 : 2.0); + const double argument = (r - width) / spread; + const double radial = r * r * std::erfc(argument); + const double phase = g * r; + const double bessel = phase > 1.0e-12 ? std::sin(phase) / phase : 1.0; + value += weight * radial * bessel; + norm += weight * radial; + } + return value / norm; +} + +// Kernel entry of G = (fx, fy, fz) in units of 2 pi/a, from whichever rank holds it. +double kernel_entry(const ModulePW::PW_Basis& basis, const std::vector& kernel, double fx) +{ + double values[2] = {0.0, 0.0}; + for (int ig = 0; ig < basis.npw; ++ig) + { + const ModuleBase::Vector3& g = basis.gdirect[ig]; + if (std::abs(g.x - fx) < 1.0e-12 && std::abs(g.y) < 1.0e-12 && std::abs(g.z) < 1.0e-12) + { + values[0] += kernel[ig]; + values[1] += 1.0; + } + } + Parallel_Reduce::reduce_pool(values, 2); + EXPECT_EQ(values[1], 1.0); + return values[0]; +} + +// Largest deviation of the kernel from the isolated-probe transform over the +// G vectors with |G| <= maximum_g (Bohr^-1). +double oblique_probe_error(const ModuleBase::Matrix3& lattice, + double scale, + const ModuleSccs::SolventAwareParameters& parameters, + double maximum_g, + double& compared) +{ + ModulePW::PW_Basis basis("cpu", "double"); + setup_basis(basis, scale, lattice, 120.0); + const std::vector kernel = ModuleSccs::sa_probe_kernel(basis, lattice, scale, parameters); + const double width = parameters.solvent_radius * parameters.radial_scale; + const double tpiba = ModuleBase::TWO_PI / scale; + double error = 0.0; + compared = 0.0; + for (int ig = 0; ig < basis.npw; ++ig) + { + const double g = tpiba * std::sqrt(basis.gg[ig]); + if (g > maximum_g) + { + continue; + } + const double expected = isolated_probe_transform(g, width, parameters.radial_spread); + const double deviation = std::abs(kernel[ig] - expected); + error = std::max(error, deviation); + compared += 1.0; + } + Parallel_Reduce::reduce_pool(compared); + return pool_max(basis, error); +} +} // namespace + +// In a cell wider than the probe the kernel is the transform of Environ's +// minimum-image probe, normalized on the grid, and it keeps constants. +TEST_F(SccsSolventAware, WideCellProbeMatchesEnvironMinimumImageProbe) +{ + // At 120 Ry the grid error of the sampled erfc edge is below 1e-10. + const ModuleSccs::SolventAwareParameters probe = probe_parameters(3.0); + set_up_cube(24.0, 120.0, probe); + const double width = parameters.solvent_radius * parameters.radial_scale; + const double spread = parameters.radial_spread; + // Environ function_erfc cuts the probe at width + 5 spread. + const double environ_cutoff = width + 5.0 * spread; + std::vector environ_probe(basis.nrxx, 0.0); + double integral = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double distance = displacement(ir).norm(); + if (distance <= environ_cutoff) + { + const double argument = (distance - width) / spread; + environ_probe[ir] = std::erfc(argument); + } + integral += environ_probe[ir]; + } + Parallel_Reduce::reduce_pool(integral); + const double volume = length * length * length; + integral *= volume / basis.nxyz; + std::vector> environ_g(basis.npw); + basis.real2recip(environ_probe.data(), environ_g.data()); + double difference = 0.0; + for (int ig = 0; ig < basis.npw; ++ig) + { + const std::complex expected = environ_g[ig] * (volume / integral); + const double deviation = std::abs(kernel[ig] - expected); + difference = std::max(difference, deviation); + } + EXPECT_LT(pool_max(basis, difference), 1.0e-13); + const double kernel_at_zero = kernel_entry(basis, kernel, 0.0); + EXPECT_NEAR(kernel_at_zero, 1.0, 1.0e-12); + const double first_shell = ModuleBase::TWO_PI / length; + const double kernel_at_shell = kernel_entry(basis, kernel, 1.0); + const double isolated_at_shell = isolated_probe_transform(first_shell, width, spread); + EXPECT_NEAR(kernel_at_shell, isolated_at_shell, 1.0e-8); + + // Any constant: the unit-integral probe keeps it. + const double constant_value = 0.37; + const std::vector constant(basis.nrxx, constant_value); + const std::vector convolved = ModuleSccs::convolve_probe(kernel, basis, constant); + double constant_error = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double deviation = std::abs(convolved[ir] - constant_value); + constant_error = std::max(constant_error, deviation); + } + EXPECT_LT(pool_max(basis, constant_error), 1.0e-12); +} + +// Image counts must follow the cell heights, not the lattice-vector lengths: +// a hexagonal cell narrower than the probe, and a cell sheared to 10 degrees +// whose 2.1 Bohr height needs images n = +-4 along a 12 Bohr vector. Counts +// from the lengths miss probe weight (kernel error 1.9e-2, against 8e-15). +TEST_F(SccsSolventAware, ObliqueCellProbeSumsPeriodicImages) +{ + const double sqrt3_half = 0.5 * std::sqrt(3.0); + const ModuleBase::Matrix3 hexagonal(1.0, 0.0, 0.0, + -0.5, sqrt3_half, 0.0, + 0.0, 0.0, 1.2); + const double shear_angle = 10.0 * ModuleBase::PI / 180.0; + const double shear_cos = std::cos(shear_angle); + const double shear_sin = std::sin(shear_angle); + const ModuleBase::Matrix3 sheared(1.0, 0.0, 0.0, + shear_cos, shear_sin, 0.0, + 0.0, 0.0, 1.0); + const ModuleSccs::SolventAwareParameters probe = probe_parameters(3.0); + double hexagonal_count = 0.0; + const double hexagonal_error = oblique_probe_error(hexagonal, 9.0, probe, 1.5, hexagonal_count); + double sheared_count = 0.0; + const double sheared_error = oblique_probe_error(sheared, 12.0, probe, 1.5, sheared_count); + EXPECT_GT(hexagonal_count, 20.0); + EXPECT_GT(sheared_count, 20.0); + EXPECT_LT(hexagonal_error, 1.0e-6); + EXPECT_LT(sheared_error, 1.0e-6); +} + +// A thick shell around a small void: the void cannot hold the probe and is +// filled, while the open solvent outside the shell stays solvent. +TEST_F(SccsSolventAware, FillsAnEnclosedVoidAndKeepsTheOpenSolvent) +{ + const ModuleSccs::SolventAwareParameters probe = probe_parameters(3.0); + set_up_cube(30.0, 20.0, probe); + const std::vector outer = sphere(origin, 8.0, 0.5); + const std::vector inner = sphere(origin, 1.5, 0.5); + std::vector shell(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + shell[ir] = outer[ir] - inner[ir]; + } + const ModuleSccs::SolventAwareBoundary filled = fill(shell); + // center local value, center filled value, filled volume + double sums[3] = {0.0, 0.0, 0.0}; + double outside_filled = 0.0; + const double rounding = 1.0e-15; + const double upper_bound = 1.0 + rounding; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double lower_bound = shell[ir] - rounding; + EXPECT_GE(filled.boundary[ir], lower_bound); + EXPECT_LE(filled.boundary[ir], upper_bound); + const double distance = displacement(ir).norm(); + if (distance < 1.0e-12) + { + sums[0] += shell[ir]; + sums[1] += filled.boundary[ir]; + } + if (distance > 10.5) + { + outside_filled = std::max(outside_filled, filled.boundary[ir]); + } + sums[2] += (filled.boundary[ir] - shell[ir]) * length * length * length / basis.nxyz; + } + Parallel_Reduce::reduce_pool(sums, 3); + EXPECT_LT(sums[0], 1.0e-4); + EXPECT_GT(sums[1], 0.999); + EXPECT_LT(pool_max(basis, outside_filled), 1.0e-6); + // The void holds about 4 pi 1.5^3 / 3 = 14 Bohr^3. + EXPECT_GT(sums[2], 5.0); + EXPECT_LT(sums[2], 40.0); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_thread_sum.cpp b/source/source_hamilt/module_sccs/test/test_sccs_thread_sum.cpp new file mode 100644 index 00000000000..ec5d9ba8574 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_thread_sum.cpp @@ -0,0 +1,98 @@ +#include "../sccs_thread_sum.h" + +#include +#include +#include +#include + +#include + +#ifdef _OPENMP +#include +#endif + +namespace +{ +// Sets the OpenMP team size for one test and turns off dynamic teams, so that +// the runtime cannot shrink the team. Both settings are restored afterwards, so +// the other tests of the binary keep OMP_NUM_THREADS and OMP_DYNAMIC. +class ThreadCount +{ +public: + explicit ThreadCount(int count) + { +#ifdef _OPENMP + previous_ = omp_get_max_threads(); + previous_dynamic_ = omp_get_dynamic(); + omp_set_dynamic(0); + omp_set_num_threads(count); +#else + (void)count; +#endif + } + + ~ThreadCount() + { +#ifdef _OPENMP + omp_set_num_threads(previous_); + omp_set_dynamic(previous_dynamic_); +#endif + } + +private: + int previous_ = 1; + int previous_dynamic_ = 0; +}; +} // namespace + +// Three threads over 1001 points: the static shares are uneven. Integer terms +// sum exactly, 1/(i + 1) matches the serial loop to rounding, and a second run +// with the same team gives the same bits. +TEST(SccsThreadSum, UnevenSharesMatchTheSerialSumAndRepeat) +{ + const ThreadCount team(3); + const std::size_t size = 1001; + std::vector owner(size, 0); + const auto add_terms = [&owner](std::size_t i, std::array& sum) { +#ifdef _OPENMP + owner[i] = omp_get_thread_num(); +#endif + const double index = static_cast(i); + sum[0] += index; + sum[1] += 1.0 / (index + 1.0); + }; + const std::array first = ModuleSccs::thread_sums<2>(size, add_terms); + const std::array second = ModuleSccs::thread_sums<2>(size, add_terms); + double serial = 0.0; + for (std::size_t i = 0; i < size; ++i) + { + const double index = static_cast(i); + serial += 1.0 / (index + 1.0); + } + const double point_count = static_cast(size); + const double integer_sum = 0.5 * point_count * (point_count - 1.0); + EXPECT_EQ(first[0], integer_sum); + EXPECT_NEAR(first[1], serial, 1e-13); + EXPECT_EQ(first[0], second[0]); + EXPECT_EQ(first[1], second[1]); +#ifdef _OPENMP + // The partial sums of more than one thread were combined. + const std::set threads(owner.begin(), owner.end()); + EXPECT_GT(threads.size(), 1u); +#endif +} + +// Four threads over two points leave two threads without points, and an empty +// range gives zero. +TEST(SccsThreadSum, ThreadsWithoutPointsAddZero) +{ + const ThreadCount team(4); + const auto add_point = [](std::size_t i, std::array& sum) { + const double index = static_cast(i); + sum[0] += 0.25 + index; + }; + const std::array pair = ModuleSccs::thread_sums<1>(2, add_point); + const std::array empty = ModuleSccs::thread_sums<1>(0, add_point); + EXPECT_EQ(pair[0], 1.5); + EXPECT_EQ(empty[0], 0.0); +} diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 8ade361f51d..892d9f02dd0 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -627,6 +627,11 @@ struct Input_para 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 + double sccs_solvent_radius = 0.0; // Bohr: solvent-aware cavity, 0 is off + double sccs_radial_scale = 2.0; // solvent-aware probe radius over sccs_solvent_radius + double sccs_radial_spread = 0.5; // Bohr: erfc spread of the solvent-aware probe + double sccs_filling_thr = 0.825; // probe solute fraction that fills a point + double sccs_filling_spread = 0.02; // erfc width of the solvent-aware filling // ============== #Parameters (14.vdW Correction) =========================== // ========================================================== diff --git a/source/source_io/module_parameter/read_inp_out.cpp b/source/source_io/module_parameter/read_inp_out.cpp index ea6361e129f..c0a1202ad64 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -178,8 +178,9 @@ In molecular dynamics calculations, the output frequency is controlled by out_fr 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 out_freq_ion = 0 they are written at every step. In out_app_flag 1, fixed names such as `sol_eps.cube` hold one complete frame per written step: ionic step 0 truncates the file and later steps append. In modes 0 and 2, positive out_freq_ion adds the geometry index, e.g. `sol_epsg1.cube`; frequency zero overwrites the fixed name at each step. Solvent files stay in OUT.{suffix} in all modes. Nothing is written while SCCS waits for its delayed start (sccs_start_drho). + * `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 the solvent-aware cavity (sccs_solvent_radius > 0) it is the filled boundary. + * `sol_cavity_local.cube` and `sol_filled_fraction.cube`: only with the solvent-aware cavity, the boundary before the filling and the solute fraction of the probe sphere around each point; points where this fraction exceeds sccs_filling_thr are filled. + With out_freq_ion > 0 the files are written every out_freq_ion ionic steps; with out_freq_ion = 0 they are written at every step. In out_app_flag 1, fixed names such as `sol_eps.cube` hold one complete frame per written step: ionic step 0 truncates the file and later steps append. In out_app_flag 2, positive out_freq_ion adds the geometry index, e.g. `sol_epsg1.cube`. In mode 0, and in mode 2 with out_freq_ion = 0, each written step overwrites the fixed name. Solvent files stay in OUT.{suffix} in all modes. 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"; diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index 3d48b6c504a..ec6d0294702 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -132,7 +132,7 @@ void ReadInput::item_sccs() item.annotation = "SCCS switching-derivative low-pass slope"; item.category = "Implicit solvation model"; item.type = "Real"; - item.description = "Low-pass filter of the SCCS switching-function derivatives, as Environ deriv_lowpass_p1 with deriv_method fft. Only with assume_isolated pcc_0d or pcc_2d. The default -1 turns it off and reproduces Environ deriv_method fft (continuum cavity potential). With lowpass disabled (the default), analytical forces may differ from finite differences of the self-consistent energy. For geometry optimization with PCC, consider enabling lowpass and check force accuracy against finite differences. With lowpass disabled, the cavity potential uses the FFT gradient of the PCC-corrected potential, which oscillates around the potential step at the cell boundary half a cell from the system center; keep the dielectric transition region several bohr away from that boundary. Typical values are sccs_lowpass_p1 10 and sccs_lowpass_p2 5. The filter changes the model energy, so compare only results computed with the same values."; + item.description = "Low-pass filter of the SCCS switching-function derivatives, as Environ deriv_lowpass_p1 with deriv_method fft. Only with assume_isolated pcc_0d or pcc_2d. The default -1 turns it off, which reproduces Environ deriv_method fft (continuum cavity potential) for the plain cavity; the solvent-aware cavity then differentiates the filling analytically on the FFT derivatives of the local cavity. With lowpass disabled (the default), analytical forces may differ from finite differences of the self-consistent energy. For geometry optimization with PCC, consider enabling lowpass and check force accuracy against finite differences. With lowpass disabled, the cavity potential uses the FFT gradient of the PCC-corrected potential, which oscillates around the potential step at the cell boundary half a cell from the system center; keep the dielectric transition region several bohr away from that boundary. Typical values are sccs_lowpass_p1 10 and sccs_lowpass_p2 5. The filter changes the model energy, so compare only results computed with the same values."; item.default_value = "-1"; item.unit = ""; item.set_availability("imp_sol==2"); @@ -183,6 +183,66 @@ void ReadInput::item_sccs() sync_doublevec(input.sccs_corespread, para.input.sccs_corespread.size(), 0.0); this->add_item(item); } + { + Input_Item item("sccs_solvent_radius"); + item.annotation = "solvent radius of the solvent-aware SCCS cavity"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Solvent radius in bohr of the solvent-aware SCCS cavity, as Environ solvent_radius (Andreussi et al., J. Chem. Theory Comput. 15, 1996 (2019)). A positive value fills the parts of the cavity that a probe sphere of radius sccs_radial_scale * sccs_solvent_radius cannot reach, for example small voids between atoms. Unlike Environ, the probe sums all periodic images within its cutoff, so cells narrower than the probe are not truncated. With PCC, see also sccs_lowpass_p1. Non-negative; default 0 (off)."; + item.default_value = "0.0"; + item.unit = "bohr"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_solvent_radius); + this->add_item(item); + } + { + Input_Item item("sccs_radial_scale"); + item.annotation = "solvent-aware probe radius over sccs_solvent_radius"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Probe radius of the solvent-aware SCCS cavity in units of sccs_solvent_radius, as Environ radial_scale; at least 1, default 2. Used only when sccs_solvent_radius is positive."; + item.default_value = "2.0"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_radial_scale); + this->add_item(item); + } + { + Input_Item item("sccs_radial_spread"); + item.annotation = "erfc spread of the solvent-aware probe"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "erfc spread of the solvent-aware probe sphere, as Environ radial_spread; positive, default 0.5 bohr. Used only when sccs_solvent_radius is positive."; + item.default_value = "0.5"; + item.unit = "bohr"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_radial_spread); + this->add_item(item); + } + { + Input_Item item("sccs_filling_thr"); + item.annotation = "probe solute fraction that fills a point"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Solute fraction of the solvent-aware probe sphere above which a point is filled, as Environ filling_threshold; between 0 and 1 (exclusive), default 0.825. Used only when sccs_solvent_radius is positive."; + item.default_value = "0.825"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_filling_thr); + this->add_item(item); + } + { + Input_Item item("sccs_filling_spread"); + item.annotation = "erfc width of the solvent-aware filling"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Width of the erfc step of the solvent-aware filling in the probe solute fraction, as Environ filling_spread; positive, default 0.02. Used only when sccs_solvent_radius is positive. As in Environ, open solvent is filled by erfc(sccs_filling_thr / sccs_filling_spread) / 2, which lowers the bulk permittivity: below 1e-17 up to a spread of sccs_filling_thr / 6, but 1e-2 at 0.5."; + item.default_value = "0.02"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_filling_spread); + this->add_item(item); + } { Input_Item item("sccs_start_drho"); item.annotation = "SCCS delayed-start density threshold"; @@ -212,7 +272,7 @@ void ReadInput::item_sccs() 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.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 (with the solvent-aware filled volume), 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); diff --git a/source/source_io/module_parameter/sccs_inp_check.cpp b/source/source_io/module_parameter/sccs_inp_check.cpp index 670fa3634e6..0ad26cad2ec 100644 --- a/source/source_io/module_parameter/sccs_inp_check.cpp +++ b/source/source_io/module_parameter/sccs_inp_check.cpp @@ -121,7 +121,12 @@ std::string nonfinite_error(const Input_para& input) {"sccs_surface_eta", input.sccs_surface_eta}, {"sccs_lowpass_p1", input.sccs_lowpass_p1}, {"sccs_lowpass_p2", input.sccs_lowpass_p2}, - {"sccs_start_drho", input.sccs_start_drho}}; + {"sccs_start_drho", input.sccs_start_drho}, + {"sccs_solvent_radius", input.sccs_solvent_radius}, + {"sccs_radial_scale", input.sccs_radial_scale}, + {"sccs_radial_spread", input.sccs_radial_spread}, + {"sccs_filling_thr", input.sccs_filling_thr}, + {"sccs_filling_spread", input.sccs_filling_spread}}; for (const std::pair& value : values) { if (!std::isfinite(value.second)) @@ -175,6 +180,35 @@ std::string model_error(const Input_para& input) return error; } +// Solvent-aware probe and filling; a zero radius turns the filling off. +std::string sa_error(const Input_para& input) +{ + std::string error; + if (input.sccs_solvent_radius < 0.0) + { + error = "sccs_solvent_radius must be non-negative"; + } + else if (input.sccs_radial_scale < 1.0) + { + error = "sccs_radial_scale must be at least 1"; + } + else if (input.sccs_radial_spread <= 0.0) + { + error = "sccs_radial_spread must be positive"; + } + // The probe solute fraction never exceeds one, so a threshold of one or + // more would fill nothing. + else if (input.sccs_filling_thr <= 0.0 || input.sccs_filling_thr >= 1.0) + { + error = "sccs_filling_thr must lie between 0 and 1"; + } + else if (input.sccs_filling_spread <= 0.0) + { + error = "sccs_filling_spread must be positive"; + } + return error; +} + // sqrt-CG solver controls and the delayed start. std::string solver_error(const Input_para& input) { @@ -230,6 +264,10 @@ bool check_sccs_inp(const Input_para& input, std::string& error) error = model_error(input); } if (error.empty()) + { + error = sa_error(input); + } + if (error.empty()) { error = solver_error(input); } 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 18cbabb6076..9f11cb1142b 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -110,6 +110,39 @@ TEST(ReadInpSccs, FullCavityWidthsAndModeValidation) EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); } +TEST(ReadInpSccs, SolventAwareFillingUsesEnvironRanges) +{ + Input_para input; + input.imp_sol = 2; + std::string error; + EXPECT_DOUBLE_EQ(input.sccs_solvent_radius, 0.0); + EXPECT_DOUBLE_EQ(input.sccs_radial_scale, 2.0); + EXPECT_DOUBLE_EQ(input.sccs_radial_spread, 0.5); + EXPECT_DOUBLE_EQ(input.sccs_filling_thr, 0.825); + EXPECT_DOUBLE_EQ(input.sccs_filling_spread, 0.02); + input.sccs_solvent_radius = 3.0; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)) << error; + const Input_para valid = input; + input.sccs_solvent_radius = -1.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + EXPECT_NE(error.find("sccs_solvent_radius"), std::string::npos); + input = valid; + input.sccs_radial_scale = 0.9; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input = valid; + input.sccs_radial_spread = 0.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input = valid; + input.sccs_filling_thr = 0.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + // The probe solute fraction never exceeds one. + input.sccs_filling_thr = 1.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input = valid; + input.sccs_filling_spread = 0.0; + EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); +} + TEST(ReadInpSccs, DelayedStartRequiresRoomBeforeScfConvergenceAndIterationLimit) { Input_para input; @@ -202,7 +235,9 @@ 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_start_drho}; + &input.sccs_lowpass_p2, &input.sccs_start_drho, &input.sccs_solvent_radius, + &input.sccs_radial_scale, &input.sccs_radial_spread, &input.sccs_filling_thr, + &input.sccs_filling_spread}; const double infinity = std::numeric_limits::infinity(); const double nonfinite[] = {std::nan(""), infinity, -infinity}; for (double* value : reals) diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 1b8261b7a8a..c135aaa9fca 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -134,3 +134,5 @@ scf_sccs_pcc_0d scf_sccs_pcc_0d_lowpass scf_sccs_full relax_sccs +scf_sccs_sa +scf_sccs_sa_pcc_2d_lowpass diff --git a/tests/01_PW/scf_sccs_sa/INPUT b/tests/01_PW/scf_sccs_sa/INPUT new file mode 100644 index 00000000000..d90d0109ae5 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/INPUT @@ -0,0 +1,25 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated none +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +out_sol 1 +sccs_solvent_radius 2.0 +sccs_filling_thr 0.3 +sccs_debug 2 diff --git a/tests/01_PW/scf_sccs_sa/KPT b/tests/01_PW/scf_sccs_sa/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_sa/README b/tests/01_PW/scf_sccs_sa/README new file mode 100644 index 00000000000..bc7eb2ce0f7 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/README @@ -0,0 +1 @@ +PW neutral water with periodic SCCS and the solvent-aware cavity (solvent radius 2 bohr, filling threshold 0.3); checks energy, forces, solvation terms, the filled volume (sccs_debug 2) and the out_sol cavity integral. diff --git a/tests/01_PW/scf_sccs_sa/STRU b/tests/01_PW/scf_sccs_sa/STRU new file mode 100644 index 00000000000..ed974dd43a5 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +12.0 0.0 0.0 +0.0 12.0 0.0 +0.0 0.0 12.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +7.433 7.111 6.0 0 0 0 +4.567 7.111 6.0 0 0 0 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/scf_sccs_sa/result.ref b/tests/01_PW/scf_sccs_sa/result.ref new file mode 100644 index 00000000000..8de04421bff --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/result.ref @@ -0,0 +1,8 @@ +etotref -441.5798082776961451 +etotperatomref -147.1932694259 +totalforceref 16.894298 +esolelref -0.1881753362 +esolcavref 0.0796469683 +sccs_filled_volumeref 29.529984 +sol_cavity_integralref 275.4640801 +totaltimeref 1.86 diff --git a/tests/01_PW/scf_sccs_sa/threshold b/tests/01_PW/scf_sccs_sa/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/INPUT b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/INPUT new file mode 100644 index 00000000000..ffbef3030ca --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/INPUT @@ -0,0 +1,27 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated pcc_2d +pcc_2d_axis 1 +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +sccs_solvent_radius 2.0 +sccs_filling_thr 0.3 +sccs_debug 2 +sccs_lowpass_p1 10 +sccs_lowpass_p2 5 diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/KPT b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/README b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/README new file mode 100644 index 00000000000..173347614cb --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/README @@ -0,0 +1 @@ +PW neutral water with SCCS, 2D PCC, lowpass and the solvent-aware cavity (solvent radius 2 bohr, filling threshold 0.3); also checks the filled volume (sccs_debug 2). diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/STRU b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/STRU new file mode 100644 index 00000000000..a9048229a55 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +10.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 10.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +6.433 11.111 5.0 0 0 0 +3.567 11.111 5.0 0 0 0 +O +0.0 +1 +5.0 10.0 5.0 0 0 0 diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/result.ref b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/result.ref new file mode 100644 index 00000000000..8ef5cae479c --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/result.ref @@ -0,0 +1,8 @@ +etotref -441.9910851843064279 +etotperatomref -147.3303617281 +totalforceref 16.437145 +esolelref -0.2181773800 +esolcavref 0.0784755279 +sccs_filled_volumeref 63.325291 +epccref 0.0349785671 +totaltimeref 3.06 diff --git a/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/threshold b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/scf_sccs_sa_pcc_2d_lowpass/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/02_NAO_Gamma/CASES_CPU.txt b/tests/02_NAO_Gamma/CASES_CPU.txt index ce8e970bdb8..ee66cdc229b 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -41,6 +41,7 @@ test.sum scf_pcc0d_water scf_pcc2d_water scf_sccs_pcc2d +scf_sccs_sa_pcc2d scf_sccs_water scf_sccs_water_delay relax_sccs_pcc0d_h3o diff --git a/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/INPUT b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/INPUT new file mode 100644 index 00000000000..102f5a2d20e --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/INPUT @@ -0,0 +1,31 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +calculation scf +basis_type lcao +gamma_only 1 +nbands 8 +ecutwfc 40 +scf_nmax 80 +scf_thr 1.0e-10 +cal_force 1 +symmetry 0 + +imp_sol 2 +assume_isolated pcc_2d +pcc_2d_axis 1 +sccs_preset custom +sccs_epsilon 1.1 +sccs_rho_min 1.0e-4 +sccs_rho_max 5.0e-3 +sccs_gamma 0.0 +sccs_pressure 0.0 +sccs_maxiter 100 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 +sccs_solvent_radius 2.0 +sccs_filling_thr 0.3 +sccs_debug 2 diff --git a/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/KPT b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/README b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/README new file mode 100644 index 00000000000..502e3f4ef4b --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/README @@ -0,0 +1 @@ +gamma-only LCAO neutral water with SCCS (epsilon 1.1), 2D PCC and the solvent-aware cavity (solvent radius 2 bohr, filling threshold 0.3); also checks the filled volume (sccs_debug 2). diff --git a/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/STRU b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/STRU new file mode 100644 index 00000000000..409f8328a37 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/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 +10.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 10.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +6.433 11.111 5.0 0 0 0 +3.567 11.111 5.0 0 0 0 + +O +0.0 +1 +5.0 10.0 5.0 0 0 0 diff --git a/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/result.ref b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/result.ref new file mode 100644 index 00000000000..642a2762aeb --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_sa_pcc2d/result.ref @@ -0,0 +1,8 @@ +etotref -466.3753101788931 +etotperatomref -155.4584367263 +totalforceref 3.090095 +esolelref -0.0148188993 +esolcavref 0.0000000000 +sccs_filled_volumeref 30.742821 +epccref 0.0516211963 +totaltimeref 1.46 diff --git a/tests/integrate/validation_tools/props_basic.sh b/tests/integrate/validation_tools/props_basic.sh index 85495f51908..2814c09db2a 100644 --- a/tests/integrate/validation_tools/props_basic.sh +++ b/tests/integrate/validation_tools/props_basic.sh @@ -99,12 +99,34 @@ run_basic_props_post_ml() { 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. + # Solvent-aware SCCS (imp_sol 2, sccs_solvent_radius > 0). With sccs_debug 2 + # the running log must report the filled volume int(s_sa - s) dV of the last + # evaluation, which is compared. + solvent_aware=0 + if [ "$imp_sol" = 2 ] && awk -v radius="$sccs_solvent_radius" 'BEGIN {exit !(radius + 0 > 0)}'; then + solvent_aware=1 + fi + if [ "$solvent_aware" = 1 ] && [ "$sccs_debug" = 2 ]; then + filled_line=$(grep "V_FILLED/Bohr3" "$running_path" | tail -1) + filled_volume=$(echo "$filled_line" | awk '{for (i = 1; i < NF; ++i) if ($i == "V_FILLED/Bohr3") print $(i + 1)}') + if [ -z "$filled_volume" ]; then + echo "Error: solvent-aware SCCS with sccs_debug 2 but no V_FILLED/Bohr3 in $running_path" + exit 1 + fi + echo "sccs_filled_volumeref $filled_volume" >>$props_result_file + fi + + # out_sol 1: the cubes must exist (the cases set imp_sol 2, or imp_sol 1 with + # eb_k > 1, so sol_cavity.cube is written; the solvent-aware cavity adds two). + # The cavity volume int s dV of sol_cavity.cube is compared; the other cubes + # are only checked for presence: the sol_eps.cube integral is dominated by the + # bulk permittivity times the volume, and eps follows from the cavity. if [[ "$out_sol" == "1" ]]; then - for cube in sol_eps.cube sol_cavity.cube; do + sol_cubes="sol_eps.cube sol_cavity.cube" + if [ "$solvent_aware" = 1 ]; then + sol_cubes="$sol_cubes sol_cavity_local.cube sol_filled_fraction.cube" + fi + for cube in $sol_cubes; do if [ ! -f "OUT.autotest/$cube" ]; then echo "Error: out_sol 1 but OUT.autotest/$cube is missing" exit 1 diff --git a/tests/integrate/validation_tools/props_common.sh b/tests/integrate/validation_tools/props_common.sh index 72855acf1bc..805ea7f9def 100644 --- a/tests/integrate/validation_tools/props_common.sh +++ b/tests/integrate/validation_tools/props_common.sh @@ -123,6 +123,8 @@ props_init() { 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") + sccs_solvent_radius=$(get_input_key_value "sccs_solvent_radius" "INPUT") + sccs_debug=$(get_input_key_value "sccs_debug" "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")