Skip to content

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

Closed
LKFEIYI wants to merge 47 commits into
deepmodeling:developfrom
LKFEIYI:codex/sccs-basic
Closed

LKFEIYI wants to merge 47 commits into
deepmodeling:developfrom
LKFEIYI:codex/sccs-basic

Conversation

@LKFEIYI

@LKFEIYI LKFEIYI commented Oct 9, 2026

Copy link
Copy Markdown

Reminder

  • I have read AGENTS.md and docs/developers_guide/agent_governance.md.
  • I have linked an issue or explained why this PR does not need one.
  • I have added adequate unit tests and/or case tests, or explained why not.
  • I have listed the exact verification commands run and their results.
  • I have described user-visible behavior changes, including INPUT parameter changes.
  • I have explained core-module impact for ESolver, HSolver, ElecState, Hamilt, Operator, Psi, or other source/ changes.
  • I have requested any needed governance exception below.

Linked Issue

Second part of #8052, split according to the review feedback. Stacked on #8075 (PCC): until #8075 is merged, the diff below also contains its commits.

Follow-up PRs, each stacked on the previous one: PCC coupling → lowpass and full cavity → SCF integration → solvent-aware cavity.

What's changed?

imp_sol 2 selects the self-consistent continuum solvation model (SCCS; Andreussi, Dabo and Marzari, J. Chem. Phys. 136, 064102 (2012)) for neutral, periodic cells, with calculation scf, energies and forces, PW and LCAO, nspin 1 and 2. The original ABACUS solvent model stays imp_sol 1, and true/false are still read as 1/0.

  • Cavity and dielectric. ε(ρ) uses the Andreussi switching function between sccs_rho_min and sccs_rho_max. sccs_preset provides Environ's vacuum, water-neutral, water-cation and water-anion parameters; custom uses sccs_epsilon, sccs_rho_min/max, sccs_gamma and sccs_pressure.
  • Electrostatics. The generalized Poisson equation is solved with Environ's square-root preconditioned CG. The solute charge is the electron density plus Gaussian ions with a fixed spread of 0.5 bohr.
  • Non-electrostatic terms. Surface tension γS and pressure PV.
  • Forces. The ionic reaction force of the Gaussian ions.

New code: source/source_hamilt/module_sccs/ (cavity, presets, periodic Coulomb operator, Gaussian ions and forces, sqrt-CG response, energy functional), the PotSccs potential component in source_estate/module_pot/, and read_inp_sccs.{h,cpp} for the INPUT keys.

Deliberate differences from Environ (details and evidence in the attached table):

  • The surface potential is the derivative of the regularized discrete surface, -div(grad s / |grad s|_eta), instead of Environ's unregularized curvature formula. With a density saddle on a grid point (He2), QE-Environ's surface potential reaches 100-3000 times its median and 4 of 6 runs do not converge in 200 steps; ABACUS converges in 5.
  • The sqrt-CG stops when both the RMS (normalized by the grid size) and the maximum of the charge residual meet sccs_tol_rms and sccs_tol_max. Environ stops on the unnormalized sum(r^2), so its accuracy depends on the grid size.
  • A previous sqrt-CG solution is used as the initial guess only if it lowers the initial residual.

Governance Notes

  • INPUT/docs changes:
    • imp_sol becomes an Integer (0/1/2); true/false remain valid.
    • New keys: sccs_preset, sccs_epsilon, sccs_rho_min, sccs_rho_max, sccs_gamma, sccs_pressure, sccs_maxiter, sccs_tol_rms, sccs_tol_max, sccs_surface_eta. They are listed right after the original solvent keys.
    • imp_sol 2 is rejected at input for anything other than CPU KS-DFT scf with PW/LCAO, nspin 1/2, assume_isolated none, and without stress, electric/gate fields, DFT-1/2, DeePKS or dm_to_rho. Later PRs may extend this list.
  • Core module impact:
    • PotBase gains three virtual hooks with empty defaults: get_solvation_energy, add_solvation_force and solvent_electrostatic_potential. Potential forwards them to its components. No existing component overrides them, so other calculations are unchanged.
    • Potential::pot_register / get_pot_type take const Input_para& explicitly; the old one-argument pot_register forwards PARAM.inp. get_pot_type now reads gga_grad and dip_cor_flag from that argument instead of PARAM.inp.
    • ElecState::cal_energies takes imp_sol as int; imp_sol 2 reads the energies from the potential.
    • PW (force_pw.cpp) and LCAO (force_stress_lcao.cpp) add the solvation force through Potential::add_solvation_force; imp_sol 1 keeps its own path.

