Skip to content

Feature: periodic self-consistent continuum solvation (SCCS) Step1: SCCS basic function - #8162

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

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

Conversation

@LKFEIYI

@LKFEIYI LKFEIYI commented Oct 17, 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.

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 (in qe-environ, this value is seldom changed, so I just fixed this value here).
  • Non-electrostatic terms. Surface tension γS and pressure PV.

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:

  • 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. So this approach seems better.
  • 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 parameters: 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.
  • 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.

@mohanchen mohanchen added the Feature Discussed The features will be discussed first but will not be implemented soon label Oct 17, 2026
@mohanchen

Copy link
Copy Markdown
Collaborator

@ABACUS-PR-BOT Could you review the PR?

@ABACUS-PR-BOT ABACUS-PR-BOT left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Independently reviewed statically at 7cdc76650742af3d40120d8969465e46bce44172. The numerical finding is an analytic finite-grid limitation; molecular impact has not been measured. No builds, runtime/numerical jobs, or CLI/doc-generation checks were run; the required CI connector was unavailable. The trusted governance check exited 0 with warnings; git diff --check exited 2 for two added blank EOF lines.

Non-blocking (verification record): Please add the exact commands, configuration/executable identity and observed results for the SCCS unit/MPI/integration checks and any Environ comparison. The PR reminder is checked but the body has no verification section. Also record documentation-generation verification or why it was not run; those generated files are intentionally untracked.

{
const ModuleBase::Vector3<double>& gradient = polarization.gradient[i];
const double gradient_square = gradient * gradient;
response.cavity_potential[i] = -response.depsilon_drho[i] * gradient_square / (8.0 * ModuleBase::PI);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[P2] Numerical limitation: this continuum cavity derivative is not generally the derivative of the implemented finite-grid reaction energy. At uniform density inside the switching interval, the latter gives epsilon_prime/(16*pi) * (2*phi*L(phi)-L(phi^2)); on an eight-point periodic grid, a neutral cos(3x) potential and zero-mean cos(2x) density perturbation distinguish the two because phi^2 aliases. Please add a varying-dielectric density finite-difference check with grid refinement and either make the discretization variational or bound/document the approximation; tighter CG tolerances alone do not remove it. This is an analytic counterexample, not a measured error in the water case.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Confirmed. This is the same continuum derivative QE-Environ uses (−ε'|∇φ|²/8π with FFT gradients). I have tried other methods, including introducing the exact derivative, but the tests showed that some new numerical problems might be involved. And this problem has to be resolved in the low-pass parts (same as QE).

Comment thread source/source_io/module_parameter/read_inp_model.cpp Outdated
Comment thread source/source_lcao/hamilt_lcao.cpp
Comment thread source/source_lcao/force_stress_lcao.cpp Outdated
response.dsolute_drho.resize(size);
response.epsilon.resize(size);
response.depsilon_drho.resize(size);
response.grad_log_epsilon.resize(size);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: grad_log_epsilon is allocated and filled but has no reader in this revision; the coefficient uses local_gradient directly. Deferring this field would save three doubles per local FFT point and its write pass on each SCCS update. This cleanup can be a follow-up.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The field is consumed by the SCCS diagnostics in the stacked follow-up PR (SCF integration, start SCCS at a desired drho), so I will keep it here to avoid churning the field across the stack.

Comment thread source/source_io/module_parameter/read_inp_sccs.cpp
@mohanchen
mohanchen self-requested a review October 17, 2026 22:31
@LKFEIYI

LKFEIYI commented Oct 18, 2026

Copy link
Copy Markdown
Author

The SCCS in the current PR is not fully ready, as it requires PCC for charged systems and a low-pass filter for cavity forces in the following PRs.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Feature Discussed The features will be discussed first but will not be implemented soon

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants