[LKFEIYI]

Linked Issue

Fourth part of #8052. It follows #8075 (PCC), #8162 (SCCS Step 1: periodic SCCS) and #8194 (SCCS Step 2: SCCS with PCC 0D/2D).

Remaining follow-up PRs, each stacked on the previous one:

1. SCF integration (delayed start, fixed-cell relax, out_sol).
2. The solvent-aware cavity.

This PR also answers the open review thread on #8194 about the accuracy of SCCS + PCC forces (see D1 and the force table below).

Summary

This PR adds two optional cavity models from QE-Environ, both off by default:

• A derivative lowpass filter, for PCC systems. It makes the solvation forces consistent with the energy.
• A full cavity. It adds core Gaussians to the cavity density, which removes the dielectric pockets that some pseudopotentials produce inside atoms. (How to find: when using APNS PP, some cases cannot converge)

Both follow Environ's definitions, so the energies agree with QE-Environ: to within 0.6 meV with the filter, and to 0.01 meV in full mode. Where ABACUS deliberately differs, the forces become derivatives of the energy, and the atom selection becomes robust. The differences and their advantages are listed in D1–D2.

What's changed?

• sccs_lowpass_p1 / sccs_lowpass_p2.
  • Corresponds to Environ's deriv_lowpass_p1/p2, and is available with assume_isolated pcc_0d / pcc_2d only.
  • Before every FFT derivative of the switching function s, its Fourier coefficients are multiplied by f(G) = erfc(p1 G²/G²cut − p2)/2. Gcut is the density cutoff, as in Environ core_fft_lowpass.f90.
  • The cavity potential is the exact derivative of the filtered, discrete sqrt-CG energy: dE/ds = L/2 (q v + lapl b + L div(b grad s)), with L = ln ε_bulk and b = ε v²/8π.
• sccs_solvent_mode full with sccs_corespread.
  • Corresponds to Environ's solvent_mode='full'.
  • One Gaussian per atom (charge zv, width sccs_corespread, default 0.5 bohr) is added to the density that defines the cavity, and only to that density. Electrostatics are unchanged.

Differences from QE-Environ

D1. With the lowpass, the cavity potential is the derivative of the energy

What Environ does.

• With deriv_lowpass_p1/p2 > 0, setup.f90 replaces the whole FFT derivative core with core_fft_lowpass. Every gradient, Laplacian and Hessian is filtered, including the gradient of the potential and the ionic force (force_fft_lowpass).
• The cavity potential calc_dedielectric_dboundary stays the continuum formula, applied to the filtered ∇̃φ.
• The energy is ½∫ρφ (calc_eelectrostatic).

What we do. It filters only the derivatives of s, which define the discrete operator, and differentiates that energy exactly.

• No adjoint solve is needed: the operator is symmetric, so dE = −vᵀ dA v / 2.
• The extra cost is four forward FFTs and one inverse FFT per SCCS evaluation, small next to the sqrt-CG solve.

Advantages.

1. Forces are the derivative of the energy. That's what relaxation, NEB and MD need. The table shows the solvation part of analytic − FD, in eV/Å:

   System	ABACUS lowpass	QE lowpass
   OH⁻ (anion, small cavity)	−1.0e-3 (step → 0)	+1.0e-2
   C₂H₅S⁻ (anion, full)	−2.7e-4	+7.5e-4
   H₃O⁺ (cation), PCC0D	+9e-4	– (total force: 3.6e-2)

An analytic gradient of the solved potential, ∇v = (∇w − w∇ln ε/2)/√ε, avoids the FFT ringing at the PCC cell boundary, but it is not the gradient the discrete FFT energy implies: the SCF potential is then not variational (energy ~2 meV too high, H₃O⁺ FD slope off by 0.045 eV/Å from QE, analytic − FD 2.3–2.6e-2 eV/Å), so we keep Environ’s FFT gradient and obtain consistent forces from the exact derivative of the lowpass-filtered energy instead.

D2. The lowpass is accepted only with pcc_0d / pcc_2d

• Periodic cells use the density chain rule from Step 1 (Environ deriv_method chain). It takes no FFT derivative of s, so the filter would act on nothing.
• The periodic forces are already consistent: −1.1e-3 eV/Å for H₂O.

Validation

Codes. ABACUS: this layer's code, measured in the stacked SCCS build. The cases use no feature of the later PRs. QE 7.6 + Environ 3.1.1. Both with the same pseudopotentials and geometries, 100/400 Ry, PCC0D unless stated otherwise.

FD protocol. A zero-sum displacement (two atoms moved ±0.01 Å in opposite directions along y) cancels the net-force subtraction that both codes apply to printed forces.

Force consistency (eV/Å). The numbers are the solvation part of analytic minus finite difference: (analytic − FD) of the solvent run minus that of vacuum + PCC.

System (preset, cavity)	ABACUS default	ABACUS lowpass 10/5	QE default	QE lowpass 10/5
H₃O⁺ (cation, electronic), PCC0D, 20 bohr	+1.7e-2	+9e-4	–	–
H₃O⁺ (cation, electronic), PCC2D, 12×20×12 bohr	−4.6e-3	−2.0e-3	–	–
OH⁻ (anion, electronic), 20 bohr	+2.1e-1	+1.4e-2 → −1.0e-3 (step → 0)	+1.4e-1	+1.0e-2
C₂H₅S⁻ (anion, full), 30 bohr	−2.0e-2	−2.7e-4	−2.2e-2	+7.5e-4