LKFEIYI and others added 30 commits October 4, 2026 09:47
Build: cmake --build build_pcc0d --target cell MODULE_CELL_cell_geometry MODULE_CELL_cell_tools --parallel 4 passed. CTest: both focused targets passed (10 tests), OMP_NUM_THREADS=1 outside sandbox. C++11 syntax checks passed for both changed production sources. Governance: no blockers; geometry header needs Vector3 and array for value members, string for validation diagnostics, and vector for input containers. No INPUT or user-visible behavior changes in this foundation step.
Place charge moments and PCC formulas beside makov_payne; geometry and atom extraction remain in source_cell. No INPUT behavior changes yet. Built both focused targets and passed 11 tests outside sandbox with OMP_NUM_THREADS=1, including energy/potential and ionic-force finite differences. C++11 syntax checks passed. Governance has no blockers; includes define the value moment type, diagnostic string and kernel interface.
Register standalone PCC for CPU KS-DFT PW/LCAO, add its energy once, and expose force and electrostatic output from the same component result. Generated INPUT YAML/Markdown from the new CLI. Full abacus_basic_para build and focused potential, energy, INPUT and grid tests passed. Neutral water, charged H3O+, LCAO and two-rank MPI SCF runs converged; old H3O+ energy/forces match at printed precision. C++11 syntax and staged governance checks passed without blockers. Remaining PARAM read is the existing PW registration dependency (net zero); header includes provide value member types, base class and standard containers/diagnostics.
Keep potential update orchestration distinct from geometry/atom preparation and distributed density integration. Rebuilt abacus_basic_para, MODULE_ESTATE_pot_pcc and MODULE_ESTATE_potentials_new; both CTest targets passed outside sandbox with OMP_NUM_THREADS=1. Quality score improves from 82 to 90. No INPUT behavior changes.
Register neutral water and charged H3O+ PW cases plus neutral water LCAO.
Collect E_pcc in eV alongside existing total-energy and force properties.
The new references use two MPI ranks; existing references are unchanged.

Verification:
- cmake --build build_pcc0d --target abacus_basic_para --parallel 2: passed.
- OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 Autotest.sh, -n 2 -o 1:
  PW 220_PW_PCC0D_WATER / 221_PW_PCC0D_H3O and LCAO
  scf_pcc0d_water / existing scf_solvation passed default thresholds.
- Nine targeted CTest executables passed, including INPUT and PotPcc.
- Staged governance and git diff --cached --check: passed.
Drop the newly added 220_PW_PCC0D_WATER case and CASES_CPU entry.
This case took about 48 seconds and overlaps neutral LCAO plus analytic
and host coverage. Retain charged H3O+ PW and neutral water LCAO so both
solver paths remain covered. Existing tests and thresholds are unchanged.

Verification (OMP_NUM_THREADS=1, OPENBLAS_NUM_THREADS=1, outside sandbox):
- Autotest.sh -n 2 -o 1 -j 1 -r '^221_PW_PCC0D_H3O$': passed,
  10.826 s; energy, force and E_pcc checked at default thresholds.
- Autotest.sh -n 2 -o 1 -j 1 -r '^scf_pcc0d_water$': passed,
  9.858 s; energy, force and E_pcc checked at default thresholds.
- Nine targeted CTest executables: passed, 0.21 s total.
- Staged governance and git diff --cached --check: passed.
Only test files and registration changed; no executable rebuild required.
Build: abacus_basic_para, MODULE_CELL_cell_geometry and MODULE_ESTATE_pcc_2d passed. Both focused CTest targets and explicit C++11 syntax checks passed. This stage adds no INPUT behavior; user documentation follows in the integration stage. pcc_2d.h includes the charge-moment API used by its public kernel declarations; cell geometry is forward declared.
Keep dimension and open axis immutable at construction. Reuse charge moments,
pool reductions, energy, force and electrostatic output consumers. Require
Gamma-only open-direction sampling and symmetry that preserves the slab.
Add pcc_2d_axis INPUT metadata and regenerate both parameter documents.
Report the actual failure stage and reject missing density storage via the
standard ABACUS WARNING_QUIT path. INPUT errors identify conflicting settings.

Verification:
- Full executable and affected unit targets built successfully.
- Ten targeted CTest targets passed in 0.22 s (OMP_NUM_THREADS=1).
- Explicit C++11 syntax checks passed for geometry, kernel and PotPcc.
- PW charged and LCAO neutral 2D SCFs converged on two MPI ranks;
  comparison with the saved original executable agrees within 2.7e-10 eV
  energy and 2.5e-5 eV/Angstrom force components (old run used one rank).
- Metadata generation and pcc_2d_axis CLI help passed.
- Staged governance and diff checks passed without blocking findings.

