I am a bit perplexed to see substantive differences in forces (as taken from OUTCAR with ASE) computed when I take AIMD frames and run static electronic relaxations on each one (no additional ionic steps, NSW=0, EDIFFG does not matter). All conditions are the same: EDIFF=1e-4, ALGO=VeryFast. All other parameters are identical. Electronic convergence is reached in both cases. The calculation is collinear spin-polarized, but all magnetic moments are rigorously zero (only Ca, O, and H atoms are in the calculation). VASP 6.3.0, r2SCAN+rVV10, 2x2x2 k-points, again everything identical.
This occurs independently of the thermostat (Langevin / npt and Nosé / nvt). Changing POTIM for the electronic relaxation does not matter. Changing the algorithm introduces only minute differences and does not change the overall picture.
Here is an example plot of the difference in forces (30 frames, 66 atoms each, NVT / Nosé):
Here is another example (also 30 frames, 66 atoms, but NPT / Langevin and double the POTIM):
I can probably live with the bulk of it when the differences are below 5 meV/angstrom. But the outliers towards 20 meV/angstrom and above are scaring me. Context: the forces are fed to training ML potentials.
What could be the reason(s) for this discrepancy in forces?