• Without the filter, the error is in both codes. The continuum cavity potential misses the energy slope in both, and most for a small anion.
• The energies are the same. The FD slopes of the two codes agree to 6e-4 (OH⁻) and 7e-6 eV/Å (C₂H₅S⁻).

Lowpass energies against QE-Environ.

System	ΔG ABACUS (eV)	ΔG QE (eV)	ΔΔG (meV)
H₃O⁺, PCC0D	−4.65902	−4.65848	−0.54
H₃O⁺, PCC2D	−1.61126	−1.61067	−0.59
OH⁻, PCC0D	−4.20709	−4.20685	−0.25
C₂H₅S⁻ (full), PCC0D	−3.18695	−3.18693	−0.01

The ABACUS value is the lower one in every case (D1). The filter changes the model energy itself, by 14–16 meV for H₃O⁺ and 6 meV for OH⁻. It therefore stays off by default, because the water presets were fitted without it.

How the comparison was run.

• Each case directory holds an ABACUS INPUT/STRU/KPT and a QE pw.in + environ.in, both written by one generator script.
• OMP_NUM_THREADS=1:
  • ABACUS: mpirun -np 8 abacus_std_para
  • QE: mpirun -np 8 pw.x -environ -in pw.in
• Energies are taken from !FINAL_ETOT_IS (ABACUS) and ! total energy (QE). Forces come from TOTAL-FORCE and "Forces acting on atoms".

D3. The lowpass is accepted only with pcc_0d / pcc_2d

• Periodic cells use the density chain rule from Step 1 (Environ deriv_method chain). It takes no FFT derivative of s, so the filter would act on nothing.
• The periodic forces are already consistent: −1.1e-3 eV/Å for H₂O (Feature: periodic self-consistent continuum solvation (SCCS) Step2: SCCS with PCC 0D/2D #8194).

The discretization error of chain at the cavity edge is proportional to the potential there. This is why it works well for neutral periodic systems (it converges faster in neutral systems without PCC) and fails for charged systems with PCC.

Situation	chain	fft
Neutral, periodic	Good: energy and forces equivalent to fft, faster surface convergence	Good, slightly slower surface convergence
Charged, PCC (open boundary)	Poor: force error about 8×10⁻² eV/Å, does not converge with the grid; energy convergence is actually better	Force error about 1.7×10⁻² eV/Å; about 10⁻³ with the lowpass exact derivative
Exact derivative / lowpass	Not possible	Possible (ABACUS's D1)
Surface term at a density saddle (unregularized)	Diverges	Diverges
SCF, charged + SA	Slow (127 steps in QE)	Faster (31 steps in QE)
SA: solvent-aware (in the following PRs)

Our choice

Path	Dielectric	Surface
Periodic, no SA	chain	spectral gradient
Periodic, with SA	chain	exact discrete chain surface
PCC	fft; exact derivative when the lowpass is on	spectral gradient; exact discrete chain surface with SA

D4. Which atoms get a core Gaussian in full mode

What Environ does. iontype.f90 compares the species label with 'H' and sets the width to 1e-10 for that species. Every other species gets the given width, or 0.5 if the given value is ≤ 0.

What we do.

• With a single sccs_corespread value, atoms whose atomic number equals their pseudopotential valence charge are skipped. Those atoms have no core electrons.
• The atomic number comes from the pseudopotential element.
• A list of nat values gives per-atom widths, and a value ≤ 0 disables that atom.
• If the element is not recognized, the single width is applied, with a warning.

Governance Notes

• INPUT/docs changes.
  • New keys:
    • sccs_lowpass_p1, sccs_lowpass_p2: both positive or both non-positive; PCC only.
    • sccs_solvent_mode: electronic or full.
    • sccs_corespread: one value or nat values. A per-atom list lowers symmetry 1 to 0 (reset hook of symmetry in read_inp_sys.cpp).
  • docs/parameters.yaml and input-main.md are generated in CI from the C++ help, so they are not edited here.
• Core module impact.
  • source_cell: unitcell::AtomData gains atomic_number, filled by the new pseudo_atom_z(const Atom&) (0 if the element is not recognized).
  • PotSccs gains check_sccs_cell. It stops the run if the sccs_corespread count is neither 1 nor nat. It warns when an unknown element prevents the automatic core exclusion.
  • read_inp_sccs.cpp is split: the parse and validation helpers move to sccs_inp_check.cpp.
  • No changes to ESolver, HSolver, Operator or Psi.
• Exceptions requested: none. The governance check reports 0 errors and 5 warnings:
  • 4 header-dependency warnings: vector3.h in sccs_cav_deriv.h, and <vector> in sccs_lowpass.h, sccs_parameters.h and read_inp_sccs.h. Each declares a member or parameter of that type.
  • 1 documentation-sync warning: the docs are generated in CI, see above.

[ABACUS-PR-BOT]

Static review of 5f305c8f7feda99e5eecb03f739658a7362c68dc: one correctness issue and two non-blocking follow-ups below. No builds, runtime tests, MPI runs or CI diagnostics were executed.

INPUT metadata is updated; the YAML/Markdown references are intentionally untracked and regenerated by the existing documentation build, so no tracked generated-file update is expected. Documentation generation was not run.

[Critsium-xy]

Static review (no build/run).

INPUT docs: docs/parameters.yaml and input-main.md are generated from the Input_Item descriptions, so no manual update is needed. Please state this in the PR body to resolve the governance checker's INPUT-doc warning.

[Critsium-xy]

Thanks a lot for your contribution!