Global references: added 2, removed 6, net -4. The existing potential factory
reads the axis from PARAM and passes it explicitly into PotPcc; PW setup now
shares one const INPUT reference for existing component selection. PotPcc and
numerical helpers add no globals. The new header include owns Pcc2dParameters
by value and thus needs the complete definition.
Add only two small vacuum PCC 2D integration cases: charged H3O+ PW,
open axis 2, and neutral water LCAO, open axis 1. Generate references
with two MPI ranks and one OpenMP thread. Keep all old references and
thresholds unchanged.

Verification with OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1:
- PW Autotest.sh -n 2 -o 1 -j 2 selecting both PCC0D/PCC2D H3O+: passed.
- LCAO Autotest.sh -n 2 -o 1 -j 2 selecting both PCC0D/PCC2D water: passed.
- New tests compare total energy, force magnitude and E_pcc.
- Staged governance and git diff --cached --check: passed.
- Two-rank runtime rejection of tilted open vector and non-Gamma open
  sampling passed; kpar=2 invalid sampling also exits 1 without hanging.
Only calculation inputs and registration changed; no rebuild needed.
Restore assume_isolated and pcc_2d_axis descriptions from abacus_sccs,
removing only solvent clauses. Regenerate parameter YAML and input-main.
Remove the PCC solvent-combination test, solvent wording in the added
energy test name/readmes, and explicit imp_sol=0 from four PCC cases.
Keep the INPUT mutual-exclusion check as requested; physics is unchanged.

Verification:
- Full executable, INPUT test and energy test rebuilt successfully.
- Both affected CTest targets passed in 0.18 s with OMP_NUM_THREADS=1.
- CLI parameter help and --check-input passed for all four PCC cases.
- Script comparison confirms exact original descriptions except the
  removed solvent clauses and confirms the mutual-exclusion check remains.
- Generated docs contain 567 parameters.
- Staged governance and git diff --cached --check passed.
Compute the explicit ionic derivative at fixed electron density in Ha/Bohr using the existing PW FFT and guarded reductions. Verify analytic Fourier forces, gauge/translation invariance, six central energy differences in a nonuniform cavity, vacuum zero force, and distributed validation with output preservation. Build sccs and MODULE_HAMILT_sccs_poisson; 16 cases pass on both one and two ranks. INPUT behavior is unchanged: force/stress remain rejected until full SCF-force validation and integration.
Dispatch solvation forces through PotBase, recover the reaction field in Ha and convert Gaussian ionic forces to Ry/Bohr once. Keep legacy solvent dispatch separate in PW/LCAO. Enable SCF cal_force for imp_sol=2 while stress, relaxation and PCC remain rejected; regenerate INPUT docs. Build abacus_basic_para and focused targets; eight CTest targets and both short integration cases pass. H2 full-SCF force finite differences at 100 Ry differ by 8.21e-4 eV/Angstrom in PW and 9.02e-6 in LCAO; PW at 200 Ry differs by 5.57e-4. New/old force differences are at most 3.4e-9 eV/Angstrom. Vacuum and legacy solvent force/energy results match exactly. Preserve the inherited finite-cutoff cavity-derivative discrepancy without loosening references.
Move polarization solver parameters into the configuration header and use forward declarations for referenced SCCS types. Include complete types directly in implementations and tests.

Validation: abacus_basic_para and MODULE_HAMILT_sccs_poisson rebuilt; SCCS serial and MPI CTest entries passed (2/2); four affected headers passed standalone C++11 syntax checks; staged governance and diff whitespace checks passed with review warnings only.

Header warnings cover standard string/vector declarations needed by the APIs and owned containers. No INPUT or numerical behavior changes; user documentation updates are not required.
- Document the PCC 0D/2D kernels, their references (Andreussi and
  Marzari 2014, Dabo et al. 2008, Rutter 2021 for the charged-slab
  constant) and the ENVIRON implementation they follow in pot_pcc.h;
  cite Rutter 2021 next to the 2D gauge in pcc_2d.h.
- Define the Potential PCC accessors after the constructor/destructor,
  declare them with the other public getters, and keep the
  pcc_component() helper with the private member functions.
- Drop the stubbed-energy PCC test from elecstate_energy_test.cpp.
- Add physical checks: periodic Hartree energy + E_pcc recovers the
  isolated energy of a charged Gaussian ion and of a polar Gaussian pair
  (0D), and of a charged layer and a polar neutral slab (2D),
  independently of the cell size.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0141raGGLbA6qqg3srvRcPnz
