Sensitivity of SCCS convergence to RHO_MIN / RHO_MAX and Poisson solver choice

[Mridhula]

Hello CP2K community,

I am having difficulty identifying stable and physically reasonable values for RHO_MIN and RHO_MAX in the SCCS implicit solvation model, and I would greatly appreciate any guidance.

My system is a highly negatively charged molecular system (total charge -6) treated under 0D (non-periodic) conditions using the MT Poisson solver, with BLYP as an initial functional. For context, the system is a 3GC nucleotide with total charge -6; the neutralized version includes six Na⁺ counterions. Due to the large net charge, the SCCS response appears to be extremely sensitive to the cavity definition.

What I tried initially

After reading a Google Groups discussion reporting SCCS parameters for OH⁻, where a relatively narrow cavity with higher density thresholds was used, I explored RHO_MIN / RHO_MAX values scaled around the OH⁻ reference point (RHO_MIN = 0.0024, RHO_MAX = 0.0155), approximately preserving their ratio. In practice, this corresponded to scanning values in the ranges RHO_MIN ≈ 0.0018-0.0042 and RHO_MAX ≈ 0.0116-0.0271 (a.u.).

However, these choices consistently led to severe convergence issues:

• the SCCS cycle did not converge even after 2000-3000 iterations,
• large oscillations and occasional energy blow-ups were observed,
• SCF iterations became extremely slow (up to ~200 s per iteration), even with BLYP.

To improve stability, I also tried adjusting SCCS-related settings, including:

• decreasing the SCCS MIXING parameter to 0.02,
• setting EPS_SCCS = 1.0E-06 so that SCCS iterations begin only after SCF convergence.

I then hypothesized that, due to the large -6 charge, placing the solvation cavity too close to the solute leads to an overly reactive solvent response that destabilizes the SCF/SCCS coupling. Based on this, I explored cavity definitions at lower electronic densities, shifting the dielectric transition further into the low-density tail of the electronic density. Specifically, I tested RHO_MIN / RHO_MAX values in the ranges 3.0 × 10⁻⁴-1.0 × 10⁻³ and 1.0 × 10⁻³-3.0 × 10⁻³ (a.u.), respectively.

While these values improved stability in some cases, I still encountered SCCS convergence issues. To check whether this behavior was specific to the highly charged system, I also neutralized the system by adding counterions (6 Na⁺), but the neutral system exhibited similar convergence problems.

Additional observation

I also noticed a strong sensitivity to the Poisson solver. For identical SCCS parameters, the number of SCCS iterations and overall stability differed significantly when switching between the MT and analytic Poisson solvers, with the analytic solver sometimes behaving more robustly for very low-density cavity definitions. Orbital transformation also seemed much more robust than using standard diagonalisation.

   

The first image shows the SCCS results obtained with the analytic Poisson solver, whereas the second image shows the corresponding results using the MT Poisson solver; in both cases, orbital transformation was employed.

The above image corresponds to a simulation using standard diagonalization together with the MT Poisson solver. As shown in the screenshots, there is significant sensitivity to these solver and SCF choices even when all other settings are kept the same. I have also attached the input file corresponding to this setup using the analytic Poisson solver.

Questions

I would appreciate any advice on the following points:

1. For highly charged molecular systems, is it generally advisable to shift the SCCS dielectric transition further into the low-density tail, even if this leads to a relatively loose cavity?
2. Are there recommended ballpark ranges for RHO_MIN / RHO_MAX for large anionic systems beyond small test cases like OH⁻?
3. Are there additional diagnostics or best practices to assess whether a given cavity definition is physically reasonable and numerically stable?

Thank you very much for your time and help.