diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 4e5456899d3..535691df663 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -20,6 +20,7 @@ - [bndpar](#bndpar) - [latname](#latname) - [assume\_isolated](#assume_isolated) + - [pcc\_2d\_axis](#pcc_2d_axis) - [init\_wfc](#init_wfc) - [init\_chg](#init_chg) - [init\_vel](#init_vel) @@ -536,6 +537,16 @@ - [tau](#tau) - [sigma\_k](#sigma_k) - [nc\_k](#nc_k) + - [sccs\_preset](#sccs_preset) + - [sccs\_epsilon](#sccs_epsilon) + - [sccs\_rho\_min](#sccs_rho_min) + - [sccs\_rho\_max](#sccs_rho_max) + - [sccs\_gamma](#sccs_gamma) + - [sccs\_pressure](#sccs_pressure) + - [sccs\_maxiter](#sccs_maxiter) + - [sccs\_tol\_rms](#sccs_tol_rms) + - [sccs\_tol\_max](#sccs_tol_max) + - [sccs\_surface\_eta](#sccs_surface_eta) - [Quasiatomic Orbital (QO) analysis](#quasiatomic-orbital-qo-analysis) - [qo\_switch](#qo_switch) - [qo\_basis](#qo_basis) @@ -769,11 +780,22 @@ Available options are: - none: regular periodic calculation without isolated-system correction. + - pcc_0d: self-consistent point-counter-charge (PCC) correction for a molecule in a cubic cell. The multipoles are taken about the mass-weighted ionic center. + - 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). - Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995). + 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. + + 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). - **Default**: none +### pcc_2d_axis + +- **Type**: Integer +- **Availability**: *[`assume_isolated`](#assume_isolated)==pcc_2d* +- **Description**: Index of the lattice vector along which assume_isolated=pcc_2d is open: 0, 1 or 2 for the first, second or third vector of LATTICE_VECTORS. The selected vector must be perpendicular to the other two, which span the periodic plane; the k-point sampling along it must be Gamma only. +- **Default**: 2 + ### init_wfc - **Type**: Vector of string @@ -4915,14 +4937,14 @@ ### imp_sol -- **Type**: Boolean -- **Description**: Calculate implicit solvation correction -- **Default**: False +- **Type**: Integer +- **Description**: Select 0 for vacuum, 1 for the original ABACUS implicit solvation model, or 2 for SCCS. Stress, external fields and other correction models are not supported. +- **Default**: 0 ### eb_k - **Type**: Real -- **Availability**: *[`imp_sol`](#imp_sol)==true* +- **Availability**: *[`imp_sol`](#imp_sol)==1* - **Description**: The relative permittivity of the bulk solvent, 80 for water - **Default**: 80 @@ -4944,6 +4966,83 @@ - **Description**: The value of the electron density at which the dielectric cavity forms - **Default**: 0.00037 +### sccs_preset + +- **Type**: String +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: SCCS parameter preset. Allowed values: custom, vacuum, water-neutral, water-cation, water-anion. Non-custom presets override sccs_epsilon, sccs_rho_min, sccs_rho_max, sccs_gamma and sccs_pressure; solver controls and surface regularization remain user-controlled. +- **Default**: custom + +### sccs_epsilon + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Bulk dielectric constant >= 1 for the custom preset. The default 1 (no dielectric) follows the Environ env_static_permittivity default for environ_type input. Water presets use 78.3 and vacuum uses 1. +- **Default**: 1.0 + +### sccs_rho_min + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Lower electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 1e-4; water-cation: 2e-4; water-anion: 2.4e-3. +- **Default**: 1e-4 +- **Unit**: e/bohr^3 + +### sccs_rho_max + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Upper electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 5e-3; water-cation: 3.5e-3; water-anion: 1.55e-2. +- **Default**: 5e-3 +- **Unit**: e/bohr^3 + +### sccs_gamma + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Surface coefficient for the custom preset. Water-neutral: 47.9; water-cation: 5; water-anion/vacuum: 0. +- **Default**: 0.0 +- **Unit**: dyn/cm + +### sccs_pressure + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Volume coefficient for the custom preset. Water-neutral: -0.36; water-cation: 0.125; water-anion: 0.45; vacuum: 0. +- **Default**: 0.0 +- **Unit**: GPa + +### sccs_maxiter + +- **Type**: Integer +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Positive maximum number of SCCS sqrt-CG iterations. Failure to satisfy both residual tolerances terminates the calculation. +- **Default**: 200 + +### sccs_tol_rms + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver. +- **Default**: 1e-10 +- **Unit**: e/bohr^3 + +### sccs_tol_max + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver. +- **Default**: 1e-8 +- **Unit**: e/bohr^3 + +### sccs_surface_eta + +- **Type**: Real +- **Availability**: *[`imp_sol`](#imp_sol)==2* +- **Description**: Positive regularization of the SCCS surface gradient norm. +- **Default**: 1e-8 +- **Unit**: bohr^-1 + [back to top](#full-list-of-input-keywords) ## Quasiatomic Orbital (QO) analysis diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 813d25e3b0a..dfc30eddd93 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -178,12 +178,24 @@ parameters: Available options are: * none: regular periodic calculation without isolated-system correction. + * pcc_0d: self-consistent point-counter-charge (PCC) correction for a molecule in a cubic cell. The multipoles are taken about the mass-weighted ionic center. + * 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). - Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995). + 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. + + 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). default_value: none unit: "" availability: "" + - name: pcc_2d_axis + category: System variables + type: Integer + description: | + Index of the lattice vector along which assume_isolated=pcc_2d is open: 0, 1 or 2 for the first, second or third vector of LATTICE_VECTORS. The selected vector must be perpendicular to the other two, which span the periodic plane; the k-point sampling along it must be Gamma only. + default_value: "2" + unit: "" + availability: assume_isolated==pcc_2d - name: init_wfc category: System variables type: Vector of string @@ -4341,10 +4353,10 @@ parameters: availability: "" - name: imp_sol category: Implicit solvation model - type: Boolean + type: Integer description: | - Calculate implicit solvation correction - default_value: "False" + Select 0 for vacuum, 1 for the original ABACUS implicit solvation model, or 2 for SCCS. Stress, external fields and other correction models are not supported. + default_value: "0" unit: "" availability: "" - name: eb_k @@ -4354,7 +4366,7 @@ parameters: The relative permittivity of the bulk solvent, 80 for water default_value: "80" unit: "" - availability: imp_sol==true + availability: imp_sol==1 - name: tau category: Implicit solvation model type: Real @@ -4379,6 +4391,86 @@ parameters: default_value: "0.00037" unit: "" availability: "" + - name: sccs_preset + category: Implicit solvation model + type: String + description: | + SCCS parameter preset. Allowed values: custom, vacuum, water-neutral, water-cation, water-anion. Non-custom presets override sccs_epsilon, sccs_rho_min, sccs_rho_max, sccs_gamma and sccs_pressure; solver controls and surface regularization remain user-controlled. + default_value: custom + unit: "" + availability: imp_sol==2 + - name: sccs_epsilon + category: Implicit solvation model + type: Real + description: | + Bulk dielectric constant >= 1 for the custom preset. The default 1 (no dielectric) follows the Environ env_static_permittivity default for environ_type input. Water presets use 78.3 and vacuum uses 1. + default_value: "1.0" + unit: "" + availability: imp_sol==2 + - name: sccs_rho_min + category: Implicit solvation model + type: Real + description: | + Lower electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 1e-4; water-cation: 2e-4; water-anion: 2.4e-3. + default_value: "1e-4" + unit: e/bohr^3 + availability: imp_sol==2 + - name: sccs_rho_max + category: Implicit solvation model + type: Real + description: | + Upper electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 5e-3; water-cation: 3.5e-3; water-anion: 1.55e-2. + default_value: "5e-3" + unit: e/bohr^3 + availability: imp_sol==2 + - name: sccs_gamma + category: Implicit solvation model + type: Real + description: | + Surface coefficient for the custom preset. Water-neutral: 47.9; water-cation: 5; water-anion/vacuum: 0. + default_value: "0.0" + unit: dyn/cm + availability: imp_sol==2 + - name: sccs_pressure + category: Implicit solvation model + type: Real + description: | + Volume coefficient for the custom preset. Water-neutral: -0.36; water-cation: 0.125; water-anion: 0.45; vacuum: 0. + default_value: "0.0" + unit: GPa + availability: imp_sol==2 + - name: sccs_maxiter + category: Implicit solvation model + type: Integer + description: | + Positive maximum number of SCCS sqrt-CG iterations. Failure to satisfy both residual tolerances terminates the calculation. + default_value: "200" + unit: "" + availability: imp_sol==2 + - name: sccs_tol_rms + category: Implicit solvation model + type: Real + description: | + Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver. + default_value: "1e-10" + unit: e/bohr^3 + availability: imp_sol==2 + - name: sccs_tol_max + category: Implicit solvation model + type: Real + description: | + Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver. + default_value: "1e-8" + unit: e/bohr^3 + availability: imp_sol==2 + - name: sccs_surface_eta + category: Implicit solvation model + type: Real + description: | + Positive regularization of the SCCS surface gradient norm. + default_value: "1e-8" + unit: bohr^-1 + availability: imp_sol==2 - name: vdw_method category: vdW correction type: String diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index cdf3db5fee0..66f9380c0d8 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -729,6 +729,7 @@ target_link_libraries( io_basic io_input surchem + sccs neighbor_search domain_decomposition neighbor diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 62df335c194..4bd5169a9c3 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -49,6 +49,7 @@ VPATH=./src_global:\ ./source_hamilt:\ ./source_hamilt/module_ewald:\ ./source_hamilt/module_surchem:\ +./source_hamilt/module_sccs:\ ./source_hamilt/module_xc:\ ./source_hamilt/module_xc/kernels:\ ./source_pw/module_proj:\ @@ -103,6 +104,7 @@ ${OBJS_PSI_INITIALIZER}\ ${OBJS_PW}\ ${OBJS_RELAXATION}\ ${OBJS_SURCHEM}\ +${OBJS_SCCS}\ ${OBJS_SYMMETRY}\ ${OBJS_XC}\ ${OBJS_IO}\ @@ -221,6 +223,7 @@ OBJS_CELL=atom_pseudo.o\ reciprocal_grid.o\ cell_index.o\ cell_tools.o\ + cell_geometry.o\ record_adj.o\ check_atomic_stru.o\ update_cell.o\ @@ -271,6 +274,11 @@ OBJS_ELECSTAT=elecstate.o\ elecstate_pw_sdft.o\ elecstate_pw_cal_tau.o\ makov_payne.o\ + charge_moments.o\ + pcc_0d.o\ + pcc_2d.o\ + pot_pcc.o\ + pot_sccs.o\ elecstate_op.o\ efield.o\ gatefield.o\ @@ -562,6 +570,7 @@ OBJS_PSI_INITIALIZER=psi_base.o\ OBJS_PW=fft_bundle.o\ fft_cpu.o\ pw_basis.o\ + pw_grid_geometry.o\ pw_basis_k.o\ pw_basis_sup.o\ pw_distributeg.o\ @@ -611,6 +620,14 @@ OBJS_SURCHEM=surchem.o\ minimize_cg.o\ sol_force.o\ +OBJS_SCCS=sccs_cavity.o\ + sccs_functional.o\ + sccs_parameters.o\ + sccs_ionic_charge.o\ + sccs_ionic_force.o\ + sccs_pw_coulomb.o\ + sccs_response.o + OBJS_SYMMETRY=symm_other.o\ symm_rot_spin.o\ symm_basic.o\ @@ -737,6 +754,7 @@ OBJS_IO=module_parameter/input_conv.o\ module_parameter/read_inp_bse.o\ module_parameter/read_inp_deepks.o\ module_parameter/read_inp_model.o\ + module_parameter/read_inp_sccs.o\ module_parameter/read_inp_postproc.o\ module_parameter/read_inp_exx_dftu.o\ module_parameter/read_inp_dfpt.o\ diff --git a/source/source_basis/module_pw/CMakeLists.txt b/source/source_basis/module_pw/CMakeLists.txt index d2abf09b418..12d037dab73 100644 --- a/source/source_basis/module_pw/CMakeLists.txt +++ b/source/source_basis/module_pw/CMakeLists.txt @@ -22,6 +22,7 @@ endif() list(APPEND objects pw_basis.cpp + pw_grid_geometry.cpp pw_basis_k.cpp pw_basis_sup.cpp pw_distributeg.cpp diff --git a/source/source_basis/module_pw/pw_grid_geometry.cpp b/source/source_basis/module_pw/pw_grid_geometry.cpp new file mode 100644 index 00000000000..e9a9cfc2a4e --- /dev/null +++ b/source/source_basis/module_pw/pw_grid_geometry.cpp @@ -0,0 +1,28 @@ +#include "pw_grid_geometry.h" + +#include "pw_basis.h" +#include "source_base/matrix3.h" + +namespace ModulePW +{ +void grid_positions(const PW_Basis& basis, + const ModuleBase::Matrix3& lattice, + const double lattice_scale, + std::vector>& positions) +{ + const ModuleBase::Vector3 a1(lattice.e11, lattice.e12, lattice.e13); + const ModuleBase::Vector3 a2(lattice.e21, lattice.e22, lattice.e23); + const ModuleBase::Vector3 a3(lattice.e31, lattice.e32, lattice.e33); + positions.resize(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const int iy = ir / basis.nplane - ix * basis.ny; + const int iz = ir % basis.nplane + basis.startz_current; + const double fx = static_cast(ix) / basis.nx; + const double fy = static_cast(iy) / basis.ny; + const double fz = static_cast(iz) / basis.nz; + positions[ir] = (a1 * fx + a2 * fy + a3 * fz) * lattice_scale; + } +} +} // namespace ModulePW diff --git a/source/source_basis/module_pw/pw_grid_geometry.h b/source/source_basis/module_pw/pw_grid_geometry.h new file mode 100644 index 00000000000..c1b46e641e1 --- /dev/null +++ b/source/source_basis/module_pw/pw_grid_geometry.h @@ -0,0 +1,25 @@ +#ifndef PW_GRID_GEOMETRY_H +#define PW_GRID_GEOMETRY_H + +#include "source_base/vector3.h" + +#include + +namespace ModuleBase +{ +class Matrix3; +} + +namespace ModulePW +{ +class PW_Basis; + +/// Cartesian coordinates in Bohr of this rank's integer FFT grid nodes. +/// Supports an empty local slab. +void grid_positions(const PW_Basis& basis, + const ModuleBase::Matrix3& lattice, + double lattice_scale, + std::vector>& positions); +} // namespace ModulePW + +#endif diff --git a/source/source_basis/module_pw/test_serial/CMakeLists.txt b/source/source_basis/module_pw/test_serial/CMakeLists.txt index 3ef55984a4d..0ab7d3216a9 100644 --- a/source/source_basis/module_pw/test_serial/CMakeLists.txt +++ b/source/source_basis/module_pw/test_serial/CMakeLists.txt @@ -12,6 +12,7 @@ add_library( ../../../source_base/module_fft/fft_bundle.cpp ../../../source_base/module_fft/fft_cpu.cpp ../pw_basis.cpp + ../pw_grid_geometry.cpp ../pw_basis_k.cpp ../pw_basis_sup.cpp ../pw_distributeg.cpp @@ -23,6 +24,12 @@ add_library( ../pw_transform_k.cpp ) +AddTest( + TARGET MODULE_PW_grid_geometry + LIBS parameter planewave_serial device base + SOURCES test_pw_grid_geometry.cpp +) + AddTest( TARGET MODULE_PW_basis_pw_serial LIBS parameter planewave_serial device base diff --git a/source/source_basis/module_pw/test_serial/test_pw_grid_geometry.cpp b/source/source_basis/module_pw/test_serial/test_pw_grid_geometry.cpp new file mode 100644 index 00000000000..a0aa25ec9ac --- /dev/null +++ b/source/source_basis/module_pw/test_serial/test_pw_grid_geometry.cpp @@ -0,0 +1,43 @@ +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_base/matrix3.h" + +#include + +TEST(PwGridGeometry, UsesGlobalSlabOffsetAndRotatedLattice) +{ + ModulePW::PW_Basis basis("cpu", "double"); + basis.nx = 2; + basis.ny = 2; + basis.nz = 4; + basis.nplane = 2; + basis.startz_current = 2; + basis.nrxx = 8; + const ModuleBase::Matrix3 lattice(0.0, 1.0, 0.0, + -2.0, 0.0, 0.0, + 0.0, 0.0, 3.0); + std::vector> positions; + ModulePW::grid_positions(basis, lattice, 10.0, positions); + ASSERT_EQ(positions.size(), 8u); + EXPECT_DOUBLE_EQ(positions[0].x, 0.0); + EXPECT_DOUBLE_EQ(positions[0].y, 0.0); + EXPECT_DOUBLE_EQ(positions[0].z, 15.0); + EXPECT_DOUBLE_EQ(positions[7].x, -10.0); + EXPECT_DOUBLE_EQ(positions[7].y, 5.0); + EXPECT_DOUBLE_EQ(positions[7].z, 22.5); +} + +TEST(PwGridGeometry, PermitsEmptySlabs) +{ + ModulePW::PW_Basis basis("cpu", "double"); + basis.nx = 2; + basis.ny = 2; + basis.nz = 4; + basis.nplane = 0; + basis.startz_current = 4; + basis.nrxx = 0; + const ModuleBase::Matrix3 lattice; + std::vector> positions; + ModulePW::grid_positions(basis, lattice, 10.0, positions); + EXPECT_TRUE(positions.empty()); +} diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index b3e574e7144..17a6d6481a5 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -33,6 +33,7 @@ add_library( parallel_kpoints.cpp cell_index.cpp cell_tools.cpp + cell_geometry.cpp record_adj.cpp check_atomic_stru.cpp update_cell.cpp diff --git a/source/source_cell/cell_geometry.cpp b/source/source_cell/cell_geometry.cpp new file mode 100644 index 00000000000..cbe439ccadf --- /dev/null +++ b/source/source_cell/cell_geometry.cpp @@ -0,0 +1,137 @@ +#include "cell_geometry.h" + +#include "source_base/matrix3.h" + +#include + +namespace unitcell +{ + +bool make_orthogonal_cell(const ModuleBase::Matrix3& lattice, + const double lattice_scale, + const double relative_tolerance, + OrthogonalCell& cell) +{ + const ModuleBase::Vector3 a(lattice.e11, lattice.e12, lattice.e13); + const ModuleBase::Vector3 b(lattice.e21, lattice.e22, lattice.e23); + const ModuleBase::Vector3 c(lattice.e31, lattice.e32, lattice.e33); + const std::array, 3> vectors = {{a, b, c}}; + for (int axis = 0; axis < 3; ++axis) + { + const ModuleBase::Vector3 physical_vector = vectors[axis] * lattice_scale; + const double length = physical_vector.norm(); + cell.lengths[axis] = length; + cell.axes[axis] = physical_vector / length; + } + const ModuleBase::Vector3 vector_sum = a + b + c; + const double half_scale = 0.5 * lattice_scale; + cell.origin = vector_sum * half_scale; + for (int left = 0; left < 3; ++left) + { + for (int right = left + 1; right < 3; ++right) + { + const double projection = cell.axes[left] * cell.axes[right]; + if (std::abs(projection) > relative_tolerance) + { + return false; + } + } + } + return true; +} + +ModuleBase::Vector3 relative_position(const ModuleBase::Vector3& position, + const OrthogonalCell& cell) +{ + const ModuleBase::Vector3 displacement = position - cell.origin; + ModuleBase::Vector3 relative; + for (int axis = 0; axis < 3; ++axis) + { + const double projection = displacement * cell.axes[axis]; + const double length = cell.lengths[axis]; + const double fractional_image = projection / length + 0.5; + const double image = std::floor(fractional_image); + const double wrapped = projection - length * image; + relative += cell.axes[axis] * wrapped; + } + return relative; +} + +ModuleBase::Vector3 weighted_center(const std::vector>& positions, + const std::vector& weights, + const OrthogonalCell& cell) +{ + OrthogonalCell reference = cell; + reference.origin = positions.front(); + double total_weight = 0.0; + ModuleBase::Vector3 weighted_displacement; + for (std::size_t index = 0; index < positions.size(); ++index) + { + const ModuleBase::Vector3 displacement = relative_position(positions[index], reference); + weighted_displacement += displacement * weights[index]; + total_weight += weights[index]; + } + const ModuleBase::Vector3 mean_displacement = weighted_displacement / total_weight; + const ModuleBase::Vector3 unwrapped_center = positions.front() + mean_displacement; + const ModuleBase::Vector3 wrapped_center = relative_position(unwrapped_center, cell); + return cell.origin + wrapped_center; +} + +bool make_slab_cell(const ModuleBase::Matrix3& lattice, + const double lattice_scale, + const int open_axis, + const double relative_tolerance, + SlabCell& cell) +{ + const ModuleBase::Vector3 a(lattice.e11, lattice.e12, lattice.e13); + const ModuleBase::Vector3 b(lattice.e21, lattice.e22, lattice.e23); + const ModuleBase::Vector3 c(lattice.e31, lattice.e32, lattice.e33); + const std::array, 3> vectors = {{a, b, c}}; + const int first_axis = (open_axis + 1) % 3; + const int second_axis = (open_axis + 2) % 3; + const ModuleBase::Vector3 open = vectors[open_axis] * lattice_scale; + const ModuleBase::Vector3 first = vectors[first_axis] * lattice_scale; + const ModuleBase::Vector3 second = vectors[second_axis] * lattice_scale; + const double length = open.norm(); + const ModuleBase::Vector3 plane_normal = first ^ second; + const double area = plane_normal.norm(); + const double orientation = (plane_normal * open < 0.0) ? -1.0 : 1.0; + const double normal_scale = orientation / area; + const ModuleBase::Vector3 normal = plane_normal * normal_scale; + cell.normal = normal; + cell.length = length; + cell.area = area; + cell.origin = 0.5 * length; + const ModuleBase::Vector3 alignment = open ^ normal; + return alignment.norm() <= relative_tolerance * length; +} + +double relative_coordinate(const ModuleBase::Vector3& position, const SlabCell& cell) +{ + const double displacement = position * cell.normal - cell.origin; + const double image_coordinate = displacement / cell.length + 0.5; + const double image = std::floor(image_coordinate); + return displacement - cell.length * image; +} + +double weighted_center(const std::vector>& positions, + const std::vector& weights, + const SlabCell& cell) +{ + SlabCell reference = cell; + reference.origin = positions.front() * cell.normal; + double total_weight = 0.0; + double weighted_displacement = 0.0; + for (std::size_t index = 0; index < positions.size(); ++index) + { + const double displacement = relative_coordinate(positions[index], reference); + weighted_displacement += weights[index] * displacement; + total_weight += weights[index]; + } + const double unwrapped = reference.origin + weighted_displacement / total_weight; + const double image_coordinate = unwrapped / cell.length; + const double image = std::floor(image_coordinate); + return unwrapped - image * cell.length; +} + +} // namespace unitcell diff --git a/source/source_cell/cell_geometry.h b/source/source_cell/cell_geometry.h new file mode 100644 index 00000000000..639bfad09cf --- /dev/null +++ b/source/source_cell/cell_geometry.h @@ -0,0 +1,76 @@ +#ifndef CELL_GEOMETRY_H +#define CELL_GEOMETRY_H + +#include "source_base/vector3.h" + +#include +#include + +namespace ModuleBase +{ +class Matrix3; +} + +namespace unitcell +{ + +/// Orthogonal cell in Cartesian coordinates; lengths and origin are in Bohr. +/// Construct with make_orthogonal_cell before using relative_position. +struct OrthogonalCell +{ + std::array, 3> axes; + std::array lengths; + ModuleBase::Vector3 origin; +}; + +/// Accept rectangular cells, including rotated cells. Return false if the +/// lattice vectors are not mutually perpendicular within relative_tolerance. +bool make_orthogonal_cell(const ModuleBase::Matrix3& lattice, + double lattice_scale, + double relative_tolerance, + OrthogonalCell& cell); + +/// Minimum-image displacement from origin, with projections in [-L/2, L/2). +/// Requires a cell constructed above and a finite Cartesian position. +ModuleBase::Vector3 relative_position(const ModuleBase::Vector3& position, + const OrthogonalCell& cell); + +/// Unwrap positions about the first position, take their positive-weight center, +/// then wrap the center about cell.origin. Intended for localized systems whose +/// extent about the first position is less than half a cell in each direction. +/// Requires a cell constructed by make_orthogonal_cell, at least one position +/// and positive weights. +ModuleBase::Vector3 weighted_center(const std::vector>& positions, + const std::vector& weights, + const OrthogonalCell& cell); + +/// A periodic plane and its perpendicular open lattice vector, in Bohr. +/// The two in-plane lattice vectors need not be orthogonal. +struct SlabCell +{ + ModuleBase::Vector3 normal; + double length = 0.0; + double area = 0.0; + double origin = 0.0; +}; + +/// Return false if the open lattice vector is not perpendicular to the +/// periodic plane within relative_tolerance. +bool make_slab_cell(const ModuleBase::Matrix3& lattice, + double lattice_scale, + int open_axis, + double relative_tolerance, + SlabCell& cell); + +/// Normal displacement in [-L/2, L/2); requires a validated cell and position. +double relative_coordinate(const ModuleBase::Vector3& position, const SlabCell& cell); + +/// Mass/weight center along the normal, unwrapped about the first position. +/// Requires a localized slab of extent less than half the normal period. +double weighted_center(const std::vector>& positions, + const std::vector& weights, + const SlabCell& cell); + +} // namespace unitcell + +#endif diff --git a/source/source_cell/cell_tools.cpp b/source/source_cell/cell_tools.cpp index 6c380927cb4..624f04145fd 100644 --- a/source/source_cell/cell_tools.cpp +++ b/source/source_cell/cell_tools.cpp @@ -6,6 +6,29 @@ namespace unitcell { + std::vector get_atom_data(const Atom* atoms, const int ntype, const double lattice_scale) + { + int nat = 0; + for (int it = 0; it < ntype; ++it) + { + nat += atoms[it].na; + } + std::vector data; + data.reserve(nat); + for (int it = 0; it < ntype; ++it) + { + for (int ia = 0; ia < atoms[it].na; ++ia) + { + AtomData atom; + atom.position = atoms[it].tau[ia] * lattice_scale; + atom.mass = atoms[it].mass; + atom.valence_charge = atoms[it].ncpp.zv; + data.push_back(atom); + } + } + return data; + } + std::vector get_atomLabels(const Atom* atoms, const int ntype) { std::vector atomLabels(ntype); diff --git a/source/source_cell/cell_tools.h b/source/source_cell/cell_tools.h index 226416166af..037d152a282 100644 --- a/source/source_cell/cell_tools.h +++ b/source/source_cell/cell_tools.h @@ -16,6 +16,18 @@ */ namespace unitcell { + /// Cartesian atom information in type-major order, matching ABACUS atom indices. + struct AtomData + { + ModuleBase::Vector3 position; ///< Bohr + double mass = 0.0; + double valence_charge = 0.0; ///< Positive ionic charge in units of e + }; + + /// Extract positions, masses and pseudopotential valence charges. + /// Requires initialized Atom storage and a positive lattice scale in Bohr. + std::vector get_atom_data(const Atom* atoms, int ntype, double lattice_scale); + /// @brief Get atom labels for each atom type. /// @param atoms atom pointer [in] /// @param ntype number of atom types [in] diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index e0db610c5bf..b58299a9636 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -51,6 +51,18 @@ list(APPEND cell_simple_srcs add_library(cell_info OBJECT ${cell_simple_srcs}) +AddTest( + TARGET MODULE_CELL_cell_geometry + LIBS base device + SOURCES test_cell_geometry.cpp ../cell_geometry.cpp +) + +AddTest( + TARGET MODULE_CELL_cell_tools + LIBS base device + SOURCES test_cell_tools.cpp ../cell_tools.cpp ../atom_spec.cpp ../atom_pseudo.cpp ../pseudo.cpp +) + if(NOT ENABLE_MPI) return() endif() diff --git a/source/source_cell/test/test_cell_geometry.cpp b/source/source_cell/test/test_cell_geometry.cpp new file mode 100644 index 00000000000..652f8bd57f1 --- /dev/null +++ b/source/source_cell/test/test_cell_geometry.cpp @@ -0,0 +1,140 @@ +#include "source_cell/cell_geometry.h" +#include "source_base/matrix3.h" + +#include + +#include + +namespace +{ +class CellGeometryTest : public testing::Test +{ + protected: + void SetUp() override + { + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + ASSERT_TRUE(unitcell::make_orthogonal_cell(lattice, 10.0, 1.0e-10, cell)); + } + + unitcell::OrthogonalCell cell; +}; + +void expect_vector(const ModuleBase::Vector3& actual, + const ModuleBase::Vector3& expected) +{ + EXPECT_NEAR(actual.x, expected.x, 1.0e-13); + EXPECT_NEAR(actual.y, expected.y, 1.0e-13); + EXPECT_NEAR(actual.z, expected.z, 1.0e-13); +} +} // namespace + +TEST_F(CellGeometryTest, BuildsRotatedRectangularCellWrapsImagesAndRejectsSkew) +{ + const ModuleBase::Matrix3 lattice(0.0, 1.0, 0.0, + -2.0, 0.0, 0.0, + 0.0, 0.0, 3.0); + ASSERT_TRUE(unitcell::make_orthogonal_cell(lattice, 10.0, 1.0e-10, cell)); + EXPECT_DOUBLE_EQ(cell.lengths[0], 10.0); + EXPECT_DOUBLE_EQ(cell.lengths[1], 20.0); + EXPECT_DOUBLE_EQ(cell.lengths[2], 30.0); + const ModuleBase::Vector3 expected_origin(-10.0, 5.0, 15.0); + expect_vector(cell.origin, expected_origin); + + const ModuleBase::Vector3 displacement(11.0, 6.0, 16.0); + const ModuleBase::Vector3 position = cell.origin + displacement; + const ModuleBase::Vector3 actual = unitcell::relative_position(position, cell); + const ModuleBase::Vector3 expected(-9.0, -4.0, -14.0); + expect_vector(actual, expected); + + // Several images away from the origin of a 10-Bohr cube, on the half-open interval. + const ModuleBase::Matrix3 cube(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + ASSERT_TRUE(unitcell::make_orthogonal_cell(cube, 10.0, 1.0e-10, cell)); + cell.origin = ModuleBase::Vector3(0.0, 0.0, 0.0); + const ModuleBase::Vector3 far_position(25.0, -25.0, 31.0); + const ModuleBase::Vector3 far_actual = unitcell::relative_position(far_position, cell); + const ModuleBase::Vector3 wrapped(-5.0, -5.0, 1.0); + expect_vector(far_actual, wrapped); + + const ModuleBase::Matrix3 skew(1.0, 0.0, 0.0, + 0.1, 1.0, 0.0, + 0.0, 0.0, 1.0); + EXPECT_FALSE(unitcell::make_orthogonal_cell(skew, 10.0, 1.0e-10, cell)); +} + +TEST_F(CellGeometryTest, WeightedCenterUnwrapsAcrossTheBoundaryAndFollowsTranslations) +{ + const std::vector> boundary_positions = { + ModuleBase::Vector3(9.8, 9.7, 9.6), + ModuleBase::Vector3(0.2, 0.3, 0.4)}; + const std::vector boundary_weights = {1.0, 3.0}; + const ModuleBase::Vector3 center = unitcell::weighted_center(boundary_positions, boundary_weights, cell); + const ModuleBase::Vector3 expected_center(0.1, 0.15, 0.2); + expect_vector(center, expected_center); + + const std::vector> original = { + ModuleBase::Vector3(9.8, 9.6, 0.1), + ModuleBase::Vector3(0.3, 0.4, 9.7)}; + const std::vector weights = {2.0, 1.0}; + unitcell::OrthogonalCell first = cell; + first.origin = unitcell::weighted_center(original, weights, cell); + + std::vector> translated = original; + for (auto& position : translated) + { + const double translated_x = position.x + 1.4; + const double translated_y = position.y + 1.4; + const double translated_z = position.z + 1.4; + position.x = std::fmod(translated_x, 10.0); + position.y = std::fmod(translated_y, 10.0); + position.z = std::fmod(translated_z, 10.0); + } + unitcell::OrthogonalCell second = cell; + second.origin = unitcell::weighted_center(translated, weights, cell); + for (std::size_t index = 0; index < original.size(); ++index) + { + const ModuleBase::Vector3 expected = unitcell::relative_position(original[index], first); + const ModuleBase::Vector3 actual = unitcell::relative_position(translated[index], second); + expect_vector(actual, expected); + } +} + +TEST(SlabGeometry, SupportsSkewPlanesRotatedNormalsAndEveryAxisAndRejectsTilt) +{ + const ModuleBase::Matrix3 lattice(0.0, 2.0, 0.0, + 0.0, 1.0, 3.0, + 8.0, 0.0, 0.0); + unitcell::SlabCell cell; + ASSERT_TRUE(unitcell::make_slab_cell(lattice, 1.0, 2, 1.0e-6, cell)); + EXPECT_DOUBLE_EQ(cell.length, 8.0); + EXPECT_DOUBLE_EQ(cell.area, 6.0); + EXPECT_DOUBLE_EQ(cell.normal.x, 1.0); + const ModuleBase::Vector3 point(8.1, 100.0, -200.0); + EXPECT_NEAR(unitcell::relative_coordinate(point, cell), -3.9, 1.0e-14); + const std::vector> positions = { + ModuleBase::Vector3(7.8, 0.0, 0.0), + ModuleBase::Vector3(0.2, 40.0, 50.0)}; + const std::vector weights = {1.0, 3.0}; + const double center = unitcell::weighted_center(positions, weights, cell); + EXPECT_NEAR(center, 0.1, 1.0e-14); + + const ModuleBase::Matrix3 rectangular(2.0, 0.0, 0.0, + 0.0, 3.0, 0.0, + 0.0, 0.0, 8.0); + const double lengths[3] = {2.0, 3.0, 8.0}; + for (int axis = 0; axis < 3; ++axis) + { + ASSERT_TRUE(unitcell::make_slab_cell(rectangular, 1.0, axis, 1.0e-6, cell)); + EXPECT_DOUBLE_EQ(cell.length, lengths[axis]); + const double expected_area = 48.0 / lengths[axis]; + EXPECT_DOUBLE_EQ(cell.area, expected_area); + } + const ModuleBase::Matrix3 tilted(2.0, 0.0, 0.0, + 0.0, 3.0, 0.0, + 0.1, 0.0, 8.0); + EXPECT_FALSE(unitcell::make_slab_cell(tilted, 1.0, 2, 1.0e-6, cell)); +} + diff --git a/source/source_cell/test/test_cell_tools.cpp b/source/source_cell/test/test_cell_tools.cpp new file mode 100644 index 00000000000..d3b9af4b966 --- /dev/null +++ b/source/source_cell/test/test_cell_tools.cpp @@ -0,0 +1,42 @@ +#include "source_cell/cell_tools.h" + +#include + +TEST(CellTools, ExtractsPositionsMassesAndValenceChargesInAtomOrder) +{ + Atom atoms[2]; + atoms[0].na = 2; + atoms[0].mass = 12.0; + atoms[0].ncpp.zv = 4.0; + atoms[0].tau = {ModuleBase::Vector3(0.1, 0.2, 0.3), + ModuleBase::Vector3(0.4, 0.5, 0.6)}; + atoms[1].na = 1; + atoms[1].mass = 1.0; + atoms[1].ncpp.zv = 1.0; + atoms[1].tau = {ModuleBase::Vector3(0.7, 0.8, 0.9)}; + + const std::vector data = unitcell::get_atom_data(atoms, 2, 10.0); + ASSERT_EQ(data.size(), 3u); + EXPECT_DOUBLE_EQ(data[0].position.x, 1.0); + EXPECT_DOUBLE_EQ(data[0].position.y, 2.0); + EXPECT_DOUBLE_EQ(data[0].position.z, 3.0); + EXPECT_DOUBLE_EQ(data[1].position.x, 4.0); + EXPECT_DOUBLE_EQ(data[1].position.y, 5.0); + EXPECT_DOUBLE_EQ(data[1].position.z, 6.0); + EXPECT_DOUBLE_EQ(data[2].position.x, 7.0); + EXPECT_DOUBLE_EQ(data[2].position.y, 8.0); + EXPECT_DOUBLE_EQ(data[2].position.z, 9.0); + EXPECT_DOUBLE_EQ(data[0].mass, 12.0); + EXPECT_DOUBLE_EQ(data[1].mass, 12.0); + EXPECT_DOUBLE_EQ(data[2].mass, 1.0); + EXPECT_DOUBLE_EQ(data[0].valence_charge, 4.0); + EXPECT_DOUBLE_EQ(data[1].valence_charge, 4.0); + EXPECT_DOUBLE_EQ(data[2].valence_charge, 1.0); + EXPECT_DOUBLE_EQ(atoms[0].tau[0].x, 0.1); +} + +TEST(CellTools, AcceptsNoAtomTypes) +{ + const std::vector data = unitcell::get_atom_data(nullptr, 0, 10.0); + EXPECT_TRUE(data.empty()); +} diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index b5dd498251c..6be568facfa 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -11,12 +11,17 @@ list(APPEND objects elecstate_pw_sdft.cpp elecstate_pw_cal_tau.cpp makov_payne.cpp + charge_moments.cpp + pcc_0d.cpp + pcc_2d.cpp module_pot/gatefield.cpp module_pot/efield.cpp module_pot/h_hartree_pw.cpp module_pot/h_tddft_pw.cpp module_pot/pot_xc.cpp module_pot/pot_local.cpp + module_pot/pot_pcc.cpp + module_pot/pot_sccs.cpp module_pot/potential_new.cpp module_pot/potential_types.cpp module_pot/pot_surchem.cpp diff --git a/source/source_estate/charge_moments.cpp b/source/source_estate/charge_moments.cpp new file mode 100644 index 00000000000..ed905fa7d5c --- /dev/null +++ b/source/source_estate/charge_moments.cpp @@ -0,0 +1,32 @@ +#include "charge_moments.h" + +namespace elecstate +{ + +ChargeMoments charge_moments(const double* charges, + const ModuleBase::Vector3* relative_positions, + const int count, + const double weight) +{ + ChargeMoments moments; + for (int index = 0; index < count; ++index) + { + const double charge = charges[index] * weight; + const ModuleBase::Vector3& relative = relative_positions[index]; + moments.charge += charge; + moments.dipole += relative * charge; + moments.second_moment += charge * relative.norm2(); + } + return moments; +} + +ChargeMoments add_charge_moments(const ChargeMoments& left, const ChargeMoments& right) +{ + ChargeMoments result; + result.charge = left.charge + right.charge; + result.dipole = left.dipole + right.dipole; + result.second_moment = left.second_moment + right.second_moment; + return result; +} + +} // namespace elecstate diff --git a/source/source_estate/charge_moments.h b/source/source_estate/charge_moments.h new file mode 100644 index 00000000000..7dcf20cf319 --- /dev/null +++ b/source/source_estate/charge_moments.h @@ -0,0 +1,29 @@ +#ifndef CHARGE_MOMENTS_H +#define CHARGE_MOMENTS_H + +#include "source_base/vector3.h" + +namespace elecstate +{ + +/// Moments about one common origin. Charges carry their physical sign. +struct ChargeMoments +{ + double charge = 0.0; + ModuleBase::Vector3 dipole; + double second_moment = 0.0; ///< Trace of the second moment, in e Bohr^2 +}; + +/// Local moments only; no MPI or cell/model dependencies. Use weight=1 for +/// point charges and the volume element for a signed grid density. +/// Zero samples are permitted. +ChargeMoments charge_moments(const double* charges, + const ModuleBase::Vector3* relative_positions, + int count, + double weight); + +ChargeMoments add_charge_moments(const ChargeMoments& left, const ChargeMoments& right); + +} // namespace elecstate + +#endif diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 68d637a4d4f..ba96ff14576 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -157,7 +157,7 @@ class ElecState /// @param dft_plus_u 0 disables the DFT+U energy term /// @param assume_isolated "makov-payne" adds the isolated-system correction void cal_energies(const int type, - const bool imp_sol, + const int imp_sol, const bool sc_mag_switch, const int dft_plus_u, const std::string& assume_isolated); diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp index b417b910091..af77832329c 100644 --- a/source/source_estate/elecstate_energy.cpp +++ b/source/source_estate/elecstate_energy.cpp @@ -316,7 +316,7 @@ void ElecState::cal_converged() * @param type: 2 means Kohn-Sham functional; */ void ElecState::cal_energies(const int type, - const bool imp_sol, + const int imp_sol, const bool sc_mag_switch, const int dft_plus_u, const std::string& assume_isolated) @@ -331,11 +331,20 @@ void ElecState::cal_energies(const int type, this->f_en.gatefield = get_etot_gatefield(); //! energy from implicit solvation model - if (imp_sol) + if (imp_sol == 1) { this->f_en.esol_el = get_solvent_model_Ael(); this->f_en.esol_cav = get_solvent_model_Acav(); } + else if (imp_sol == 2) + { + this->pot->get_solvation_energy(this->f_en.esol_el, this->f_en.esol_cav); + } + else + { + this->f_en.esol_el = 0.0; + this->f_en.esol_cav = 0.0; + } //! spin constrained energy if (sc_mag_switch) @@ -374,6 +383,10 @@ void ElecState::cal_energies(const int type, } this->f_en.correction_el = makov_payne_correction(*ucell, *this->charge, v_elecstat_ptr).total; } + else if (assume_isolated == "pcc_0d" || assume_isolated == "pcc_2d") + { + this->f_en.correction_el = this->pot->pcc_energy_rydberg(); + } else { this->f_en.correction_el = 0.0; diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index 05566d819ab..4b978c65f70 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -287,6 +287,12 @@ void print_etot(const Magnetism& magnet, energies_Ry.push_back(elec.f_en.esol_cav); } + if (inp.assume_isolated == "pcc_0d" || inp.assume_isolated == "pcc_2d") + { + titles.push_back("E_pcc"); + energies_Ry.push_back(elec.f_en.correction_el); + } + //! electric field energy if (inp.efield_flag) { diff --git a/source/source_estate/module_pot/pot_base.h b/source/source_estate/module_pot/pot_base.h index aaa578bc5bf..658962f276e 100644 --- a/source/source_estate/module_pot/pot_base.h +++ b/source/source_estate/module_pot/pot_base.h @@ -28,7 +28,18 @@ class PotBase virtual void cal_fixed_v(double* vl_pseudo){} - virtual double get_energy() const { return 0.0; } + virtual double get_energy() const { return 0.0; } + + // Solvation contributions in Ry, added once after the band-energy potential subtraction. + virtual void get_solvation_energy(double& electrostatic, double& non_electrostatic) const + { + electrostatic = 0.0; + non_electrostatic = 0.0; + } + // Explicit ionic solvation force in Ry/Bohr, added to existing contributions. + virtual void add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} + // Electronic electrostatic correction only; excludes cavity/surface derivatives. + virtual const std::vector* solvent_electrostatic_potential() const { return nullptr; } bool fixed_mode = 0; bool dynamic_mode = 0; diff --git a/source/source_estate/module_pot/pot_pcc.cpp b/source/source_estate/module_pot/pot_pcc.cpp new file mode 100644 index 00000000000..271547e4fe3 --- /dev/null +++ b/source/source_estate/module_pot/pot_pcc.cpp @@ -0,0 +1,343 @@ +#include "pot_pcc.h" + +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_base/tool_quit.h" +#include "source_basis/module_pw/pw_grid_geometry.h" +#include "source_cell/cell_tools.h" + +#include +#include +#include + +namespace +{ +bool fractional_translation(const ModuleBase::Matrix3& rotation, + const ModuleBase::Vector3& translation) +{ + const double tolerance = 1.0e-6; + const double diagonal_error_x = rotation.e11 - 1.0; + const double diagonal_error_y = rotation.e22 - 1.0; + const double diagonal_error_z = rotation.e33 - 1.0; + const double identity_error = std::abs(diagonal_error_x) + std::abs(diagonal_error_y) + + std::abs(diagonal_error_z) + std::abs(rotation.e12) + + std::abs(rotation.e13) + std::abs(rotation.e21) + + std::abs(rotation.e23) + std::abs(rotation.e31) + + std::abs(rotation.e32); + if (identity_error > tolerance) + { + return false; + } + const double dx = translation.x - std::round(translation.x); + const double dy = translation.y - std::round(translation.y); + const double dz = translation.z - std::round(translation.z); + return std::abs(dx) > tolerance || std::abs(dy) > tolerance || std::abs(dz) > tolerance; +} + +bool primitive_symmetry(const ModuleSymmetry::Symmetry& symmetry) +{ + for (int operation = 0; operation < symmetry.nrotk; ++operation) + { + const bool forbidden = fractional_translation(symmetry.gmatrix[operation], symmetry.gtrans[operation]); + if (forbidden) + { + return false; + } + } + for (int operation = 0; operation < symmetry.nrotk_anti; ++operation) + { + const bool forbidden = fractional_translation(symmetry.gmatrix_anti[operation], + symmetry.gtrans_anti[operation]); + if (forbidden) + { + return false; + } + } + const ModuleBase::Matrix3 identity; + for (const auto& translation : symmetry.ptrans) + { + const bool forbidden = fractional_translation(identity, translation); + if (forbidden) + { + return false; + } + } + return true; +} +bool slab_operation_valid(const ModuleBase::Matrix3& rotation, + const ModuleBase::Vector3& translation, + const int axis) +{ + const double entries[3][3] = {{rotation.e11, rotation.e12, rotation.e13}, + {rotation.e21, rotation.e22, rotation.e23}, + {rotation.e31, rotation.e32, rotation.e33}}; + const double translations[3] = {translation.x, translation.y, translation.z}; + for (int other = 0; other < 3; ++other) + { + if (other != axis && (std::abs(entries[axis][other]) > 1.0e-6 + || std::abs(entries[other][axis]) > 1.0e-6)) + { + return false; + } + } + const double fraction = translations[axis] - std::round(translations[axis]); + return entries[axis][axis] < 0.0 || std::abs(fraction) <= 1.0e-6; +} + +bool slab_symmetry_valid(const ModuleSymmetry::Symmetry& symmetry, const int axis) +{ + for (int operation = 0; operation < symmetry.nrotk; ++operation) + { + const bool valid = slab_operation_valid(symmetry.gmatrix[operation], symmetry.gtrans[operation], axis); + if (!valid) { return false; } + } + for (int operation = 0; operation < symmetry.nrotk_anti; ++operation) + { + const bool valid = slab_operation_valid(symmetry.gmatrix_anti[operation], symmetry.gtrans_anti[operation], axis); + if (!valid) { return false; } + } + const ModuleBase::Matrix3 identity; + for (const auto& translation : symmetry.ptrans) + { + const bool valid = slab_operation_valid(identity, translation, axis); + if (!valid) { return false; } + } + return true; +} +} // namespace + +namespace elecstate +{ + +PotPcc::PotPcc(const ModulePW::PW_Basis* basis, const double electron_count) + : dimension_(Dimension::molecule), open_axis_(2), electron_count_(electron_count) +{ + this->rho_basis_ = basis; + this->dynamic_mode = true; +} + +PotPcc::PotPcc(const ModulePW::PW_Basis* basis, + const Dimension dimension, + const int open_axis, + const double electron_count) + : dimension_(dimension), open_axis_(open_axis), electron_count_(electron_count) +{ + this->rho_basis_ = basis; + this->dynamic_mode = true; +} + +void PotPcc::validate_kpoints(const std::vector>& points, + const int count, + const int open_axis) +{ + double maximum = 0.0; + for (int ik = 0; ik < count; ++ik) + { + const double components[3] = {points[ik].x, points[ik].y, points[ik].z}; + const double absolute = std::abs(components[open_axis]); + maximum = std::max(maximum, absolute); + } +#ifdef __MPI + // Each pool holds its own k points; all pools must stop together. + Parallel_Reduce::reduce_max(maximum); +#endif + if (maximum > 1.0e-12) + { + ModuleBase::WARNING_QUIT("PotPcc::validate_kpoints", "pcc_2d requires Gamma-only sampling along pcc_2d_axis"); + } +} + +ModuleBase::Vector3 PotPcc::relative_position(const ModuleBase::Vector3& position) const +{ + if (dimension_ == Dimension::molecule) + { + return unitcell::relative_position(position, geometry_); + } + const double coordinate = unitcell::relative_coordinate(position, slab_); + return slab_.normal * coordinate; +} + +double PotPcc::correction_energy() const +{ + if (dimension_ == Dimension::molecule) { return pcc_0d_energy(moments_, parameters_); } + return pcc_2d_energy(moments_, slab_parameters_); +} + +double PotPcc::correction_potential(const ModuleBase::Vector3& position) const +{ + if (dimension_ == Dimension::molecule) { return pcc_0d_potential(moments_, position, parameters_); } + const double coordinate = position * slab_.normal; + return pcc_2d_potential(moments_, coordinate, slab_.normal, slab_parameters_); +} + +ModuleBase::Vector3 PotPcc::correction_force(const double charge, + const ModuleBase::Vector3& position) const +{ + if (dimension_ == Dimension::molecule) { return pcc_0d_force(moments_, charge, position, parameters_); } + const double coordinate = position * slab_.normal; + return pcc_2d_force(moments_, charge, coordinate, slab_.normal, slab_parameters_); +} + +void PotPcc::prepare_ions(const UnitCell& cell, ChargeMoments& ionic_moments) +{ + ModuleBase::timer::start("PotPcc", "prepare_ions"); + // The cell, symmetry and atoms are replicated, so every rank stops together. + if (dimension_ == Dimension::molecule) + { + const bool orthogonal = unitcell::make_orthogonal_cell(cell.latvec, cell.lat0, 1.0e-10, geometry_); + const bool cubic = orthogonal && make_pcc_0d_parameters(geometry_, 1.0e-10, parameters_); + if (!cubic) + { + ModuleBase::WARNING_QUIT("PotPcc::prepare_ions", "PCC 0D requires an equal-edge cubic cell"); + } + if (!primitive_symmetry(cell.symm)) + { + ModuleBase::WARNING_QUIT("PotPcc::prepare_ions", + "PCC 0D requires a primitive cell without fractional-translation symmetry"); + } + } + else + { + const bool perpendicular = unitcell::make_slab_cell(cell.latvec, cell.lat0, open_axis_, 1.0e-6, slab_); + if (!perpendicular) + { + ModuleBase::WARNING_QUIT("PotPcc::prepare_ions", + "PCC 2D requires the open lattice vector to be perpendicular to the periodic plane"); + } + slab_parameters_ = make_pcc_2d_parameters(slab_); + if (!slab_symmetry_valid(cell.symm, open_axis_)) + { + ModuleBase::WARNING_QUIT("PotPcc::prepare_ions", + "PCC 2D symmetry must preserve the open axis without fractional translations along it"); + } + } + + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + const int atom_count = static_cast(atoms.size()); + std::vector masses(atom_count); + ionic_positions_.resize(atom_count); + ionic_charges_.resize(atom_count); + for (int atom = 0; atom < atom_count; ++atom) + { + ionic_positions_[atom] = atoms[atom].position; + ionic_charges_[atom] = atoms[atom].valence_charge; + masses[atom] = atoms[atom].mass; + } + if (dimension_ == Dimension::molecule) + { + geometry_.origin = unitcell::weighted_center(ionic_positions_, masses, geometry_); + } + else + { + slab_.origin = unitcell::weighted_center(ionic_positions_, masses, slab_); + } + std::vector> relative_ions(atom_count); + for (int atom = 0; atom < atom_count; ++atom) + { + relative_ions[atom] = this->relative_position(ionic_positions_[atom]); + } + const double* ionic_charge_data = ionic_charges_.data(); + const ModuleBase::Vector3* ionic_position_data = relative_ions.data(); + ionic_moments = charge_moments(ionic_charge_data, ionic_position_data, atom_count, 1.0); + + ModuleBase::timer::end("PotPcc", "prepare_ions"); +} + +ChargeMoments PotPcc::collect_electrons( + const Charge& charge, + const UnitCell& cell, + std::vector>& positions) const +{ + ModuleBase::timer::start("PotPcc", "collect_electrons"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + std::vector electronic_charge(basis.nrxx, 0.0); + ModulePW::grid_positions(basis, cell.latvec, cell.lat0, positions); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + for (int spin = 0; spin < charge.nspin; ++spin) + { + electronic_charge[ir] -= charge.rho[spin][ir]; + } + positions[ir] = this->relative_position(positions[ir]); + } + const double volume_element = cell.omega / basis.nxyz; + const double* electronic_data = electronic_charge.data(); + const ModuleBase::Vector3* position_data = positions.data(); + ChargeMoments electronic_moments = charge_moments(electronic_data, position_data, basis.nrxx, volume_element); + double reduced[5] = {electronic_moments.charge, electronic_moments.dipole.x, + electronic_moments.dipole.y, electronic_moments.dipole.z, + electronic_moments.second_moment}; +#ifdef __MPI + Parallel_Reduce::reduce_pool(reduced, 5); +#endif + electronic_moments.charge = reduced[0]; + electronic_moments.dipole = ModuleBase::Vector3(reduced[1], reduced[2], reduced[3]); + electronic_moments.second_moment = reduced[4]; + ModuleBase::timer::end("PotPcc", "collect_electrons"); + return electronic_moments; +} + +void PotPcc::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) +{ + ModuleBase::timer::start("PotPcc", "cal_v_eff"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + ChargeMoments ionic_moments; + this->prepare_ions(*cell, ionic_moments); + std::vector> positions; + const ChargeMoments electronic_moments = this->collect_electrons(*charge, *cell, positions); + // Ions are replicated on every rank and must be added after the reduction. + moments_ = add_charge_moments(ionic_moments, electronic_moments); + // The grid density need not integrate to nelec exactly, e.g. after a very + // tight diagonalization; the correction uses the grid charge and reports + // the mismatch. The moments are pool-reduced, so every rank agrees. + const double expected_charge = ionic_moments.charge - electron_count_; + const double charge_error = moments_.charge - expected_charge; + if (std::abs(charge_error) > 1.0e-6) + { + std::ostringstream message; + message << "PCC net charge differs from the expected value by " << charge_error + << " e; the correction uses the grid charge"; + ModuleBase::WARNING("PotPcc::cal_v_eff", message.str()); + } + energy_rydberg_ = 2.0 * this->correction_energy(); + electron_potential_.resize(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double positive_potential = this->correction_potential(positions[ir]); + electron_potential_[ir] = -2.0 * positive_potential; + } + for (int spin = 0; spin < charge->nspin; ++spin) + { + for (int ir = 0; ir < basis.nrxx; ++ir) + { + potential(spin, ir) += electron_potential_[ir]; + } + } + ModuleBase::timer::end("PotPcc", "cal_v_eff"); +} + +double PotPcc::get_energy() const +{ + return energy_rydberg_; +} + +const std::vector& PotPcc::electron_potential() const +{ + return electron_potential_; +} + +void PotPcc::add_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + ModuleBase::timer::start("PotPcc", "add_force"); + for (int atom = 0; atom < cell.nat; ++atom) + { + const ModuleBase::Vector3 relative = this->relative_position(ionic_positions_[atom]); + const ModuleBase::Vector3 correction = this->correction_force(ionic_charges_[atom], relative); + force(atom, 0) += 2.0 * correction.x; + force(atom, 1) += 2.0 * correction.y; + force(atom, 2) += 2.0 * correction.z; + } + ModuleBase::timer::end("PotPcc", "add_force"); +} + +} // namespace elecstate diff --git a/source/source_estate/module_pot/pot_pcc.h b/source/source_estate/module_pot/pot_pcc.h new file mode 100644 index 00000000000..c567cfa5400 --- /dev/null +++ b/source/source_estate/module_pot/pot_pcc.h @@ -0,0 +1,92 @@ +#ifndef POT_PCC_H +#define POT_PCC_H + +#include "pot_base.h" +#include "source_cell/cell_geometry.h" +#include "source_estate/pcc_0d.h" +#include "source_estate/pcc_2d.h" + +namespace elecstate +{ + +/** + * @brief Self-consistent vacuum point-counter-charge (PCC) correction. + * + * Owns the result of its last density update. The correction removes the + * spurious interaction of the total charge with its periodic images by + * expanding that charge to second order in its multipoles about the + * mass-weighted ionic center (charge q, dipole d, second moment Q). The + * charge seen here is ionic + electronic only; any solvent polarization charge + * must be corrected by its own component. + * In Hartree atomic units, for a unit positive test charge at r (or at the + * normal coordinate z for slabs): + * + * 0D (cubic cell of edge L, volume V = L^3, Madelung constant alpha of the + * simple cubic lattice): + * v(r) = alpha q / L - 2 pi / (3 V) * (q r^2 - 2 d.r + Q) + * + * 2D (periodic area A, normal period L): + * v(z) = -pi q L / (3 A) - 2 pi / (A L) * (q z^2 - 2 d_z z + Q_zz) + * + * The correction energy is E = 1/2 * integral rho_total(r) v(r) dr and the + * forces follow from the gradient of v at the ionic positions. The kernels are + * implemented in pcc_0d.cpp and pcc_2d.cpp; see pcc_2d.h for the gauge of the + * charged-slab constant. + * + * References: + * - O. Andreussi and N. Marzari, Phys. Rev. B 90, 245101 (2014) + * (parabolic 0D and 2D point-counter-charge corrections). + * - I. Dabo, B. Kozinsky, N. E. Singh-Miller and N. Marzari, + * Phys. Rev. B 77, 115139 (2008) (density-countercharge corrections). + * - M. J. Rutter, Electron. Struct. 3, 015002 (2021) (charged slabs: the + * zero-mean periodic potential of a charged plane fixes the 2D constant). + * + * The implementation follows the PCC corrections of the ENVIRON library + * (www.quantum-environ.org). + */ +class PotPcc : public PotBase +{ + public: + enum class Dimension { molecule, slab }; + + /// electron_count is the expected number of electrons (INPUT nelec); a + /// grid charge that differs from it by more than 1e-6 e is reported in the + /// warning log, and the correction uses the grid charge. + PotPcc(const ModulePW::PW_Basis* basis, double electron_count); + + PotPcc(const ModulePW::PW_Basis* basis, Dimension dimension, int open_axis, double electron_count); + static void validate_kpoints(const std::vector>& points, int count, int open_axis); + + void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override; + double get_energy() const override; + const std::vector& electron_potential() const; + void add_force(const UnitCell& cell, ModuleBase::matrix& force) const; + + private: + void prepare_ions(const UnitCell& cell, ChargeMoments& ionic_moments); + ChargeMoments collect_electrons(const Charge& charge, + const UnitCell& cell, + std::vector>& positions) const; + + ModuleBase::Vector3 relative_position(const ModuleBase::Vector3& position) const; + double correction_energy() const; + double correction_potential(const ModuleBase::Vector3& position) const; + ModuleBase::Vector3 correction_force(double charge, const ModuleBase::Vector3& position) const; + + const Dimension dimension_; + const int open_axis_; + const double electron_count_; + unitcell::SlabCell slab_; + Pcc2dParameters slab_parameters_; + unitcell::OrthogonalCell geometry_; + Pcc0dParameters parameters_; + ChargeMoments moments_; + std::vector> ionic_positions_; + std::vector ionic_charges_; + std::vector electron_potential_; + double energy_rydberg_ = 0.0; +}; + +} // namespace elecstate + +#endif diff --git a/source/source_estate/module_pot/pot_sccs.cpp b/source/source_estate/module_pot/pot_sccs.cpp new file mode 100644 index 00000000000..7d5dd323b19 --- /dev/null +++ b/source/source_estate/module_pot/pot_sccs.cpp @@ -0,0 +1,157 @@ +#include "pot_sccs.h" + +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_base/tool_quit.h" +#include "source_cell/cell_tools.h" +#include "source_hamilt/module_sccs/sccs_functional.h" +#include "source_hamilt/module_sccs/sccs_ionic_charge.h" +#include "source_hamilt/module_sccs/sccs_ionic_force.h" +#include "source_hamilt/module_sccs/sccs_response.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include +#include + +namespace elecstate +{ +void make_sccs_config_from_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver) +{ + const ModuleSccs::Preset preset = ModuleSccs::parse_preset(input.sccs_preset); + if (preset == ModuleSccs::Preset::Custom) + { + config = ModuleSccs::SccsConfig(); + config.cavity.density_min = input.sccs_rho_min; + config.cavity.density_max = input.sccs_rho_max; + config.cavity.epsilon_bulk = input.sccs_epsilon; + config.surface_tension = ModuleSccs::dyn_per_cm_to_hartree_per_bohr2(input.sccs_gamma); + config.pressure = ModuleSccs::gpa_to_hartree_per_bohr3(input.sccs_pressure); + } + else + { + config = ModuleSccs::make_sccs_config(preset); + } + config.surface_regularization = input.sccs_surface_eta; + solver.max_iterations = input.sccs_maxiter; + solver.tolerance_rms = input.sccs_tol_rms; + solver.tolerance_max = input.sccs_tol_max; +} + +std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count) +{ + double ionic_charge = 0.0; + for (int it = 0; it < cell.ntype; ++it) + { + ionic_charge += cell.atoms[it].ncpp.zv * cell.atoms[it].na; + } + const double net_charge = ionic_charge - electron_count; + if (config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= 1e-6) { return std::string(); } + return "charged SCCS with periodic boundaries: the periodic Poisson solver drops the G = 0 component " + "of the net charge, so the energy depends on the cell size."; +} + +PotSccs::PotSccs(const ModulePW::PW_Basis* basis, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver, + double electron_count) + : config_(config), solver_(solver), electron_count_(electron_count) +{ + this->rho_basis_ = basis; + this->dynamic_mode = true; +} + +void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) +{ + ModuleBase::timer::start("PotSccs", "cal_v_eff"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + const std::vector atoms = unitcell::get_atom_data(cell->atoms, cell->ntype, cell->lat0); + std::vector ions; + ModuleSccs::gaussian_ionic_density(atoms, basis, cell->tpiba, ModuleSccs::gaussian_ion_spread, ions); + std::vector density(basis.nrxx, 0.0); + std::vector solute_charge(basis.nrxx); + double grid_charges[2] = {0.0, 0.0}; // electrons, ions + for (int ir = 0; ir < basis.nrxx; ++ir) + { + for (int spin = 0; spin < charge->nspin; ++spin) { density[ir] += charge->rho[spin][ir]; } + solute_charge[ir] = ions[ir] - density[ir]; + grid_charges[0] += density[ir]; + grid_charges[1] += ions[ir]; + } + // Pool-reduced, so every rank takes the same decision. + Parallel_Reduce::reduce_pool(grid_charges, 2); + const double dv = basis.omega / basis.nxyz; + const double electron_error = grid_charges[0] * dv - electron_count_; + double valence_charge = 0.0; + for (const unitcell::AtomData& atom : atoms) { valence_charge += atom.valence_charge; } + const double ionic_error = grid_charges[1] * dv - valence_charge; + if (std::abs(ionic_error) > 1e-6) + { + ModuleBase::WARNING_QUIT("PotSccs::cal_v_eff", "SCCS ionic density normalization does not match the valence charge"); + } + if (std::abs(electron_error) > 1e-6) + { + std::ostringstream message; + message << "SCCS grid electron count differs from the expected value by " << electron_error + << " e; the solute charge uses the grid density"; + ModuleBase::WARNING("PotSccs::cal_v_eff", message.str()); + } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis, + cell->tpiba, response); + ModuleSccs::FunctionalResult functional; + ModuleSccs::evaluate_functional(solute_charge, response, config_, basis, cell->tpiba, functional); + electrostatic_rydberg_ = 2.0 * functional.reaction_energy; + non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy); + electrostatic_potential_.resize(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + electrostatic_potential_[ir] = -2.0 * functional.reaction_potential[ir]; + const double value = 2.0 * functional.electron_potential[ir]; + for (int spin = 0; spin < charge->nspin; ++spin) { potential(spin, ir) += value; } + } + restart_potential_ = std::move(response.restart_potential); + ModuleBase::timer::end("PotSccs", "cal_v_eff"); +} + +void PotSccs::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + ModuleBase::timer::start("PotSccs", "add_solvation_force"); + const ModulePW::PW_Basis& basis = *this->rho_basis_; + const std::vector atoms = unitcell::get_atom_data(cell.atoms, cell.ntype, cell.lat0); + std::vector reaction(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + reaction[ir] = -0.5 * electrostatic_potential_[ir]; + } + std::vector> ionic_force; + ModuleSccs::gaussian_ionic_force(atoms, reaction, basis, cell.tpiba, ModuleSccs::gaussian_ion_spread, + ionic_force); + for (int ia = 0; ia < cell.nat; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + force(ia, axis) += 2.0 * ionic_force[ia][axis]; + } + } + ModuleBase::timer::end("PotSccs", "add_solvation_force"); +} + +double PotSccs::get_energy() const +{ + return electrostatic_rydberg_ + non_electrostatic_rydberg_; +} + +void PotSccs::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = electrostatic_rydberg_; + non_electrostatic = non_electrostatic_rydberg_; +} + +const std::vector* PotSccs::solvent_electrostatic_potential() const +{ + return &electrostatic_potential_; +} +} // namespace elecstate diff --git a/source/source_estate/module_pot/pot_sccs.h b/source/source_estate/module_pot/pot_sccs.h new file mode 100644 index 00000000000..2a2d02e3dcf --- /dev/null +++ b/source/source_estate/module_pot/pot_sccs.h @@ -0,0 +1,48 @@ +#ifndef POT_SCCS_H +#define POT_SCCS_H + +#include "pot_base.h" +#include "source_hamilt/module_sccs/sccs_parameters.h" + +#include + +struct Input_para; +namespace elecstate +{ +// Map INPUT values, already validated by the INPUT reader, onto the solver. +void make_sccs_config_from_input(const Input_para& input, + ModuleSccs::SccsConfig& config, + ModuleSccs::PolarizationSolverParameters& solver); + +// Returns a warning (empty if none) for a charged cell in a dielectric: the +// periodic Poisson solver drops the G = 0 component of the net charge, so the +// energy depends on the cell size. electron_count is INPUT nelec. +std::string check_sccs_charge(const ModuleSccs::SccsConfig& config, const UnitCell& cell, double electron_count); + +class PotSccs : public PotBase +{ +public: + // electron_count is the expected number of electrons (INPUT nelec); a grid + // density that differs from it by more than 1e-6 e is reported in the + // warning log, and the solute charge uses the grid density. + PotSccs(const ModulePW::PW_Basis* basis, + const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver, + double electron_count); + void cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::matrix& potential) override; + double get_energy() const override; + void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; + void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; + const std::vector* solvent_electrostatic_potential() const override; + +private: + const ModuleSccs::SccsConfig config_; + const ModuleSccs::PolarizationSolverParameters solver_; + const double electron_count_; + double electrostatic_rydberg_ = 0.0; + double non_electrostatic_rydberg_ = 0.0; + std::vector electrostatic_potential_; + std::vector restart_potential_; +}; +} // namespace elecstate +#endif diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index 150dfe893a1..51dbced99f1 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -1,4 +1,5 @@ #include "potential_new.h" +#include "pot_pcc.h" #include "pot_ml_exx.h" #include "source_base/global_function.h" @@ -14,7 +15,6 @@ namespace elecstate { - Potential::Potential(const ModulePW::PW_Basis* rho_basis_in, const ModulePW::PW_Basis* rho_basis_smooth_in, const UnitCell* ucell_in, @@ -62,6 +62,11 @@ Potential::~Potential() } void Potential::pot_register(const std::vector& components_list) +{ + this->pot_register(components_list, PARAM.inp); +} + +void Potential::pot_register(const std::vector& components_list, const Input_para& input) { ModuleBase::TITLE("Potential", "pot_register"); // delete old components first. @@ -80,7 +85,7 @@ void Potential::pot_register(const std::vector& components_list) //--------------------------- for (auto comp: components_list) { - PotBase* tmp = this->get_pot_type(comp); + PotBase* tmp = this->get_pot_type(comp, input); this->components.push_back(tmp); } @@ -356,4 +361,78 @@ double Potential::get_ml_exx_energy() const #endif } +const PotPcc* Potential::pcc_component() const +{ + for (const PotBase* component : this->components) + { + const PotPcc* pcc = dynamic_cast(component); + if (pcc != nullptr) + { + return pcc; + } + } + return nullptr; +} + +double Potential::pcc_energy_rydberg() const +{ + const PotPcc* pcc = this->pcc_component(); + if (pcc == nullptr) + { + return 0.0; + } + return pcc->get_energy(); +} + +const std::vector* Potential::pcc_electron_potential() const +{ + const PotPcc* pcc = this->pcc_component(); + if (pcc == nullptr) + { + return nullptr; + } + return &pcc->electron_potential(); +} + +void Potential::add_pcc_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + const PotPcc* pcc = this->pcc_component(); + if (pcc != nullptr) + { + pcc->add_force(cell, force); + } +} + +void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ + for (const PotBase* component : this->components) + { + component->add_solvation_force(cell, force); + } +} + +void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = 0.0; + non_electrostatic = 0.0; + for (const PotBase* component : this->components) + { + double el = 0.0; + double nonel = 0.0; + component->get_solvation_energy(el, nonel); + electrostatic += el; + non_electrostatic += nonel; + } +} + +const std::vector* Potential::solvent_electrostatic_potential() const +{ + for (const PotBase* component : this->components) + { + const std::vector* value = component->solvent_electrostatic_potential(); + if (value != nullptr) { return value; } + } + return nullptr; +} + } // namespace elecstate diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index 77ab063c99f..9bf4683a139 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -10,8 +10,11 @@ #include #include +struct Input_para; + namespace elecstate { +class PotPcc; class TDFieldManager; /** @@ -72,12 +75,17 @@ class Potential : public PotBase void init_pot(const Charge* const chg); // initialize potential components before SCF void pot_register(const std::vector& components_list); + void pot_register(const std::vector& components_list, const Input_para& input); // update potential from current charge void update_from_charge(const Charge* const chg, const UnitCell* const ucell); // interface for SCF-converged, etxc vtxc for Energy, vnew for force_scc void get_vnew(const Charge* chg, ModuleBase::matrix& vnew); - PotBase* get_pot_type(const std::string& pot_type); + PotBase* get_pot_type(const std::string& pot_type, const Input_para& input); + + void get_solvation_energy(double& electrostatic, double& non_electrostatic) const override; + void add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const override; + const std::vector* solvent_electrostatic_potential() const override; /** * @brief Inject the shared RT-TDDFT field state before potential setup. @@ -211,6 +219,12 @@ class Potential : public PotBase /// @return E_ML-EXX double get_ml_exx_energy() const; + /// Results of the registered point-counter-charge (PCC) component, used by + /// the energy, force and output stages; neutral values when PCC is off. + double pcc_energy_rydberg() const; + const std::vector* pcc_electron_potential() const; + void add_pcc_force(const UnitCell& cell, ModuleBase::matrix& force) const; + private: void cal_v_eff(const Charge* const chg, const UnitCell* const ucell, ModuleBase::matrix& v_eff) override; void cal_fixed_v(double* vl_pseudo) override; @@ -219,6 +233,10 @@ class Potential : public PotBase void allocate(); + // the registered PCC component, or nullptr when PCC is off + const PotPcc* pcc_component() const; + PotBase* make_sccs_potential(const Input_para& input); + std::vector v_eff_fixed; ModuleBase::matrix v_eff; diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 8b2f2b5e033..6890d083bdc 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -2,6 +2,8 @@ #include "efield.h" #include "gatefield.h" #include "pot_local.h" +#include "pot_pcc.h" +#include "pot_sccs.h" #include "pot_sep.h" #include "pot_surchem.h" #include "pot_xc.h" @@ -20,7 +22,20 @@ namespace elecstate { -PotBase* Potential::get_pot_type(const std::string& pot_type) +PotBase* Potential::make_sccs_potential(const Input_para& input) +{ + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + make_sccs_config_from_input(input, config, solver); + if (this->ucell_ != nullptr) + { + const std::string warning = check_sccs_charge(config, *this->ucell_, input.nelec); + if (!warning.empty()) { ModuleBase::WARNING("Potential::make_sccs_potential", warning); } + } + return new PotSccs(this->rho_basis_, config, solver, input.nelec); +} + +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para& input) { ModuleBase::TITLE("Potential", "get_pot_type"); if (pot_type == "local") @@ -33,7 +48,21 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) } else if (pot_type == "xc") { - return new PotXC(PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); + return new PotXC(PARAM.globalv.domag, PARAM.globalv.domag_z, input.gga_grad, this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); + } + else if (pot_type == "pcc" || pot_type == "pcc_2d") + { + const double electron_count = input.nelec; + if (pot_type == "pcc") + { + return new PotPcc(this->rho_basis_, electron_count); + } + const int open_axis = input.pcc_2d_axis; + return new PotPcc(this->rho_basis_, PotPcc::Dimension::slab, open_axis, electron_count); + } + else if (pot_type == "sccs") + { + return this->make_sccs_potential(input); } else if (pot_type == "surchem") { @@ -41,7 +70,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) } else if (pot_type == "efield") { - return new PotEfield(this->rho_basis_, this->ucell_, this->solvent_, PARAM.inp.dip_cor_flag); + return new PotEfield(this->rho_basis_, this->ucell_, this->solvent_, input.dip_cor_flag); } else if (pot_type == "gatefield") { diff --git a/source/source_estate/pcc_0d.cpp b/source/source_estate/pcc_0d.cpp new file mode 100644 index 00000000000..132ac499a43 --- /dev/null +++ b/source/source_estate/pcc_0d.cpp @@ -0,0 +1,80 @@ +#include "pcc_0d.h" + +#include "source_base/constants.h" +#include "source_cell/cell_geometry.h" + +#include + +namespace elecstate +{ + +bool make_pcc_0d_parameters(const unitcell::OrthogonalCell& cell, + const double relative_tolerance, + Pcc0dParameters& parameters) +{ + const double length = cell.lengths[0]; + parameters.length = length; + for (int axis = 1; axis < 3; ++axis) + { + const double difference = cell.lengths[axis] - length; + if (std::abs(difference) > relative_tolerance * length) + { + return false; + } + } + return true; +} + +double pcc_0d_potential(const ChargeMoments& moments, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters) +{ + const double length = parameters.length; + const double volume = length * length * length; + const double dipole_projection = moments.dipole * position; + const double parabolic = moments.charge * position.norm2() + - 2.0 * dipole_projection + moments.second_moment; + return parameters.madelung * moments.charge / length + - 2.0 * ModuleBase::PI * parabolic / (3.0 * volume); +} + +ModuleBase::Vector3 pcc_0d_gradient(const ChargeMoments& moments, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters) +{ + const double length = parameters.length; + const double volume = length * length * length; + const double factor = -4.0 * ModuleBase::PI / (3.0 * volume); + const ModuleBase::Vector3 field = position * moments.charge - moments.dipole; + return field * factor; +} + +double pcc_0d_bilinear_energy(const ChargeMoments& left, + const ChargeMoments& right, + const Pcc0dParameters& parameters) +{ + const double length = parameters.length; + const double volume = length * length * length; + const double monopole = parameters.madelung * left.charge * right.charge / length; + const double dipole_product = left.dipole * right.dipole; + const double multipole = left.second_moment * right.charge + + left.charge * right.second_moment - 2.0 * dipole_product; + return monopole - 2.0 * ModuleBase::PI * multipole / (3.0 * volume); +} + +double pcc_0d_energy(const ChargeMoments& moments, const Pcc0dParameters& parameters) +{ + const double bilinear = pcc_0d_bilinear_energy(moments, moments, parameters); + return 0.5 * bilinear; +} + +ModuleBase::Vector3 pcc_0d_force(const ChargeMoments& moments, + const double ionic_charge, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters) +{ + const ModuleBase::Vector3 gradient = pcc_0d_gradient(moments, position, parameters); + return gradient * (-ionic_charge); +} + +} // namespace elecstate diff --git a/source/source_estate/pcc_0d.h b/source/source_estate/pcc_0d.h new file mode 100644 index 00000000000..8931c6ac73e --- /dev/null +++ b/source/source_estate/pcc_0d.h @@ -0,0 +1,51 @@ +#ifndef PCC_0D_H +#define PCC_0D_H + +#include "charge_moments.h" + +namespace unitcell +{ +struct OrthogonalCell; +} + +namespace elecstate +{ + +struct Pcc0dParameters +{ + double length = 0.0; ///< Cubic cell edge in Bohr + double madelung = 2.837297479480619; +}; + +/// Return false unless the OrthogonalCell has equal edges within +/// relative_tolerance. No cell geometry is stored here. +bool make_pcc_0d_parameters(const unitcell::OrthogonalCell& cell, + double relative_tolerance, + Pcc0dParameters& parameters); + +/// Numerical kernels require validated parameters and moments about a common +/// origin. Coordinates are relative to that origin. Energies/potentials are in +/// Hartree; gradients/forces are in Hartree/Bohr. The potential is for unit +/// positive charge; electron potential conversion belongs to the host adapter. +double pcc_0d_potential(const ChargeMoments& moments, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters); + +ModuleBase::Vector3 pcc_0d_gradient(const ChargeMoments& moments, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters); + +double pcc_0d_bilinear_energy(const ChargeMoments& left, + const ChargeMoments& right, + const Pcc0dParameters& parameters); + +double pcc_0d_energy(const ChargeMoments& moments, const Pcc0dParameters& parameters); + +ModuleBase::Vector3 pcc_0d_force(const ChargeMoments& moments, + double ionic_charge, + const ModuleBase::Vector3& position, + const Pcc0dParameters& parameters); + +} // namespace elecstate + +#endif diff --git a/source/source_estate/pcc_2d.cpp b/source/source_estate/pcc_2d.cpp new file mode 100644 index 00000000000..63fbd9289bd --- /dev/null +++ b/source/source_estate/pcc_2d.cpp @@ -0,0 +1,59 @@ +#include "pcc_2d.h" + +#include "source_base/constants.h" +#include "source_cell/cell_geometry.h" + +#include + +namespace elecstate +{ +Pcc2dParameters make_pcc_2d_parameters(const unitcell::SlabCell& cell) +{ + Pcc2dParameters parameters; + parameters.area = cell.area; + parameters.length = cell.length; + return parameters; +} + +double pcc_2d_potential(const ChargeMoments& moments, + const double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters) +{ + const double dipole = moments.dipole * normal; + const double parabolic = moments.charge * coordinate * coordinate + - 2.0 * dipole * coordinate + moments.second_moment; + const double constant = -ModuleBase::PI * parameters.length / (3.0 * parameters.area); + const double factor = 2.0 * ModuleBase::PI / (parameters.area * parameters.length); + return constant * moments.charge - factor * parabolic; +} + +ModuleBase::Vector3 pcc_2d_gradient(const ChargeMoments& moments, + const double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters) +{ + const double dipole = moments.dipole * normal; + const double factor = -4.0 * ModuleBase::PI / (parameters.area * parameters.length); + const double derivative = factor * (moments.charge * coordinate - dipole); + return normal * derivative; +} + +double pcc_2d_energy(const ChargeMoments& moments, const Pcc2dParameters& parameters) +{ + const double constant = -ModuleBase::PI * parameters.length / (3.0 * parameters.area); + const double factor = 2.0 * ModuleBase::PI / (parameters.area * parameters.length); + const double multipole = 2.0 * (moments.charge * moments.second_moment - moments.dipole.norm2()); + return 0.5 * (constant * moments.charge * moments.charge - factor * multipole); +} + +ModuleBase::Vector3 pcc_2d_force(const ChargeMoments& moments, + const double ionic_charge, + const double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters) +{ + const ModuleBase::Vector3 gradient = pcc_2d_gradient(moments, coordinate, normal, parameters); + return gradient * (-ionic_charge); +} +} // namespace elecstate diff --git a/source/source_estate/pcc_2d.h b/source/source_estate/pcc_2d.h new file mode 100644 index 00000000000..1c1afb12b98 --- /dev/null +++ b/source/source_estate/pcc_2d.h @@ -0,0 +1,45 @@ +#ifndef PCC_2D_H +#define PCC_2D_H + +#include "charge_moments.h" + +namespace unitcell +{ +struct SlabCell; +} + +namespace elecstate +{ +struct Pcc2dParameters +{ + double area = 0.0; ///< Periodic area in Bohr^2 + double length = 0.0; ///< Normal period in Bohr +}; + +Pcc2dParameters make_pcc_2d_parameters(const unitcell::SlabCell& cell); + +/// Hartree kernels. Parameters must be validated; moments use normal-projected +/// positions r = normal * u, so second_moment is the normal second moment. +/// The open planar kernel is zero on the source plane; its periodic counterpart +/// has zero mean. This gauge fixes the charged-slab constant to -pi*L/(3*A); +/// see M. J. Rutter, Electron. Struct. 3, 015002 (2021). +double pcc_2d_potential(const ChargeMoments& moments, + double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters); + +ModuleBase::Vector3 pcc_2d_gradient(const ChargeMoments& moments, + double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters); + +double pcc_2d_energy(const ChargeMoments& moments, const Pcc2dParameters& parameters); + +ModuleBase::Vector3 pcc_2d_force(const ChargeMoments& moments, + double ionic_charge, + double coordinate, + const ModuleBase::Vector3& normal, + const Pcc2dParameters& parameters); +} // namespace elecstate + +#endif diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index d3aecdbef57..a39b330c356 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -7,8 +7,38 @@ abacus_disable_feature_definitions(__ROCM) abacus_disable_feature_definitions(__MLALGO) abacus_disable_feature_definitions(_OPENMP) +AddTest( + TARGET MODULE_ESTATE_charge_moments + LIBS base device + SOURCES test_charge_moments.cpp ../charge_moments.cpp +) + +AddTest( + TARGET MODULE_ESTATE_pcc_0d + LIBS base device + SOURCES test_pcc_0d.cpp ../pcc_0d.cpp ../charge_moments.cpp ../../source_cell/cell_geometry.cpp +) + +AddTest( + TARGET MODULE_ESTATE_pcc_2d + LIBS base device + SOURCES test_pcc_2d.cpp ../pcc_2d.cpp ../charge_moments.cpp ../../source_cell/cell_geometry.cpp +) + if (ENABLE_MPI) +set(pcc_test_sources + ../module_pot/pot_pcc.cpp ../pcc_0d.cpp ../pcc_2d.cpp ../charge_moments.cpp + ../../source_cell/cell_geometry.cpp ../../source_cell/cell_tools.cpp + ../../source_cell/atom_spec.cpp ../../source_cell/atom_pseudo.cpp ../../source_cell/pseudo.cpp +) + +AddTest( + TARGET MODULE_ESTATE_pot_pcc + LIBS parameter base device planewave_serial + SOURCES test_pot_pcc.cpp ${pcc_test_sources} +) + # Copy at configure time so a plain `make` + `ctest` run finds the data. # install() only runs during `cmake --install`, which local test runs skip. file(COPY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) @@ -104,6 +134,7 @@ AddTest( LIBS parameter base device planewave_serial SOURCES potential_new_test.cpp td_field_test.cpp ../module_pot/potential_new.cpp ../module_pot/td_field.cpp ../module_pot/td_field_profiles.cpp ../module_pot/td_field_manager.cpp + ${pcc_test_sources} ) endif() diff --git a/source/source_estate/test/elecstate_base_test.cpp b/source/source_estate/test/elecstate_base_test.cpp index 317ef2805a2..8bcae6b32de 100644 --- a/source/source_estate/test/elecstate_base_test.cpp +++ b/source/source_estate/test/elecstate_base_test.cpp @@ -37,6 +37,10 @@ int main(int argc, char** argv) // Mock functions for testing elecstate.cpp namespace elecstate { +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/elecstate_energy_test.cpp b/source/source_estate/test/elecstate_energy_test.cpp index 136b1d08f78..0ccadc08dbf 100644 --- a/source/source_estate/test/elecstate_energy_test.cpp +++ b/source/source_estate/test/elecstate_energy_test.cpp @@ -12,6 +12,22 @@ int XC_Functional::func_type = 1; bool XC_Functional::ked_flag = false; namespace elecstate { +Potential::~Potential() {} +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::cal_v_eff(const Charge*, const UnitCell*, ModuleBase::matrix&) {} +void Potential::cal_fixed_v(double*) {} +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + +void Potential::get_solvation_energy(double& el, double& cav) const +{ + el = 0.8; + cav = 0.9; +} + +double Potential::pcc_energy_rydberg() const +{ + return 0.0; +} void Potential::get_vnew(Charge const*, ModuleBase::matrix&) { return; @@ -76,7 +92,7 @@ class ElecStateEnergyTest : public ::testing::Test /// explicitly, so the fixture owns them instead of staging them in the /// global parameter singleton. The values mirror the Input_para defaults, /// except sc_mag_switch, which the original fixture turned on. - bool imp_sol = false; + int imp_sol = 0; bool sc_mag_switch = true; int dft_plus_u = 0; std::string assume_isolated = "none"; @@ -267,3 +283,21 @@ TEST_F(ElecStateEnergyTest, CalBandgapUpDwBoundaryConditions) EXPECT_DOUBLE_EQ(elecstate->bandgap_up, 5.0); EXPECT_DOUBLE_EQ(elecstate->bandgap_dw, 5.0); } + +TEST_F(ElecStateEnergyTest, AddsSccsEnergyOnceAndResetsVacuum) +{ + elecstate::Potential potential; + elecstate->pot = &potential; + elecstate->f_en.deband = 0.1; + imp_sol = 2; + elecstate->cal_energies(2, imp_sol, sc_mag_switch, dft_plus_u, assume_isolated); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_el, 0.8); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_cav, 0.9); + EXPECT_DOUBLE_EQ(elecstate->f_en.etot, 2.4); + imp_sol = 0; + elecstate->cal_energies(2, imp_sol, sc_mag_switch, dft_plus_u, assume_isolated); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_el, 0.0); + EXPECT_DOUBLE_EQ(elecstate->f_en.esol_cav, 0.0); + EXPECT_DOUBLE_EQ(elecstate->f_en.etot, 0.7); + elecstate->pot = nullptr; +} diff --git a/source/source_estate/test/elecstate_pw_test.cpp b/source/source_estate/test/elecstate_pw_test.cpp index 1c47c8356b8..fc74f035816 100644 --- a/source/source_estate/test/elecstate_pw_test.cpp +++ b/source/source_estate/test/elecstate_pw_test.cpp @@ -26,6 +26,10 @@ class TestParameters int XC_Functional::func_type = 1; namespace elecstate { +void Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* Potential::solvent_electrostatic_potential() const { return nullptr; } + void Potential::init_pot(Charge const*) { } diff --git a/source/source_estate/test/potential_new_test.cpp b/source/source_estate/test/potential_new_test.cpp index 6dc01b9f321..026b2222d70 100644 --- a/source/source_estate/test/potential_new_test.cpp +++ b/source/source_estate/test/potential_new_test.cpp @@ -130,6 +130,15 @@ class MockPotComponent : public PotBase } } + void add_solvation_force(const UnitCell&, ModuleBase::matrix& force) const override + { + if (type_ != "solvation") { return; } + for (int ia = 0; ia < force.nr; ++ia) + { + for (int axis = 0; axis < force.nc; ++axis) { force(ia, axis) += 1.25; } + } + } + static void reset() { created = 0; @@ -153,7 +162,7 @@ int MockPotComponent::destroyed = 0; int MockPotComponent::fixed_calls = 0; int MockPotComponent::dynamic_calls = 0; -PotBase* Potential::get_pot_type(const std::string& pot_type) +PotBase* Potential::get_pot_type(const std::string& pot_type, const Input_para&) { const int grid_size = this->get_rho_basis() == nullptr ? 0 : this->get_rho_basis()->nrxx; return new MockPotComponent(pot_type, grid_size); @@ -411,3 +420,26 @@ TEST_F(PotentialNewTest, DifferentBasisObjectsRequireMatchingGammaOnly) EXPECT_EXIT(potential->update_from_charge(&charge, ucell.get()), ::testing::ExitedWithCode(1), ""); } + +TEST_F(PotentialNewTest, SolvationForcesUseRegisteredComponentsOnce) +{ + setup_smooth_basis(); + create_potential(smooth_basis.get(), smooth_basis.get()); + potential->pot_register({"fixed", "solvation"}); + ModuleBase::matrix force(2, 3); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { force(ia, axis) = 7.0; } + } + potential->add_solvation_force(*ucell, force); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } + } + potential->pot_register({"fixed"}); + potential->add_solvation_force(*ucell, force); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { EXPECT_DOUBLE_EQ(force(ia, axis), 8.25); } + } +} diff --git a/source/source_estate/test/test_charge_moments.cpp b/source/source_estate/test/test_charge_moments.cpp new file mode 100644 index 00000000000..77818118b93 --- /dev/null +++ b/source/source_estate/test/test_charge_moments.cpp @@ -0,0 +1,36 @@ +#include "source_estate/charge_moments.h" + +#include + +TEST(ChargeMoments, IntegratesAddsAndAcceptsAnEmptyLocalGrid) +{ + const double charges[2] = {2.0, -1.0}; + const ModuleBase::Vector3 positions[2] = { + ModuleBase::Vector3(1.0, 0.0, 0.0), + ModuleBase::Vector3(0.0, 2.0, 0.0)}; + const elecstate::ChargeMoments moments = elecstate::charge_moments(charges, positions, 2, 0.5); + EXPECT_DOUBLE_EQ(moments.charge, 0.5); + EXPECT_DOUBLE_EQ(moments.dipole.x, 1.0); + EXPECT_DOUBLE_EQ(moments.dipole.y, -1.0); + EXPECT_DOUBLE_EQ(moments.dipole.z, 0.0); + EXPECT_DOUBLE_EQ(moments.second_moment, -1.0); + + elecstate::ChargeMoments ions; + ions.charge = 2.0; + ions.dipole.x = 3.0; + ions.second_moment = 4.0; + elecstate::ChargeMoments electrons; + electrons.charge = -1.0; + electrons.dipole.y = -2.0; + electrons.second_moment = -3.0; + const elecstate::ChargeMoments total = elecstate::add_charge_moments(ions, electrons); + EXPECT_DOUBLE_EQ(total.charge, 1.0); + EXPECT_DOUBLE_EQ(total.dipole.x, 3.0); + EXPECT_DOUBLE_EQ(total.dipole.y, -2.0); + EXPECT_DOUBLE_EQ(total.second_moment, 1.0); + + electrons = elecstate::charge_moments(nullptr, nullptr, 0, 1.0); + EXPECT_DOUBLE_EQ(electrons.charge, 0.0); + EXPECT_DOUBLE_EQ(electrons.second_moment, 0.0); +} + diff --git a/source/source_estate/test/test_pcc_0d.cpp b/source/source_estate/test/test_pcc_0d.cpp new file mode 100644 index 00000000000..ac73b764eb7 --- /dev/null +++ b/source/source_estate/test/test_pcc_0d.cpp @@ -0,0 +1,199 @@ +#include "source_estate/pcc_0d.h" +#include "source_base/constants.h" +#include "source_base/matrix3.h" +#include "source_cell/cell_geometry.h" + +#include + +#include + +namespace +{ +/// Hartree energy per cell of periodic Gaussian charges (width sigma) under the +/// usual G = 0 exclusion. One charge at the origin, or a +q/-q pair separated +/// by `separation` along z when `dipole_pair` is true. +double periodic_gaussian_energy(const double length, + const double sigma, + const double charge, + const double separation, + const bool dipole_pair) +{ + const double step = 2.0 * ModuleBase::PI / length; + const int max_index = static_cast(8.0 / (sigma * step)) + 1; + const double volume = length * length * length; + double sum = 0.0; + for (int i = -max_index; i <= max_index; ++i) + { + for (int j = -max_index; j <= max_index; ++j) + { + for (int k = -max_index; k <= max_index; ++k) + { + if (i == 0 && j == 0 && k == 0) + { + continue; + } + const double gx = step * i; + const double gy = step * j; + const double gz = step * k; + const double g2 = gx * gx + gy * gy + gz * gz; + double structure = 1.0; + if (dipole_pair) + { + const double phase = std::sin(0.5 * gz * separation); + structure = 4.0 * phase * phase; + } + sum += charge * charge * structure * std::exp(-g2 * sigma * sigma) / g2; + } + } + } + return 2.0 * ModuleBase::PI * sum / volume; +} + +class Pcc0dTest : public testing::Test +{ + protected: + void SetUp() override + { + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + ASSERT_TRUE(unitcell::make_orthogonal_cell(lattice, 10.0, 1.0e-10, cell)); + ASSERT_TRUE(elecstate::make_pcc_0d_parameters(cell, 1.0e-10, parameters)); + moments.charge = 1.2; + moments.dipole = ModuleBase::Vector3(0.1, -0.2, 0.3); + moments.second_moment = 0.8; + } + + unitcell::OrthogonalCell cell; + elecstate::Pcc0dParameters parameters; + elecstate::ChargeMoments moments; +}; +} // namespace + +TEST_F(Pcc0dTest, SeparatesCubicRestrictionFromOrthogonalGeometry) +{ + EXPECT_DOUBLE_EQ(parameters.length, 10.0); + const ModuleBase::Matrix3 rectangular(1.0, 0.0, 0.0, + 0.0, 2.0, 0.0, + 0.0, 0.0, 1.0); + ASSERT_TRUE(unitcell::make_orthogonal_cell(rectangular, 10.0, 1.0e-10, cell)); + EXPECT_FALSE(elecstate::make_pcc_0d_parameters(cell, 1.0e-10, parameters)); +} + +// The gradient, the ionic force and the electron potential are derivatives of +// the potential and of the energy with respect to a point charge. +TEST_F(Pcc0dTest, DerivativesMatchFiniteDifferences) +{ + const double step = 1.0e-5; + const ModuleBase::Vector3 position(0.4, -0.5, 0.6); + const ModuleBase::Vector3 gradient = elecstate::pcc_0d_gradient(moments, position, parameters); + const double charges[2] = {1.3, -0.4}; + ModuleBase::Vector3 positions[2] = {position, ModuleBase::Vector3(-0.5, 0.6, -0.7)}; + const elecstate::ChargeMoments ions = elecstate::charge_moments(charges, positions, 2, 1.0); + const ModuleBase::Vector3 force = elecstate::pcc_0d_force(ions, charges[0], position, parameters); + for (int axis = 0; axis < 3; ++axis) + { + ModuleBase::Vector3 plus = position; + ModuleBase::Vector3 minus = position; + plus[axis] += step; + minus[axis] -= step; + const double plus_potential = elecstate::pcc_0d_potential(moments, plus, parameters); + const double minus_potential = elecstate::pcc_0d_potential(moments, minus, parameters); + const double potential_slope = (plus_potential - minus_potential) / (2.0 * step); + EXPECT_NEAR(potential_slope, gradient[axis], 1.0e-10); + positions[0] = plus; + const elecstate::ChargeMoments plus_ions = elecstate::charge_moments(charges, positions, 2, 1.0); + positions[0] = minus; + const elecstate::ChargeMoments minus_ions = elecstate::charge_moments(charges, positions, 2, 1.0); + const double plus_energy = elecstate::pcc_0d_energy(plus_ions, parameters); + const double minus_energy = elecstate::pcc_0d_energy(minus_ions, parameters); + const double numerical_force = -(plus_energy - minus_energy) / (2.0 * step); + EXPECT_NEAR(numerical_force, force[axis], 1.0e-10); + } + // Adding electrons -step at position changes the energy by the electron potential (Ry). + const double positive_potential = elecstate::pcc_0d_potential(moments, position, parameters); + const double electron_potential_ry = -2.0 * positive_potential; + elecstate::ChargeMoments plus = moments; + plus.charge -= step; + plus.dipole -= position * step; + plus.second_moment -= step * position.norm2(); + elecstate::ChargeMoments minus = moments; + minus.charge += step; + minus.dipole += position * step; + minus.second_moment += step * position.norm2(); + const double plus_energy = 2.0 * elecstate::pcc_0d_energy(plus, parameters); + const double minus_energy = 2.0 * elecstate::pcc_0d_energy(minus, parameters); + const double energy_slope = (plus_energy - minus_energy) / (2.0 * step); + EXPECT_NEAR(energy_slope, electron_potential_ry, 1.0e-10); +} + +TEST_F(Pcc0dTest, EnergyIsIndependentOfMultipoleOrigin) +{ + const double charges[2] = {1.2, -0.3}; + ModuleBase::Vector3 positions[2] = { + ModuleBase::Vector3(0.3, -0.2, 0.1), + ModuleBase::Vector3(-0.5, 0.6, -0.7)}; + moments = elecstate::charge_moments(charges, positions, 2, 1.0); + const double original_energy = elecstate::pcc_0d_energy(moments, parameters); + const ModuleBase::Vector3 shift(0.4, -0.3, 0.2); + positions[0] -= shift; + positions[1] -= shift; + moments = elecstate::charge_moments(charges, positions, 2, 1.0); + const double shifted_energy = elecstate::pcc_0d_energy(moments, parameters); + EXPECT_NEAR(original_energy, shifted_energy, 1.0e-14); +} + +// Physical check of the model: a charged Gaussian ion (H3O+-like) in cubic +// cells of growing size. The periodic Hartree energy drifts as 1/L, while +// periodic + PCC reproduces the isolated self-energy q^2 / (2 sigma sqrt(pi)). +TEST_F(Pcc0dTest, RecoversIsolatedEnergyOfChargedGaussianForAnyCellSize) +{ + const double sigma = 1.0; + const double charge = 1.0; + const double isolated = charge * charge / (2.0 * sigma * std::sqrt(ModuleBase::PI)); + const double lengths[3] = {16.0, 20.0, 24.0}; + for (int index = 0; index < 3; ++index) + { + const double length = lengths[index]; + parameters.length = length; + elecstate::ChargeMoments gaussian; + gaussian.charge = charge; + gaussian.second_moment = 3.0 * sigma * sigma * charge; + const double periodic = periodic_gaussian_energy(length, sigma, charge, 0.0, false); + const double correction = elecstate::pcc_0d_energy(gaussian, parameters); + const double periodic_error = std::abs(periodic - isolated); + EXPECT_GT(periodic_error, 1.0e-2); + EXPECT_NEAR(periodic + correction, isolated, 1.0e-10); + } +} + +// Neutral polar molecule modelled as a +q/-q Gaussian pair. PCC removes the +// leading 1/L^3 dipole interaction; the residual error comes from higher +// multipoles and must be far smaller and decay faster than the uncorrected one. +TEST_F(Pcc0dTest, RemovesDipoleImageInteractionOfNeutralGaussianPair) +{ + const double sigma = 1.0; + const double charge = 1.0; + const double separation = 2.0; + const double self_energy = charge * charge / (2.0 * sigma * std::sqrt(ModuleBase::PI)); + const double attraction = charge * charge * std::erf(0.5 * separation / sigma) / separation; + const double isolated = 2.0 * self_energy - attraction; + elecstate::ChargeMoments pair; + pair.dipole.z = charge * separation; + const double lengths[2] = {16.0, 24.0}; + double corrected_errors[2] = {0.0, 0.0}; + for (int index = 0; index < 2; ++index) + { + const double length = lengths[index]; + parameters.length = length; + const double periodic = periodic_gaussian_energy(length, sigma, charge, separation, true); + const double correction = elecstate::pcc_0d_energy(pair, parameters); + const double periodic_error = std::abs(periodic - isolated); + const double corrected_error = std::abs(periodic + correction - isolated); + EXPECT_LT(corrected_error, 0.05 * periodic_error); + corrected_errors[index] = corrected_error; + } + const double ratio = lengths[1] / lengths[0]; + const double dipole_decay = std::pow(ratio, 3.0); + EXPECT_GT(corrected_errors[0] / corrected_errors[1], dipole_decay); +} diff --git a/source/source_estate/test/test_pcc_2d.cpp b/source/source_estate/test/test_pcc_2d.cpp new file mode 100644 index 00000000000..0065d5aefac --- /dev/null +++ b/source/source_estate/test/test_pcc_2d.cpp @@ -0,0 +1,139 @@ +#include "source_estate/pcc_2d.h" +#include "source_base/constants.h" +#include "source_cell/cell_geometry.h" + +#include + +#include + +namespace +{ +/// Hartree energy per cell of periodic Gaussian sheets (surface charge +/// `density`, normal width sigma) with the G = 0 term excluded. One sheet at +/// the origin, or a +density/-density pair separated by `separation` along the +/// normal when `dipole_pair` is true. Only G parallel to the normal contributes. +double periodic_sheet_energy(const double area, + const double length, + const double sigma, + const double density, + const double separation, + const bool dipole_pair) +{ + const double step = 2.0 * ModuleBase::PI / length; + const int max_index = static_cast(8.0 / (sigma * step)) + 1; + double sum = 0.0; + for (int k = 1; k <= max_index; ++k) + { + const double g = step * k; + double structure = 1.0; + if (dipole_pair) + { + const double phase = std::sin(0.5 * g * separation); + structure = 4.0 * phase * phase; + } + sum += density * density * structure * std::exp(-g * g * sigma * sigma) / (g * g); + } + return 4.0 * ModuleBase::PI * area * sum / length; +} +} // namespace + +// The ionic force is normal and the electron potential (Ry) is the energy +// derivative with respect to a point charge. +TEST(Pcc2d, ForceAndElectronPotentialMatchEnergyDerivatives) +{ + elecstate::Pcc2dParameters parameters; + parameters.area = 6.0; + parameters.length = 8.0; + const double step = 1.0e-5; + const ModuleBase::Vector3 normal(0.0, 1.0, 0.0); + const double charges[2] = {1.3, -0.4}; + ModuleBase::Vector3 positions[2] = { + ModuleBase::Vector3(0.0, 0.2, 0.0), + ModuleBase::Vector3(0.0, -0.5, 0.0)}; + const elecstate::ChargeMoments moments = elecstate::charge_moments(charges, positions, 2, 1.0); + const double coordinate = positions[0].y; + const ModuleBase::Vector3 force + = elecstate::pcc_2d_force(moments, charges[0], coordinate, normal, parameters); + EXPECT_DOUBLE_EQ(force.x, 0.0); + EXPECT_DOUBLE_EQ(force.z, 0.0); + positions[0].y = coordinate + step; + const elecstate::ChargeMoments plus_moments = elecstate::charge_moments(charges, positions, 2, 1.0); + positions[0].y = coordinate - step; + const elecstate::ChargeMoments minus_moments = elecstate::charge_moments(charges, positions, 2, 1.0); + const double plus = elecstate::pcc_2d_energy(plus_moments, parameters); + const double minus = elecstate::pcc_2d_energy(minus_moments, parameters); + const double numerical_force = -(plus - minus) / (2.0 * step); + EXPECT_NEAR(force.y, numerical_force, 1.0e-10); + + const double positive_potential = elecstate::pcc_2d_potential(moments, coordinate, normal, parameters); + const double potential_ry = -2.0 * positive_potential; + const ModuleBase::Vector3 dipole_step = normal * (coordinate * step); + const double second_step = coordinate * coordinate * step; + elecstate::ChargeMoments removed = moments; + removed.charge -= step; + removed.dipole -= dipole_step; + removed.second_moment -= second_step; + elecstate::ChargeMoments added = moments; + added.charge += step; + added.dipole += dipole_step; + added.second_moment += second_step; + const double removed_energy = 2.0 * elecstate::pcc_2d_energy(removed, parameters); + const double added_energy = 2.0 * elecstate::pcc_2d_energy(added, parameters); + const double energy_slope = (removed_energy - added_energy) / (2.0 * step); + EXPECT_NEAR(energy_slope, potential_ry, 1.0e-10); +} + +// Physical check of the charged-slab gauge: a charged Gaussian layer in cells +// of growing vacuum. periodic + PCC must give the open-boundary energy with the +// potential referenced to zero on the plane of the charge, -2 sqrt(pi) A s^2 sigma, +// independently of the vacuum size. +TEST(Pcc2d, RecoversOpenBoundaryEnergyOfChargedLayerForAnyVacuum) +{ + const double area = 50.0; + const double sigma = 1.0; + const double density = 0.1; + const double charge = area * density; + const double isolated = -2.0 * std::sqrt(ModuleBase::PI) * area * density * density * sigma; + const double lengths[3] = {20.0, 30.0, 40.0}; + for (int index = 0; index < 3; ++index) + { + elecstate::Pcc2dParameters parameters; + parameters.area = area; + parameters.length = lengths[index]; + elecstate::ChargeMoments layer; + layer.charge = charge; + layer.second_moment = charge * sigma * sigma; + const double periodic = periodic_sheet_energy(area, lengths[index], sigma, density, 0.0, false); + const double correction = elecstate::pcc_2d_energy(layer, parameters); + EXPECT_NEAR(periodic + correction, isolated, 1.0e-10); + } +} + +// Polar neutral slab modelled as two opposite Gaussian layers. Periodic images +// impose a spurious zero-average field that changes the energy by ~1/L; PCC +// restores the isolated dipole-layer energy for every vacuum size. +TEST(Pcc2d, RecoversIsolatedEnergyOfPolarNeutralSlabForAnyVacuum) +{ + const double area = 50.0; + const double sigma = 1.0; + const double density = 0.1; + const double separation = 3.0; + const double same_layer = 2.0 * sigma / std::sqrt(ModuleBase::PI); + const double gap = 0.5 * separation / sigma; + const double cross_layer = separation * std::erf(gap) + same_layer * std::exp(-gap * gap); + const double isolated = -2.0 * ModuleBase::PI * area * density * density * (same_layer - cross_layer); + const double lengths[3] = {20.0, 30.0, 40.0}; + for (int index = 0; index < 3; ++index) + { + elecstate::Pcc2dParameters parameters; + parameters.area = area; + parameters.length = lengths[index]; + elecstate::ChargeMoments slab; + slab.dipole.z = area * density * separation; + const double periodic = periodic_sheet_energy(area, lengths[index], sigma, density, separation, true); + const double correction = elecstate::pcc_2d_energy(slab, parameters); + const double periodic_error = std::abs(periodic - isolated); + EXPECT_GT(periodic_error, 0.1); + EXPECT_NEAR(periodic + correction, isolated, 1.0e-10); + } +} diff --git a/source/source_estate/test/test_pot_pcc.cpp b/source/source_estate/test/test_pot_pcc.cpp new file mode 100644 index 00000000000..f4cadd7819c --- /dev/null +++ b/source/source_estate/test/test_pot_pcc.cpp @@ -0,0 +1,213 @@ +#include "source_estate/module_pot/pot_pcc.h" +#include "source_base/constants.h" +#include "source_base/global_variable.h" + +#include + +#include +#include +#include + +// These fixtures use public cell/charge storage without running INPUT setup. +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +SepPot::SepPot() {} +SepPot::~SepPot() {} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} +Charge::Charge() {} +Charge::~Charge() {} + +namespace +{ +class PotPccTest : public testing::Test +{ + protected: + void SetUp() override + { + basis.nx = 2; + basis.ny = 2; + basis.nz = 2; + basis.nplane = 2; + basis.startz_current = 0; + basis.nrxx = 8; + basis.nxyz = 8; + cell.lat0 = 10.0; + cell.omega = 1000.0; + cell.latvec = ModuleBase::Matrix3(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + cell.ntype = 1; + cell.nat = 1; + atom.na = 1; + atom.mass = 1.0; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.45, 0.45, 0.45)}; + cell.atoms = &atom; + density.assign(8, 0.0); + density[7] = 1.0 / 125.0; + channels[0] = density.data(); + charge.rho = channels; + charge.nspin = 1; + } + + ModulePW::PW_Basis basis{"cpu", "double"}; + UnitCell cell; + Atom atom; + Charge charge; + std::vector density; + double* channels[2] = {nullptr, nullptr}; +}; +} // namespace + +TEST_F(PotPccTest, NeutralDipoleProducesConsistentEnergyPotentialAndForce) +{ + elecstate::PotPcc correction(&basis, 1.0); + ModuleBase::matrix potential(1, 8); + correction.cal_v_eff(&charge, &cell, potential); + const double expected_energy = 4.0 * ModuleBase::PI * 0.75 / 3000.0; + EXPECT_NEAR(correction.get_energy(), expected_energy, 1.0e-14); + EXPECT_NEAR(potential(0, 7), expected_energy, 1.0e-14); + const std::vector& stored = correction.electron_potential(); + ASSERT_EQ(stored.size(), 8u); + EXPECT_DOUBLE_EQ(stored[7], potential(0, 7)); + ModuleBase::matrix force(1, 3); + correction.add_force(cell, force); + const double expected_force = 8.0 * ModuleBase::PI * 0.5 / 3000.0; + EXPECT_NEAR(force(0, 0), expected_force, 1.0e-14); + EXPECT_NEAR(force(0, 1), expected_force, 1.0e-14); + EXPECT_NEAR(force(0, 2), expected_force, 1.0e-14); +} + +TEST_F(PotPccTest, RecomputesDensityKeepsInstancesIndependentAndSharesSpinChannels) +{ + elecstate::PotPcc first(&basis, 1.0); + elecstate::PotPcc second(&basis, 1.0); + ModuleBase::matrix potential(1, 8); + first.cal_v_eff(&charge, &cell, potential); + const double neutral_energy = first.get_energy(); + density[7] = 2.0 / 125.0; + potential.zero_out(); + second.cal_v_eff(&charge, &cell, potential); + const double charged_energy = second.get_energy(); + const double expected_charged = 2.837297479480619 / 10.0 + 4.0 * ModuleBase::PI * 1.5 / 3000.0; + EXPECT_NEAR(charged_energy, expected_charged, 1.0e-14); + EXPECT_DOUBLE_EQ(first.get_energy(), neutral_energy); + potential.zero_out(); + first.cal_v_eff(&charge, &cell, potential); + EXPECT_NEAR(first.get_energy(), charged_energy, 1.0e-14); + + // Both spin channels get the same correction on top of their potentials. + std::vector down(8, 0.0); + density[7] = 0.3 / 125.0; + down[7] = 0.7 / 125.0; + channels[1] = down.data(); + charge.nspin = 2; + elecstate::PotPcc correction(&basis, 1.0); + ModuleBase::matrix spin_potential(2, 8); + spin_potential(0, 7) = 2.0; + spin_potential(1, 7) = 3.0; + correction.cal_v_eff(&charge, &cell, spin_potential); + const double energy = correction.get_energy(); + const double first_expected = 2.0 + energy; + const double second_expected = 3.0 + energy; + EXPECT_NEAR(spin_potential(0, 7), first_expected, 1.0e-14); + EXPECT_NEAR(spin_potential(1, 7), second_expected, 1.0e-14); +} + +TEST_F(PotPccTest, ForceMatchesFixedDensityEnergyDerivativeAfterMovingIons) +{ + elecstate::PotPcc correction(&basis, 1.0); + ModuleBase::matrix potential(1, 8); + correction.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix force(1, 3); + correction.add_force(cell, force); + const double step_bohr = 1.0e-4; + const double step_lat0 = step_bohr / cell.lat0; + atom.tau[0].x += step_lat0; + potential.zero_out(); + correction.cal_v_eff(&charge, &cell, potential); + const double plus_energy = correction.get_energy(); + atom.tau[0].x -= 2.0 * step_lat0; + potential.zero_out(); + correction.cal_v_eff(&charge, &cell, potential); + const double minus_energy = correction.get_energy(); + const double numerical_force = -(plus_energy - minus_energy) / (2.0 * step_bohr); + EXPECT_NEAR(force(0, 0), numerical_force, 1.0e-10); +} + +TEST_F(PotPccTest, SlabAxisControlsPotentialEnergyForceAndIgnoresInPlaneCoordinates) +{ + const double expected_energy = 4.0 * ModuleBase::PI * 0.25 / 1000.0; + const double expected_force = 8.0 * ModuleBase::PI * 0.5 / 1000.0; + for (int axis = 0; axis < 3; ++axis) + { + elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, axis, 1.0); + ModuleBase::matrix potential(1, 8); + correction.cal_v_eff(&charge, &cell, potential); + EXPECT_NEAR(correction.get_energy(), expected_energy, 1.0e-14); + EXPECT_NEAR(potential(0, 7), expected_energy, 1.0e-14); + ModuleBase::matrix force(1, 3); + correction.add_force(cell, force); + for (int component = 0; component < 3; ++component) + { + const double expected = (component == axis) ? expected_force : 0.0; + EXPECT_NEAR(force(0, component), expected, 1.0e-14); + } + } +} + +TEST_F(PotPccTest, RejectsUnsupportedCellsKpointsAndSymmetry) +{ + const std::vector> points = {ModuleBase::Vector3(0.25, 0.0, 0.0)}; + elecstate::PotPcc::validate_kpoints(points, 1, 2); + EXPECT_EXIT(elecstate::PotPcc::validate_kpoints(points, 1, 0), testing::ExitedWithCode(1), ""); + elecstate::PotPcc correction(&basis, elecstate::PotPcc::Dimension::slab, 2, 1.0); + cell.symm.ptrans = {ModuleBase::Vector3(0.5, 0.0, 0.0)}; + ModuleBase::matrix potential(1, 8); + correction.cal_v_eff(&charge, &cell, potential); + cell.symm.ptrans[0].z = 0.5; + EXPECT_EXIT(correction.cal_v_eff(&charge, &cell, potential), testing::ExitedWithCode(1), ""); + cell.symm.ptrans.clear(); + cell.symm.nrotk = 1; + cell.symm.gmatrix[0] = ModuleBase::Matrix3(0.0, 0.0, 1.0, + 0.0, 1.0, 0.0, + 1.0, 0.0, 0.0); + EXPECT_EXIT(correction.cal_v_eff(&charge, &cell, potential), testing::ExitedWithCode(1), ""); + cell.symm.nrotk = 0; + cell.latvec.e22 = 2.0; + elecstate::PotPcc molecule(&basis, 1.0); + EXPECT_EXIT(molecule.cal_v_eff(&charge, &cell, potential), testing::ExitedWithCode(1), ""); + cell.latvec.e22 = 1.0; + cell.latvec.e31 = 0.1; + elecstate::PotPcc tilted_slab(&basis, elecstate::PotPcc::Dimension::slab, 2, 1.0); + EXPECT_EXIT(tilted_slab.cal_v_eff(&charge, &cell, potential), testing::ExitedWithCode(1), ""); +} + +TEST_F(PotPccTest, ReportsGridChargeMismatchAndUsesGridCharge) +{ + const std::string log_name = "pot_pcc_charge_warning.log"; + std::ofstream& warning_log = GlobalV::ofs_warning; + warning_log.open(log_name.c_str()); + ModuleBase::matrix matched_potential(1, 8); + elecstate::PotPcc matched(&basis, 1.0); + matched.cal_v_eff(&charge, &cell, matched_potential); + warning_log.flush(); + std::ifstream matched_log(log_name.c_str()); + std::stringstream matched_text; + matched_text << matched_log.rdbuf(); + EXPECT_EQ(matched_text.str().find("PCC net charge"), std::string::npos); + + ModuleBase::matrix mismatched_potential(1, 8); + elecstate::PotPcc mismatched(&basis, 2.0); + mismatched.cal_v_eff(&charge, &cell, mismatched_potential); + warning_log.close(); + std::ifstream mismatched_log(log_name.c_str()); + std::stringstream mismatched_text; + mismatched_text << mismatched_log.rdbuf(); + EXPECT_NE(mismatched_text.str().find("PCC net charge differs from the expected value"), std::string::npos); + EXPECT_DOUBLE_EQ(mismatched.get_energy(), matched.get_energy()); + std::remove(log_name.c_str()); +} diff --git a/source/source_estate/test/test_pot_sccs.cpp b/source/source_estate/test/test_pot_sccs.cpp new file mode 100644 index 00000000000..4bf0aa26977 --- /dev/null +++ b/source/source_estate/test/test_pot_sccs.cpp @@ -0,0 +1,244 @@ +#include "source_estate/module_pot/pot_sccs.h" +#include "source_io/module_parameter/input_parameter.h" +#include "source_hamilt/module_sccs/test/sccs_test.h" +#include "source_base/global_variable.h" + +#include +#include +#include +#include +#include + +// Public storage fixtures: no INPUT initialization or privately owned atom maps. +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +SepPot::SepPot() {} +SepPot::~SepPot() {} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} +Charge::Charge() {} +Charge::~Charge() {} + +class PotSccsTest : public SccsTest::PwTest +{ +protected: + // Two atoms in a uniform electron density: the solvation force is the + // finite difference of the Rydberg energy and is added once to an existing + // force. + void check_solvation_force(const ModuleSccs::SccsConfig& config, + const ModuleSccs::PolarizationSolverParameters& solver) + { + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 2; + cell.atoms = &atom; + atom.na = 2; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.21, 0.32, 0.43), + ModuleBase::Vector3(0.64, 0.51, 0.27)}; + const double density_value = 2.0 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + elecstate::PotSccs component(&basis, config, solver, 2.0); + ModuleBase::matrix potential(1, basis.nrxx); + component.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix force(2, 3); + component.add_solvation_force(cell, force); + const double step = 1e-4; + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + const double original = atom.tau[ia][axis]; + atom.tau[ia][axis] = original + step / length; + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + const double positive = component.get_energy(); + atom.tau[ia][axis] = original - step / length; + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + const double negative = component.get_energy(); + atom.tau[ia][axis] = original; + const double finite_difference = -(positive - negative) / (2.0 * step); + EXPECT_NEAR(force(ia, axis), finite_difference, 1e-8); + } + } + potential.zero_out(); + component.cal_v_eff(&charge, &cell, potential); + ModuleBase::matrix accumulated(2, 3); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) { accumulated(ia, axis) = 7.0; } + } + component.add_solvation_force(cell, accumulated); + for (int ia = 0; ia < 2; ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + const double expected_force = 7.0 + force(ia, axis); + EXPECT_NEAR(accumulated(ia, axis), expected_force, 1e-10); + } + } + } +}; + +TEST_F(PotSccsTest, IndependentInstancesAndSpinChannels) +{ + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + const double density_value = 1.0 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data(), density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + Input_para input; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + elecstate::PotSccs first(&basis, config, solver, 1.0); + config.cavity.epsilon_bulk = 1.0; + elecstate::PotSccs vacuum(&basis, config, solver, 1.0); + EXPECT_DOUBLE_EQ(first.get_energy(), 0.0); + ModuleBase::matrix potential(1, basis.nrxx); + first.cal_v_eff(&charge, &cell, potential); + const double first_energy = first.get_energy(); + EXPECT_LT(first_energy, 0.0); + double el = 0.0; + double cav = 0.0; + first.get_solvation_energy(el, cav); + EXPECT_DOUBLE_EQ(el, first_energy); + EXPECT_DOUBLE_EQ(cav, 0.0); + ModuleBase::matrix zero(1, basis.nrxx); + vacuum.cal_v_eff(&charge, &cell, zero); + EXPECT_NEAR(vacuum.get_energy(), 0.0, 1e-12); + EXPECT_DOUBLE_EQ(first.get_energy(), first_energy); + for (int ir = 0; ir < basis.nrxx; ++ir) { EXPECT_NEAR(zero(0, ir), 0.0, 1e-11); } + for (double& value : density) { value *= 0.5; } + charge.nspin = 2; + ModuleBase::matrix spin_potential(2, basis.nrxx); + first.cal_v_eff(&charge, &cell, spin_potential); + EXPECT_NEAR(first.get_energy(), first_energy, 1e-12); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(spin_potential(0, ir), potential(0, ir), 1e-11); + EXPECT_DOUBLE_EQ(spin_potential(0, ir), spin_potential(1, ir)); + } +} + +TEST_F(PotSccsTest, InputMapsOntoSccsConfig) +{ + // The custom defaults follow Environ environ_type input: + // env_static_permittivity 1, rhomin 1e-4, rhomax 5e-3, no surface tension + // or pressure. + Input_para input; + ModuleSccs::SccsConfig defaults; + ModuleSccs::PolarizationSolverParameters default_solver; + elecstate::make_sccs_config_from_input(input, defaults, default_solver); + EXPECT_EQ(input.sccs_preset, "custom"); + EXPECT_DOUBLE_EQ(defaults.cavity.epsilon_bulk, 1.0); + EXPECT_DOUBLE_EQ(defaults.cavity.density_min, 1e-4); + EXPECT_DOUBLE_EQ(defaults.cavity.density_max, 5e-3); + EXPECT_DOUBLE_EQ(defaults.surface_tension, 0.0); + EXPECT_DOUBLE_EQ(defaults.pressure, 0.0); + EXPECT_EQ(default_solver.max_iterations, 200); + input.sccs_preset = "water-neutral"; + input.sccs_epsilon = 2.0; // replaced by the preset + input.sccs_maxiter = 42; + input.sccs_surface_eta = 1e-6; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); + EXPECT_DOUBLE_EQ(config.cavity.density_min, 1e-4); + EXPECT_LT(config.pressure, 0.0); + EXPECT_DOUBLE_EQ(config.surface_regularization, 1e-6); + EXPECT_EQ(solver.max_iterations, 42); +} + +TEST_F(PotSccsTest, IonicForceAddsRydbergDerivativeOnce) +{ + Input_para input; + input.sccs_epsilon = 5.0; + input.sccs_tol_rms = 1e-13; + input.sccs_tol_max = 1e-12; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + check_solvation_force(config, solver); +} + +TEST_F(PotSccsTest, ChargedDielectricCellWarnsAndRuns) +{ + UnitCell cell; + Atom atom; + cell.lat0 = length; + cell.tpiba = tpiba; + cell.omega = basis.omega; + cell.ntype = 1; + cell.nat = 1; + cell.atoms = &atom; + atom.na = 1; + atom.ncpp.zv = 1.0; + atom.tau = {ModuleBase::Vector3(0.25, 0.25, 0.25)}; + Input_para input; + ModuleSccs::SccsConfig config; + ModuleSccs::PolarizationSolverParameters solver; + elecstate::make_sccs_config_from_input(input, config, solver); + config.cavity.epsilon_bulk = 5.0; + EXPECT_TRUE(elecstate::check_sccs_charge(config, cell, 1.0).empty()); + EXPECT_NE(elecstate::check_sccs_charge(config, cell, 0.5).find("G = 0"), std::string::npos); + ModuleSccs::SccsConfig vacuum_config = config; + vacuum_config.cavity.epsilon_bulk = 1.0; + EXPECT_TRUE(elecstate::check_sccs_charge(vacuum_config, cell, 0.5).empty()); + + const double density_value = 0.5 / basis.omega; + std::vector density(basis.nrxx, density_value); + double* channels[] = {density.data()}; + Charge charge; + charge.nspin = 1; + charge.rho = channels; + // One log per rank: the MPI variant runs this test on several ranks at once. + const std::string log_name = "pot_sccs_charge_warning_" + std::to_string(basis.poolrank) + ".log"; + std::ofstream& warning_log = GlobalV::ofs_warning; + warning_log.open(log_name.c_str()); + elecstate::PotSccs matched(&basis, config, solver, 0.5); + ModuleBase::matrix matched_potential(1, basis.nrxx); + matched.cal_v_eff(&charge, &cell, matched_potential); + warning_log.flush(); + std::ifstream matched_log(log_name.c_str()); + std::stringstream matched_text; + matched_text << matched_log.rdbuf(); + EXPECT_EQ(matched_text.str().find("SCCS grid electron count"), std::string::npos); + EXPECT_TRUE(std::isfinite(matched.get_energy())); + + elecstate::PotSccs mismatched(&basis, config, solver, 1.0); + ModuleBase::matrix mismatched_potential(1, basis.nrxx); + mismatched.cal_v_eff(&charge, &cell, mismatched_potential); + warning_log.close(); + std::ifstream mismatched_log(log_name.c_str()); + std::stringstream mismatched_text; + mismatched_text << mismatched_log.rdbuf(); + EXPECT_NE(mismatched_text.str().find("SCCS grid electron count differs"), std::string::npos); + EXPECT_DOUBLE_EQ(mismatched.get_energy(), matched.get_energy()); + std::remove(log_name.c_str()); +} diff --git a/source/source_estate/write_elecstat_pot.cpp b/source/source_estate/write_elecstat_pot.cpp index 15a5363ed62..0872d9e132c 100644 --- a/source/source_estate/write_elecstat_pot.cpp +++ b/source/source_estate/write_elecstat_pot.cpp @@ -22,6 +22,8 @@ void write_elecstat_pot( const UnitCell* ucell, const double* v_eff, const surchem& solvent, + const std::vector* pcc_potential, + const std::vector* solvent_potential, const int precision) { ModuleBase::TITLE("ModuleIO", "write_elecstat_pot"); @@ -29,10 +31,15 @@ void write_elecstat_pot( std::vector v_elecstat(rho_basis->nrxx, 0.0); + if (solvent_potential != nullptr && solvent_potential->size() != static_cast(rho_basis->nrxx)) + { + ModuleBase::WARNING_QUIT("write_elecstat_pot", "Solvent potential does not match the output grid"); + } + const int nspin = PARAM.inp.nspin; const int efield = PARAM.inp.efield_flag; const int dip_corr = PARAM.inp.dip_cor_flag; - const bool imp_sol = PARAM.inp.imp_sol; + const int imp_sol = PARAM.inp.imp_sol; //========================================== // Hartree potential @@ -66,7 +73,15 @@ void write_elecstat_pot( { v_elecstat[ir] += v_efield(0, ir); } - if(imp_sol == true) + if (pcc_potential != nullptr) + { + v_elecstat[ir] += (*pcc_potential)[ir]; + } + if (solvent_potential != nullptr) + { + v_elecstat[ir] += (*solvent_potential)[ir]; + } + if (imp_sol == 1) { v_elecstat[ir] += solvent.delta_phi[ir]; } diff --git a/source/source_estate/write_elecstat_pot.h b/source/source_estate/write_elecstat_pot.h index bee575b95c3..f21239935b4 100644 --- a/source/source_estate/write_elecstat_pot.h +++ b/source/source_estate/write_elecstat_pot.h @@ -33,6 +33,8 @@ void write_elecstat_pot( const UnitCell* ucell_, const double* v_eff_fixed, const surchem& solvent, + const std::vector* pcc_potential, + const std::vector* solvent_potential, const int precision); } // namespace ModuleIO diff --git a/source/source_hamilt/CMakeLists.txt b/source/source_hamilt/CMakeLists.txt index 90fe2257c77..d0cdba7a1ef 100644 --- a/source/source_hamilt/CMakeLists.txt +++ b/source/source_hamilt/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(module_vdw) add_subdirectory(module_surchem) +add_subdirectory(module_sccs) add_subdirectory(module_xc) add_subdirectory(module_hcontainer) diff --git a/source/source_hamilt/module_sccs/CMakeLists.txt b/source/source_hamilt/module_sccs/CMakeLists.txt new file mode 100644 index 00000000000..f76b46ec716 --- /dev/null +++ b/source/source_hamilt/module_sccs/CMakeLists.txt @@ -0,0 +1,16 @@ +add_library(sccs OBJECT + sccs_cavity.cpp + sccs_parameters.cpp + sccs_functional.cpp + sccs_ionic_charge.cpp + sccs_ionic_force.cpp + sccs_pw_coulomb.cpp + sccs_response.cpp) + +if(ENABLE_COVERAGE) + add_coverage(sccs) +endif() + +if(BUILD_TESTING) + add_subdirectory(test) +endif() diff --git a/source/source_hamilt/module_sccs/sccs_cavity.cpp b/source/source_hamilt/module_sccs/sccs_cavity.cpp new file mode 100644 index 00000000000..94f69ed0abc --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cavity.cpp @@ -0,0 +1,39 @@ +#include "sccs_cavity.h" +#include "source_base/constants.h" + +#include + +namespace ModuleSccs +{ +CavityPoint evaluate_cavity(double density, const CavityParameters& parameters) +{ + CavityPoint point; + if (density <= parameters.density_min) + { + point.epsilon = parameters.epsilon_bulk; + } + else if (density >= parameters.density_max) + { + point.solute = 1.0; + point.epsilon = 1.0; + } + else + { + const double density_ratio = parameters.density_max / parameters.density_min; + const double log_width = std::log(density_ratio); + const double relative_density = parameters.density_max / density; + const double log_relative_density = std::log(relative_density); + const double x = log_relative_density / log_width; + const double angle = ModuleBase::TWO_PI * x; + const double solvent = x - std::sin(angle) / ModuleBase::TWO_PI; + const double dsolvent_drho = -(1.0 - std::cos(angle)) / (log_width * density); + const double log_epsilon = std::log(parameters.epsilon_bulk); + const double log_epsilon_point = log_epsilon * solvent; + point.solute = 1.0 - solvent; + point.dsolute_drho = -dsolvent_drho; + point.epsilon = std::exp(log_epsilon_point); + point.depsilon_drho = point.epsilon * log_epsilon * dsolvent_drho; + } + return point; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_cavity.h b/source/source_hamilt/module_sccs/sccs_cavity.h new file mode 100644 index 00000000000..af3e4eeb9d9 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_cavity.h @@ -0,0 +1,27 @@ +#ifndef SCCS_CAVITY_H +#define SCCS_CAVITY_H + +namespace ModuleSccs +{ +// Densities are in electrons/Bohr^3. The cavity uses the electronic density. +struct CavityParameters +{ + double density_min = 0.0; + double density_max = 0.0; + double epsilon_bulk = 0.0; +}; + +struct CavityPoint +{ + double solute = 0.0; + double dsolute_drho = 0.0; + double epsilon = 0.0; + double depsilon_drho = 0.0; +}; + +// Parameters are validated by the INPUT reader. Negative Fourier ringing +// uses the bulk limit. +CavityPoint evaluate_cavity(double density, const CavityParameters& parameters); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_functional.cpp b/source/source_hamilt/module_sccs/sccs_functional.cpp new file mode 100644 index 00000000000..7547d7dae7b --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_functional.cpp @@ -0,0 +1,68 @@ +#include "sccs_functional.h" +#include "sccs_parameters.h" +#include "sccs_response.h" +#include "sccs_pw_coulomb.h" + +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include + +namespace ModuleSccs +{ +void evaluate_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result) +{ + const std::size_t size = charge.size(); + const double dv = basis.omega / basis.nxyz; + PeriodicCoulombOperator coulomb(basis, tpiba); + std::vector vacuum; + coulomb.apply_potential(charge, vacuum); + FunctionalResult candidate; + candidate.reaction_potential.resize(size); + candidate.electron_potential.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const double reaction = response.polarization.potential[i] - vacuum[i]; + candidate.reaction_potential[i] = reaction; + candidate.reaction_energy += 0.5 * charge[i] * reaction * dv; + candidate.electron_potential[i] = -reaction + response.cavity_potential[i]; + candidate.volume += response.solute[i] * dv; + } + Parallel_Reduce::reduce_pool(candidate.reaction_energy); + Parallel_Reduce::reduce_pool(candidate.volume); + candidate.volume_energy = config.pressure * candidate.volume; + + // Spectral surface derivative: dE/dn = (p - gamma div(grad s/|grad s|)) s'. + std::vector> solute_g(basis.npw); + basis.real2recip(response.solute.data(), solute_g.data()); + std::vector> unit_gradient(size); + XC_Functional::grad_rho(solute_g.data(), unit_gradient.data(), &basis, tpiba); + const double eta = config.surface_regularization; + for (std::size_t i = 0; i < size; ++i) + { + const double norm_squared = unit_gradient[i] * unit_gradient[i] + eta * eta; + const double norm = std::sqrt(norm_squared); + candidate.surface += (norm - eta) * dv; + unit_gradient[i] /= norm; + } + Parallel_Reduce::reduce_pool(candidate.surface); + candidate.surface_energy = config.surface_tension * candidate.surface; + std::vector divergence(size); + XC_Functional::grad_dot(unit_gradient.data(), divergence.data(), &basis, tpiba); + for (std::size_t i = 0; i < size; ++i) + { + const double nonel = (config.pressure - config.surface_tension * divergence[i]) + * response.dsolute_drho[i]; + candidate.electron_potential[i] += nonel; + } + result = std::move(candidate); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_functional.h b/source/source_hamilt/module_sccs/sccs_functional.h new file mode 100644 index 00000000000..ed880ef939d --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_functional.h @@ -0,0 +1,36 @@ +#ifndef SCCS_FUNCTIONAL_H +#define SCCS_FUNCTIONAL_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct SccsConfig; +struct SccsResponse; + +struct FunctionalResult +{ + double reaction_energy = 0.0; // Hartree + double surface_energy = 0.0; + double volume_energy = 0.0; + double surface = 0.0; // Bohr^2 + double volume = 0.0; // Bohr^3 + std::vector reaction_potential; // dielectric minus vacuum, Hartree + std::vector electron_potential; // derivative of all three energy terms +}; + +// Collective over the PW pool; a non-finite energy stops the run. +void evaluate_functional(const std::vector& charge, + const SccsResponse& response, + const SccsConfig& config, + const ModulePW::PW_Basis& basis, + double tpiba, + FunctionalResult& result); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp new file mode 100644 index 00000000000..a1015081797 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_charge.cpp @@ -0,0 +1,33 @@ +#include "sccs_ionic_charge.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/cell_tools.h" + +#include +#include + +namespace ModuleSccs +{ +void gaussian_ionic_density(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density) +{ + std::vector> charge(basis.npw, 0.0); + const double tpiba2 = tpiba * tpiba; + for (const unitcell::AtomData& atom : atoms) + { + for (int ig = 0; ig < basis.npw; ++ig) + { + const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; + const double phase = tpiba * (basis.gcar[ig] * atom.position); + const std::complex phase_factor = ModuleBase::NEG_IMAG_UNIT * phase; + charge[ig] += atom.valence_charge / basis.omega * std::exp(exponent) * std::exp(phase_factor); + } + } + density.resize(basis.nrxx); + basis.recip2real(charge.data(), density.data()); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ionic_charge.h b/source/source_hamilt/module_sccs/sccs_ionic_charge.h new file mode 100644 index 00000000000..38cec3a8c9b --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_charge.h @@ -0,0 +1,29 @@ +#ifndef SCCS_IONIC_CHARGE_H +#define SCCS_IONIC_CHARGE_H + +#include + +namespace ModulePW +{ +class PW_Basis; +} +namespace unitcell +{ +struct AtomData; +} + +namespace ModuleSccs +{ +// rho ~ exp(-r^2/spread^2), in Bohr; preserves the original SCCS ionic spread. +const double gaussian_ion_spread = 0.5; + +// AtomData is extracted by source_cell; positions are Cartesian Bohr. +// Collective over the PW pool. +void gaussian_ionic_density(const std::vector& atoms, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector& density); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/sccs_ionic_force.cpp new file mode 100644 index 00000000000..0d6ee4f6726 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_force.cpp @@ -0,0 +1,45 @@ +#include "sccs_ionic_force.h" + +#include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/cell_tools.h" + +#include +#include + +namespace ModuleSccs +{ +void gaussian_ionic_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces) +{ + std::vector> potential_g(basis.npw); + basis.real2recip(reaction_potential.data(), potential_g.data()); + forces.assign(atoms.size(), ModuleBase::Vector3(0.0, 0.0, 0.0)); + const double tpiba2 = tpiba * tpiba; + // F = -integral(v_reaction d rho_ion/d R), using ABACUS's normalized FFT. + // The cell volume cancels the Z/omega in the Gaussian charge coefficients. + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + const unitcell::AtomData& atom = atoms[ia]; + for (int ig = 0; ig < basis.npw; ++ig) + { + const double exponent = -0.25 * spread * spread * tpiba2 * basis.gg[ig]; + const double gaussian = std::exp(exponent); + const double phase = tpiba * (basis.gcar[ig] * atom.position); + const std::complex phase_argument = ModuleBase::NEG_IMAG_UNIT * phase; + const std::complex phase_factor = std::exp(phase_argument); + const std::complex weighted = std::conj(potential_g[ig]) * phase_factor; + const double factor = -atom.valence_charge * tpiba * gaussian * weighted.imag(); + forces[ia] += basis.gcar[ig] * factor; + } + Parallel_Reduce::reduce_pool(forces[ia].x); + Parallel_Reduce::reduce_pool(forces[ia].y); + Parallel_Reduce::reduce_pool(forces[ia].z); + } +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_ionic_force.h b/source/source_hamilt/module_sccs/sccs_ionic_force.h new file mode 100644 index 00000000000..150781b03d8 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_ionic_force.h @@ -0,0 +1,30 @@ +#ifndef SCCS_IONIC_FORCE_H +#define SCCS_IONIC_FORCE_H + +#include + +namespace ModuleBase +{ +template class Vector3; +} +namespace ModulePW +{ +class PW_Basis; +} +namespace unitcell +{ +struct AtomData; +} +namespace ModuleSccs +{ +// Explicit Gaussian-ion derivative of the reaction energy at fixed electron +// density. reaction_potential is in Ha; force is in Ha/Bohr, atom-major. +// Collective over the PW pool. +void gaussian_ionic_force(const std::vector& atoms, + const std::vector& reaction_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + double spread, + std::vector>& forces); +} // namespace ModuleSccs +#endif diff --git a/source/source_hamilt/module_sccs/sccs_parameters.cpp b/source/source_hamilt/module_sccs/sccs_parameters.cpp new file mode 100644 index 00000000000..0a8712174f6 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_parameters.cpp @@ -0,0 +1,78 @@ +#include "sccs_parameters.h" + +#include "source_base/tool_quit.h" + + +namespace ModuleSccs +{ +namespace +{ +// CODATA 2018 SI values retain the original SCCS preset conversion. +// ModuleBase's older SI values differ, so using them would change these presets. +const double hartree_joule = 4.3597447222071e-18; +const double bohr_metre = 5.29177210903e-11; +} + +double dyn_per_cm_to_hartree_per_bohr2(double value) +{ + const double joule_per_square_metre = value * 1e-3; + return joule_per_square_metre * bohr_metre * bohr_metre / hartree_joule; +} + +double gpa_to_hartree_per_bohr3(double value) +{ + const double pascal = value * 1e9; + return pascal * bohr_metre * bohr_metre * bohr_metre / hartree_joule; +} + +Preset parse_preset(const std::string& name) +{ + if (name == "custom") { return Preset::Custom; } + if (name == "vacuum") { return Preset::Vacuum; } + if (name == "water-neutral") { return Preset::WaterNeutral; } + if (name == "water-cation") { return Preset::WaterCation; } + if (name == "water-anion") { return Preset::WaterAnion; } + const std::string message = "Unknown SCCS preset: " + name; + ModuleBase::WARNING_QUIT("ModuleSccs::parse_preset", message); +} + +SccsConfig make_sccs_config(Preset preset) +{ + if (preset == Preset::Custom) + { + ModuleBase::WARNING_QUIT("ModuleSccs::make_sccs_config", + "SCCS custom configurations must be supplied explicitly"); + } + SccsConfig config; + config.cavity.epsilon_bulk = 78.3; + config.surface_regularization = 1e-8; + if (preset == Preset::Vacuum || preset == Preset::WaterNeutral) + { + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + if (preset == Preset::Vacuum) + { + config.cavity.epsilon_bulk = 1.0; + } + else + { + config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(47.9); + config.pressure = gpa_to_hartree_per_bohr3(-0.36); + } + } + else if (preset == Preset::WaterCation) + { + config.cavity.density_min = 2e-4; + config.cavity.density_max = 3.5e-3; + config.surface_tension = dyn_per_cm_to_hartree_per_bohr2(5.0); + config.pressure = gpa_to_hartree_per_bohr3(0.125); + } + else + { + config.cavity.density_min = 2.4e-3; + config.cavity.density_max = 1.55e-2; + config.pressure = gpa_to_hartree_per_bohr3(0.45); + } + return config; +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_parameters.h b/source/source_hamilt/module_sccs/sccs_parameters.h new file mode 100644 index 00000000000..707c5c0c386 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_parameters.h @@ -0,0 +1,36 @@ +#ifndef SCCS_PARAMETERS_H +#define SCCS_PARAMETERS_H + +#include "sccs_cavity.h" + +#include + +namespace ModuleSccs +{ +enum class Preset { Custom, Vacuum, WaterNeutral, WaterCation, WaterAnion }; + +struct SccsConfig +{ + CavityParameters cavity; + double surface_tension = 0.0; // Hartree/Bohr^2 + double pressure = 0.0; // Hartree/Bohr^3 + double surface_regularization = 0.0; +}; + +struct PolarizationSolverParameters +{ + int max_iterations = 200; + double tolerance_rms = 1e-10; + double tolerance_max = 1e-8; +}; + +// Exact, lower-case names as accepted by the INPUT reader; an unknown name +// stops the run. +Preset parse_preset(const std::string& name); +// Published presets only; custom configurations are supplied by callers. +SccsConfig make_sccs_config(Preset preset); +double dyn_per_cm_to_hartree_per_bohr2(double value); +double gpa_to_hartree_per_bohr3(double value); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp new file mode 100644 index 00000000000..62fbebabd20 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.cpp @@ -0,0 +1,32 @@ +#include "sccs_pw_coulomb.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" + +namespace ModuleSccs +{ +PeriodicCoulombOperator::PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba) + : basis_(basis), tpiba_(tpiba) +{ +} + +void PeriodicCoulombOperator::apply_potential(const std::vector& charge, std::vector& potential) +{ + reciprocal_work_.resize(basis_.npw); + basis_.real2recip(charge.data(), reciprocal_work_.data()); + const double tpiba2 = tpiba_ * tpiba_; + for (int ig = 0; ig < basis_.npw; ++ig) + { + if (ig == basis_.ig_gge0 || basis_.gg[ig] == 0.0) + { + reciprocal_work_[ig] = 0.0; + } + else + { + reciprocal_work_[ig] = ModuleBase::FOUR_PI * reciprocal_work_[ig] / (tpiba2 * basis_.gg[ig]); + } + } + potential.resize(basis_.nrxx); + basis_.recip2real(reciprocal_work_.data(), potential.data()); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_pw_coulomb.h b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h new file mode 100644 index 00000000000..41245a45608 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_pw_coulomb.h @@ -0,0 +1,31 @@ +#ifndef SCCS_PW_COULOMB_H +#define SCCS_PW_COULOMB_H + +#include +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +// Potential in Hartree for a signed charge density in e/Bohr^3. Calls are +// collective over the ABACUS pool used to initialize basis. G=0 is removed: +// this is the periodic Coulomb inverse with a uniform background. +class PeriodicCoulombOperator +{ +public: + PeriodicCoulombOperator(const ModulePW::PW_Basis& basis, double tpiba); + void apply_potential(const std::vector& charge, std::vector& potential); + +private: + const ModulePW::PW_Basis& basis_; + const double tpiba_; + // Scratch only, overwritten on each call; PW_Basis FFTs are also non-reentrant. + std::vector> reciprocal_work_; +}; +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/sccs_response.cpp b/source/source_hamilt/module_sccs/sccs_response.cpp new file mode 100644 index 00000000000..ea0a81c414f --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.cpp @@ -0,0 +1,316 @@ +#include "sccs_response.h" +#include "sccs_parameters.h" +#include "sccs_pw_coulomb.h" + +#include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_quit.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_xc/xc_functional.h" + +#include +#include +#include +#include +#include + +namespace ModuleSccs +{ +namespace +{ +void prepare_cavity(const std::vector& density, + const CavityParameters& cavity, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& response, + std::vector& coefficient) +{ + const std::size_t size = density.size(); + response.solute.resize(size); + response.dsolute_drho.resize(size); + response.epsilon.resize(size); + response.depsilon_drho.resize(size); + response.grad_log_epsilon.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const CavityPoint point = evaluate_cavity(density[i], cavity); + response.solute[i] = point.solute; + response.dsolute_drho[i] = point.dsolute_drho; + response.epsilon[i] = point.epsilon; + response.depsilon_drho[i] = point.depsilon_drho; + } + + std::vector> density_g(basis.npw); + basis.real2recip(density.data(), density_g.data()); + std::vector> gradient(size); + std::vector laplacian(size); + // These XC helpers use G in units of tpiba and ABACUS's normalized FFT. + XC_Functional::grad_rho(density_g.data(), gradient.data(), &basis, tpiba); + XC_Functional::laplacian_rho(density_g.data(), laplacian.data(), &basis, tpiba); + const double density_ratio = cavity.density_max / cavity.density_min; + const double width = std::log(density_ratio); + const double log_bulk = std::log(cavity.epsilon_bulk); + coefficient.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + double second_log = 0.0; + if (density[i] > cavity.density_min && density[i] < cavity.density_max) + { + const double local_ratio = cavity.density_max / density[i]; + const double x = std::log(local_ratio) / width; + const double angle = ModuleBase::TWO_PI * x; + second_log = log_bulk * (1.0 - std::cos(angle) + ModuleBase::TWO_PI * std::sin(angle) / width) + / (width * density[i] * density[i]); + } + const double first_log = response.depsilon_drho[i] / response.epsilon[i]; + const ModuleBase::Vector3& local_gradient = gradient[i]; + const double gradient_square = local_gradient * local_gradient; + response.grad_log_epsilon[i] = local_gradient * first_log; + const double lap_log = first_log * laplacian[i] + second_log * gradient_square; + coefficient[i] = response.epsilon[i] * (0.5 * lap_log + 0.25 * first_log * first_log * gradient_square) + / ModuleBase::FOUR_PI; + } +} + +// Pool-reduced norms, so every rank takes the same convergence decision. A +// non-finite residual fails convergence and stops at the next CG step. +void residual_norms(const std::vector& values, + const ModulePW::PW_Basis& basis, + double& rms, + double& maximum) +{ + double square = 0.0; + maximum = 0.0; + for (double value : values) + { + square += value * value; + const double magnitude = std::abs(value); + maximum = std::max(maximum, magnitude); + } + Parallel_Reduce::reduce_pool(square); + Parallel_Reduce::reduce_max_pool(basis.poolnproc, maximum); + const double mean_square = square / basis.nxyz; + rms = std::sqrt(mean_square); +} + +bool converged(const PolarizationResult& result, const PolarizationSolverParameters& solver) +{ + return result.residual_rms <= solver.tolerance_rms && result.residual_max <= solver.tolerance_max; +} + +double grid_dot(const std::vector& left, + const std::vector& right, + const ModulePW::PW_Basis& basis) +{ + double value = 0.0; + for (std::size_t i = 0; i < left.size(); ++i) + { + value += left[i] * right[i]; + } + Parallel_Reduce::reduce_pool(value); + return value * basis.omega / basis.nxyz; +} + +// P r = eps^-1/2 G eps^-1/2 r; only FFT scratch survives an application. +void apply_preconditioner(const std::vector& rhs, + const std::vector& invsqrt, + PeriodicCoulombOperator& coulomb, + std::vector& weighted, + std::vector& value) +{ + const std::size_t size = rhs.size(); + weighted.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + weighted[i] = rhs[i] * invsqrt[i]; + } + coulomb.apply_potential(weighted, value); + for (std::size_t i = 0; i < size; ++i) + { + value[i] *= invsqrt[i]; + } +} + +// The preconditioned sqrt-CG system of one response. +struct SqrtCgOperator +{ + const std::vector& coefficient; + const std::vector& invsqrt; + const ModulePW::PW_Basis& basis; + PeriodicCoulombOperator& coulomb; +}; + +// Potential and charge residual of the sqrt-CG, with preconditioner scratch. +struct SqrtCgState +{ + std::vector potential; + std::vector residual; + std::vector weighted; + std::vector preconditioned; +}; + +// Replace the cold start by the old fixed point when that lowers the RMS +// charge residual. +void try_warm_start(const SqrtCgOperator& system, + const std::vector& charge, + const std::vector& initial_potential, + SqrtCgState& state, + PolarizationResult& polarization) +{ + const std::size_t size = charge.size(); + const std::vector& coefficient = system.coefficient; + std::vector& z = state.preconditioned; + std::vector guess_residual(size); + for (std::size_t i = 0; i < size; ++i) + { + guess_residual[i] = charge[i] - coefficient[i] * initial_potential[i]; + } + apply_preconditioner(guess_residual, system.invsqrt, system.coulomb, state.weighted, z); + for (std::size_t i = 0; i < size; ++i) + { + guess_residual[i] = coefficient[i] * (initial_potential[i] - z[i]); + } + double guess_rms = 0.0; + double guess_max = 0.0; + residual_norms(guess_residual, system.basis, guess_rms, guess_max); + if (guess_rms < polarization.residual_rms) + { + state.potential.swap(z); + state.residual.swap(guess_residual); + polarization.warm_started = true; + polarization.residual_rms = guess_rms; + polarization.residual_max = guess_max; + } +} + +// Preconditioned CG until both residual tolerances hold; a breakdown or a +// missed tolerance stops the run with WARNING_QUIT. +void iterate_sqrt_cg(const SqrtCgOperator& system, + const PolarizationSolverParameters& solver, + SqrtCgState& state, + PolarizationResult& polarization) +{ + const ModulePW::PW_Basis& basis = system.basis; + const std::vector& coefficient = system.coefficient; + std::vector& potential = state.potential; + std::vector& residual = state.residual; + std::vector& z = state.preconditioned; + const std::size_t size = residual.size(); + std::vector direction(size, 0.0); + std::vector image(size, 0.0); + double old_rz = 0.0; + for (int iteration = 1; !converged(polarization, solver) && iteration <= solver.max_iterations; ++iteration) + { + apply_preconditioner(residual, system.invsqrt, system.coulomb, state.weighted, z); + const double rz = grid_dot(residual, z, basis); + if (!std::isfinite(rz) || std::abs(rz) < 1e-30) + { + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", + "SCCS sqrt-CG has a null or nonfinite preconditioned residual"); + } + const double beta = std::abs(old_rz) > 1e-30 ? rz / old_rz : 0.0; + old_rz = rz; + for (std::size_t i = 0; i < size; ++i) + { + direction[i] = z[i] + beta * direction[i]; + image[i] = coefficient[i] * z[i] + residual[i] + beta * image[i]; + } + const double curvature = grid_dot(direction, image, basis); + if (!std::isfinite(curvature) || curvature == 0.0) + { + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", "SCCS sqrt-CG has invalid curvature"); + } + const double alpha = rz / curvature; + for (std::size_t i = 0; i < size; ++i) + { + potential[i] += alpha * direction[i]; + residual[i] -= alpha * image[i]; + } + polarization.iterations = iteration; + residual_norms(residual, basis, polarization.residual_rms, polarization.residual_max); + } + if (!converged(polarization, solver)) + { + std::ostringstream message; + message << "SCCS sqrt-CG did not reach both residual tolerances within sccs_maxiter = " + << solver.max_iterations << " iterations (RMS " << polarization.residual_rms << ", MAX " + << polarization.residual_max << ")"; + ModuleBase::WARNING_QUIT("ModuleSccs::solve_sccs_response", message.str()); + } +} + +// Keep the unshifted solution for warm starts, then store the zero-mean +// potential, its gradient and the cavity potential. +void finish_response(const ModulePW::PW_Basis& basis, + double tpiba, + std::vector& potential, + SccsResponse& response) +{ + PolarizationResult& polarization = response.polarization; + const std::size_t size = potential.size(); + response.restart_potential = potential; + double mean = 0.0; + for (double value : potential) + { + mean += value; + } + Parallel_Reduce::reduce_pool(mean); + mean /= basis.nxyz; + for (double& value : potential) + { + value -= mean; + } + std::vector> potential_g(basis.npw); + basis.real2recip(potential.data(), potential_g.data()); + polarization.gradient.resize(size); + XC_Functional::grad_rho(potential_g.data(), polarization.gradient.data(), &basis, tpiba); + polarization.potential.swap(potential); + response.cavity_potential.resize(size); + for (std::size_t i = 0; i < size; ++i) + { + const ModuleBase::Vector3& gradient = polarization.gradient[i]; + const double gradient_square = gradient * gradient; + response.cavity_potential[i] = -response.depsilon_drho[i] * gradient_square / (8.0 * ModuleBase::PI); + } +} +} // namespace + +void solve_sccs_response(const std::vector& density, + const std::vector& charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result) +{ + // Determine start mode collectively, including ranks with no real-space points. + double initial_count = initial_potential.size(); + Parallel_Reduce::reduce_pool(initial_count); + const bool warm_start = initial_count != 0.0; + + SccsResponse candidate; + std::vector coefficient; + prepare_cavity(density, cavity, basis, tpiba, candidate, coefficient); + const std::size_t size = density.size(); + std::vector invsqrt(size); + for (std::size_t i = 0; i < size; ++i) + { + invsqrt[i] = 1.0 / std::sqrt(candidate.epsilon[i]); + } + PeriodicCoulombOperator coulomb(basis, tpiba); + const SqrtCgOperator system = {coefficient, invsqrt, basis, coulomb}; + SqrtCgState state; + state.residual = charge; + state.potential.assign(size, 0.0); + PolarizationResult& polarization = candidate.polarization; + residual_norms(state.residual, basis, polarization.residual_rms, polarization.residual_max); + if (!converged(polarization, solver) && warm_start) + { + try_warm_start(system, charge, initial_potential, state, polarization); + } + iterate_sqrt_cg(system, solver, state, polarization); + finish_response(basis, tpiba, state.potential, candidate); + result = std::move(candidate); +} +} // namespace ModuleSccs diff --git a/source/source_hamilt/module_sccs/sccs_response.h b/source/source_hamilt/module_sccs/sccs_response.h new file mode 100644 index 00000000000..bc81ad622b9 --- /dev/null +++ b/source/source_hamilt/module_sccs/sccs_response.h @@ -0,0 +1,55 @@ +#ifndef SCCS_RESPONSE_H +#define SCCS_RESPONSE_H + +#include "sccs_cavity.h" +#include "source_base/vector3.h" + +#include + +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleSccs +{ +struct PolarizationSolverParameters; + +struct PolarizationResult +{ + int iterations = 0; + double residual_rms = 0.0; + double residual_max = 0.0; + bool warm_started = false; + std::vector potential; // Hartree, shifted to periodic zero mean + std::vector> gradient; // Hartree/Bohr +}; + +struct SccsResponse +{ + std::vector solute; + std::vector dsolute_drho; + std::vector epsilon; + std::vector depsilon_drho; + std::vector> grad_log_epsilon; + PolarizationResult polarization; + std::vector cavity_potential; // -epsilon' |grad v|^2/(8 pi), Hartree + std::vector restart_potential; // unshifted solution for explicit warm starts +}; + +// Original chain-derivative sqrt-CG periodic response; charge is ions minus electrons. +// All ranks in the PW pool must call together with the same configuration. +// initial_potential is empty on all ranks for a cold start, otherwise a local grid. +// A solver breakdown or missed tolerance stops the run with WARNING_QUIT. +// No state or INPUT globals are read or retained. +void solve_sccs_response(const std::vector& cavity_density, + const std::vector& solute_charge, + const CavityParameters& cavity, + const PolarizationSolverParameters& solver, + const std::vector& initial_potential, + const ModulePW::PW_Basis& basis, + double tpiba, + SccsResponse& result); +} // namespace ModuleSccs + +#endif diff --git a/source/source_hamilt/module_sccs/test/CMakeLists.txt b/source/source_hamilt/module_sccs/test/CMakeLists.txt new file mode 100644 index 00000000000..722c01ebd01 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/CMakeLists.txt @@ -0,0 +1,27 @@ +AddTest(TARGET MODULE_HAMILT_sccs_cavity + SOURCES test_sccs_cavity.cpp ../sccs_cavity.cpp) +AddTest(TARGET MODULE_HAMILT_sccs_parameters + LIBS base device + SOURCES test_sccs_parameters.cpp ../sccs_parameters.cpp ../sccs_cavity.cpp) + +AddTest(TARGET MODULE_HAMILT_sccs_poisson + LIBS parameter planewave device base + SOURCES test_sccs_pw_coulomb.cpp test_sccs_ionic_charge.cpp test_sccs_ionic_force.cpp test_sccs_response.cpp test_sccs_functional.cpp + ../../../source_estate/test/test_pot_sccs.cpp + ../../../source_estate/module_pot/pot_sccs.cpp + ../../../source_cell/cell_tools.cpp ../../../source_cell/atom_spec.cpp + ../../../source_cell/atom_pseudo.cpp ../../../source_cell/pseudo.cpp + ../sccs_functional.cpp + ../sccs_parameters.cpp + ../sccs_ionic_charge.cpp ../sccs_ionic_force.cpp ../sccs_pw_coulomb.cpp ../sccs_response.cpp ../sccs_cavity.cpp + ../../module_xc/xc_grad_utils.cpp) +set_tests_properties(MODULE_HAMILT_sccs_poisson PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60) + +if(ENABLE_MPI) + add_test(NAME MODULE_HAMILT_sccs_poisson_mpi + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS}) + set_tests_properties(MODULE_HAMILT_sccs_poisson_mpi PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" TIMEOUT 60 PROCESSORS 2) +endif() diff --git a/source/source_hamilt/module_sccs/test/sccs_test.h b/source/source_hamilt/module_sccs/test/sccs_test.h new file mode 100644 index 00000000000..6a1570cca4f --- /dev/null +++ b/source/source_hamilt/module_sccs/test/sccs_test.h @@ -0,0 +1,60 @@ +#ifndef SCCS_TEST_H +#define SCCS_TEST_H + +#include "source_base/constants.h" +#include "source_base/matrix3.h" +#include "source_base/parallel_comm.h" +#include "source_basis/module_pw/pw_basis.h" +#include + +#include + +namespace SccsTest +{ +extern int pool_size; +extern int pool_rank; + +class PwTest : public testing::Test +{ +protected: + PwTest() : basis("cpu", "double") {} + + void SetUp() override + { +#ifdef __MPI + basis.initmpi(pool_size, pool_rank, POOL_WORLD); +#endif + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + basis.initgrids(length, lattice, 20.0); + basis.initparameters(false, 20.0, 1, false); + basis.setuptransform(); + basis.collect_local_pw(); + tpiba = ModuleBase::TWO_PI / length; + } + + // A cosine along x or y, in the PW local real-grid ordering. + std::vector cosine_mode(int direction) const + { + std::vector values(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int ix = ir / (basis.ny * basis.nplane); + const int iy = ir / basis.nplane - ix * basis.ny; + const int coordinate = direction == 0 ? ix : iy; + const int count = direction == 0 ? basis.nx : basis.ny; + const double angle = ModuleBase::TWO_PI * coordinate / count; + values[ir] = std::cos(angle); + } + return values; + } + + const double length = 10.0; + double tpiba = 0.0; + ModulePW::PW_Basis basis; + std::string error; +}; +} // namespace SccsTest + +#endif diff --git a/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp new file mode 100644 index 00000000000..c4d8fc15c4d --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_cavity.cpp @@ -0,0 +1,110 @@ +#include "../sccs_cavity.h" +#include + +#include + +namespace +{ +class SccsCavityTest : public ::testing::Test +{ +protected: + ModuleSccs::CavityParameters parameters; + + void SetUp() override + { + parameters.density_min = 1e-4; + parameters.density_max = 5e-3; + parameters.epsilon_bulk = 78.3; + } +}; + +// Bulk and solute limits, continuity at both thresholds, and no dielectric +// response in vacuum. +TEST_F(SccsCavityTest, LimitsThresholdsAndVacuum) +{ + const double bulk_densities[] = {-1e-6, 0.0, 1e-4}; + for (double density : bulk_densities) + { + ModuleSccs::CavityPoint point; + point = ModuleSccs::evaluate_cavity(density, parameters); + EXPECT_DOUBLE_EQ(point.solute, 0.0); + EXPECT_DOUBLE_EQ(point.epsilon, 78.3); + EXPECT_DOUBLE_EQ(point.dsolute_drho, 0.0); + EXPECT_DOUBLE_EQ(point.depsilon_drho, 0.0); + } + const double solute_densities[] = {5e-3, 0.1}; + for (double density : solute_densities) + { + ModuleSccs::CavityPoint point; + point = ModuleSccs::evaluate_cavity(density, parameters); + EXPECT_DOUBLE_EQ(point.solute, 1.0); + EXPECT_DOUBLE_EQ(point.epsilon, 1.0); + EXPECT_DOUBLE_EQ(point.dsolute_drho, 0.0); + EXPECT_DOUBLE_EQ(point.depsilon_drho, 0.0); + } + const double near_min = parameters.density_min * (1.0 + 1e-8); + const double near_max = parameters.density_max * (1.0 - 1e-8); + ModuleSccs::CavityPoint lower; + ModuleSccs::CavityPoint upper; + lower = ModuleSccs::evaluate_cavity(near_min, parameters); + upper = ModuleSccs::evaluate_cavity(near_max, parameters); + EXPECT_NEAR(lower.solute, 0.0, 1e-14); + EXPECT_NEAR(upper.solute, 1.0, 1e-14); + EXPECT_NEAR(lower.epsilon, 78.3, 1e-12); + EXPECT_NEAR(upper.epsilon, 1.0, 1e-12); + EXPECT_NEAR(lower.dsolute_drho, 0.0, 1e-10); + EXPECT_NEAR(upper.dsolute_drho, 0.0, 1e-10); + EXPECT_NEAR(lower.depsilon_drho, 0.0, 1e-8); + EXPECT_NEAR(upper.depsilon_drho, 0.0, 1e-8); + parameters.epsilon_bulk = 1.0; + const ModuleSccs::CavityPoint vacuum = ModuleSccs::evaluate_cavity(1e-3, parameters); + EXPECT_GT(vacuum.solute, 0.0); + EXPECT_LT(vacuum.solute, 1.0); + EXPECT_DOUBLE_EQ(vacuum.epsilon, 1.0); + EXPECT_DOUBLE_EQ(vacuum.depsilon_drho, 0.0); +} + +// Analytic values at the logarithmic midpoint and finite differences across +// the transition. +TEST_F(SccsCavityTest, MidpointAndDerivativesAcrossTransition) +{ + const double product = parameters.density_min * parameters.density_max; + const double midpoint_density = std::sqrt(product); + const double ratio = parameters.density_max / parameters.density_min; + const double width = std::log(ratio); + const double epsilon = std::sqrt(parameters.epsilon_bulk); + const double log_epsilon = std::log(parameters.epsilon_bulk); + const double dsolute = 2.0 / (width * midpoint_density); + const double depsilon = -epsilon * log_epsilon * dsolute; + const ModuleSccs::CavityPoint midpoint = ModuleSccs::evaluate_cavity(midpoint_density, parameters); + EXPECT_NEAR(midpoint.solute, 0.5, 1e-14); + EXPECT_NEAR(midpoint.epsilon, epsilon, 1e-13); + EXPECT_NEAR(midpoint.dsolute_drho, dsolute, 1e-10); + EXPECT_NEAR(midpoint.depsilon_drho, depsilon, 1e-9); + + const double fractions[] = {0.2, 0.5, 0.8}; + for (double fraction : fractions) + { + const double scale = std::pow(ratio, fraction); + const double density = parameters.density_min * scale; + const double step = density * 1e-5; + const double lower_density = density - step; + const double upper_density = density + step; + ModuleSccs::CavityPoint point; + ModuleSccs::CavityPoint lower; + ModuleSccs::CavityPoint upper; + point = ModuleSccs::evaluate_cavity(density, parameters); + lower = ModuleSccs::evaluate_cavity(lower_density, parameters); + upper = ModuleSccs::evaluate_cavity(upper_density, parameters); + const double solute_fd = (upper.solute - lower.solute) / (2.0 * step); + const double epsilon_fd = (upper.epsilon - lower.epsilon) / (2.0 * step); + const double solute_ratio = solute_fd / point.dsolute_drho; + const double epsilon_ratio = epsilon_fd / point.depsilon_drho; + EXPECT_NEAR(solute_ratio, 1.0, 1e-8); + EXPECT_NEAR(epsilon_ratio, 1.0, 1e-8); + EXPECT_GT(point.dsolute_drho, 0.0); + EXPECT_LT(point.depsilon_drho, 0.0); + } +} + +} // namespace diff --git a/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp new file mode 100644 index 00000000000..c25cde778d4 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_functional.cpp @@ -0,0 +1,84 @@ +#include "sccs_test.h" +#include "../sccs_functional.h" +#include "../sccs_parameters.h" +#include "../sccs_response.h" +#include "source_base/parallel_reduce.h" + +using SccsFunctionalTest = SccsTest::PwTest; + +TEST_F(SccsFunctionalTest, UniformDielectricEnergyAndElectronDerivative) +{ + const std::vector mode = cosine_mode(0); + std::vector charge = mode; + for (double& value : charge) + { + value *= 1e-3; + } + ModuleSccs::SccsConfig config; + config.cavity.density_min = 1e-4; + config.cavity.density_max = 5e-3; + config.cavity.epsilon_bulk = 5.0; + config.surface_regularization = 1e-8; + ModuleSccs::PolarizationSolverParameters solver; + const std::vector density(basis.nrxx, 0.0); + const std::vector cold; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold, + basis, tpiba, response); + ModuleSccs::FunctionalResult result; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, result); + const double kernel = ModuleBase::FOUR_PI / (tpiba * tpiba); + const double expected = 0.25 * basis.omega * 1e-6 * kernel * (0.2 - 1.0); + EXPECT_NEAR(result.reaction_energy, expected, 1e-13); + EXPECT_DOUBLE_EQ(result.surface_energy, 0.0); + EXPECT_DOUBLE_EQ(result.volume_energy, 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double potential = 0.8 * kernel * charge[ir]; + EXPECT_NEAR(result.electron_potential[ir], potential, 1e-13); + } +} + +TEST_F(SccsFunctionalTest, NonElectrostaticDiscreteEnergyDerivative) +{ + ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + config.surface_tension = 0.01; + config.pressure = -0.001; + config.surface_regularization = 0.02; + const std::vector mode = cosine_mode(0); + ModuleSccs::SccsResponse response; + response.solute.resize(basis.nrxx); + response.dsolute_drho.assign(basis.nrxx, 1.0); + response.polarization.potential.assign(basis.nrxx, 0.0); + response.cavity_potential.assign(basis.nrxx, 0.0); + const std::vector charge(basis.nrxx, 0.0); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + response.solute[ir] = 0.5 + 0.2 * mode[ir]; + } + ModuleSccs::FunctionalResult center; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, center); + // Avoid cancellation of grid-integrated energies at excessively small steps. + const double step = 1e-5; + double predicted = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + predicted += center.electron_potential[ir] * mode[ir]; + response.solute[ir] += step * mode[ir]; + } + Parallel_Reduce::reduce_pool(predicted); + predicted *= basis.omega / basis.nxyz; + ModuleSccs::FunctionalResult plus; + ModuleSccs::FunctionalResult minus; + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, plus); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + response.solute[ir] -= 2.0 * step * mode[ir]; + } + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, minus); + const double fd = (plus.surface_energy + plus.volume_energy - minus.surface_energy - minus.volume_energy) + / (2.0 * step); + EXPECT_NEAR(fd, predicted, 1e-8); + const double expected_volume = 0.5 * basis.omega; + EXPECT_NEAR(center.volume, expected_volume, 1e-10); +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp new file mode 100644 index 00000000000..27c2b795746 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_ionic_charge.cpp @@ -0,0 +1,53 @@ +#include "sccs_test.h" +#include "../sccs_ionic_charge.h" + +#include "source_base/parallel_reduce.h" +#include "source_cell/cell_tools.h" + +using SccsIonicChargeTest = SccsTest::PwTest; + +TEST_F(SccsIonicChargeTest, ValenceNormalizationAndFourierPhase) +{ + std::vector atoms(2); + atoms[0].position.x = length / 4.0; + atoms[0].valence_charge = 8.0; + atoms[1].valence_charge = 1.0; + const double spread = ModuleSccs::gaussian_ion_spread; + std::vector density; + ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, density); + double integral = 0.0; + for (double value : density) + { + integral += value; + } + Parallel_Reduce::reduce_pool(integral); + integral *= basis.omega / basis.nxyz; + EXPECT_NEAR(integral, 9.0, 1e-12); + + std::vector> charge_g(basis.npw); + basis.real2recip(density.data(), charge_g.data()); + const double exponent = -spread * spread * tpiba * tpiba / 4.0; + const double amplitude = 8.0 / basis.omega * std::exp(exponent); + const double second_amplitude = 1.0 / basis.omega * std::exp(exponent); + double selected = 0.0; + for (int ig = 0; ig < basis.npw; ++ig) + { + const ModuleBase::Vector3& g = basis.gdirect[ig]; + if (g.x == 1.0 && g.y == 0.0 && g.z == 0.0) + { + selected = 1.0; + EXPECT_NEAR(charge_g[ig].real(), second_amplitude, 1e-14); + const double expected_imaginary = -amplitude; + EXPECT_NEAR(charge_g[ig].imag(), expected_imaginary, 1e-14); + } + } + Parallel_Reduce::reduce_pool(selected); + EXPECT_DOUBLE_EQ(selected, 1.0); + atoms[0].position.x += length; + std::vector translated; + ModuleSccs::gaussian_ionic_density(atoms, basis, tpiba, spread, translated); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(translated[ir], density[ir], 1e-13); + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp b/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp new file mode 100644 index 00000000000..9b482dec854 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_ionic_force.cpp @@ -0,0 +1,97 @@ +#include "sccs_test.h" +#include "../sccs_ionic_force.h" +#include "../sccs_ionic_charge.h" +#include "../sccs_response.h" +#include "../sccs_parameters.h" +#include "../sccs_functional.h" + +#include "source_cell/cell_tools.h" + +using SccsIonicForceTest = SccsTest::PwTest; + +TEST_F(SccsIonicForceTest, AnalyticFourierForceAndGaugeInvariance) +{ + std::vector atoms(1); + atoms[0].valence_charge = 2.0; + atoms[0].position.x = length / 4.0; + std::vector potential = cosine_mode(0); + const double spread = ModuleSccs::gaussian_ion_spread; + std::vector> forces; + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + const double exponent = -0.25 * spread * spread * tpiba * tpiba; + const double expected = 2.0 * tpiba * std::exp(exponent); + EXPECT_NEAR(forces[0].x, expected, 1e-12); + EXPECT_NEAR(forces[0].y, 0.0, 1e-12); + EXPECT_NEAR(forces[0].z, 0.0, 1e-12); + atoms[0].position.x += length; + for (double& value : potential) { value += 4.0; } + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + EXPECT_NEAR(forces[0].x, expected, 1e-12); + potential.assign(basis.nrxx, 3.0); + ModuleSccs::gaussian_ionic_force(atoms, potential, basis, tpiba, spread, forces); + EXPECT_NEAR(forces[0].x, 0.0, 1e-12); + EXPECT_NEAR(forces[0].y, 0.0, 1e-12); + EXPECT_NEAR(forces[0].z, 0.0, 1e-12); +} + +TEST_F(SccsIonicForceTest, ReactionEnergyFiniteDifferenceAtFixedNonuniformCavity) +{ + std::vector atoms(2); + atoms[0].valence_charge = 1.0; + atoms[0].position = ModuleBase::Vector3(2.1, 3.2, 4.3); + atoms[1].valence_charge = 1.0; + atoms[1].position = ModuleBase::Vector3(6.4, 5.1, 2.7); + std::vector density = cosine_mode(0); + for (double& value : density) { value = (2.0 / basis.omega) * (1.0 + 0.1 * value); } + ModuleSccs::SccsConfig config; + config.cavity.density_min = 1e-3; + config.cavity.density_max = 3e-3; + config.cavity.epsilon_bulk = 5.0; + config.surface_regularization = 1e-8; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-13; + solver.tolerance_max = 1e-12; + const double spread = ModuleSccs::gaussian_ion_spread; + const std::vector cold_start; + auto evaluate = [&](const std::vector& displaced, + ModuleSccs::FunctionalResult& functional) { + std::vector charge; + ModuleSccs::gaussian_ionic_density(displaced, basis, tpiba, spread, charge); + for (int ir = 0; ir < basis.nrxx; ++ir) { charge[ir] -= density[ir]; } + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, config.cavity, solver, cold_start, basis, tpiba, response); + ModuleSccs::evaluate_functional(charge, response, config, basis, tpiba, functional); + }; + ModuleSccs::FunctionalResult baseline; + evaluate(atoms, baseline); + std::vector> forces; + ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); + const double step = 1e-4; + for (std::size_t ia = 0; ia < atoms.size(); ++ia) + { + for (int axis = 0; axis < 3; ++axis) + { + std::vector plus = atoms; + std::vector minus = atoms; + plus[ia].position[axis] += step; + minus[ia].position[axis] -= step; + ModuleSccs::FunctionalResult positive; + ModuleSccs::FunctionalResult negative; + evaluate(plus, positive); + evaluate(minus, negative); + const double finite_difference = -(positive.reaction_energy - negative.reaction_energy) / (2.0 * step); + EXPECT_NEAR(forces[ia][axis], finite_difference, 1e-8); + } + } + config.cavity.epsilon_bulk = 1.0; + evaluate(atoms, baseline); + ModuleSccs::gaussian_ionic_force(atoms, baseline.reaction_potential, basis, + tpiba, spread, forces); + for (const ModuleBase::Vector3& force : forces) + { + EXPECT_NEAR(force.x, 0.0, 1e-12); + EXPECT_NEAR(force.y, 0.0, 1e-12); + EXPECT_NEAR(force.z, 0.0, 1e-12); + } +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp new file mode 100644 index 00000000000..aeedf8f4447 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_parameters.cpp @@ -0,0 +1,43 @@ +#include "../sccs_parameters.h" +#include + +TEST(SccsParameters, ExactPresetNamesOnly) +{ + EXPECT_EQ(ModuleSccs::parse_preset("water-neutral"), ModuleSccs::Preset::WaterNeutral); + EXPECT_EQ(ModuleSccs::parse_preset("water-cation"), ModuleSccs::Preset::WaterCation); + EXPECT_EQ(ModuleSccs::parse_preset("water-anion"), ModuleSccs::Preset::WaterAnion); + EXPECT_EQ(ModuleSccs::parse_preset("vacuum"), ModuleSccs::Preset::Vacuum); + EXPECT_EQ(ModuleSccs::parse_preset("custom"), ModuleSccs::Preset::Custom); + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleSccs::parse_preset("unknown"), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(ModuleSccs::parse_preset("Water-Neutral"), ::testing::ExitedWithCode(1), ""); + testing::internal::GetCapturedStdout(); +} + +// The water presets in Hartree atomic units, which also fixes the dyn/cm and +// GPa conversions, and the vacuum preset. +TEST(SccsParameters, PresetValues) +{ + const ModuleSccs::Preset presets[] = {ModuleSccs::Preset::WaterNeutral, + ModuleSccs::Preset::WaterCation, + ModuleSccs::Preset::WaterAnion}; + const double minima[] = {1e-4, 2e-4, 2.4e-3}; + const double maxima[] = {5e-3, 3.5e-3, 1.55e-2}; + const double tensions[] = {3.076640259577260e-5, 3.211524279308205e-6, 0.0}; + const double pressures[] = {-1.223615131827540e-5, 4.248663652178958e-6, 1.529518914784425e-5}; + for (int index = 0; index < 3; ++index) + { + const ModuleSccs::SccsConfig config = ModuleSccs::make_sccs_config(presets[index]); + EXPECT_DOUBLE_EQ(config.cavity.density_min, minima[index]); + EXPECT_DOUBLE_EQ(config.cavity.density_max, maxima[index]); + EXPECT_DOUBLE_EQ(config.cavity.epsilon_bulk, 78.3); + EXPECT_DOUBLE_EQ(config.surface_regularization, 1e-8); + EXPECT_NEAR(config.surface_tension, tensions[index], 1e-18); + EXPECT_NEAR(config.pressure, pressures[index], 1e-18); + } + const ModuleSccs::SccsConfig vacuum = ModuleSccs::make_sccs_config(ModuleSccs::Preset::Vacuum); + EXPECT_DOUBLE_EQ(vacuum.cavity.epsilon_bulk, 1.0); + EXPECT_DOUBLE_EQ(vacuum.surface_tension, 0.0); + EXPECT_DOUBLE_EQ(vacuum.pressure, 0.0); +} + diff --git a/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp b/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp new file mode 100644 index 00000000000..a36b7247fef --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_pw_coulomb.cpp @@ -0,0 +1,68 @@ +#include "sccs_test.h" +#include "../sccs_pw_coulomb.h" + +#include "source_base/parallel_global.h" + +namespace SccsTest +{ +int pool_size = 1; +int pool_rank = 0; +} + +using SccsPwCoulombTest = SccsTest::PwTest; + +TEST_F(SccsPwCoulombTest, FourierModeAndConstantBackground) +{ + const std::vector mode = cosine_mode(0); + std::vector charge = mode; + for (double& value : charge) + { + value += 0.3; + } + ModuleSccs::PeriodicCoulombOperator coulomb(basis, tpiba); + std::vector potential; + coulomb.apply_potential(charge, potential); + const double kernel = ModuleBase::FOUR_PI / (tpiba * tpiba); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = kernel * mode[ir]; + EXPECT_NEAR(potential[ir], expected, 1e-11); + } + charge.assign(basis.nrxx, 1.0); + coulomb.apply_potential(charge, potential); + for (double value : potential) + { + EXPECT_NEAR(value, 0.0, 1e-12); + } +} + +int main(int argc, char** argv) +{ + int threads = 1; + Parallel_Global::read_pal_param(argc, argv, SccsTest::pool_size, threads, SccsTest::pool_rank); + int band_size = 1; + int band_rank = 0; + int band_group = 0; + int pool = 0; +#ifdef __MPI + // These tests create only pool/band groups, not diagonalization/grid groups. + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; + Parallel_Global::init_pools(SccsTest::pool_size, + SccsTest::pool_rank, + 1, + 1, + band_size, + band_rank, + band_group, + SccsTest::pool_size, + SccsTest::pool_rank, + pool); +#endif + testing::InitGoogleTest(&argc, argv); + const int status = RUN_ALL_TESTS(); +#ifdef __MPI + Parallel_Global::finalize_mpi(); +#endif + return status; +} diff --git a/source/source_hamilt/module_sccs/test/test_sccs_response.cpp b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp new file mode 100644 index 00000000000..1ff9242bc88 --- /dev/null +++ b/source/source_hamilt/module_sccs/test/test_sccs_response.cpp @@ -0,0 +1,149 @@ +#include "sccs_test.h" +#include "../sccs_response.h" +#include "../sccs_parameters.h" + +#include "source_base/parallel_reduce.h" + +#include + +using SccsResponseTest = SccsTest::PwTest; + +TEST_F(SccsResponseTest, VacuumAndUniformDielectricAnalyticLimits) +{ + const double dielectrics[] = {1.0, 5.0}; + const std::vector charge = cosine_mode(0); + const std::vector density(basis.nrxx, 0.0); + const std::vector cold_start; + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + for (double epsilon : dielectrics) + { + cavity.epsilon_bulk = epsilon; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); + const double kernel = ModuleBase::FOUR_PI / (epsilon * tpiba * tpiba); + EXPECT_EQ(response.polarization.iterations, 1); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double expected = kernel * charge[ir]; + EXPECT_NEAR(response.polarization.potential[ir], expected, 1e-11); + EXPECT_DOUBLE_EQ(response.cavity_potential[ir], 0.0); + } + } +} + +TEST_F(SccsResponseTest, ManufacturedNonuniformDielectricSolution) +{ + const std::vector mode_x = cosine_mode(0); + const std::vector mode_y = cosine_mode(1); + ModuleSccs::CavityParameters cavity; + cavity.density_min = 1e-4; + cavity.density_max = 5e-3; + cavity.epsilon_bulk = 78.3; + std::vector density(basis.nrxx); + std::vector charge(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 1e-3 + 2e-5 * mode_x[ir]; + const ModuleSccs::CavityPoint point = ModuleSccs::evaluate_cavity(density[ir], cavity); + // v=cos(ky), eps=eps(x): -div(eps grad v)/(4 pi)=eps k^2 cos(ky)/(4 pi). + charge[ir] = point.epsilon * tpiba * tpiba * mode_y[ir] / ModuleBase::FOUR_PI; + } + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-12; + solver.tolerance_max = 1e-12; + const std::vector cold_start; + ModuleSccs::SccsResponse response; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, response); + double maximum_error = 0.0; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const double difference = response.polarization.potential[ir] - mode_y[ir]; + const double magnitude = std::abs(difference); + maximum_error = std::max(maximum_error, magnitude); + EXPECT_NEAR(response.polarization.gradient[ir].x, 0.0, 1e-9); + EXPECT_NEAR(response.polarization.gradient[ir].z, 0.0, 1e-9); + } + Parallel_Reduce::reduce_max_pool(basis.poolnproc, maximum_error); + EXPECT_LT(maximum_error, 1e-9); + EXPECT_GT(response.polarization.iterations, 1); +} + +TEST_F(SccsResponseTest, ConvergenceFailureBothTolerancesAndExplicitWarmStart) +{ + const std::vector mode_x = cosine_mode(0); + std::vector charge = cosine_mode(1); + std::vector density(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + density[ir] = 0.009 + 0.008 * mode_x[ir]; + charge[ir] *= 1e-3; + } + ModuleSccs::CavityParameters cavity; + cavity.density_min = 0.0024; + cavity.density_max = 0.0155; + cavity.epsilon_bulk = 78.3; + ModuleSccs::PolarizationSolverParameters solver; + solver.tolerance_rms = 1e-11; + solver.tolerance_max = 1e-10; + const std::vector cold_start; + ModuleSccs::SccsResponse response; + // A missed tolerance stops the run; a forked death test cannot share the + // pool collectives of a multi-rank run. + if (SccsTest::pool_size == 1) + { + solver.max_iterations = 1; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, + response), + ::testing::ExitedWithCode(1), ""); + testing::internal::GetCapturedStdout(); + } + solver.max_iterations = 200; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, basis, tpiba, response); + EXPECT_LE(response.polarization.residual_rms, solver.tolerance_rms); + EXPECT_LE(response.polarization.residual_max, solver.tolerance_max); + for (double& value : density) + { + value *= 1.001; + } + ModuleSccs::SccsResponse cold; + ModuleSccs::SccsResponse warm; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, cold); + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, response.restart_potential, + basis, tpiba, warm); + EXPECT_TRUE(warm.polarization.warm_started); + EXPECT_LT(warm.polarization.iterations, cold.polarization.iterations); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + EXPECT_NEAR(warm.polarization.potential[ir], cold.polarization.potential[ir], 1e-9); + } + + // Each of the two tolerances alone drives the iteration count. + solver.tolerance_rms = 1e-5; + solver.tolerance_max = 1e-5; + ModuleSccs::SccsResponse loose; + ModuleSccs::SccsResponse tight_maximum; + ModuleSccs::SccsResponse tight_rms; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, loose); + solver.tolerance_max = 1e-11; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_maximum); + solver.tolerance_rms = 1e-11; + solver.tolerance_max = 1.0; + ModuleSccs::solve_sccs_response(density, charge, cavity, solver, cold_start, + basis, tpiba, tight_rms); + EXPECT_LE(tight_maximum.polarization.residual_max, 1e-11); + EXPECT_LE(tight_rms.polarization.residual_rms, 1e-11); + EXPECT_GT(tight_maximum.polarization.iterations, loose.polarization.iterations); + EXPECT_GT(tight_rms.polarization.iterations, loose.polarization.iterations); +} + diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index 76da6b23eba..e3c52fbd540 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -86,6 +86,21 @@ void Potential::cal_fixed_v(double* vl_pseudo) { } +void Potential::get_solvation_energy(double& electrostatic, double& non_electrostatic) const +{ + electrostatic = 0.0; + non_electrostatic = 0.0; +} + +void Potential::add_solvation_force(const UnitCell& cell, ModuleBase::matrix& force) const +{ +} + +const std::vector* Potential::solvent_electrostatic_potential() const +{ + return nullptr; +} + } // namespace elecstate #include "source_cell/klist.h" diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 488235c6b1e..afa0b154b90 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -131,6 +131,7 @@ add_library( module_parameter/read_inp_dfpt.cpp module_parameter/read_inp_deepks.cpp module_parameter/read_inp_model.cpp + module_parameter/read_inp_sccs.cpp module_parameter/read_inp_postproc.cpp module_parameter/read_inp_exx_dftu.cpp module_parameter/read_inp_other.cpp diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index b057d58fe7a..12c8e7d4c82 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -152,6 +152,8 @@ void ctrl_output_fp(UnitCell& ucell, std::string fn = PARAM.globalv.global_out_dir + "potes"; fn += geom_block + ".cube"; + const std::vector* pcc_potential = pelec->pot->pcc_electron_potential(); + const std::vector* solvent_potential = pelec->pot->solvent_electrostatic_potential(); ModuleIO::write_elecstat_pot( #ifdef __MPI pw_big->bz, @@ -164,6 +166,8 @@ void ctrl_output_fp(UnitCell& ucell, &(ucell), pelec->pot->get_fixed_v(), solvent, + pcc_potential, + solvent_potential, inp.out_pot[1]); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 7deb3732dbf..54e408c593b 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -35,7 +35,8 @@ struct Input_para int kpar = 1; ///< ecch pool is for one k point int bndpar = 1; ///< parallel for stochastic/deterministic bands std::string latname = "user_defined_lattice"; ///< lattice name - std::string assume_isolated = "none"; ///< isolated-system correction: none or makov-payne + std::string assume_isolated = "none"; ///< isolated-system correction + int pcc_2d_axis = 2; ///< open lattice vector for PCC 2D double ecutwfc = 0; ///< energy cutoff for wavefunctions double ecutrho = 0; ///< energy cutoff for charge/potential @@ -571,12 +572,21 @@ struct Input_para double block_height = 0.1; ///< height of the block // implicit solvation model Menglin Sun added on 2022-04-04 - bool imp_sol = false; ///< true: implicit solvation correction; false: - ///< vacuum calculation(default) + int imp_sol = 0; ///< 0: vacuum, 1: legacy solvation, 2: SCCS double eb_k = 80; ///< the relative permittivity of the bulk solvent double tau = 1.0798e-05; ///< the effective surface tension parameter double sigma_k = 0.6; ///< the width of the diffuse cavity double nc_k = 0.00037; ///< the cut-off charge density + std::string sccs_preset = "custom"; + double sccs_epsilon = 1.0; // Environ env_static_permittivity default + double sccs_rho_min = 1e-4; + double sccs_rho_max = 5e-3; + double sccs_gamma = 0.0; // dyn/cm + double sccs_pressure = 0.0; // GPa + int sccs_maxiter = 200; + double sccs_tol_rms = 1e-10; + double sccs_tol_max = 1e-8; + double sccs_surface_eta = 1e-8; // Bohr^-1 // ============== #Parameters (14.vdW Correction) =========================== // ========================================================== diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index aa1d232feb5..b9976ecacf6 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/source/source_io/module_parameter/read_inp_model.cpp @@ -2,6 +2,7 @@ #include "source_base/tool_quit.h" #include "read_input.h" #include "read_input_tool.h" +#include "read_inp_sccs.h" namespace ModuleIO { void ReadInput::item_model() @@ -181,16 +182,25 @@ void ReadInput::item_model() this->add_item(item); } - // imlicit_solvation + // Implicit solvation: imp_sol, the original model, then SCCS. { Input_Item item("imp_sol"); - item.annotation = "calculate implicit solvation correction or not"; + item.annotation = "implicit solvent model"; item.category = "Implicit solvation model"; - item.type = "Boolean"; - item.description = "Calculate implicit solvation correction"; - item.default_value = "False"; - item.unit = ""; - read_sync_bool(input.imp_sol); + item.type = "Integer"; + item.description = "Select 0 for vacuum, 1 for the original ABACUS implicit solvation model, or 2 for SCCS. Stress, external fields and other correction models are not supported."; + item.default_value = "0"; + item.read_value = [](const Input_Item& item, Parameter& para) { + std::string error; + const bool valid = parse_solvation_model(item.str_values[0], para.input.imp_sol, error); + if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } + }; + sync_int(input.imp_sol); + item.check_value = [](const Input_Item&, const Parameter& para) { + std::string error; + const bool valid = validate_sccs_input(para.input, error); + if (!valid) { ModuleBase::WARNING_QUIT("ReadInput", error); } + }; this->add_item(item); } { @@ -201,7 +211,7 @@ void ReadInput::item_model() item.description = "The relative permittivity of the bulk solvent, 80 for water"; item.default_value = "80"; item.unit = ""; - item.set_availability("imp_sol==true"); + item.set_availability("imp_sol==1"); read_sync_double(input.eb_k); this->add_item(item); } @@ -238,6 +248,7 @@ void ReadInput::item_model() read_sync_double(input.nc_k); this->add_item(item); } + this->item_sccs(); // vdW Correction { diff --git a/source/source_io/module_parameter/read_inp_sccs.cpp b/source/source_io/module_parameter/read_inp_sccs.cpp new file mode 100644 index 00000000000..1d651c177fe --- /dev/null +++ b/source/source_io/module_parameter/read_inp_sccs.cpp @@ -0,0 +1,220 @@ +#include "read_inp_sccs.h" +#include "read_input.h" +#include "read_input_tool.h" +#include "source_base/tool_quit.h" + +#include +#include + +namespace ModuleIO +{ +bool parse_solvation_model(const std::string& value, int& model, std::string& error) +{ + error.clear(); + std::string normalized = value; + for (char& c : normalized) + { + const unsigned char byte = static_cast(c); + const int lower = std::tolower(byte); + c = static_cast(lower); + } + const std::vector yes = {"true", "1", "t", "yes", "y", "on", ".true."}; + const std::vector no = {"false", "0", "f", "no", "n", "off", ".false."}; + if (normalized == "2") { model = 2; } + else if (std::find(yes.begin(), yes.end(), normalized) != yes.end()) { model = 1; } + else if (std::find(no.begin(), no.end(), normalized) != no.end()) { model = 0; } + else + { + error = "imp_sol must be 0, 1, 2 or a legacy Boolean value"; + return false; + } + return true; +} + +namespace +{ +// Calculation settings that SCCS supports; empty when all are supported. +std::string sccs_context_error(const Input_para& input) +{ + std::string error; + if (input.assume_isolated != "none") { error = "SCCS currently requires assume_isolated none"; } + else if (input.device != "cpu" || input.esolver_type != "ksdft") { error = "SCCS requires CPU KS-DFT"; } + else if (input.basis_type != "pw" && input.basis_type != "lcao") { error = "SCCS requires basis_type pw or lcao"; } + else if (input.calculation != "scf" || input.cal_stress) + { error = "SCCS currently supports SCF energies and forces without stress"; } + else if (input.nspin != 1 && input.nspin != 2) { error = "SCCS requires nspin 1 or 2"; } + else if (input.efield_flag || input.gate_flag || input.dfthalf_type != 0 + || input.deepks_scf || input.deepks_out_labels || input.deepks_bandgap + || input.deepks_v_delta || input.dm_to_rho) + { error = "SCCS does not support electric/gate fields, DFT-1/2, DeePKS or dm_to_rho"; } + else if (!input.vl_in_h || !input.vion_in_h || !input.vh_in_h) + { error = "SCCS requires the local ionic and Hartree potentials in the Hamiltonian"; } + return error; +} + +// Cavity, surface and boundary model parameters, starting with the preset. +std::string sccs_model_error(const Input_para& input) +{ + const std::vector presets = {"custom", "vacuum", "water-neutral", "water-cation", "water-anion"}; + if (std::find(presets.begin(), presets.end(), input.sccs_preset) == presets.end()) + { + return "Unknown sccs_preset"; + } + std::string error; + if (input.sccs_epsilon < 1.0) { error = "sccs_epsilon must be at least 1"; } + else if (input.sccs_rho_min <= 0.0 || input.sccs_rho_max <= input.sccs_rho_min) + { error = "sccs_rho_min and sccs_rho_max must satisfy 0 < sccs_rho_min < sccs_rho_max"; } + else if (input.sccs_surface_eta <= 0.0) { error = "sccs_surface_eta must be positive"; } + return error; +} + +// sqrt-CG solver controls. +std::string sccs_solver_error(const Input_para& input) +{ + std::string error; + if (input.sccs_maxiter <= 0) { error = "sccs_maxiter must be positive"; } + else if (input.sccs_tol_rms <= 0.0) { error = "sccs_tol_rms must be positive"; } + else if (input.sccs_tol_max <= 0.0) { error = "sccs_tol_max must be positive"; } + return error; +} +} // namespace + +bool validate_sccs_input(const Input_para& input, std::string& error) +{ + error.clear(); + if (input.imp_sol < 0 || input.imp_sol > 2) + { + error = "imp_sol must be 0, 1 or 2"; + return false; + } + if (input.imp_sol != 2) { return true; } + error = sccs_context_error(input); + if (error.empty()) { error = sccs_model_error(input); } + if (error.empty()) { error = sccs_solver_error(input); } + return error.empty(); +} + +void ReadInput::item_sccs() +{ + { + Input_Item item("sccs_preset"); + item.annotation = "SCCS parameter preset"; + item.category = "Implicit solvation model"; + item.type = "String"; + item.description = "SCCS parameter preset. Allowed values: custom, vacuum, water-neutral, water-cation, water-anion. Non-custom presets override sccs_epsilon, sccs_rho_min, sccs_rho_max, sccs_gamma and sccs_pressure; solver controls and surface regularization remain user-controlled."; + item.default_value = "custom"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_string(input.sccs_preset); + this->add_item(item); + } + { + Input_Item item("sccs_epsilon"); + item.annotation = "Bulk dielectric constant >= 1 for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Bulk dielectric constant >= 1 for the custom preset. The default 1 (no dielectric) follows the Environ env_static_permittivity default for environ_type input. Water presets use 78.3 and vacuum uses 1."; + item.default_value = "1.0"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_epsilon); + this->add_item(item); + } + { + Input_Item item("sccs_rho_min"); + item.annotation = "Lower electronic-density cavity threshold for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Lower electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 1e-4; water-cation: 2e-4; water-anion: 2.4e-3."; + item.default_value = "1e-4"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_rho_min); + this->add_item(item); + } + { + Input_Item item("sccs_rho_max"); + item.annotation = "Upper electronic-density cavity threshold for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Upper electronic-density cavity threshold for the custom preset. Water-neutral/vacuum: 5e-3; water-cation: 3.5e-3; water-anion: 1.55e-2."; + item.default_value = "5e-3"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_rho_max); + this->add_item(item); + } + { + Input_Item item("sccs_gamma"); + item.annotation = "Surface coefficient for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Surface coefficient for the custom preset. Water-neutral: 47.9; water-cation: 5; water-anion/vacuum: 0."; + item.default_value = "0.0"; + item.unit = "dyn/cm"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_gamma); + this->add_item(item); + } + { + Input_Item item("sccs_pressure"); + item.annotation = "Volume coefficient for the custom preset"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Volume coefficient for the custom preset. Water-neutral: -0.36; water-cation: 0.125; water-anion: 0.45; vacuum: 0."; + item.default_value = "0.0"; + item.unit = "GPa"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_pressure); + this->add_item(item); + } + { + Input_Item item("sccs_maxiter"); + item.annotation = "Positive maximum number of SCCS sqrt-CG iterations"; + item.category = "Implicit solvation model"; + item.type = "Integer"; + item.description = "Positive maximum number of SCCS sqrt-CG iterations. Failure to satisfy both residual tolerances terminates the calculation."; + item.default_value = "200"; + item.unit = ""; + item.set_availability("imp_sol==2"); + read_sync_int(input.sccs_maxiter); + this->add_item(item); + } + { + Input_Item item("sccs_tol_rms"); + item.annotation = "Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive RMS charge-residual tolerance for the SCCS sqrt-CG solver."; + item.default_value = "1e-10"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_tol_rms); + this->add_item(item); + } + { + Input_Item item("sccs_tol_max"); + item.annotation = "Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive maximum charge-residual tolerance for the SCCS sqrt-CG solver."; + item.default_value = "1e-8"; + item.unit = "e/bohr^3"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_tol_max); + this->add_item(item); + } + { + Input_Item item("sccs_surface_eta"); + item.annotation = "Positive regularization of the SCCS surface gradient norm"; + item.category = "Implicit solvation model"; + item.type = "Real"; + item.description = "Positive regularization of the SCCS surface gradient norm."; + item.default_value = "1e-8"; + item.unit = "bohr^-1"; + item.set_availability("imp_sol==2"); + read_sync_double(input.sccs_surface_eta); + this->add_item(item); + } +} +} // namespace ModuleIO diff --git a/source/source_io/module_parameter/read_inp_sccs.h b/source/source_io/module_parameter/read_inp_sccs.h new file mode 100644 index 00000000000..8736ed5768e --- /dev/null +++ b/source/source_io/module_parameter/read_inp_sccs.h @@ -0,0 +1,10 @@ +#ifndef READ_INP_SCCS_H +#define READ_INP_SCCS_H +#include +struct Input_para; +namespace ModuleIO +{ +bool parse_solvation_model(const std::string& value, int& model, std::string& error); +bool validate_sccs_input(const Input_para& input, std::string& error); +} +#endif diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 97a40123eb8..a0ba7ff398a 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -506,9 +506,13 @@ Available options are: Available options are: * none: regular periodic calculation without isolated-system correction. +* pcc_0d: self-consistent point-counter-charge (PCC) correction for a molecule in a cubic cell. The multipoles are taken about the mass-weighted ionic center. +* 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). -Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995).)"; +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. + +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"; read_sync_string(input.assume_isolated); item.reset_value = [](const Input_Item& item, Parameter& para) { @@ -518,7 +522,7 @@ Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995).)"; } }; item.check_value = [](const Input_Item& item, const Parameter& para) { - const std::vector allowed = {"none", "makov-payne", "m-p", "mp"}; + const std::vector allowed = {"none", "makov-payne", "m-p", "mp", "pcc_0d", "pcc_2d"}; if (std::find(allowed.begin(), allowed.end(), para.input.assume_isolated) == allowed.end()) { ModuleBase::WARNING_QUIT("ReadInput", nofound_str(allowed, "assume_isolated")); @@ -530,6 +534,48 @@ Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995).)"; ModuleBase::WARNING_QUIT("ReadInput", "Makov-Payne correction is available only for cubic lattices: latname = sc, fcc, or bcc."); } + const Input_para& inp = para.input; + 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"; } + 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"; } + else if (inp.basis_type != "pw" && inp.basis_type != "lcao") { error = "requires basis_type pw or lcao"; } + else if (inp.calculation != "scf" && inp.calculation != "relax") { error = "requires calculation scf or fixed-cell relax"; } + else if (inp.efield_flag) { error = "cannot be combined with efield_flag"; } + else if (inp.gate_flag) { error = "cannot be combined with gate_flag"; } + else if (inp.cal_stress) { error = "does not support cal_stress"; } + else if (inp.dfthalf_type != 0) { error = "cannot be combined with dfthalf_type"; } + else if (inp.deepks_scf || inp.deepks_out_labels || inp.deepks_bandgap || inp.deepks_v_delta) + { + error = "cannot be combined with DeePKS correction or output"; + } + else if (inp.dm_to_rho) { error = "cannot be combined with dm_to_rho"; } + if (!error.empty()) + { + const std::string message = inp.assume_isolated + " " + error; + ModuleBase::WARNING_QUIT("ReadInput", message); + } + } + }; + this->add_item(item); + } + { + Input_Item item("pcc_2d_axis"); + item.annotation = "open lattice vector of assume_isolated pcc_2d"; + item.category = "System variables"; + item.type = "Integer"; + item.description = "Index of the lattice vector along which assume_isolated=pcc_2d is open: 0, 1 or 2 for the first, second or third vector of LATTICE_VECTORS. The selected vector must be perpendicular to the other two, which span the periodic plane; the k-point sampling along it must be Gamma only."; + item.default_value = "2"; + item.set_availability("assume_isolated==pcc_2d"); + read_sync_int(input.pcc_2d_axis); + item.check_value = [](const Input_Item&, const Parameter& para) { + if (para.input.pcc_2d_axis < 0 || para.input.pcc_2d_axis > 2) + { + ModuleBase::WARNING_QUIT("ReadInput", "pcc_2d_axis must be 0, 1 or 2"); + } }; this->add_item(item); } diff --git a/source/source_io/module_parameter/read_input.h b/source/source_io/module_parameter/read_input.h index 37484a21247..02d6a95d71c 100644 --- a/source/source_io/module_parameter/read_input.h +++ b/source/source_io/module_parameter/read_input.h @@ -139,6 +139,7 @@ class ReadInput void item_postprocess(); // items for some models void item_model(); + void item_sccs(); // items for exx void item_exx(); // items for dft+u diff --git a/source/source_io/test_serial/CMakeLists.txt b/source/source_io/test_serial/CMakeLists.txt index 3238fb743d4..0201fb14e01 100644 --- a/source/source_io/test_serial/CMakeLists.txt +++ b/source/source_io/test_serial/CMakeLists.txt @@ -17,6 +17,7 @@ add_library( ../module_parameter/read_inp_dfpt.cpp ../module_parameter/read_inp_deepks.cpp ../module_parameter/read_inp_model.cpp + ../module_parameter/read_inp_sccs.cpp ../module_parameter/read_inp_postproc.cpp ../module_parameter/read_inp_exx_dftu.cpp ../module_parameter/read_inp_other.cpp @@ -40,7 +41,7 @@ AddTest( AddTest( TARGET MODULE_IO_read_item_serial LIBS parameter base device io_input_serial - SOURCES read_input_item_test.cpp + SOURCES read_input_item_test.cpp test_read_inp_sccs.cpp ) AddTest( diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index ffe37e6ba80..4ace2f6f5a1 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -98,6 +98,53 @@ TEST_F(InputTest, RelaxMethod) EXPECT_EQ(find_label("relax_new", readinput.input_lists), readinput.input_lists.end()); } +TEST_F(InputTest, PccSelectionAndUnsupportedCombinations) +{ + ModuleIO::ReadInput readinput(0); + Parameter param; + Input_para& input = TestParameters::input(param); + input.device = "cpu"; + input.esolver_type = "ksdft"; + input.basis_type = "pw"; + input.calculation = "scf"; + auto it = find_label("assume_isolated", readinput.input_lists); + ASSERT_NE(it, readinput.input_lists.end()); + it->second.str_values = {"pcc_0d"}; + it->second.read_value(it->second, param); + EXPECT_EQ(input.assume_isolated, "pcc_0d"); + it->second.check_value(it->second, param); + input.basis_type = "lcao"; + input.calculation = "relax"; + it->second.check_value(it->second, param); + const Input_para valid = input; + + input.nspin = 4; + EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); + input = valid; + input.device = "gpu"; + EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); + input = valid; + input.cal_stress = true; + EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); + input = valid; + input.calculation = "cell-relax"; + EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); + input = valid; + input.assume_isolated = "pcc_2d"; + it->second.check_value(it->second, param); + auto axis = find_label("pcc_2d_axis", readinput.input_lists); + ASSERT_NE(axis, readinput.input_lists.end()); + for (int value = 0; value < 3; ++value) + { + input.pcc_2d_axis = value; + axis->second.check_value(axis->second, param); + } + input.pcc_2d_axis = 3; + EXPECT_EXIT(axis->second.check_value(axis->second, param), testing::ExitedWithCode(1), ""); + input.assume_isolated = "pcc_invalid"; + EXPECT_EXIT(it->second.check_value(it->second, param), testing::ExitedWithCode(1), ""); +} + TEST_F(InputTest, Item_test) { ModuleIO::ReadInput readinput(0); diff --git a/source/source_io/test_serial/test_read_inp_sccs.cpp b/source/source_io/test_serial/test_read_inp_sccs.cpp new file mode 100644 index 00000000000..94cb6bad718 --- /dev/null +++ b/source/source_io/test_serial/test_read_inp_sccs.cpp @@ -0,0 +1,56 @@ +#include "source_io/module_parameter/read_inp_sccs.h" +#include "source_io/module_parameter/input_parameter.h" +#include + +TEST(ReadInpSccs, LegacyBooleansAndModelSelection) +{ + std::string error; + const std::string yes[] = {"true", "TRUE", "yes", "On", ".true.", "1", "t", "y"}; + const std::string no[] = {"false", "FALSE", "no", "Off", ".false.", "0", "f", "n"}; + int model = 0; + for (const std::string& value : yes) + { + ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + EXPECT_EQ(model, 1); + } + for (const std::string& value : no) + { + ASSERT_TRUE(ModuleIO::parse_solvation_model(value, model, error)); + EXPECT_EQ(model, 0); + } + ASSERT_TRUE(ModuleIO::parse_solvation_model("2", model, error)); + EXPECT_EQ(model, 2); + EXPECT_FALSE(ModuleIO::parse_solvation_model("3", model, error)); + EXPECT_EQ(model, 2); + EXPECT_FALSE(ModuleIO::parse_solvation_model("2.0", model, error)); +} + +TEST(ReadInpSccs, SupportedScopeAndNumericalValidation) +{ + Input_para input; + input.imp_sol = 2; + input.device = "cpu"; + input.basis_type = "pw"; + std::string error; + ASSERT_TRUE(ModuleIO::validate_sccs_input(input, error)) << error; + input.cal_force = true; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + input.cal_stress = true; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.cal_stress = false; + input.cal_force = false; + input.assume_isolated = "pcc_0d"; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.assume_isolated = "none"; + input.sccs_rho_min = input.sccs_rho_max; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.sccs_rho_min = 1e-4; + input.sccs_preset = "unknown"; + EXPECT_FALSE(ModuleIO::validate_sccs_input(input, error)); + input.sccs_preset = "vacuum"; + input.basis_type = "lcao"; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); + input.imp_sol = 1; + input.sccs_preset = "unused"; + EXPECT_TRUE(ModuleIO::validate_sccs_input(input, error)); +} diff --git a/source/source_lcao/force_stress_assemble.cpp b/source/source_lcao/force_stress_assemble.cpp index 6590de0c2a2..331b75d39f3 100644 --- a/source/source_lcao/force_stress_assemble.cpp +++ b/source/source_lcao/force_stress_assemble.cpp @@ -38,6 +38,10 @@ void sum_force_terms(const int iat, + parts.fcc(iat, i) // nonlinear core correction force (pw) + parts.fscc(iat, i) // self consistent corretion force (pw) + parts.fpothybrid(iat, i); // pulay force for hybrid gauge rt-tddft + if (parts.fpcc.nr != 0) + { + fcs(iat, i) += parts.fpcc(iat, i); + } // Force contribution from DFT+U, Quxin add on 20201029 if (PARAM.inp.dft_plus_u) diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index aab0dc61166..9f80519426c 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -186,6 +186,13 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // vdW force/stress and external-field forces LCAO_domain::cal_vdw_fields_fs(vdw_result, ucell, solvent, rhopw, locpp, isforce, isstress, parts, sparts); + if (isforce) + { + parts.fpcc.create(ucell.nat, 3); + pelec->pot->add_pcc_force(ucell, parts.fpcc); + if (parts.fsol.nr == 0) { parts.fsol.create(ucell.nat, 3); } + pelec->pot->add_solvation_force(ucell, parts.fsol); + } // DFT+U force/stress LCAO_domain::cal_dftu_fs(ucell, gd, pv, orb, kv, dmat, two_center_bundle, dftu, isforce, isstress, parts, sparts); diff --git a/source/source_lcao/force_stress_lcao.h b/source/source_lcao/force_stress_lcao.h index 3eeb5c1a605..d1a9b29d776 100644 --- a/source/source_lcao/force_stress_lcao.h +++ b/source/source_lcao/force_stress_lcao.h @@ -69,6 +69,7 @@ struct LCAOForceParts ModuleBase::matrix fefield_tddft; ModuleBase::matrix fgate; ModuleBase::matrix fsol; + ModuleBase::matrix fpcc; }; struct LCAOStressParts diff --git a/source/source_lcao/force_stress_terms.cpp b/source/source_lcao/force_stress_terms.cpp index 0ca1cc1db03..671d5e13628 100644 --- a/source/source_lcao/force_stress_terms.cpp +++ b/source/source_lcao/force_stress_terms.cpp @@ -92,7 +92,7 @@ void cal_external_field_forces(UnitCell& ucell, } //! atomic forces from implicit solvation model - if (PARAM.inp.imp_sol) + if (PARAM.inp.imp_sol == 1) { parts.fsol.create(ucell.nat, 3); solvent.cal_force_sol(ucell, rhopw, locpp.vloc, PARAM.inp.nspin, parts.fsol); diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index 3ea53637d95..22cfa837430 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -7,6 +7,7 @@ #include "source_lcao/setup_deepks.h" #include "source_estate/module_dm/density_matrix.h" #include "source_estate/module_pot/potential_new.h" +#include "source_estate/module_pot/pot_pcc.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include #ifdef __MLALGO @@ -86,7 +87,15 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, if (inp.vion_in_h) { pot_register_in.push_back("local"); } if (inp.vh_in_h) { pot_register_in.push_back("hartree"); } pot_register_in.push_back("xc"); - if (inp.imp_sol) { pot_register_in.push_back("surchem"); } + if (inp.imp_sol == 1) { pot_register_in.push_back("surchem"); } + if (inp.imp_sol == 2) { pot_register_in.push_back("sccs"); } + if (inp.assume_isolated == "pcc_0d") { pot_register_in.push_back("pcc"); } + if (inp.assume_isolated == "pcc_2d") + { + const int count = this->kv->get_nks(); + elecstate::PotPcc::validate_kpoints(this->kv->kvec_d, count, inp.pcc_2d_axis); + pot_register_in.push_back("pcc_2d"); + } if (inp.efield_flag) { pot_register_in.push_back("efield"); } if (inp.gate_flag) { pot_register_in.push_back("gatefield"); } if (inp.esolver_type == "tddft") { pot_register_in.push_back("tddft"); } diff --git a/source/source_lcao/hamilt_lcao_factory.cpp b/source/source_lcao/hamilt_lcao_factory.cpp index 05a7b34b42b..6280c7c3402 100644 --- a/source/source_lcao/hamilt_lcao_factory.cpp +++ b/source/source_lcao/hamilt_lcao_factory.cpp @@ -169,7 +169,7 @@ LcaoOpsBundle build_gamma_ops(const UnitCell& ucell, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, inp); // effective potential term Operator* veff = new Veff>(hsk, kv->kvec_d, pot_in, @@ -230,7 +230,7 @@ LcaoOpsBundle build_multik_ops(const UnitCell& ucell, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, inp); // Veff term ops = new Veff>(hsk, kv->kvec_d, pot_in, diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index d37e2379954..e5c93e312cf 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -137,11 +137,21 @@ void Forces::cal_force(const int nspin, } } + const int solvation_model = PARAM.inp.imp_sol; ModuleBase::matrix forcesol; - if (PARAM.inp.imp_sol) + ModuleBase::matrix forcepcc(this->nat, 3); + elec.pot->add_pcc_force(ucell, forcepcc); + if (solvation_model != 0) { forcesol.create(this->nat, 3); - solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, nspin, forcesol); + if (solvation_model == 1) + { + solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, nspin, forcesol); + } + else + { + elec.pot->add_solvation_force(ucell, forcesol); + } if (PARAM.inp.test_force) { ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); @@ -158,7 +168,7 @@ void Forces::cal_force(const int nspin, for (int ia = 0; ia < ucell.atoms[it].na; ia++) { force(iat, ipol) = forcelc(iat, ipol) + forceion(iat, ipol) + forcenl(iat, ipol) + forcecc(iat, ipol) - + forcescc(iat, ipol); + + forcescc(iat, ipol) + forcepcc(iat, ipol); if (vdw_result != nullptr) // linpz and jiyy added vdw force, modified by zhengdy { @@ -175,7 +185,7 @@ void Forces::cal_force(const int nspin, force(iat, ipol) = force(iat, ipol) + force_gate(iat, ipol); } - if (PARAM.inp.imp_sol) + if (solvation_model != 0) { force(iat, ipol) = force(iat, ipol) + forcesol(iat, ipol); } @@ -259,7 +269,7 @@ void Forces::cal_force(const int nspin, force_gate, false); } - if (PARAM.inp.imp_sol) + if (solvation_model != 0) { ModuleIO::print_force(GlobalV::ofs_running, ucell, diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index 18beccada0c..d3d284643f5 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -13,6 +13,7 @@ #include "source_base/parallel_reduce.h" #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type #include "source_io/module_parameter/parameter.h" +#include "source_estate/module_pot/pot_pcc.h" namespace hamilt { @@ -68,28 +69,40 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, } } } - if (PARAM.inp.vl_in_h) + const Input_para& input = PARAM.inp; + if (input.vl_in_h) { std::vector pot_register_in; - if (PARAM.inp.vion_in_h) + if (input.vion_in_h) { pot_register_in.push_back("local"); } - if (PARAM.inp.vh_in_h) + if (input.vh_in_h) { pot_register_in.push_back("hartree"); } // no variable can choose xc, maybe it is necessary pot_register_in.push_back("xc"); - if (PARAM.inp.imp_sol) + if (input.imp_sol == 1) { pot_register_in.push_back("surchem"); } - if (PARAM.inp.efield_flag) + if (input.imp_sol == 2) { pot_register_in.push_back("sccs"); } + if (input.assume_isolated == "pcc_0d") + { + pot_register_in.push_back("pcc"); + } + if (input.assume_isolated == "pcc_2d") + { + const int count = pkv->get_nks(); + elecstate::PotPcc::validate_kpoints(pkv->kvec_d, count, input.pcc_2d_axis); + pot_register_in.push_back("pcc_2d"); + } + if (input.efield_flag) { pot_register_in.push_back("efield"); } - if (PARAM.inp.gate_flag) + if (input.gate_flag) { pot_register_in.push_back("gatefield"); } @@ -110,7 +123,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, if (pot_register_in.size() > 0) { // register Potential by gathered operator - pot_in->pot_register(pot_register_in); + pot_in->pot_register(pot_register_in, input); Operator* veff = new Veff>(isk, pot_in->get_veff_smooth_data(), pot_in->get_veff_smooth().nr, diff --git a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp index 7023d90719b..78dfc1a44aa 100644 --- a/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp +++ b/source/source_pw/module_stodft/test/test_sto_hamilt_pw.cpp @@ -6,6 +6,9 @@ #include elecstate::Potential::~Potential(){} +void elecstate::Potential::add_solvation_force(const UnitCell&, ModuleBase::matrix&) const {} +void elecstate::Potential::get_solvation_energy(double& el, double& cav) const { el = 0.0; cav = 0.0; } +const std::vector* elecstate::Potential::solvent_electrostatic_potential() const { return nullptr; } void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase::matrix&){} void elecstate::Potential::cal_fixed_v(double*){} diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT new file mode 100644 index 00000000000..8eb0415ef14 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/INPUT @@ -0,0 +1,21 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +calculation scf +basis_type pw +init_wfc random +pw_seed 1 +nbands 8 +ecutwfc 20 +ecutrho 240 +scf_nmax 80 +scf_thr 1.0e-8 +cal_force 1 + +imp_sol 2 +assume_isolated none +sccs_preset water-neutral +sccs_maxiter 200 +sccs_tol_rms 1.0e-10 +sccs_tol_max 1.0e-8 +sccs_surface_eta 1.0e-8 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README new file mode 100644 index 00000000000..44f4976649c --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/README @@ -0,0 +1,2 @@ +PW neutral water with periodic SCCS (water-neutral preset), ecutwfc 20 and ecutrho 240. +Checks total energy, forces and electrostatic, surface and volume solvation terms. diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU new file mode 100644 index 00000000000..ed974dd43a5 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +H 1.008 H_ONCV_PBE-1.0.upf +O 15.9994 O_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +12.0 0.0 0.0 +0.0 12.0 0.0 +0.0 0.0 12.0 + +ATOMIC_POSITIONS +Cartesian + +H +0.0 +2 +7.433 7.111 6.0 0 0 0 +4.567 7.111 6.0 0 0 0 +O +0.0 +1 +6.0 6.0 6.0 0 0 0 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref new file mode 100644 index 00000000000..a1c90a54215 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/result.ref @@ -0,0 +1,6 @@ +etotref -441.6484359463214560 +etotperatomref -147.2161453154 +totalforceref 16.644169 +esolelref -0.2883000404 +esolcavref 0.0807425541 +totaltimeref 3.33 diff --git a/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold new file mode 100644 index 00000000000..a125b952d62 --- /dev/null +++ b/tests/01_PW/220_PW_SCCS_PERIODIC_WATER/threshold @@ -0,0 +1,3 @@ +threshold 0.000001 +force_threshold 0.0001 +fatal_threshold 0.01 diff --git a/tests/01_PW/221_PW_PCC0D_H3O/INPUT b/tests/01_PW/221_PW_PCC0D_H3O/INPUT new file mode 100644 index 00000000000..25328671483 --- /dev/null +++ b/tests/01_PW/221_PW_PCC0D_H3O/INPUT @@ -0,0 +1,19 @@ +INPUT_PARAMETERS +suffix autotest +calculation scf +basis_type pw +ks_solver cg +pseudo_dir ../../PP_ORB +ntype 2 +nbands 8 +ecutwfc 20 +scf_thr 1e-12 +scf_nmax 200 +mixing_type pulay +mixing_beta 0.2 +smearing_method fixed +symmetry 0 +cal_force 1 +init_wfc atomic +assume_isolated pcc_0d +nelec_delta -1 diff --git a/tests/01_PW/221_PW_PCC0D_H3O/KPT b/tests/01_PW/221_PW_PCC0D_H3O/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/221_PW_PCC0D_H3O/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/221_PW_PCC0D_H3O/README.md b/tests/01_PW/221_PW_PCC0D_H3O/README.md new file mode 100644 index 00000000000..2be075dd83d --- /dev/null +++ b/tests/01_PW/221_PW_PCC0D_H3O/README.md @@ -0,0 +1,5 @@ +Vacuum PCC 0D regression for charged H3O with `assume_isolated pcc_0d`. + +The integration harness compares total energy, ionic force magnitudes and `E_pcc` +(in eV). References use the repository pseudopotentials, `scf_thr 1e-12`, +two MPI ranks and one OpenMP thread per rank. diff --git a/tests/01_PW/221_PW_PCC0D_H3O/STRU b/tests/01_PW/221_PW_PCC0D_H3O/STRU new file mode 100644 index 00000000000..75b5d3b13e0 --- /dev/null +++ b/tests/01_PW/221_PW_PCC0D_H3O/STRU @@ -0,0 +1,26 @@ +ATOMIC_SPECIES +O 15.9994 O_ONCV_PBE-1.0.upf +H 1.008 H_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +O +0.0 +1 +10.0000000000 10.0000000000 10.0000000000 0 0 0 + +H +0.0 +3 +11.8077500000 9.6969000000 10.2918300000 0 0 0 +8.9439100000 8.5571900000 10.4978700000 0 0 0 +9.6429200000 10.5858300000 8.2760700000 0 0 0 diff --git a/tests/01_PW/221_PW_PCC0D_H3O/result.ref b/tests/01_PW/221_PW_PCC0D_H3O/result.ref new file mode 100644 index 00000000000..2b81bdee45c --- /dev/null +++ b/tests/01_PW/221_PW_PCC0D_H3O/result.ref @@ -0,0 +1,5 @@ +etotref -448.4562759354126911 +etotperatomref -112.1140689839 +totalforceref 20.652937 +epccref 2.0044964276 +totaltimeref 1.80 diff --git a/tests/01_PW/222_PW_PCC2D_H3O/INPUT b/tests/01_PW/222_PW_PCC2D_H3O/INPUT new file mode 100644 index 00000000000..d2fbefd04c6 --- /dev/null +++ b/tests/01_PW/222_PW_PCC2D_H3O/INPUT @@ -0,0 +1,20 @@ +INPUT_PARAMETERS +suffix autotest +calculation scf +basis_type pw +ks_solver cg +pseudo_dir ../../PP_ORB +ntype 2 +nbands 8 +ecutwfc 20 +scf_thr 1e-12 +scf_nmax 200 +mixing_type pulay +mixing_beta 0.2 +smearing_method fixed +symmetry 0 +cal_force 1 +init_wfc atomic +assume_isolated pcc_2d +nelec_delta -1 +pcc_2d_axis 2 diff --git a/tests/01_PW/222_PW_PCC2D_H3O/KPT b/tests/01_PW/222_PW_PCC2D_H3O/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/01_PW/222_PW_PCC2D_H3O/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/01_PW/222_PW_PCC2D_H3O/README.md b/tests/01_PW/222_PW_PCC2D_H3O/README.md new file mode 100644 index 00000000000..635d18e9523 --- /dev/null +++ b/tests/01_PW/222_PW_PCC2D_H3O/README.md @@ -0,0 +1,6 @@ +Vacuum PCC 2D regression for a charged H3O periodic layer with +`assume_isolated pcc_2d` and open axis 2 (`pcc_2d_axis 2`). + +The integration harness compares total energy, ionic force magnitudes and `E_pcc` +(in eV). References use the repository pseudopotentials, `scf_thr 1e-12`, +two MPI ranks and one OpenMP thread per rank. diff --git a/tests/01_PW/222_PW_PCC2D_H3O/STRU b/tests/01_PW/222_PW_PCC2D_H3O/STRU new file mode 100644 index 00000000000..75b5d3b13e0 --- /dev/null +++ b/tests/01_PW/222_PW_PCC2D_H3O/STRU @@ -0,0 +1,26 @@ +ATOMIC_SPECIES +O 15.9994 O_ONCV_PBE-1.0.upf +H 1.008 H_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +O +0.0 +1 +10.0000000000 10.0000000000 10.0000000000 0 0 0 + +H +0.0 +3 +11.8077500000 9.6969000000 10.2918300000 0 0 0 +8.9439100000 8.5571900000 10.4978700000 0 0 0 +9.6429200000 10.5858300000 8.2760700000 0 0 0 diff --git a/tests/01_PW/222_PW_PCC2D_H3O/result.ref b/tests/01_PW/222_PW_PCC2D_H3O/result.ref new file mode 100644 index 00000000000..448dacbcd64 --- /dev/null +++ b/tests/01_PW/222_PW_PCC2D_H3O/result.ref @@ -0,0 +1,5 @@ +etotref -451.0967403187477203 +etotperatomref -112.7741850797 +totalforceref 20.659011 +epccref -0.6360805933 +totaltimeref 1.91 diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index 6f2f88dedda..9dd0c9730af 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -123,3 +123,6 @@ scf_dftu_spin2 scf_dftu_spin4 scf_deltaspin2 scf_deltaspin4 +221_PW_PCC0D_H3O +222_PW_PCC2D_H3O +220_PW_SCCS_PERIODIC_WATER diff --git a/tests/02_NAO_Gamma/CASES_CPU.txt b/tests/02_NAO_Gamma/CASES_CPU.txt index de87db6ecd4..9ba3fa93de9 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -38,3 +38,5 @@ scf_u_spin2_old #scf_u_yukawa scf_upf100 test.sum +scf_pcc0d_water +scf_pcc2d_water diff --git a/tests/02_NAO_Gamma/scf_pcc0d_water/INPUT b/tests/02_NAO_Gamma/scf_pcc0d_water/INPUT new file mode 100644 index 00000000000..600f548e27b --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc0d_water/INPUT @@ -0,0 +1,20 @@ +INPUT_PARAMETERS +suffix autotest +calculation scf +basis_type lcao +ks_solver genelpa +pseudo_dir ../../PP_ORB +ntype 2 +nbands 8 +ecutwfc 20 +scf_thr 1e-7 +scf_nmax 200 +mixing_type pulay +mixing_beta 0.2 +smearing_method fixed +symmetry 0 +cal_force 1 +init_wfc atomic +assume_isolated pcc_0d +orbital_dir ../../PP_ORB +gamma_only 1 diff --git a/tests/02_NAO_Gamma/scf_pcc0d_water/KPT b/tests/02_NAO_Gamma/scf_pcc0d_water/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc0d_water/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/02_NAO_Gamma/scf_pcc0d_water/README.md b/tests/02_NAO_Gamma/scf_pcc0d_water/README.md new file mode 100644 index 00000000000..218f62738eb --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc0d_water/README.md @@ -0,0 +1,7 @@ +Vacuum PCC 0D regression with `assume_isolated pcc_0d`. + +The integration harness compares total energy, ionic force magnitudes and `E_pcc` +(in eV). References were generated with two MPI ranks, one OpenMP thread and +the repository pseudopotentials/orbitals. Numerical kernels have independent +analytic and finite-difference unit tests; PW neutral/charged results were also +compared against the original full PCC implementation. diff --git a/tests/02_NAO_Gamma/scf_pcc0d_water/STRU b/tests/02_NAO_Gamma/scf_pcc0d_water/STRU new file mode 100644 index 00000000000..331f96c9743 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc0d_water/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +O 15.999 O_ONCV_PBE-1.0.upf +H 1.008 H_ONCV_PBE-1.0.upf +NUMERICAL_ORBITAL +O_gga_7au_100Ry_2s2p1d.orb +H_gga_8au_100Ry_2s1p.orb + +LATTICE_CONSTANT +1.8897261258369282 +LATTICE_VECTORS +12.0 0 0 +0 12.0 0 +0 0 12.0 +ATOMIC_POSITIONS +Cartesian +O +0 +1 +6.0 6.0 6.0 1 1 1 +H +0 +2 +6.75716000000000 6.0 6.58626 1 1 1 +5.24284 6.0 6.58626 1 1 1 diff --git a/tests/02_NAO_Gamma/scf_pcc0d_water/result.ref b/tests/02_NAO_Gamma/scf_pcc0d_water/result.ref new file mode 100644 index 00000000000..cb1d9c0a986 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc0d_water/result.ref @@ -0,0 +1,5 @@ +etotref -466.5443393467753 +etotperatomref -155.5147797823 +totalforceref 4.003482 +epccref 0.0033441637 +totaltimeref 1.44 diff --git a/tests/02_NAO_Gamma/scf_pcc2d_water/INPUT b/tests/02_NAO_Gamma/scf_pcc2d_water/INPUT new file mode 100644 index 00000000000..59ee882807c --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc2d_water/INPUT @@ -0,0 +1,21 @@ +INPUT_PARAMETERS +suffix autotest +calculation scf +basis_type lcao +ks_solver genelpa +pseudo_dir ../../PP_ORB +ntype 2 +nbands 8 +ecutwfc 20 +scf_thr 1e-7 +scf_nmax 200 +mixing_type pulay +mixing_beta 0.2 +smearing_method fixed +symmetry 0 +cal_force 1 +init_wfc atomic +assume_isolated pcc_2d +orbital_dir ../../PP_ORB +gamma_only 1 +pcc_2d_axis 1 diff --git a/tests/02_NAO_Gamma/scf_pcc2d_water/KPT b/tests/02_NAO_Gamma/scf_pcc2d_water/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc2d_water/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/02_NAO_Gamma/scf_pcc2d_water/README.md b/tests/02_NAO_Gamma/scf_pcc2d_water/README.md new file mode 100644 index 00000000000..7924bf4081c --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc2d_water/README.md @@ -0,0 +1,6 @@ +Standalone vacuum PCC 2D regression. Compare total energy, ionic forces +and E_pcc against a two-rank reference with one OpenMP thread. The PW case +checks a charged periodic layer with open axis 2; the LCAO case checks a +neutral water layer with open axis 1. Analytic sheet/dipole, normal-force and +electronic-potential derivative tests provide independent kernel references. +The saved original PCC executable was also used for numerical comparison. diff --git a/tests/02_NAO_Gamma/scf_pcc2d_water/STRU b/tests/02_NAO_Gamma/scf_pcc2d_water/STRU new file mode 100644 index 00000000000..05933dd2b32 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc2d_water/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +O 15.999 O_ONCV_PBE-1.0.upf +H 1.008 H_ONCV_PBE-1.0.upf +NUMERICAL_ORBITAL +O_gga_7au_100Ry_2s2p1d.orb +H_gga_8au_100Ry_2s1p.orb + +LATTICE_CONSTANT +1.8897261258369282 +LATTICE_VECTORS +12.0 0 0 +0 12.0 0 +0 0 12.0 +ATOMIC_POSITIONS +Cartesian +O +0 +1 +6.0 6.0 6.0 1 1 1 +H +0 +2 +6.75716000000000 6.58626 6.0 1 1 1 +5.24284 6.58626 6.0 1 1 1 diff --git a/tests/02_NAO_Gamma/scf_pcc2d_water/result.ref b/tests/02_NAO_Gamma/scf_pcc2d_water/result.ref new file mode 100644 index 00000000000..11f6923b1f4 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_pcc2d_water/result.ref @@ -0,0 +1,5 @@ +etotref -466.5376820791657 +etotperatomref -155.5125606931 +totalforceref 3.987612 +epccref 0.0099396021 +totaltimeref 1.45 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 94b4f526119..f4f80ac113b 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -104,6 +104,7 @@ out_dm=$(get_input_key_value "out_dm" "INPUT") out_mul=$(get_input_key_value "out_mul" "INPUT") gamma_only=$(get_input_key_value "gamma_only" "INPUT") imp_sol=$(get_input_key_value "imp_sol" "INPUT") +assume_isolated=$(get_input_key_value "assume_isolated" "INPUT") run_rpa=$(get_input_key_value "rpa" "INPUT") out_pot=$(get_input_key_value "out_pot" "INPUT") out_elf=$(get_input_key_value "out_elf" "INPUT") @@ -810,13 +811,21 @@ fi #-------------------------------------------- # implicit solvation model #-------------------------------------------- -if ! test -z "$imp_sol" && [ $imp_sol == 1 ]; then +if ! test -z "$imp_sol" && { [ "$imp_sol" = 1 ] || [ "$imp_sol" = 2 ]; }; then esol_el=`grep E_sol_el $running_path | awk '{print $3}'` esol_cav=`grep E_sol_cav $running_path | awk '{print $3}'` echo "esolelref $esol_el" >>$1 echo "esolcavref $esol_cav" >>$1 fi +#-------------------------------------------- +# vacuum point-countercharge correction +#-------------------------------------------- +if [[ "$assume_isolated" == "pcc_0d" || "$assume_isolated" == "pcc_2d" ]]; then + epcc=$(grep E_pcc "$running_path" | tail -1 | awk '{print $3}') + echo "epccref $epcc" >> "$1" +fi + #-------------------------------------------- # random phase approximation #--------------------------------------------