Feature: self-consistent continuum solvation (SCCS) Step1: SCCS base

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

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

source/CMakeLists.txt

@@ -729,6 +729,7 @@ target_link_libraries(
     io_basic
     io_input
     surchem
+    sccs
     neighbor_search
     domain_decomposition
     neighbor

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\

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

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<ModuleBase::Vector3<double>>& positions)
+{
+    const ModuleBase::Vector3<double> a1(lattice.e11, lattice.e12, lattice.e13);
+    const ModuleBase::Vector3<double> a2(lattice.e21, lattice.e22, lattice.e23);
+    const ModuleBase::Vector3<double> 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<double>(ix) / basis.nx;
+        const double fy = static_cast<double>(iy) / basis.ny;
+        const double fz = static_cast<double>(iz) / basis.nz;
+        positions[ir] = (a1 * fx + a2 * fy + a3 * fz) * lattice_scale;
+    }
+}
+} // namespace ModulePW