diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 8c66bde4d2..ced2b16b82 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -284,6 +284,8 @@ OBJS_ELECSTAT=elecstate.o\ pcc_2d.o\ pot_pcc.o\ pot_sccs.o\ + sccs_pcc_0d.o\ + sccs_pcc_2d.o\ elecstate_op.o\ efield.o\ gatefield.o\ @@ -631,10 +633,14 @@ OBJS_SURCHEM=surchem.o\ sol_force.o\ OBJS_SCCS=sccs_cavity.o\ + sccs_cav_deriv.o\ sccs_functional.o\ sccs_parameters.o\ sccs_ion_charge.o\ sccs_ion_force.o\ + sccs_coul_0d.o\ + sccs_coul_2d.o\ + sccs_moments.o\ sccs_pw_coulomb.o\ sccs_response.o diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index bd3264574e..b78f842266 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -22,6 +22,8 @@ list(APPEND objects module_pot/pot_local.cpp module_pot/pot_pcc.cpp module_pot/pot_sccs.cpp + module_pot/sccs_pcc_0d.cpp + module_pot/sccs_pcc_2d.cpp module_pot/potential_new.cpp module_pot/potential_types.cpp module_pot/pot_surchem.cpp diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp index 6624e17624..cfdf978761 100644 --- a/source/source_estate/module_pot/pot_sccs.cpp +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -1,4 +1,6 @@ #include "pot_sccs.h" +#include "sccs_pcc_0d.h" +#include "sccs_pcc_2d.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" @@ -8,6 +10,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_pw_coulomb.h" #include "source_io/module_parameter/input_parameter.h" #include @@ -16,6 +19,13 @@ namespace elecstate { +namespace +{ +// Charge tolerance in e, the same as the PotPcc electron-count check: far +// above the grid round-off (~1e-12 e) and far below any physical net charge. +const double charge_tolerance = 1e-6; +} // namespace + void sccs_cfg_input(const Input_para& input, ModuleSccs::SccsConfig& config, ModuleSccs::PolarizationSolverParameters& solver) @@ -35,6 +45,19 @@ void sccs_cfg_input(const Input_para& input, config = ModuleSccs::preset_config(preset); } config.surface_regularization = input.sccs_surface_eta; + if (input.assume_isolated == "pcc_0d") + { + config.boundary = ModuleSccs::Boundary::Pcc0d; + } + else if (input.assume_isolated == "pcc_2d") + { + config.boundary = ModuleSccs::Boundary::Pcc2d; + config.pcc_2d_axis = input.pcc_2d_axis; + } + else + { + config.boundary = ModuleSccs::Boundary::Periodic; + } solver.max_iterations = input.sccs_maxiter; solver.tolerance_rms = input.sccs_tol_rms; solver.tolerance_max = input.sccs_tol_max; @@ -48,9 +71,14 @@ std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& c ionic_charge += cell.atoms[it].ncpp.zv * cell.atoms[it].na; } const double net_charge = ionic_charge - electron_count; - if (config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= 1e-6) { return std::string(); } + const bool periodic = config.boundary == ModuleSccs::Boundary::Periodic; + if (!periodic || config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= charge_tolerance) + { + return std::string(); + } return "charged SCCS with periodic boundaries: the periodic Poisson solver drops the G = 0 component " - "of the net charge, so the energy depends on the cell size."; + "of the net charge, so the energy depends on the cell size. Use assume_isolated pcc_0d or pcc_2d " + "for converged charged energies."; } PotSccs::PotSccs(const ModulePW::PW_Basis* basis, @@ -87,22 +115,36 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase:: double valence_charge = 0.0; for (const unitcell::AtomData& atom : atoms) { valence_charge += atom.valence_charge; } const double ionic_error = grid_charges[1] * dv - valence_charge; - if (std::abs(ionic_error) > 1e-6) + if (std::abs(ionic_error) > charge_tolerance) { ModuleBase::WARNING_QUIT("PotSccs::cal_v_eff", "SCCS ionic density normalization does not match the valence charge"); } - if (std::abs(electron_error) > 1e-6) + // With PCC the PCC component reports the same mismatch. + if (config_.boundary == ModuleSccs::Boundary::Periodic && std::abs(electron_error) > charge_tolerance) { std::ostringstream message; message << "SCCS grid electron count differs from the expected value by " << electron_error << " e; the solute charge uses the grid density"; ModuleBase::WARNING("PotSccs::cal_v_eff", message.str()); } + std::unique_ptr coulomb; + if (config_.boundary == ModuleSccs::Boundary::Pcc0d) + { + pcc_0d_operator(*cell, basis, atoms, coulomb); + } + else if (config_.boundary == ModuleSccs::Boundary::Pcc2d) + { + pcc_2d_operator(*cell, basis, atoms, config_.pcc_2d_axis, coulomb); + } + else + { + coulomb.reset(new ModuleSccs::PeriodicCoulombOperator(basis, cell->tpiba)); + } ModuleSccs::SccsResponse response; ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, - cell->tpiba, response); + cell->tpiba, *coulomb, response); ModuleSccs::FunctionalResult functional; - ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, functional); + ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, *coulomb, functional); electrostatic_rydberg_ = 2.0 * functional.reaction_energy; non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); electrostatic_potential_.resize(basis.nrxx); diff --git a/source/source_estate/module_pot/sccs_pcc_0d.cpp b/source/source_estate/module_pot/sccs_pcc_0d.cpp new file mode 100644 index 0000000000..f2d1de6f3e --- /dev/null +++ b/source/source_estate/module_pot/sccs_pcc_0d.cpp @@ -0,0 +1,58 @@ +#include "sccs_pcc_0d.h" + +#include "source_base/tool_quit.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_geometry.h" +#include "source_cell/cell_tools.h" +#include "source_cell/unitcell.h" +#include "source_hamilt/module_sccs/sccs_coul_0d.h" + +namespace elecstate +{ +void pcc_0d_geometry(const UnitCell& cell, + const std::vector& atoms, + unitcell::OrthogonalCell& geometry, + Pcc0dParameters& parameters) +{ + // Same orthogonality and equal-edge tolerance as PotPcc. + const double cell_tolerance = 1e-10; + const bool orthogonal = unitcell::make_orthogonal_cell(cell.latvec, cell.lat0, cell_tolerance, geometry); + const bool cubic = orthogonal && make_pcc_0d_parameters(geometry, cell_tolerance, parameters); + if (!cubic) + { + ModuleBase::WARNING_QUIT("pcc_0d_geometry", "PCC 0D requires an equal-edge cubic cell"); + } + std::vector> positions; + std::vector masses; + for (const unitcell::AtomData& atom : atoms) + { + positions.push_back(atom.position); + masses.push_back(atom.mass); + } + geometry.origin = unitcell::weighted_center(positions, masses, geometry); +} + +void pcc_0d_relative(const unitcell::OrthogonalCell& geometry, + std::vector>& positions) +{ + for (ModuleBase::Vector3& position : positions) + { + position = unitcell::relative_position(position, geometry); + } +} + +void pcc_0d_operator(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + std::unique_ptr& coulomb) +{ + unitcell::OrthogonalCell geometry; + Pcc0dParameters parameters; + pcc_0d_geometry(cell, atoms, geometry, parameters); + std::vector> relative_positions; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, relative_positions); + pcc_0d_relative(geometry, relative_positions); + coulomb.reset(new ModuleSccs::Pcc0dCoulombOperator(basis, cell.tpiba, relative_positions, parameters)); +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_pcc_0d.h b/source/source_estate/module_pot/sccs_pcc_0d.h new file mode 100644 index 0000000000..97a6b8cd01 --- /dev/null +++ b/source/source_estate/module_pot/sccs_pcc_0d.h @@ -0,0 +1,48 @@ +#ifndef SCCS_PCC_0D_ADAPTER_H +#define SCCS_PCC_0D_ADAPTER_H + +#include +#include + +class UnitCell; +namespace ModuleBase +{ +template class Vector3; +} +namespace ModulePW +{ +class PW_Basis; +} +namespace ModuleSccs +{ +class CoulombOperator; +} +namespace unitcell +{ +struct AtomData; +struct OrthogonalCell; +} +namespace elecstate +{ +struct Pcc0dParameters; + +// Cubic PCC cell with the mass-weighted ionic origin, as PotPcc; a cell that +// is not an equal-edge cube stops the run. +void pcc_0d_geometry(const UnitCell& cell, + const std::vector& atoms, + unitcell::OrthogonalCell& geometry, + Pcc0dParameters& parameters); + +// Minimum-image displacements of Cartesian positions from the origin, in place. +void pcc_0d_relative(const unitcell::OrthogonalCell& geometry, + std::vector>& positions); + +// Assemble an SCCS grid-charge operator using the existing cell geometry and +// mass-weighted ionic origin. An unsupported cell stops the run. +void pcc_0d_operator(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + std::unique_ptr& coulomb); +} // namespace elecstate + +#endif diff --git a/source/source_estate/module_pot/sccs_pcc_2d.cpp b/source/source_estate/module_pot/sccs_pcc_2d.cpp new file mode 100644 index 0000000000..44a403fefa --- /dev/null +++ b/source/source_estate/module_pot/sccs_pcc_2d.cpp @@ -0,0 +1,65 @@ +#include "sccs_pcc_2d.h" + +#include "source_base/tool_quit.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_geometry.h" +#include "source_cell/cell_tools.h" +#include "source_cell/unitcell.h" +#include "source_hamilt/module_sccs/sccs_coul_2d.h" + +namespace elecstate +{ +void pcc_2d_geometry(const UnitCell& cell, + const std::vector& atoms, + int open_axis, + unitcell::SlabCell& geometry, + Pcc2dParameters& parameters) +{ + // Same alignment tolerance as PotPcc, loose enough for lattice vectors + // typed with rounded digits. + const double alignment_tolerance = 1e-6; + const bool perpendicular = unitcell::make_slab_cell(cell.latvec, cell.lat0, open_axis, alignment_tolerance, + geometry); + if (!perpendicular) + { + ModuleBase::WARNING_QUIT("pcc_2d_geometry", + "PCC 2D requires the open lattice vector to be perpendicular to the periodic plane"); + } + parameters = make_pcc_2d_parameters(geometry); + std::vector> positions; + std::vector masses; + for (const unitcell::AtomData& atom : atoms) + { + positions.push_back(atom.position); + masses.push_back(atom.mass); + } + geometry.origin = unitcell::weighted_center(positions, masses, geometry); +} + +void pcc_2d_project(const unitcell::SlabCell& geometry, + std::vector>& positions) +{ + for (ModuleBase::Vector3& position : positions) + { + const double coordinate = unitcell::relative_coordinate(position, geometry); + position = geometry.normal * coordinate; + } +} + +void pcc_2d_operator(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + int open_axis, + std::unique_ptr& coulomb) +{ + unitcell::SlabCell geometry; + Pcc2dParameters parameters; + pcc_2d_geometry(cell, atoms, open_axis, geometry, parameters); + std::vector> projected_positions; + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, projected_positions); + pcc_2d_project(geometry, projected_positions); + coulomb.reset(new ModuleSccs::Pcc2dCoulombOperator(basis, cell.tpiba, projected_positions, geometry.normal, + parameters)); +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/sccs_pcc_2d.h b/source/source_estate/module_pot/sccs_pcc_2d.h new file mode 100644 index 0000000000..dfb2a4f6f8 --- /dev/null +++ b/source/source_estate/module_pot/sccs_pcc_2d.h @@ -0,0 +1,51 @@ +#ifndef SCCS_PCC_2D_ADAPTER_H +#define SCCS_PCC_2D_ADAPTER_H + +#include +#include + +class UnitCell; +namespace ModuleBase +{ +template class Vector3; +} +namespace ModulePW +{ +class PW_Basis; +} +namespace ModuleSccs +{ +class CoulombOperator; +} +namespace unitcell +{ +struct AtomData; +struct SlabCell; +} +namespace elecstate +{ +struct Pcc2dParameters; + +// Slab PCC cell with the mass-weighted ionic origin along the normal, with +// the PotPcc alignment tolerance; a tilted open axis stops the run. +void pcc_2d_geometry(const UnitCell& cell, + const std::vector& atoms, + int open_axis, + unitcell::SlabCell& geometry, + Pcc2dParameters& parameters); + +// Replace Cartesian positions, in place, by the normal times their +// minimum-image coordinate from the slab origin. +void pcc_2d_project(const unitcell::SlabCell& geometry, + std::vector>& positions); + +// Assemble an SCCS grid-charge operator using the existing cell geometry and +// mass-weighted ionic origin. An unsupported cell stops the run. +void pcc_2d_operator(const UnitCell& cell, + const ModulePW::PW_Basis& basis, + const std::vector& atoms, + int open_axis, + std::unique_ptr& coulomb); +} // namespace elecstate + +#endif diff --git a/source/source_estate/test/sccs_cell_test.h b/source/source_estate/test/sccs_cell_test.h new file mode 100644 index 0000000000..9436d5bba5 --- /dev/null +++ b/source/source_estate/test/sccs_cell_test.h @@ -0,0 +1,44 @@ +#ifndef SCCS_CELL_TEST_H +#define SCCS_CELL_TEST_H + +#include "source_cell/unitcell.h" +#include "source_estate/module_charge/charge.h" +#include "source_hamilt/module_sccs/test/sccs_test.h" + +#include + +namespace SccsTest +{ +// One unit-charge atom at fractional (1/4, 1/4, 1/4) with 0.6 uniform +// electrons in one spin channel, on the cell of the PW basis. +class CellTest : public PwTest +{ +protected: + void set_up_cell() + { + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.mass = 1.0; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 0.6 / basis.omega; + density.assign(basis.nrxx, density_value); + channel = density.data(); + charge.nspin = 1; + charge.rho = &channel; + } + + UnitCell cell; + Atom atom; + Charge charge; + std::vector density; + double* channel = nullptr; +}; +} // namespace SccsTest + +#endif diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp index 9ed4b9c89f..523d28d16f 100644 --- a/source/source_estate/test/test_pot_sccs.cpp +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -173,6 +173,18 @@ TEST_F(PotSccsTest, InputMapsOntoSccsConfig) EXPECT_LT(config.pressure, 0.0); EXPECT_DOUBLE_EQ(config.surface_regularization, 1e-6); EXPECT_EQ(solver.max_iterations, 42); + EXPECT_EQ(config.boundary, ModuleSccs::Boundary::Periodic); + input.assume_isolated = "pcc_0d"; + elecstate::sccs_cfg_input(input, config, solver); + EXPECT_EQ(config.boundary, ModuleSccs::Boundary::Pcc0d); + input.assume_isolated = "pcc_2d"; + for (int axis = 0; axis < 3; ++axis) + { + input.pcc_2d_axis = axis; + elecstate::sccs_cfg_input(input, config, solver); + EXPECT_EQ(config.boundary, ModuleSccs::Boundary::Pcc2d); + EXPECT_EQ(config.pcc_2d_axis, axis); + } } TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) @@ -210,6 +222,9 @@ TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) ModuleSccs::SccsConfig vacuum_config = config; vacuum_config.cavity.epsilon_bulk = 1.0; EXPECT_TRUE(elecstate::check_charge(vacuum_config, cell, 0.5).empty()); + ModuleSccs::SccsConfig pcc_config = config; + pcc_config.boundary = ModuleSccs::Boundary::Pcc0d; + EXPECT_TRUE(elecstate::check_charge(pcc_config, cell, 0.5).empty()); const double density_value = 0.5 / basis.omega; std::vector density(basis.nrxx, density_value); diff --git a/source/source_estate/test/test_sccs_pcc_0d.cpp b/source/source_estate/test/test_sccs_pcc_0d.cpp new file mode 100644 index 0000000000..d9ab205bc7 --- /dev/null +++ b/source/source_estate/test/test_sccs_pcc_0d.cpp @@ -0,0 +1,88 @@ +#include "source_estate/module_pot/pot_sccs.h" +#include "source_estate/module_pot/pot_pcc.h" +#include "source_estate/module_pot/sccs_pcc_0d.h" +#include "source_cell/cell_tools.h" +#include "source_hamilt/module_sccs/sccs_coulomb.h" +#include "source_hamilt/module_sccs/sccs_functional.h" +#include "source_hamilt/module_sccs/sccs_ion_charge.h" +#include "source_hamilt/module_sccs/sccs_ion_force.h" +#include "source_hamilt/module_sccs/sccs_response.h" +#include "sccs_cell_test.h" +#include "source_io/module_parameter/input_parameter.h" + +class SccsPcc0dTest : public SccsTest::CellTest +{ +protected: + void SetUp() override + { + SccsTest::PwTest::SetUp(); + set_up_cell(); + } +}; + +// With epsilon 1 the solvent adds nothing, so the point-ion PCC is applied +// exactly once; with a dielectric it adds only the reaction and cavity terms. +TEST_F(SccsPcc0dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) +{ + const double epsilons[] = {1.0, 5.0}; + for (double epsilon : epsilons) + { + Input_para input; + input.assume_isolated = "pcc_0d"; + input.sccs_epsilon = epsilon; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::sccs_cfg_input(input, config, solver); + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::unique_ptr coulomb; + elecstate::pcc_0d_operator(cell, basis, atoms, coulomb); + std::vector ions; + ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, + ModuleSccs::gaussian_ion_spread, ions); + std::vector solute_charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + solute_charge[ir] = ions[ir] - density[ir]; + } + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(density, solute_charge, config.cavity, solver, cold, + basis, tpiba, *coulomb, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, + *coulomb, functional); + elecstate::PotSccs solvent(&basis, config, solver, 0.6); + elecstate::PotPcc correction(&basis, 0.6); + ModuleBase::matrix potential(1, basis.nrxx); + correction.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix point_potential = potential; + const double point_energy = correction.get_energy(); + solvent.cal_v_eff(&charge, &cell, potential); + const double reaction_energy_rydberg = 2.0 * functional.reaction_energy; + EXPECT_NEAR(solvent.get_energy(), reaction_energy_rydberg, 1e-12); + const double combined_energy = solvent.get_energy() + correction.get_energy(); + const double expected_energy = point_energy + 2.0 * functional.reaction_energy; + EXPECT_NEAR(combined_energy, expected_energy, 1e-12); + std::vector> reaction_force; + ModuleSccs::gauss_ion_force(atoms, functional.reaction_potential, basis, + tpiba, ModuleSccs::gaussian_ion_spread, reaction_force); + ModuleBase::matrix point_force(cell.nat, 3); + correction.add_force(cell, point_force); + ModuleBase::matrix combined_force = point_force; + solvent.add_sol_force(cell, combined_force); + for (int axis = 0; axis < 3; ++axis) + { + const double expected = point_force(0, axis) + 2.0 * reaction_force[0][axis]; + EXPECT_NEAR(combined_force(0, axis), expected, 1e-12); + } + const std::vector* electrostatic = solvent.sol_elec_pot(); + ASSERT_NE(electrostatic, nullptr); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = point_potential(0, ir) + 2.0 * functional.electron_potential[ir]; + EXPECT_NEAR(potential(0, ir), expected, 1e-12); + const double expected_potential = -2.0 * functional.reaction_potential[ir]; + EXPECT_NEAR((*electrostatic)[ir], expected_potential, 1e-12); + } + } +} diff --git a/source/source_estate/test/test_sccs_pcc_2d.cpp b/source/source_estate/test/test_sccs_pcc_2d.cpp new file mode 100644 index 0000000000..6c0f86ef3e --- /dev/null +++ b/source/source_estate/test/test_sccs_pcc_2d.cpp @@ -0,0 +1,110 @@ +#include "source_estate/module_pot/pot_sccs.h" +#include "source_estate/module_pot/pot_pcc.h" +#include "source_estate/module_pot/sccs_pcc_2d.h" +#include "source_cell/cell_tools.h" +#include "source_hamilt/module_sccs/sccs_coulomb.h" +#include "source_hamilt/module_sccs/sccs_functional.h" +#include "source_hamilt/module_sccs/sccs_ion_charge.h" +#include "source_hamilt/module_sccs/sccs_ion_force.h" +#include "source_hamilt/module_sccs/sccs_response.h" +#include "sccs_cell_test.h" +#include "source_io/module_parameter/input_parameter.h" + + +class SccsPcc2dTest : public SccsTest::CellTest +{ +protected: + void SetUp() override + { + cell.latvec.e33 = 1.5; + set_up_basis(cell.latvec); + set_up_cell(); + } +}; + +// With epsilon 1 the solvent adds nothing, so the point-ion PCC is applied +// exactly once; with a dielectric it adds only the reaction and cavity terms. +TEST_F(SccsPcc2dTest, SolventAddsOnlyReactionAndCavityTermsToThePointIonPcc) +{ + const double epsilons[] = {1.0, 5.0}; + for (double epsilon : epsilons) + { + Input_para input; + input.assume_isolated = "pcc_2d"; + input.sccs_epsilon = epsilon; + input.sccs_rho_min = 1e-3; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::sccs_cfg_input(input, config, solver); + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::unique_ptr coulomb; + elecstate::pcc_2d_operator(cell, basis, atoms, input.pcc_2d_axis, coulomb); + std::vector ions; + ModuleSccs::gauss_ion_rho(atoms, basis, tpiba, + ModuleSccs::gaussian_ion_spread, ions); + std::vector solute_charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + solute_charge[ir] = ions[ir] - density[ir]; + } + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(density, solute_charge, config.cavity, solver, cold, + basis, tpiba, *coulomb, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::eval_functional(solute_charge, response, config, basis, tpiba, + *coulomb, functional); + elecstate::PotSccs solvent(&basis, config, solver, 0.6); + elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, input.pcc_2d_axis, 0.6); + ModuleBase::matrix potential(1, basis.nrxx); + correction.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix point_potential = potential; + const double point_energy = correction.get_energy(); + solvent.cal_v_eff(&charge, &cell, potential); + const double reaction_energy_rydberg = 2.0 * functional.reaction_energy; + EXPECT_NEAR(solvent.get_energy(), reaction_energy_rydberg, 1e-12); + const double combined_energy = solvent.get_energy() + correction.get_energy(); + const double expected_energy = point_energy + 2.0 * functional.reaction_energy; + EXPECT_NEAR(combined_energy, expected_energy, 1e-12); + std::vector> reaction_force; + ModuleSccs::gauss_ion_force(atoms, functional.reaction_potential, basis, + tpiba, ModuleSccs::gaussian_ion_spread, reaction_force); + ModuleBase::matrix point_force(cell.nat, 3); + correction.add_force(cell, point_force); + ModuleBase::matrix combined_force = point_force; + solvent.add_sol_force(cell, combined_force); + for (int axis = 0; axis < 3; ++axis) + { + const double expected = point_force(0, axis) + 2.0 * reaction_force[0][axis]; + EXPECT_NEAR(combined_force(0, axis), expected, 1e-12); + } + const std::vector* electrostatic = solvent.sol_elec_pot(); + ASSERT_NE(electrostatic, nullptr); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = point_potential(0, ir) + 2.0 * functional.electron_potential[ir]; + EXPECT_NEAR(potential(0, ir), expected, 1e-12); + const double expected_potential = -2.0 * functional.reaction_potential[ir]; + EXPECT_NEAR((*electrostatic)[ir], expected_potential, 1e-12); + } + } +} + +TEST_F(SccsPcc2dTest, AcceptsThePccSlabAlignmentTolerance) +{ + // A tilt of 1e-8 relative to the open vector, as from rounded lattice + // digits, passes both PCC and the SCCS operator; a real tilt stops both. + cell.latvec.e31 = 1.5e-8; + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::unique_ptr coulomb; + elecstate::pcc_2d_operator(cell, basis, atoms, 2, coulomb); + EXPECT_NE(coulomb.get(), nullptr); + elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, 2, 0.6); + ModuleBase::matrix potential(1, basis.nrxx); + correction.cal_v_eff(&charge, &cell, potential); + cell.latvec.e31 = 1.5e-4; + std::unique_ptr tilted; + EXPECT_EXIT(elecstate::pcc_2d_operator(cell, basis, atoms, 2, tilted), + testing::ExitedWithCode(1), + ""); +} diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt index f0186d6722..d6f210a273 100644 --- a/source/source_hamilt/module_sccs/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -1,9 +1,13 @@ add_library(sccs OBJECT sccs_cavity.cpp + sccs_cav_deriv.cpp sccs_parameters.cpp sccs_functional.cpp sccs_ion_charge.cpp sccs_ion_force.cpp + sccs_coul_0d.cpp + sccs_coul_2d.cpp + sccs_moments.cpp sccs_pw_coulomb.cpp sccs_response.cpp) diff --git a/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp b/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp new file mode 100644 index 0000000000..ddfbad0246 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cav_deriv.cpp @@ -0,0 +1,98 @@ +#include "sccs_cav_deriv.h" +#include "sccs_cavity.h" +#include "sccs_response.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include + +namespace ModuleSccs +{ +namespace +{ +// 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. +void eps_of_boundary(const CavityParameters& cavity, + const std::vector>& gradient, + const std::vector& laplacian, + SccsResponse& response, + std::vector& coefficient) +{ + const std::size_t size = response.solute.size(); + const double log_bulk = std::log(cavity.epsilon_bulk); + coefficient.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const double log_epsilon = log_bulk * (1.0 - response.solute[i]); + response.epsilon[i] = std::exp(log_epsilon); + const ModuleBase::Vector3& local_gradient = gradient[i]; + const double gradient_square = local_gradient * local_gradient; + response.grad_log_epsilon[i] = local_gradient * (-log_bulk); + const double lap_log = -log_bulk * laplacian[i]; + coefficient[i] = response.epsilon[i] * (0.5 * lap_log + 0.25 * log_bulk * log_bulk * gradient_square) + / ModuleBase::FOUR_PI; + } +} +} // namespace + +void make_cav_derivs(const std::vector& density, + const CavityParameters& cavity, + const ModulePW::PW_Basis& basis, + double tpiba, + bool open_boundary, + SccsResponse& response, + std::vector& coefficient) +{ + const std::size_t size = density.size(); + response.solute.resize(size); + response.dsolute_drho.resize(size); + response.epsilon.resize(size); + response.depsilon_drho.resize(size); + response.grad_log_epsilon.resize(size); + std::vector d2solute_drho2(size); + for (std::size_t i = 0; i < size; ++i) + { + const CavityPoint point = evaluate_cavity(density[i], cavity); + response.solute[i] = point.solute; + response.dsolute_drho[i] = point.dsolute_drho; + response.epsilon[i] = point.epsilon; + response.depsilon_drho[i] = point.depsilon_drho; + d2solute_drho2[i] = point.d2solute_drho2; + } + + 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()); + XC_Functional::grad_rho(boundary_g.data(), gradient.data(), &basis, tpiba); + XC_Functional::laplacian_rho(boundary_g.data(), laplacian.data(), &basis, tpiba); + } + else + { + // 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); + 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; + } + } + 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 new file mode 100644 index 0000000000..eb2d86a752 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cav_deriv.h @@ -0,0 +1,27 @@ +#ifndef SCCS_CAV_DERIV_H +#define SCCS_CAV_DERIV_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct CavityParameters; +struct SccsResponse; + +// Collective over the PW pool. Periodic cells use density chain derivatives; +// open boundaries differentiate solute on the FFT grid. +void make_cav_derivs(const std::vector& density, + const CavityParameters& cavity, + const ModulePW::PW_Basis& basis, + double tpiba, + bool open_boundary, + SccsResponse& response, + std::vector& coefficient); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_cavity.cpp b/source/source_hamilt/module_sccs/sccs_cavity.cpp index 94f69ed0ab..08aeb09798 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.cpp +++ b/source/source_hamilt/module_sccs/sccs_cavity.cpp @@ -27,10 +27,13 @@ CavityPoint evaluate_cavity(double density, const CavityParameters& parameters) const double angle = ModuleBase::TWO_PI * x; const double solvent = x - std::sin(angle) / ModuleBase::TWO_PI; const double dsolvent_drho = -(1.0 - std::cos(angle)) / (log_width * density); + const double d2solvent_drho2 = (1.0 - std::cos(angle) + ModuleBase::TWO_PI * std::sin(angle) / log_width) + / (log_width * density * density); const double log_epsilon = std::log(parameters.epsilon_bulk); const double log_epsilon_point = log_epsilon * solvent; point.solute = 1.0 - solvent; point.dsolute_drho = -dsolvent_drho; + point.d2solute_drho2 = -d2solvent_drho2; point.epsilon = std::exp(log_epsilon_point); point.depsilon_drho = point.epsilon * log_epsilon * dsolvent_drho; } diff --git a/source/source_hamilt/module_sccs/sccs_cavity.h b/source/source_hamilt/module_sccs/sccs_cavity.h index af3e4eeb9d..ae4d987eda 100644 --- a/source/source_hamilt/module_sccs/sccs_cavity.h +++ b/source/source_hamilt/module_sccs/sccs_cavity.h @@ -15,6 +15,7 @@ struct CavityPoint { double solute = 0.0; double dsolute_drho = 0.0; + double d2solute_drho2 = 0.0; double epsilon = 0.0; double depsilon_drho = 0.0; }; diff --git a/source/source_hamilt/module_sccs/sccs_coul_0d.cpp b/source/source_hamilt/module_sccs/sccs_coul_0d.cpp new file mode 100644 index 0000000000..615c69c285 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_coul_0d.cpp @@ -0,0 +1,32 @@ +#include "sccs_coul_0d.h" + +#include "sccs_moments.h" + +#include "source_basis/module_pw/pw_basis.h" + +namespace ModuleSccs +{ +Pcc0dCoulombOperator::Pcc0dCoulombOperator( + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector>& relative_positions, + const elecstate::Pcc0dParameters& parameters) + : basis_(basis), periodic_(basis, tpiba), relative_positions_(relative_positions), parameters_(parameters) +{ +} + +bool Pcc0dCoulombOperator::has_correction() const +{ + return true; +} + +void Pcc0dCoulombOperator::apply_potential(const std::vector& charge, std::vector& potential) +{ + const elecstate::ChargeMoments moments = pool_moments(charge, relative_positions_, basis_); + periodic_.apply_potential(charge, potential); + for (int ir = 0; ir < basis_.nrxx; ++ir) + { + potential[ir] += elecstate::pcc_0d_potential(moments, relative_positions_[ir], parameters_); + } +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_coul_0d.h b/source/source_hamilt/module_sccs/sccs_coul_0d.h new file mode 100644 index 0000000000..1c90d7afd3 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_coul_0d.h @@ -0,0 +1,30 @@ +#ifndef SCCS_COUL_0D_H +#define SCCS_COUL_0D_H + +#include "sccs_pw_coulomb.h" +#include "source_estate/pcc_0d.h" + +namespace ModuleSccs +{ +// Grid positions are minimum-image displacements about a common, fixed origin +// supplied by the cell layer. Only grid-charge moments enter this operator; +// the host's point-ion PCC energy must not be added here. +class Pcc0dCoulombOperator : public CoulombOperator +{ +public: + Pcc0dCoulombOperator(const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector>& relative_positions, + const elecstate::Pcc0dParameters& parameters); + bool has_correction() const override; + void apply_potential(const std::vector& charge, std::vector& potential) override; + +private: + const ModulePW::PW_Basis& basis_; + PeriodicCoulombOperator periodic_; + const std::vector> relative_positions_; + const elecstate::Pcc0dParameters parameters_; +}; +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_coul_2d.cpp b/source/source_hamilt/module_sccs/sccs_coul_2d.cpp new file mode 100644 index 0000000000..64bd937523 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_coul_2d.cpp @@ -0,0 +1,35 @@ +#include "sccs_coul_2d.h" + +#include "sccs_moments.h" + +#include "source_basis/module_pw/pw_basis.h" + +namespace ModuleSccs +{ +Pcc2dCoulombOperator::Pcc2dCoulombOperator( + const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector>& projected_positions, + const ModuleBase::Vector3& normal, + const elecstate::Pcc2dParameters& parameters) + : basis_(basis), periodic_(basis, tpiba), projected_positions_(projected_positions), + normal_(normal), parameters_(parameters) +{ +} + +bool Pcc2dCoulombOperator::has_correction() const +{ + return true; +} + +void Pcc2dCoulombOperator::apply_potential(const std::vector& charge, std::vector& potential) +{ + const elecstate::ChargeMoments moments = pool_moments(charge, projected_positions_, basis_); + periodic_.apply_potential(charge, potential); + for (int ir = 0; ir < basis_.nrxx; ++ir) + { + const double coordinate = projected_positions_[ir] * normal_; + potential[ir] += elecstate::pcc_2d_potential(moments, coordinate, normal_, parameters_); + } +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_coul_2d.h b/source/source_hamilt/module_sccs/sccs_coul_2d.h new file mode 100644 index 0000000000..5b1f3a333c --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_coul_2d.h @@ -0,0 +1,32 @@ +#ifndef SCCS_COUL_2D_H +#define SCCS_COUL_2D_H + +#include "sccs_pw_coulomb.h" +#include "source_estate/pcc_2d.h" + +namespace ModuleSccs +{ +// Positions are minimum-image normal projections about a fixed origin supplied +// by the cell layer. Only grid-charge moments enter this operator; the host +// supplies the point-ion vacuum PCC energy separately. +class Pcc2dCoulombOperator : public CoulombOperator +{ +public: + Pcc2dCoulombOperator(const ModulePW::PW_Basis& basis, + double tpiba, + const std::vector>& projected_positions, + const ModuleBase::Vector3& normal, + const elecstate::Pcc2dParameters& parameters); + bool has_correction() const override; + void apply_potential(const std::vector& charge, std::vector& potential) override; + +private: + const ModulePW::PW_Basis& basis_; + PeriodicCoulombOperator periodic_; + const std::vector> projected_positions_; + const ModuleBase::Vector3 normal_; + const elecstate::Pcc2dParameters parameters_; +}; +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_coulomb.h b/source/source_hamilt/module_sccs/sccs_coulomb.h new file mode 100644 index 0000000000..0878021579 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_coulomb.h @@ -0,0 +1,19 @@ +#ifndef SCCS_COULOMB_H +#define SCCS_COULOMB_H + +#include + +namespace ModuleSccs +{ +// Potential in Hartree for signed charge density in e/Bohr^3. Applications +// are collective over the PW pool. +class CoulombOperator +{ +public: + virtual ~CoulombOperator() {} + virtual bool has_correction() const = 0; + virtual void apply_potential(const std::vector& charge, std::vector& potential) = 0; +}; +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp index cf7ce44ba5..19636b61ac 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -19,10 +19,21 @@ void eval_functional(const std::vector& charge, const ModulePW::PW_Basis& basis, double tpiba, FunctionalResult& result) +{ + PeriodicCoulombOperator coulomb(basis, tpiba); + eval_functional(charge, response, config, basis, tpiba, coulomb, result); +} + +void eval_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + CoulombOperator& coulomb, + FunctionalResult& result) { const std::size_t size = charge.size(); const double dv = basis.omega / basis.nxyz; - PeriodicCoulombOperator coulomb(basis, tpiba); std::vector vacuum; coulomb.apply_potential(charge, vacuum); FunctionalResult candidate; diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h index 5d6f7d04b9..9b9111abc8 100644 --- a/source/source_hamilt/module_sccs/sccs_functional.h +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -10,6 +10,8 @@ class PW_Basis; namespace ModuleSccs { +class CoulombOperator; + struct SccsConfig; struct SccsResponse; @@ -31,6 +33,15 @@ void eval_functional(const std::vector& charge, const ModulePW::PW_Basis& basis, double tpiba, FunctionalResult& result); + +// Use the same operator for the dielectric response and vacuum subtraction. +void eval_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + CoulombOperator& coulomb, + FunctionalResult& result); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/sccs_moments.cpp b/source/source_hamilt/module_sccs/sccs_moments.cpp new file mode 100644 index 0000000000..4a2a18b0a2 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_moments.cpp @@ -0,0 +1,25 @@ +#include "sccs_moments.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" + +namespace ModuleSccs +{ +elecstate::ChargeMoments pool_moments(const std::vector& charge, + const std::vector>& positions, + const ModulePW::PW_Basis& basis) +{ + const double dv = basis.omega / basis.nxyz; + const int count = charge.size(); + const double* charge_data = charge.data(); + const ModuleBase::Vector3* position_data = positions.data(); + elecstate::ChargeMoments moments = elecstate::charge_moments(charge_data, position_data, count, dv); + double values[5] = {moments.charge, moments.dipole.x, moments.dipole.y, moments.dipole.z, + moments.second_moment}; + Parallel_Reduce::reduce_pool(values, 5); + moments.charge = values[0]; + moments.dipole = ModuleBase::Vector3(values[1], values[2], values[3]); + moments.second_moment = values[4]; + return moments; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_moments.h b/source/source_hamilt/module_sccs/sccs_moments.h new file mode 100644 index 0000000000..fcc7c4d2b2 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_moments.h @@ -0,0 +1,22 @@ +#ifndef SCCS_MOMENTS_H +#define SCCS_MOMENTS_H + +#include "source_estate/charge_moments.h" + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +// Moments of a local grid charge about the common origin of positions, +// reduced over the PW pool, so every rank holds the same values. +elecstate::ChargeMoments pool_moments(const std::vector& charge, + const std::vector>& positions, + const ModulePW::PW_Basis& basis); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h index 0c01362923..bcc805af4a 100644 --- a/source/source_hamilt/module_sccs/sccs_parameters.h +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -7,11 +7,15 @@ namespace ModuleSccs { +enum class Boundary { Periodic, Pcc0d, Pcc2d }; + enum class Preset { Custom, Vacuum, WaterNeutral, WaterCation, WaterAnion }; struct SccsConfig { CavityParameters cavity; + Boundary boundary = Boundary::Periodic; + int pcc_2d_axis = 2; double surface_tension = 0.0; // Hartree/Bohr^2 double pressure = 0.0; // Hartree/Bohr^3 double surface_regularization = 0.0; diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp index 62fbebabd2..36f7e678a5 100644 --- a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp @@ -10,6 +10,11 @@ PeriodicCoulombOperator::PeriodicCoulombOperator(const ModulePW::PW_Basis& basis { } +bool PeriodicCoulombOperator::has_correction() const +{ + return false; +} + void PeriodicCoulombOperator::apply_potential(const std::vector& charge, std::vector& potential) { reciprocal_work_.resize(basis_.npw); diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h index 41245a4560..adcce673e3 100644 --- a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h @@ -1,6 +1,8 @@ #ifndef SCCS_PW_COULOMB_H #define SCCS_PW_COULOMB_H +#include "sccs_coulomb.h" + #include #include @@ -14,11 +16,12 @@ namespace ModuleSccs // Potential in Hartree for a signed charge density in e/Bohr^3. Calls are // collective over the ABACUS pool used to initialize basis. G=0 is removed: // this is the periodic Coulomb inverse with a uniform background. -class PeriodicCoulombOperator +class PeriodicCoulombOperator : public CoulombOperator { public: PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba); - void apply_potential(const std::vector& charge, std::vector& potential); + bool has_correction() const override; + void apply_potential(const std::vector& charge, std::vector& potential) override; private: const ModulePW::PW_Basis& basis_; diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp index b9d52725f6..e7e5cc6b40 100644 --- a/source/source_hamilt/module_sccs/sccs_response.cpp +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -1,4 +1,5 @@ #include "sccs_response.h" +#include "sccs_cav_deriv.h" #include "sccs_parameters.h" #include "sccs_pw_coulomb.h" @@ -18,58 +19,29 @@ namespace ModuleSccs { namespace { -void prepare_cavity(const std::vector& density, +// sqrt(eps) v = C_PCC(s), s = (q - f v)/sqrt(eps). The far-field +// screening charge is int(s)/sqrt(eps_bulk) - int(q), including its physical sign. +void finish_open_bc(const std::vector& charge, + const std::vector& coefficient, + const std::vector& invsqrt, const CavityParameters& cavity, const ModulePW::PW_Basis& basis, - double tpiba, - SccsResponse& response, - std::vector& coefficient) + SccsResponse& result) { - const std::size_t size = density.size(); - response.solute.resize(size); - response.dsolute_drho.resize(size); - response.epsilon.resize(size); - response.depsilon_drho.resize(size); - response.grad_log_epsilon.resize(size); - for (std::size_t i = 0; i < size; ++i) - { - const CavityPoint point = evaluate_cavity(density[i], cavity); - response.solute[i] = point.solute; - response.dsolute_drho[i] = point.dsolute_drho; - response.epsilon[i] = point.epsilon; - response.depsilon_drho[i] = point.depsilon_drho; - } - - std::vector> density_g(basis.npw); - basis.real2recip(density.data(), density_g.data()); - std::vector> gradient(size); - std::vector laplacian(size); - // These XC helpers use G in units of tpiba and ABACUS's normalized FFT. - XC_Functional::grad_rho(density_g.data(), gradient.data(), &basis, tpiba); - XC_Functional::laplacian_rho(density_g.data(), laplacian.data(), &basis, tpiba); - const double density_ratio = cavity.density_max / cavity.density_min; - const double width = std::log(density_ratio); - const double log_bulk = std::log(cavity.epsilon_bulk); - coefficient.resize(size); + 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) { - double second_log = 0.0; - if (density[i] > cavity.density_min && density[i] < cavity.density_max) - { - const double local_ratio = cavity.density_max / density[i]; - const double x = std::log(local_ratio) / width; - const double angle = ModuleBase::TWO_PI * x; - second_log = log_bulk * (1.0 - std::cos(angle) + ModuleBase::TWO_PI * std::sin(angle) / width) - / (width * density[i] * density[i]); - } - const double first_log = response.depsilon_drho[i] / response.epsilon[i]; - const ModuleBase::Vector3& local_gradient = gradient[i]; - const double gradient_square = local_gradient * local_gradient; - response.grad_log_epsilon[i] = local_gradient * first_log; - const double lap_log = first_log * laplacian[i] + second_log * gradient_square; - coefficient[i] = response.epsilon[i] * (0.5 * lap_log + 0.25 * first_log * first_log * gradient_square) - / ModuleBase::FOUR_PI; + source_sum += (charge[i] - coefficient[i] * potential[i]) * invsqrt[i]; + solute_sum += charge[i]; } + Parallel_Reduce::reduce_pool(source_sum); + Parallel_Reduce::reduce_pool(solute_sum); + const double volume_element = basis.omega / basis.nxyz; + const double bulk_invsqrt = 1.0 / std::sqrt(cavity.epsilon_bulk); + result.far_field_polarization_charge = (source_sum * bulk_invsqrt - solute_sum) * volume_element; } // Pool-reduced norms, so every rank takes the same convergence decision. A @@ -114,7 +86,7 @@ double grid_dot(const std::vector& left, // P r = eps^-1/2 G eps^-1/2 r; only FFT scratch survives an application. void apply_precond(const std::vector& rhs, const std::vector& invsqrt, - PeriodicCoulombOperator& coulomb, + CoulombOperator& coulomb, std::vector& weighted, std::vector& value) { @@ -137,7 +109,7 @@ struct SqrtCgOperator const std::vector& coefficient; const std::vector& invsqrt; const ModulePW::PW_Basis& basis; - PeriodicCoulombOperator& coulomb; + CoulombOperator& coulomb; }; // Potential and charge residual of the sqrt-CG, with preconditioner scratch. @@ -239,32 +211,46 @@ void iterate_sqrt_cg(const SqrtCgOperator& system, } } -// Keep the unshifted solution for warm starts, then store the zero-mean -// potential, its gradient and the cavity potential. -void finish_response(const ModulePW::PW_Basis& basis, +// Keep the unshifted solution for warm starts, fix the gauge of a periodic +// cell and store the potential, its gradient and the boundary potentials. +void finish_response(const SqrtCgOperator& system, + const std::vector& charge, + const CavityParameters& cavity, double tpiba, std::vector& potential, SccsResponse& response) { - PolarizationResult& polarization = response.polarization; + const ModulePW::PW_Basis& basis = system.basis; + const std::vector& coefficient = system.coefficient; + const std::vector& invsqrt = system.invsqrt; + const bool open_boundary = system.coulomb.has_correction(); const std::size_t size = potential.size(); + PolarizationResult& polarization = response.polarization; response.restart_potential = potential; - double mean = 0.0; - for (double value : potential) + // Preserve the physical gauge fixed by a boundary correction. + if (!open_boundary) { - mean += value; - } - Parallel_Reduce::reduce_pool(mean); - mean /= basis.nxyz; - for (double& value : potential) - { - value -= mean; + double mean = 0.0; + for (double value : potential) + { + mean += value; + } + Parallel_Reduce::reduce_pool(mean); + mean /= basis.nxyz; + for (double& value : potential) + { + value -= mean; + } } std::vector> potential_g(basis.npw); basis.real2recip(potential.data(), potential_g.data()); polarization.gradient.resize(size); XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); polarization.potential.swap(potential); + if (open_boundary) + { + finish_open_bc(charge, coefficient, invsqrt, cavity, basis, response); + } // Continuum derivative, as in Environ: exact for the discrete energy only // when |grad v|^2 is resolved by the FFT grid. response.cavity_potential.resize(size); @@ -285,22 +271,36 @@ void solve_response(const std::vector& density, const ModulePW::PW_Basis& basis, double tpiba, SccsResponse& result) +{ + PeriodicCoulombOperator coulomb(basis, tpiba); + solve_response(density, charge, cavity, solver, initial_potential, basis, tpiba, coulomb, result); +} + +void solve_response(const std::vector& density, + const std::vector& charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + CoulombOperator& coulomb, + SccsResponse& result) { // Determine start mode collectively, including ranks with no real-space points. double initial_count = initial_potential.size(); Parallel_Reduce::reduce_pool(initial_count); const bool warm_start = initial_count != 0.0; + const bool open_boundary = coulomb.has_correction(); SccsResponse candidate; std::vector coefficient; - prepare_cavity(density, cavity, basis, tpiba, candidate, coefficient); + make_cav_derivs(density, cavity, basis, tpiba, open_boundary, candidate, coefficient); const std::size_t size = density.size(); std::vector invsqrt(size); for (std::size_t i = 0; i < size; ++i) { invsqrt[i] = 1.0 / std::sqrt(candidate.epsilon[i]); } - PeriodicCoulombOperator coulomb(basis, tpiba); const SqrtCgOperator system = {coefficient, invsqrt, basis, coulomb}; SqrtCgState state; state.residual = charge; @@ -312,7 +312,8 @@ void solve_response(const std::vector& density, try_warm_start(system, charge, initial_potential, state, polarization); } iterate_sqrt_cg(system, solver, state, polarization); - finish_response(basis, tpiba, state.potential, candidate); + std::vector& potential = state.potential; + finish_response(system, charge, cavity, tpiba, potential, candidate); result = std::move(candidate); } } // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h index 71be6fcf9c..076ffbb645 100644 --- a/source/source_hamilt/module_sccs/sccs_response.h +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -13,6 +13,8 @@ class PW_Basis; namespace ModuleSccs { +class CoulombOperator; + struct PolarizationSolverParameters; struct PolarizationResult @@ -21,7 +23,7 @@ struct PolarizationResult double residual_rms = 0.0; double residual_max = 0.0; bool warm_started = false; - std::vector potential; // Hartree, shifted to periodic zero mean + std::vector potential; // Hartree; zero mean only for periodic cells std::vector> gradient; // Hartree/Bohr }; @@ -35,6 +37,7 @@ struct SccsResponse PolarizationResult polarization; std::vector cavity_potential; // -epsilon' |grad v|^2/(8 pi), Hartree std::vector restart_potential; // unshifted solution for explicit warm starts + double far_field_polarization_charge = 0.0; // open-boundary screened charge minus solute charge }; // Periodic sqrt-CG response with chain-rule cavity derivatives, as Environ's @@ -51,6 +54,17 @@ void solve_response(const std::vector& cavity_density, const ModulePW::PW_Basis& basis, double tpiba, SccsResponse& result); + +// Use the same operator for the dielectric response and vacuum subtraction. +void solve_response(const std::vector& cavity_density, + const std::vector& solute_charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + CoulombOperator& coulomb, + SccsResponse& result); } // namespace ModuleSccs #endif diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt index 6801ed8760..deb7efe830 100644 --- a/source/source_hamilt/module_sccs/test/CMakeLists.txt +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -6,11 +6,22 @@ AddTest(TARGET MODULE_HAMILT_sccs_parameters AddTest(TARGET MODULE_HAMILT_sccs_poisson LIBS parameter planewave device base - SOURCES test_sccs_pw_coulomb.cpp test_sccs_ion_charge.cpp test_sccs_ion_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + SOURCES test_sccs_cav_deriv.cpp test_sccs_coul_0d.cpp test_sccs_coul_2d.cpp test_sccs_pw_coulomb.cpp test_sccs_ion_charge.cpp test_sccs_ion_force.cpp test_sccs_response.cpp test_sccs_functional.cpp ../../../source_estate/test/test_pot_sccs.cpp + ../../../source_estate/test/test_sccs_pcc_0d.cpp + ../../../source_estate/test/test_sccs_pcc_2d.cpp ../../../source_estate/module_pot/pot_sccs.cpp - ../../../source_cell/cell_tools.cpp ../../../source_cell/atom_spec.cpp + ../../../source_estate/module_pot/sccs_pcc_0d.cpp + ../../../source_estate/module_pot/sccs_pcc_2d.cpp + ../../../source_estate/module_pot/pot_pcc.cpp + ../../../source_estate/pcc_2d.cpp + ../../../source_cell/cell_geometry.cpp ../../../source_cell/cell_tools.cpp ../../../source_cell/atom_spec.cpp ../../../source_cell/atom_pseudo.cpp ../../../source_cell/pseudo.cpp + ../../../source_estate/charge_moments.cpp ../../../source_estate/pcc_0d.cpp + ../sccs_coul_0d.cpp + ../sccs_coul_2d.cpp + ../sccs_moments.cpp + ../sccs_cav_deriv.cpp ../sccs_functional.cpp ../sccs_parameters.cpp ../sccs_ion_charge.cpp ../sccs_ion_force.cpp ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp diff --git a/source/source_hamilt/module_sccs/test/sccs_test.h b/source/source_hamilt/module_sccs/test/sccs_test.h index 6a1570cca4..c79017964e 100644 --- a/source/source_hamilt/module_sccs/test/sccs_test.h +++ b/source/source_hamilt/module_sccs/test/sccs_test.h @@ -11,6 +11,21 @@ namespace SccsTest { +// Model charges of the open-boundary tests, mean + amplitude * cosine mode. +// The mean gives a monopole and the mode has an analytic potential; the +// operators are linear, so only the shape matters, and the values are of +// valence-density size with the density kept positive. +const double charge_mean = 1e-3; +const double charge_amplitude = 2e-4; +// A second charge, not proportional to the first, for symmetry checks. +const double probe_mean = 3e-4; +const double probe_amplitude = 1e-4; +// Cavity density inside the water-neutral switching range (1e-4 to 5e-3), so +// the dielectric varies in space; the small amplitude keeps it mildly +// nonuniform and the sqrt-CG short. +const double cavity_mean = 1e-3; +const double cavity_amplitude = 2e-5; + extern int pool_size; extern int pool_rank; @@ -21,12 +36,18 @@ class PwTest : public testing::Test void SetUp() override { -#ifdef __MPI - basis.initmpi(pool_size, pool_rank, POOL_WORLD); -#endif const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); + set_up_basis(lattice); + } + + // Rows of lattice in units of length; ecut 20 Ry. + void set_up_basis(const ModuleBase::Matrix3& lattice) + { +#ifdef __MPI + basis.initmpi(pool_size, pool_rank, POOL_WORLD); +#endif basis.initgrids(length, lattice, 20.0); basis.initparameters(false, 20.0, 1, false); basis.setuptransform(); 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 new file mode 100644 index 0000000000..0c3cfd6eba --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_cav_deriv.cpp @@ -0,0 +1,84 @@ +#include "sccs_test.h" +#include "../sccs_cav_deriv.h" +#include "../sccs_cavity.h" +#include "../sccs_response.h" + +#include "source_base/parallel_reduce.h" + +using SccsCavityDerivativesTest = SccsTest::PwTest; + +TEST_F(SccsCavityDerivativesTest, OpenBoundaryCoefficientUsesSwitchingFunctionFft) +{ + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + // Manufacture a single-shell boundary rather than deriving the expected + // coefficient with the same FFT helpers as production. + const double amplitude = 0.45; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double solute = 0.5 + amplitude * mode[ir]; + double lower = cavity.density_min; + double upper = cavity.density_max; + for (int iteration = 0; iteration < 64; ++iteration) + { + const double middle = 0.5 * (lower + upper); + const ModuleSccs::CavityPoint point = ModuleSccs::evaluate_cavity(middle, cavity); + if (point.solute < solute) + { + lower = middle; + } + else + { + upper = middle; + } + } + density[ir] = 0.5 * (lower + upper); + } + const bool open_boundary = true; + ModuleSccs::SccsResponse response; + std::vector coefficient; + ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + response, coefficient); + const double log_bulk = std::log(cavity.epsilon_bulk); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const double angle = ModuleBase::TWO_PI * ix / basis.nx; + const double solute = 0.5 + amplitude * mode[ir]; + const double log_epsilon = log_bulk * (1.0 - solute); + const double epsilon = std::exp(log_epsilon); + const double gradient = -amplitude * tpiba * std::sin(angle); + const double laplacian = -amplitude * tpiba * tpiba * mode[ir]; + const double expected = epsilon * (-0.5 * log_bulk * laplacian + + 0.25 * log_bulk * log_bulk * gradient * gradient) + / ModuleBase::FOUR_PI; + EXPECT_NEAR(response.solute[ir], solute, 1e-14); + EXPECT_NEAR(coefficient[ir], expected, 1e-12); + const double expected_gradient = -log_bulk * gradient; + EXPECT_NEAR(response.grad_log_epsilon[ir].x, expected_gradient, 1e-12); + EXPECT_NEAR(response.grad_log_epsilon[ir].y, 0.0, 1e-12); + EXPECT_NEAR(response.grad_log_epsilon[ir].z, 0.0, 1e-12); + } + const bool periodic_boundary = false; + ModuleSccs::SccsResponse periodic; + std::vector periodic_coefficient; + ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, periodic_boundary, + periodic, periodic_coefficient); + double difference = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double coefficient_difference = coefficient[ir] - periodic_coefficient[ir]; + const double local_difference = std::abs(coefficient_difference); + if (local_difference > difference) + { + difference = local_difference; + } + } + Parallel_Reduce::reduce_max_pool(basis.poolnproc, difference); + EXPECT_GT(difference, 1e-8); + +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp index c4d8fc15c4..18f641d400 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp @@ -64,8 +64,8 @@ TEST_F(SccsCavityTest, LimitsThresholdsAndVacuum) EXPECT_DOUBLE_EQ(vacuum.depsilon_drho, 0.0); } -// Analytic values at the logarithmic midpoint and finite differences across -// the transition. +// Analytic values at the logarithmic midpoint and finite differences of the +// first and second derivatives across the transition. TEST_F(SccsCavityTest, MidpointAndDerivativesAcrossTransition) { const double product = parameters.density_min * parameters.density_max; @@ -98,6 +98,9 @@ TEST_F(SccsCavityTest, MidpointAndDerivativesAcrossTransition) upper = ModuleSccs::evaluate_cavity(upper_density, parameters); const double solute_fd = (upper.solute - lower.solute) / (2.0 * step); const double epsilon_fd = (upper.epsilon - lower.epsilon) / (2.0 * step); + const double second_fd = (upper.dsolute_drho - lower.dsolute_drho) / (2.0 * step); + const double second_scale = std::abs(point.d2solute_drho2); + EXPECT_NEAR(point.d2solute_drho2, second_fd, 1e-6 * second_scale); const double solute_ratio = solute_fd / point.dsolute_drho; const double epsilon_ratio = epsilon_fd / point.depsilon_drho; EXPECT_NEAR(solute_ratio, 1.0, 1e-8); 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 new file mode 100644 index 0000000000..9e3d9b1cf3 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_0d.cpp @@ -0,0 +1,220 @@ +#include "sccs_test.h" +#include "../sccs_coul_0d.h" +#include "../sccs_cav_deriv.h" +#include "../sccs_functional.h" +#include "../sccs_parameters.h" +#include "../sccs_response.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_geometry.h" + + +class SccsPcc0dCoulombTest : public SccsTest::PwTest +{ +protected: + void SetUp() override + { + SccsTest::PwTest::SetUp(); + const ModuleBase::Matrix3 lattice; + unitcell::OrthogonalCell cell; + // Same cell tolerance as PotPcc. + ASSERT_TRUE(unitcell::make_orthogonal_cell(lattice, length, 1e-10, cell)); + ASSERT_TRUE(elecstate::make_pcc_0d_parameters(cell, 1e-10, parameters)); + const double half_length = 0.5 * length; + cell.origin = ModuleBase::Vector3(half_length, half_length, half_length); + ModulePW::grid_positions(basis, lattice, length, positions); + for (ModuleBase::Vector3& position : positions) + { + position = unitcell::relative_position(position, cell); + } + } + + double integrate_product(const std::vector& left, const std::vector& right) const + { + double value = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + value += left[ir] * right[ir]; + } + Parallel_Reduce::reduce_pool(value); + return value * basis.omega / basis.nxyz; + } + + elecstate::Pcc0dParameters parameters; + std::vector> positions; +}; + +TEST_F(SccsPcc0dCoulombTest, CorrectionMatchesExistingPccEnergyAndIsSymmetric) +{ + std::vector charge = cosine_mode(0); + std::vector direction = cosine_mode(1); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + charge[ir] = SccsTest::charge_mean + SccsTest::charge_amplitude * charge[ir]; + direction[ir] = SccsTest::probe_mean + SccsTest::probe_amplitude * direction[ir]; + } + ModuleSccs::PeriodicCoulombOperator periodic(basis, tpiba); + ModuleSccs::Pcc0dCoulombOperator corrected(basis, tpiba, positions, parameters); + EXPECT_TRUE(corrected.has_correction()); + std::vector periodic_potential; + std::vector potential; + std::vector direction_potential; + periodic.apply_potential(charge, periodic_potential); + corrected.apply_potential(charge, potential); + corrected.apply_potential(direction, direction_potential); + const double charge_direction = integrate_product(charge, direction_potential); + const double direction_charge = integrate_product(direction, potential); + EXPECT_NEAR(charge_direction, direction_charge, 1e-13); + + const double dv = basis.omega / basis.nxyz; + const double* charge_data = charge.data(); + const ModuleBase::Vector3* position_data = positions.data(); + elecstate::ChargeMoments moments = elecstate::charge_moments(charge_data, position_data, basis.nrxx, dv); + Parallel_Reduce::reduce_pool(moments.charge); + Parallel_Reduce::reduce_pool(moments.dipole.x); + Parallel_Reduce::reduce_pool(moments.dipole.y); + Parallel_Reduce::reduce_pool(moments.dipole.z); + Parallel_Reduce::reduce_pool(moments.second_moment); + std::vector correction(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + correction[ir] = potential[ir] - periodic_potential[ir]; + const double expected = elecstate::pcc_0d_potential(moments, positions[ir], parameters); + EXPECT_NEAR(correction[ir], expected, 1e-14); + } + const double grid_energy = 0.5 * integrate_product(charge, correction); + const double moment_energy = elecstate::pcc_0d_energy(moments, parameters); + EXPECT_NEAR(grid_energy, moment_energy, 1e-13); +} + +TEST_F(SccsPcc0dCoulombTest, ChargedUniformDielectricPreservesGaugeAndVacuumSubtraction) +{ + std::vector charge = cosine_mode(0); + std::vector direction = cosine_mode(1); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + charge[ir] = SccsTest::charge_mean + SccsTest::charge_amplitude * charge[ir]; + // Perturbation direction for the finite-difference check, not proportional to charge. + direction[ir] = 1e-3 + 3e-4 * direction[ir]; + } + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, parameters); + const std::vector density(basis.nrxx, 0.0); + const std::vector cold; + std::vector vacuum; + coulomb.apply_potential(charge, vacuum); + const std::vector ones(basis.nrxx, 1.0); + const double vacuum_mean = integrate_product(vacuum, ones) / basis.omega; + EXPECT_GT(std::abs(vacuum_mean), 1e-3); + ModuleSccs::PolarizationSolverParameters solver; + const double dielectrics[] = {1.0, 5.0}; + for (double epsilon : dielectrics) + { + 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, + 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, + coulomb, result); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double dielectric_potential = vacuum[ir] / epsilon; + const double reaction = dielectric_potential - vacuum[ir]; + EXPECT_NEAR(response.polarization.potential[ir], dielectric_potential, 1e-13); + EXPECT_NEAR(result.reaction_potential[ir], reaction, 1e-13); + const double expected_electron_potential = -reaction; + EXPECT_NEAR(result.electron_potential[ir], expected_electron_potential, 1e-13); + } + const double expected_energy = 0.5 * (1.0 / epsilon - 1.0) * integrate_product(charge, vacuum); + EXPECT_NEAR(result.reaction_energy, expected_energy, 1e-13); + // Electron potential carries the opposite sign to signed solute charge. + const double predicted = -integrate_product(result.electron_potential, direction); + // Relative perturbation of 1e-4 x 1e-3: the central difference error is below 1e-10. + const double step = 1e-4; + const double signs[] = {-1.0, 1.0}; + double energy[2]; + int index = 0; + for (double sign : signs) + { + std::vector shifted = charge; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + shifted[ir] += sign * step * direction[ir]; + } + ModuleSccs::SccsResponse shifted_response; + ModuleSccs::FunctionalResult shifted_result; + ModuleSccs::solve_response(density, shifted, config.cavity, solver, cold, + basis, tpiba, coulomb, shifted_response); + ModuleSccs::eval_functional(shifted, shifted_response, config, basis, tpiba, + coulomb, shifted_result); + energy[index++] = shifted_result.reaction_energy; + } + const double derivative = (energy[1] - energy[0]) / (2.0 * step); + EXPECT_NEAR(derivative, predicted, 1e-10); + } +} + +TEST_F(SccsPcc0dCoulombTest, NonuniformResponseSatisfiesFixedPointAndFarFieldCharge) +{ + ModuleSccs::CavityParameters cavity; + // Water-neutral switching range. + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = SccsTest::cavity_mean + SccsTest::cavity_amplitude * mode[ir]; + charge[ir] = SccsTest::charge_mean + SccsTest::charge_amplitude * mode[ir]; + } + ModuleSccs::Pcc0dCoulombOperator coulomb(basis, tpiba, positions, parameters); + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(density, charge, cavity, solver, cold, + basis, tpiba, coulomb, response); + EXPECT_GT(response.polarization.iterations, 1); + const bool open_boundary = true; + ModuleSccs::SccsResponse derivatives; + std::vector coefficient; + ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + derivatives, coefficient); + std::vector source(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double invsqrt = 1.0 / std::sqrt(response.epsilon[ir]); + source[ir] = (charge[ir] - coefficient[ir] * response.polarization.potential[ir]) * invsqrt; + } + std::vector reconstructed; + coulomb.apply_potential(source, reconstructed); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = reconstructed[ir] / std::sqrt(response.epsilon[ir]); + EXPECT_NEAR(response.polarization.potential[ir], expected, 1e-10); + } + const std::vector ones(basis.nrxx, 1.0); + const double source_sum = integrate_product(source, ones); + const double solute_sum = integrate_product(charge, ones); + 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, + basis, tpiba, coulomb, restarted); + EXPECT_TRUE(restarted.polarization.warm_started); + EXPECT_LT(restarted.polarization.iterations, response.polarization.iterations); + EXPECT_NEAR(restarted.far_field_polarization_charge, response.far_field_polarization_charge, 1e-10); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(restarted.polarization.potential[ir], response.polarization.potential[ir], 1e-10); + } +} 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 new file mode 100644 index 0000000000..da3ed19379 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_coul_2d.cpp @@ -0,0 +1,201 @@ +#include "sccs_test.h" +#include "../sccs_coul_2d.h" +#include "../sccs_cav_deriv.h" +#include "../sccs_functional.h" +#include "../sccs_parameters.h" +#include "../sccs_response.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_geometry.h" + + +class SccsPcc2dCoulombTest : public SccsTest::PwTest +{ +protected: + void SetUp() override + { + // Oblique periodic plane, with a perpendicular open vector along x. + const ModuleBase::Matrix3 lattice(0.0, 1.0, 0.0, + 0.0, 0.2, 1.0, + 1.5, 0.0, 0.0); + set_up_basis(lattice); + unitcell::SlabCell cell; + ASSERT_TRUE(unitcell::make_slab_cell(lattice, length, 2, 1e-10, cell)); + cell.origin = 0.5 * cell.length; + normal = cell.normal; + parameters = elecstate::make_pcc_2d_parameters(cell); + ModulePW::grid_positions(basis, lattice, length, positions); + for (ModuleBase::Vector3& position : positions) + { + const double coordinate = unitcell::relative_coordinate(position, cell); + position = normal * coordinate; + } + } + + double integrate_product(const std::vector& left, const std::vector& right) const + { + double value = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + value += left[ir] * right[ir]; + } + Parallel_Reduce::reduce_pool(value); + return value * basis.omega / basis.nxyz; + } + + elecstate::Pcc2dParameters parameters; + ModuleBase::Vector3 normal; + std::vector> positions; +}; + +TEST_F(SccsPcc2dCoulombTest, ChargedPlanarModeMatchesAnalyticGauge) +{ + // One cosine mode along the open axis on a uniform part: the mode has the + // analytic potential 4 pi amplitude / k^2 and the uniform part carries the + // monopole, so the charged gauge -pi L Q / (3 A) is tested as well. + const double density_mean = SccsTest::charge_mean; + const double amplitude = SccsTest::charge_amplitude; + double monopole = 0.0; + double dipole = 0.0; + double second = 0.0; + const double plane_volume = basis.omega / basis.nz; + for (int iz = 0; iz < basis.nz; ++iz) + { + const double angle = ModuleBase::TWO_PI * iz / basis.nz; + const double weight = (density_mean + amplitude * std::cos(angle)) * plane_volume; + const double coordinate = parameters.length * iz / basis.nz - 0.5 * parameters.length; + monopole += weight; + dipole += weight * coordinate; + second += weight * coordinate * coordinate; + } + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int iz = ir % basis.nplane + basis.startz_current; + const double angle = ModuleBase::TWO_PI * iz / basis.nz; + charge[ir] = density_mean + amplitude * std::cos(angle); + } + ModuleSccs::Pcc2dCoulombOperator coulomb(basis, tpiba, positions, normal, parameters); + EXPECT_TRUE(coulomb.has_correction()); + std::vector potential; + coulomb.apply_potential(charge, potential); + const double wavevector = ModuleBase::TWO_PI / parameters.length; + const double mode_factor = ModuleBase::FOUR_PI * amplitude / (wavevector * wavevector); + const double constant = -ModuleBase::PI * parameters.length / (3.0 * parameters.area); + const double factor = -2.0 * ModuleBase::PI / basis.omega; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int iz = ir % basis.nplane + basis.startz_current; + const double angle = ModuleBase::TWO_PI * iz / basis.nz; + const double coordinate = positions[ir] * normal; + const double multipole = monopole * coordinate * coordinate - 2.0 * dipole * coordinate + second; + const double expected = mode_factor * std::cos(angle) + constant * monopole + factor * multipole; + // The operator is exact for this charge; the tolerance is round-off. + EXPECT_NEAR(potential[ir], expected, 1e-12); + } + std::vector direction = cosine_mode(0); + for (double& value : direction) + { + value = SccsTest::probe_mean + SccsTest::probe_amplitude * value; + } + std::vector direction_potential; + coulomb.apply_potential(direction, direction_potential); + const double left = integrate_product(charge, direction_potential); + const double right = integrate_product(direction, potential); + EXPECT_NEAR(left, right, 1e-13); +} + +TEST_F(SccsPcc2dCoulombTest, UniformDielectricRetainsChargedGaugeAndReactionEnergy) +{ + std::vector charge = cosine_mode(0); + for (double& value : charge) + { + value = SccsTest::charge_mean + SccsTest::charge_amplitude * value; + } + ModuleSccs::Pcc2dCoulombOperator coulomb(basis, tpiba, positions, normal, parameters); + std::vector vacuum; + coulomb.apply_potential(charge, vacuum); + const std::vector density(basis.nrxx, 0.0); + const std::vector ones(basis.nrxx, 1.0); + const std::vector cold; + ModuleSccs::PolarizationSolverParameters solver; + const double dielectrics[] = {1.0, 5.0}; + for (double epsilon : dielectrics) + { + 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, + basis, tpiba, coulomb, response); + EXPECT_EQ(response.polarization.iterations, 1); + ModuleSccs::FunctionalResult result; + ModuleSccs::eval_functional(charge, response, config, basis, tpiba, + coulomb, result); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected_potential = vacuum[ir] / epsilon; + const double expected_reaction = expected_potential - vacuum[ir]; + EXPECT_NEAR(response.polarization.potential[ir], expected_potential, 1e-13); + EXPECT_NEAR(result.reaction_potential[ir], expected_reaction, 1e-13); + const double expected_electron_potential = -expected_reaction; + EXPECT_NEAR(result.electron_potential[ir], expected_electron_potential, 1e-13); + } + const double expected_energy = 0.5 * (1.0 / epsilon - 1.0) * integrate_product(charge, vacuum); + EXPECT_NEAR(result.reaction_energy, expected_energy, 1e-13); + const double expected_screening = (1.0 / epsilon - 1.0) * integrate_product(charge, ones); + EXPECT_NEAR(response.far_field_polarization_charge, expected_screening, 1e-13); + } +} + +TEST_F(SccsPcc2dCoulombTest, NonuniformResponseSatisfiesFixedPointAndWarmStart) +{ + const std::vector mode = cosine_mode(0); + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = SccsTest::cavity_mean + SccsTest::cavity_amplitude * mode[ir]; + charge[ir] = SccsTest::charge_mean + SccsTest::charge_amplitude * mode[ir]; + } + ModuleSccs::CavityParameters cavity; + // Water-neutral switching range. + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 5.0; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + ModuleSccs::Pcc2dCoulombOperator coulomb(basis, tpiba, positions, normal, parameters); + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_response(density, charge, cavity, solver, cold, + basis, tpiba, coulomb, response); + const bool open_boundary = true; + ModuleSccs::SccsResponse derivatives; + std::vector coefficient; + ModuleSccs::make_cav_derivs(density, cavity, basis, tpiba, open_boundary, + derivatives, coefficient); + std::vector source(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double invsqrt = 1.0 / std::sqrt(response.epsilon[ir]); + source[ir] = (charge[ir] - coefficient[ir] * response.polarization.potential[ir]) * invsqrt; + } + std::vector reconstructed; + coulomb.apply_potential(source, reconstructed); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = reconstructed[ir] / std::sqrt(response.epsilon[ir]); + EXPECT_NEAR(response.polarization.potential[ir], expected, 1e-10); + } + ModuleSccs::SccsResponse restarted; + ModuleSccs::solve_response(density, charge, cavity, solver, response.restart_potential, + basis, tpiba, coulomb, restarted); + EXPECT_TRUE(restarted.polarization.warm_started); + EXPECT_LT(restarted.polarization.iterations, response.polarization.iterations); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(restarted.polarization.potential[ir], response.polarization.potential[ir], 1e-10); + } +} 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 35bb798207..a6c91db1d8 100644 --- a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -2,6 +2,7 @@ #include "../sccs_functional.h" #include "../sccs_parameters.h" #include "../sccs_response.h" +#include "../sccs_pw_coulomb.h" #include "source_base/parallel_reduce.h" using SccsFunctionalTest = SccsTest::PwTest; @@ -82,3 +83,41 @@ TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) const double expected_volume = 0.5 * basis.omega; EXPECT_NEAR(center.volume, expected_volume, 1e-10); } + +TEST_F(SccsFunctionalTest, ExplicitPeriodicOperatorPreservesResponseAndFunctional) +{ + ModuleSccs::SccsConfig config = ModuleSccs::preset_config(ModuleSccs::Preset::WaterNeutral); + const std::vector mode = cosine_mode(0); + std::vector charge(basis.nrxx); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + charge[ir] = 1e-3 * mode[ir]; + density[ir] = 1e-3 + 1e-5 * mode[ir]; + } + ModuleSccs::PolarizationSolverParameters solver; + const std::vector cold; + ModuleSccs::SccsResponse legacy; + ModuleSccs::SccsResponse explicit_response; + ModuleSccs::solve_response(density, charge, config.cavity, solver, cold, + basis, tpiba, legacy); + ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); + EXPECT_FALSE(coulomb.has_correction()); + ModuleSccs::solve_response(density, charge, config.cavity, 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); + EXPECT_EQ(legacy.cavity_potential, explicit_response.cavity_potential); + EXPECT_EQ(legacy.restart_potential, explicit_response.restart_potential); + ModuleSccs::FunctionalResult legacy_result; + ModuleSccs::FunctionalResult explicit_result; + ModuleSccs::eval_functional(charge, legacy, config, basis, tpiba, + legacy_result); + ModuleSccs::eval_functional(charge, explicit_response, config, 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); + EXPECT_DOUBLE_EQ(legacy_result.volume_energy, explicit_result.volume_energy); + EXPECT_EQ(legacy_result.reaction_potential, explicit_result.reaction_potential); + EXPECT_EQ(legacy_result.electron_potential, explicit_result.electron_potential); +} diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index d5d4e71c23..961002946c 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -193,10 +193,11 @@ void ReadInput::item_model() "* 1: the original ABACUS implicit solvation model.\n" "* 2: the self-consistent continuum solvation model (SCCS). Supported scope: CPU KS-DFT " "(esolver_type ksdft, device cpu) with basis_type pw or lcao; calculation scf energies and " - "forces, no stress (cal_stress is rejected); nspin 1 or 2; periodic cells with " - "assume_isolated none; no electric/gate field, DFT-1/2, DeePKS or dm_to_rho. SCCS is meant " - "for neutral cells: a charged cell only prints a warning and runs, and its energy then " - "depends on the cell size because the periodic Poisson solver drops the G = 0 net charge.\n\n" + "forces, no stress (cal_stress is rejected); nspin 1 or 2; assume_isolated none, pcc_0d or " + "pcc_2d; no electric/gate field, DFT-1/2, DeePKS or dm_to_rho. Use assume_isolated pcc_0d " + "or pcc_2d for charged systems: with assume_isolated none a charged cell only prints a " + "warning and runs, and its energy then depends on the cell size because the periodic " + "Poisson solver drops the G = 0 net charge.\n\n" "Legacy Boolean values (true/false) are read as 1/0."; item.default_value = "0"; item.read_value = [](const Input_Item& item, Parameter& para) { diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp index fa19b401f3..816ea7c656 100644 --- a/source/source_io/module_parameter/read_inp_sccs.cpp +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -39,18 +39,37 @@ namespace std::string context_error(const Input_para& input) { std::string error; - if (input.assume_isolated != "none") { error = "SCCS currently requires assume_isolated none"; } - else if (input.device != "cpu" || input.esolver_type != "ksdft") { error = "SCCS requires CPU KS-DFT"; } - else if (input.basis_type != "pw" && input.basis_type != "lcao") { error = "SCCS requires basis_type pw or lcao"; } + if (input.assume_isolated != "none" && input.assume_isolated != "pcc_0d" + && input.assume_isolated != "pcc_2d") + { + error = "SCCS currently supports assume_isolated none, pcc_0d or pcc_2d"; + } + else if (input.device != "cpu" || input.esolver_type != "ksdft") + { + error = "SCCS requires CPU KS-DFT"; + } + else if (input.basis_type != "pw" && input.basis_type != "lcao") + { + error = "SCCS requires basis_type pw or lcao"; + } else if (input.calculation != "scf" || input.cal_stress) - { error = "SCCS currently supports SCF energies and forces without stress"; } - else if (input.nspin != 1 && input.nspin != 2) { error = "SCCS requires nspin 1 or 2"; } + { + error = "SCCS currently supports SCF energies and forces without stress"; + } + else if (input.nspin != 1 && input.nspin != 2) + { + error = "SCCS requires nspin 1 or 2"; + } else if (input.efield_flag || input.gate_flag || input.dfthalf_type != 0 || input.deepks_scf || input.deepks_out_labels || input.deepks_bandgap || input.deepks_v_delta || input.dm_to_rho) - { error = "SCCS does not support electric/gate fields, DFT-1/2, DeePKS or dm_to_rho"; } + { + error = "SCCS does not support electric/gate fields, DFT-1/2, DeePKS or dm_to_rho"; + } else if (!input.vl_in_h || !input.vion_in_h || !input.vh_in_h) - { error = "SCCS requires the local ionic and Hartree potentials in the Hamiltonian"; } + { + error = "SCCS requires the local ionic and Hartree potentials in the Hamiltonian"; + } return error; } diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 79b933f7f3..70e69d9736 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -556,7 +556,7 @@ Available options are: * pcc_2d: self-consistent PCC correction for a slab, open along the lattice vector selected by pcc_2d_axis (default: the third) and periodic along the other two. The open vector must be perpendicular to the periodic plane (rewrite a tilted cell as the equivalent perpendicular one), and the k-point sampling along it must be Gamma only. The monopole term uses the open planar kernel that vanishes on the plane of the charge, -pi*Q*L/(3*A) with L the cell length along the open vector and A the periodic area, so the energy of a charged slab converges with the vacuum size; it is referenced to zero potential on the plane of the charge. ENVIRON uses -pi*Q/(3*L), so charged-slab energies agree with ENVIRON only for A = L^2. * makov-payne, m-p, mp: compute the Makov-Payne correction to the total energy and estimate a corrected vacuum level for eigenvalue alignment. This option is available only for cubic lattices (latname = sc, fcc, or bcc). -pcc_0d and pcc_2d contribute to the energy, the potential and the fixed-cell forces; the correction energy is printed as E_pcc. They require CPU KS-DFT (esolver_type ksdft) with basis_type pw or lcao, calculation scf or relax, and nspin 1 or 2, without efield_flag, gate_flag, cal_stress, DFT-1/2, deepks output or dm_to_rho. +pcc_0d and pcc_2d contribute to the energy, the potential and the fixed-cell forces; the correction energy is printed as E_pcc. They require CPU KS-DFT (esolver_type ksdft) with basis_type pw or lcao, calculation scf or relax, and nspin 1 or 2, without efield_flag, gate_flag, cal_stress, DFT-1/2, deepks output or dm_to_rho. With imp_sol=2, pcc_0d and pcc_2d support SCF energies and forces; imp_sol=1 is not supported. Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995); PCC: O. Andreussi and N. Marzari, Phys. Rev. B 90, 245101 (2014).)"; item.default_value = "none"; @@ -584,7 +584,7 @@ Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995); PCC: O. Andreuss if (inp.assume_isolated == "pcc_0d" || inp.assume_isolated == "pcc_2d") { std::string error; - if (inp.imp_sol) { error = "cannot be combined with imp_sol"; } + if (inp.imp_sol == 1) { error = "cannot be combined with imp_sol=1"; } else if (inp.nspin != 1 && inp.nspin != 2) { error = "requires nspin 1 or 2"; } else if (inp.device != "cpu") { error = "requires device cpu"; } else if (inp.esolver_type != "ksdft") { error = "requires esolver_type ksdft"; } diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 7fd18f1d1f..3b4365dbff 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -144,6 +144,16 @@ TEST_F(InputTest, PccSelectionAndUnsupportedCombinations) input.basis_type = "lcao"; input.calculation = "relax"; check_item(it->second, param); + input.imp_sol = 1; + EXPECT_EXIT(check_item(it->second, param), testing::ExitedWithCode(1), ""); + input.imp_sol = 2; + input.calculation = "scf"; + input.cal_force = false; + check_item(it->second, param); + input.assume_isolated = "pcc_2d"; + check_item(it->second, param); + input.imp_sol = 0; + input.assume_isolated = "pcc_0d"; const Input_para valid = input; input.nspin = 4; 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 01ea319b8b..ea18435b48 100644 --- a/source/source_io/test_serial/test_read_inp_sccs.cpp +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -42,7 +42,18 @@ TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) input.cal_stress = false; input.cal_force = false; input.assume_isolated = "pcc_0d"; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.cal_force = true; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.cal_force = false; + input.assume_isolated = "pcc_2d"; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.cal_force = true; + EXPECT_TRUE(ModuleIO::check_sccs_inp(input, error)); + input.cal_force = false; + input.calculation = "relax"; EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); + input.calculation = "scf"; input.assume_isolated = "none"; input.sccs_rho_min = input.sccs_rho_max; EXPECT_FALSE(ModuleIO::check_sccs_inp(input, error)); diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index e872199795..e395806413 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -128,3 +128,5 @@ nscf_hse_path scf_pcc_0d scf_pcc_2d scf_sccs +scf_sccs_pcc_2d +scf_sccs_pcc_0d diff --git a/tests/01_PW/scf_sccs_pcc_0d/INPUT b/tests/01_PW/scf_sccs_pcc_0d/INPUT new file mode 100644 index 0000000000..f287686f37 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/INPUT @@ -0,0 +1,21 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated pcc_0d +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 diff --git a/tests/01_PW/scf_sccs_pcc_0d/KPT b/tests/01_PW/scf_sccs_pcc_0d/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_pcc_0d/README b/tests/01_PW/scf_sccs_pcc_0d/README new file mode 100644 index 0000000000..ebc67d3501 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/README @@ -0,0 +1 @@ +PW neutral water with SCCS (water-neutral preset) and 0D PCC, ecutwfc 20 and ecutrho 240 diff --git a/tests/01_PW/scf_sccs_pcc_0d/STRU b/tests/01_PW/scf_sccs_pcc_0d/STRU new file mode 100644 index 0000000000..ed974dd43a --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +12.0 0.0 0.0 +0.0 12.0 0.0 +0.0 0.0 12.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +7.433 7.111 6.0 0 0 0 +4.567 7.111 6.0 0 0 0 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/scf_sccs_pcc_0d/result.ref b/tests/01_PW/scf_sccs_pcc_0d/result.ref new file mode 100644 index 0000000000..6b87966cdf --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/result.ref @@ -0,0 +1,7 @@ +etotref -441.6509526968036994 +etotperatomref -147.2169842323 +totalforceref 16.638619 +esolelref -0.3079260694 +esolcavref 0.0807450240 +epccref 0.0166182392 +totaltimeref 4.31 diff --git a/tests/01_PW/scf_sccs_pcc_0d/threshold b/tests/01_PW/scf_sccs_pcc_0d/threshold new file mode 100644 index 0000000000..a125b952d6 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_0d/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/scf_sccs_pcc_2d/INPUT b/tests/01_PW/scf_sccs_pcc_2d/INPUT new file mode 100644 index 0000000000..ad4803dd65 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/INPUT @@ -0,0 +1,22 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated 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 diff --git a/tests/01_PW/scf_sccs_pcc_2d/KPT b/tests/01_PW/scf_sccs_pcc_2d/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/scf_sccs_pcc_2d/README b/tests/01_PW/scf_sccs_pcc_2d/README new file mode 100644 index 0000000000..6b253dd127 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/README @@ -0,0 +1 @@ +PW neutral water with SCCS (water-neutral preset) and 2D PCC open along y, ecutwfc 20 and ecutrho 240 diff --git a/tests/01_PW/scf_sccs_pcc_2d/STRU b/tests/01_PW/scf_sccs_pcc_2d/STRU new file mode 100644 index 0000000000..a9048229a5 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/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_pcc_2d/result.ref b/tests/01_PW/scf_sccs_pcc_2d/result.ref new file mode 100644 index 0000000000..2b06e1c629 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/result.ref @@ -0,0 +1,7 @@ +etotref -441.9953561042844967 +etotperatomref -147.3317853681 +totalforceref 16.524125 +esolelref -0.2690294917 +esolcavref 0.0811367710 +epccref 0.0426195690 +totaltimeref 6.07 diff --git a/tests/01_PW/scf_sccs_pcc_2d/threshold b/tests/01_PW/scf_sccs_pcc_2d/threshold new file mode 100644 index 0000000000..a125b952d6 --- /dev/null +++ b/tests/01_PW/scf_sccs_pcc_2d/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 37a3e8404f..2bdb478b1c 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -40,4 +40,5 @@ scf_upf100 test.sum scf_pcc0d_water scf_pcc2d_water +scf_sccs_pcc2d scf_sccs_water diff --git a/tests/02_NAO_Gamma/scf_sccs_pcc2d/INPUT b/tests/02_NAO_Gamma/scf_sccs_pcc2d/INPUT new file mode 100644 index 0000000000..1931e68c62 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_pcc2d/INPUT @@ -0,0 +1,28 @@ +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 diff --git a/tests/02_NAO_Gamma/scf_sccs_pcc2d/KPT b/tests/02_NAO_Gamma/scf_sccs_pcc2d/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_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_pcc2d/README b/tests/02_NAO_Gamma/scf_sccs_pcc2d/README new file mode 100644 index 0000000000..4bdd18f5d2 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_pcc2d/README @@ -0,0 +1 @@ +gamma-only LCAO neutral water with SCCS (epsilon 1.1) and 2D PCC open along y diff --git a/tests/02_NAO_Gamma/scf_sccs_pcc2d/STRU b/tests/02_NAO_Gamma/scf_sccs_pcc2d/STRU new file mode 100644 index 0000000000..409f8328a3 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_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_pcc2d/result.ref b/tests/02_NAO_Gamma/scf_sccs_pcc2d/result.ref new file mode 100644 index 0000000000..018e0b0287 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_sccs_pcc2d/result.ref @@ -0,0 +1,7 @@ +etotref -466.3788933482496 +etotperatomref -155.4596311161 +totalforceref 3.106598 +esolelref -0.0184368798 +esolcavref 0.0000000000 +epccref 0.0517652102 +totaltimeref 4.05 diff --git a/tests/03_NAO_multik/CASES_CPU.txt b/tests/03_NAO_multik/CASES_CPU.txt index c1a91cb6ce..3744c0193e 100644 --- a/tests/03_NAO_multik/CASES_CPU.txt +++ b/tests/03_NAO_multik/CASES_CPU.txt @@ -71,3 +71,4 @@ scf_deltaspin2 scf_deltaspin4 scf_nspin4_nsoc_hr scf_nspin4_soc_hr +scf_sccs_pcc2d diff --git a/tests/03_NAO_multik/scf_sccs_pcc2d/INPUT b/tests/03_NAO_multik/scf_sccs_pcc2d/INPUT new file mode 100644 index 0000000000..47c69bf2d8 --- /dev/null +++ b/tests/03_NAO_multik/scf_sccs_pcc2d/INPUT @@ -0,0 +1,29 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +calculation scf +basis_type lcao +gamma_only 0 +kpar 2 +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 diff --git a/tests/03_NAO_multik/scf_sccs_pcc2d/KPT b/tests/03_NAO_multik/scf_sccs_pcc2d/KPT new file mode 100644 index 0000000000..51b8d09335 --- /dev/null +++ b/tests/03_NAO_multik/scf_sccs_pcc2d/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 1 2 0 0 0 diff --git a/tests/03_NAO_multik/scf_sccs_pcc2d/README b/tests/03_NAO_multik/scf_sccs_pcc2d/README new file mode 100644 index 0000000000..adc9ef6e42 --- /dev/null +++ b/tests/03_NAO_multik/scf_sccs_pcc2d/README @@ -0,0 +1 @@ +LCAO neutral water with SCCS and 2D PCC, 2x1x2 k mesh and kpar 2 diff --git a/tests/03_NAO_multik/scf_sccs_pcc2d/STRU b/tests/03_NAO_multik/scf_sccs_pcc2d/STRU new file mode 100644 index 0000000000..409f8328a3 --- /dev/null +++ b/tests/03_NAO_multik/scf_sccs_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/03_NAO_multik/scf_sccs_pcc2d/result.ref b/tests/03_NAO_multik/scf_sccs_pcc2d/result.ref new file mode 100644 index 0000000000..9ae78e2297 --- /dev/null +++ b/tests/03_NAO_multik/scf_sccs_pcc2d/result.ref @@ -0,0 +1,7 @@ +etotref -466.427275825557 +etotperatomref -155.4757586085 +totalforceref 3.272637 +esolelref -0.0188949736 +esolcavref 0.0000000000 +epccref 0.0519353453 +totaltimeref 3.90 diff --git a/tests/integrate/validation_tools/props_basic.sh b/tests/integrate/validation_tools/props_basic.sh index 2959876110..ea6ce37248 100644 --- a/tests/integrate/validation_tools/props_basic.sh +++ b/tests/integrate/validation_tools/props_basic.sh @@ -95,8 +95,8 @@ run_basic_props_post_ml(){ # implicit solvation model #-------------------------------------------- if ! test -z "$imp_sol" && { [ "$imp_sol" = 1 ] || [ "$imp_sol" = 2 ]; }; then - esol_el=`grep E_sol_el $running_path | awk '{print $3}'` - esol_cav=`grep E_sol_cav $running_path | awk '{print $3}'` + esol_el=`grep E_sol_el $running_path | tail -1 | awk '{print $3}'` + esol_cav=`grep E_sol_cav $running_path | tail -1 | awk '{print $3}'` echo "esolelref $esol_el" >>$props_result_file echo "esolcavref $esol_cav" >>$props_result_file fi