<div>Dear developers,</div><div><br /></div><div>With r2SCAN and the default XC_DERIV (PW), the analytic forces of CP2K 2025.2 are not the derivative of the energy at a cutoff of 500 Ry. I could only test 2025.2, the version installed on our HPC system.<br /><b><br /></b></div><div><b>Minimal example (attached)</b>: one fluoroacetonitrile molecule (6 atoms) in a 12.77 A cubic cell, GTH-SCAN, DZVP-MOLOPT-SR-GTH. RUN_TYPE DEBUG with DX = 0.001 bohr, x component of the force on F:</div> numerical analytical difference<br /> (meV/A) (meV/A) (meV/A)<br /> r2SCAN, XC_DERIV PW (default) -256 -1178 -922<br /> r2SCAN, XC_DERIV SPLINE3 -200 -200 -0.2<br /> rSCAN, XC_DERIV PW -427 -427 -0.3<br /><div> r2SCAN, PW, CUTOFF 1500 Ry 28 28 -0.1</div><br /><div>In the first case the x-forces summed over all atoms give -1236 meV/A, against -323 from finite differences. In the other three cases every component checked agrees within 0.3 meV/A.<br /><br /></div><div><b>Suspected cause</b>. In src/xc/xc_libxc.F, tau is replaced by MAX(tau, sigma/(8 rho)) before the libxc call, and vrho, vsigma and vtau are then used as returned. Where the replacement is active the energy is e(rho, sigma, tau = sigma/(8 rho)), but the potential lacks the terms that come from the dependence of the replaced tau on rho and sigma. The same lines are in 2026.1 and 2026.2 versions. <br /><b><br /></b></div><div><b>Test of this reading.</b> I preloaded a libxc 7.0.0 with one change (patch attached): at grid points where tau arrives equal to sigma/(8 rho), it returns vrho - vtau*sigma/(8 rho^2), vsigma + vtau/(8 rho) and vtau = 0, i.e. the derivatives of the energy that is actually evaluated. With nothing else changed, the net force on the molecule follows the slope of the energy again when the molecule is shifted across the grid (ratio 1.01 instead of 4.50; table.txt in the attachment).</div><div><br /></div><div><b>Other settings.</b> On the same molecule, the ratio of net force to energy slope under a rigid shift (4.50 with the settings above) is unchanged with NGRIDS 1, REL_CUTOFF 200, EPS_DEFAULT 1e-14 or EPS_SCF 1e-9. It is 1.00 with XC_DERIV SPLINE2, SPLINE3, NN10_SMOOTH or SPLINE2_SMOOTH, but 3.84 with XC_SMOOTH_RHO NN10 and the default derivative. With the default derivative it is 1.37, 0.16, 1.83 and 1.15 at 600, 700, 800 and 1000 Ry, and 1.00 from 1200 Ry. In the liquid, larger basis sets (DZVP-MOLOPT-GTH, TZV2P-MOLOPT-GTH) do not reduce the error on F.<br /><br /></div><div>The same library counts the affected points. With the default derivative at 500 Ry they are 0.6% of the grid points with density between 1e-4 and 1e-2 e/bohr^3; with XC_DERIV SPLINE3 at 500 Ry, or the default at 1500 Ry, there are none. rSCAN has the same affected points, but its vtau is zero there, which is why it stays consistent.</div><div><br /></div><div><b>Why it matters</b>: in a 664-atom liquid containing this molecule, the forces on F at 500 Ry are off by 0.6-1.0 eV/A per component (against 1500 Ry) with the default, and by 0.10-0.14 eV/A with SPLINE3.<br /><br /></div><div><b>Questions</b>:</div><div>1. Does PR #6171 on master cover this case? I have not been able to test master.</div><div>2. Would you recommend a non-PW XC_DERIV for meta-GGAs? Even with a consistent potential, the PW gradient still exceeds the bound at those points. </div><div><br /></div><div>I am happy to run further tests.<br /><br /></div><div>Best, </div><div>Carlos</div>
<p></p>
-- <br />
You received this message because you are subscribed to the Google Groups "cp2k" group.<br />
To unsubscribe from this group and stop receiving emails from it, send an email to <a href="mailto:cp2k+unsubscribe@googlegroups.com">cp2k+unsubscribe@googlegroups.com</a>.<br />
To view this discussion visit <a href="https://groups.google.com/d/msgid/cp2k/5d6c25d8-0e80-4a60-ae22-961775f5b616n%40googlegroups.com?utm_medium=email&utm_source=footer">https://groups.google.com/d/msgid/cp2k/5d6c25d8-0e80-4a60-ae22-961775f5b616n%40googlegroups.com</a>.<br />