INPUT values are validated only by the INPUT reader. Internal invariants
(cavity evaluation, CG breakdown or non-convergence, non-finite energies,
Gaussian ionic normalization, non-neutral periodic cells) stop the run with
ModuleBase::WARNING_QUIT instead of returning bool plus an error string.
The CG non-convergence message reports sccs_maxiter and the final RMS/MAX
residuals.

Also add the cheap periodic water case 220_PW_SCCS_PERIODIC_WATER
(ecutwfc 20, ecutrho 240), which reproduces the original implementation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BKqk2SNWXoZn2PRRGRzGca
The geometry, grid and moment helpers only receive values that ABACUS
has already built from STRU and INPUT, so their finite/size/storage
checks and bool+error returns could never fire.

- grid_positions and charge_moments return their result directly.
- weighted_center returns the center; make_pcc_2d_parameters returns
  the parameters.
- make_orthogonal_cell, make_slab_cell and make_pcc_0d_parameters keep a
  bool for the structure conditions a user can violate (perpendicular,
  cubic); PotPcc turns them into WARNING_QUIT messages.
- PotPcc no longer reduces replicated structure checks over the pool
  (every rank sees the same cell and stops together), and drops the
  storage, nonfinite-result and stale-result checks. validate_kpoints
  keeps its cross-pool reduction, since each pool holds its own k points.
- write_elecstat_pot drops the PCC potential size check.

Tests for the removed failure paths are removed; a death test covers
the non-cubic 0D and tilted 2D cells.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BKqk2SNWXoZn2PRRGRzGca
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BKqk2SNWXoZn2PRRGRzGca

# Conflicts:
#	source/source_estate/module_pot/potential_new.cpp
#	source/source_estate/module_pot/potential_new.h
claude and others added 17 commits October 6, 2026 17:58
- The Gaussian ionic density is built in reciprocal space, so its grid
  integral equals the valence charge to rounding; drop the
  normalization check that could never fire.
- A non-finite residual already fails convergence and stops at the next
  sqrt-CG step; drop the separate residual, net-charge and functional
  energy finiteness checks.
- validate_sccs_input reports which SCCS parameter is out of range and
  no longer checks finiteness, matching other INPUT items.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BKqk2SNWXoZn2PRRGRzGca
PotPcc takes the expected electron count (INPUT nelec) and writes a
warning when the PCC net charge differs from it by more than 1e-6 e,
e.g. after a very tight diagonalization. The correction keeps using the
grid charge, as before.
The periodic Poisson solver drops the G = 0 component of the net charge.
In a dielectric this leaves a cell-size error with no closed-form
correction, so a charged cell with epsilon above 1 now runs with a
warning instead of stopping; vacuum (epsilon 1) runs silently.

PotSccs takes the expected electron count (INPUT nelec). It stops if
the Gaussian ionic density does not integrate to the valence charge,
and warns when the grid electron count differs from nelec by more than
1e-6 e; the solute charge keeps using the grid density.
# Conflicts:
#	source/source_estate/module_pot/potential_types.cpp
The custom preset is the counterpart of Environ environ_type input, whose
env_static_permittivity defaults to 1 (no dielectric). sccs_epsilon
defaulted to 78.3 instead. The vacuum and water presets already match
Environ environ_type_setup and are unchanged.
solve_sccs_response now prepares the cavity and calls try_warm_start,
iterate_sqrt_cg and finish_response, which share the operator data in
SqrtCgOperator and the potential, residual and scratch in SqrtCgState.
The periodic water case is bit-identical.
validate_sccs_input runs the calculation-context, model and solver checks
in sccs_context_error, sccs_model_error and sccs_solver_error, so each
later key joins a short list of related checks.
Drop the symmetry and origin details from assume_isolated pcc_0d and
pcc_2d and the LCAO load-balance note from pcc_2d_axis; the run still
stops with a clear message when the cell or symmetry is unsupported.
… model

imp_sol keeps only the model choice and the unsupported features. It is
registered again where the original solvation keys are, and the SCCS keys
follow eb_k, tau, sigma_k and nc_k in INPUT and the documentation.
The exact monopole and dipole values are implied by the isolated-energy
tests, and the bilinear kernel is an internal step of pcc_0d_energy. The
finite-difference checks of each kernel, the orthogonal and slab geometry
checks, the PotPcc input rejections and the charge-moment checks now share
one test each. Every property stays checked; 31 tests become 18.
The cavity limits, thresholds and vacuum share one test, as do the
midpoint values and the finite-difference derivatives; the preset names,
the preset values (which fix the unit conversions) and the solver
tolerance checks are merged likewise. The custom-preset defaults join the
INPUT mapping test, and the solvation-force finite difference becomes a
fixture method for every cavity model.
@LKFEIYI LKFEIYI closed this Oct 9, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants