Quick follow up relate to this, @hess. Am I correct in the code analysis that when computing the dispersion correction with potental-shift, you compute the shift value, and then integrate it out from r=0 to r_cutoff assuming that g(r)=1 inside the cutoff? Then for computing the dispersion correction outside the cutoff, it’s the same as an abrupt cutoff/potential switching (with potential switching adding in the contribution in the switch region). I’m looking at lines ~415-500 of dispersioncorrection.cpp.
The assumption that g(r)=1 inside the cutoff, if i understand the code correctly, is a pretty big approximation, though the contributions from overestimating and underestimating g(r) in different regions actually probably do cancel out a decent amount.
This does bring me back, I remember my first contribution to GROMACS as a graduate student was generalizing the dispersion cutoff for potential-switching . . .
For the shift value it doesn’t matter what g(r) inside the cut-off is, as long as it is 1 on average. This will be the case for any dense system, but I realize now that this is not the case if you have a system with partial vacuum or two phases.
I’m not entirely certain about this. The average energy of a pairwise function can be written as an integral of the pair potential = Nrho/2 \int 4pi r^2 U(r) g(r) dr. What I’m seeing (which I could be misreading) is that the code is calculating \int 4pi r^2 (E_c), where E_c is the energy shift due to both repulsion and dispersion. The true correction would be \int 4pi r^2 (U(r)_shift + E_c) g(r) \approx \int 4pi r^2 (U(r)_shift) g(r) dr + \int 4pi r^2 E(c) g(r) dr, and then the g(r)=1 approximation is being applied for the correction. Since g(r) dips up and down, then it would mostly balance out. For it to exactly balance out in homogeneous systems, I think you are saying that we have to have \int_0^r_c 4pi r^2 g(r) = N_inside rc, and that as long as r_c is out at the point the g(r) \approx = 1, then N_inside/(4 pi r_c^2) = 1?
What we have seen is that this correction results a relatively small change in /N with cutoff distance (0.1-0.2 kJ) for the shift cutoff modifier, compared to a completely flat variation with the switch cutoff modifier, so that supports that assuming that \int 4pi^2 g(r)/V_c \approx 1. This compares to the free energies, where shift and switch differ by a couple kJ/mol for similar changes in cutoff (switch being cutoff independent, shift not cutoff dependent), and we are still trying to understand why it does this if the dispersion correction (+ shift correction) for the shift modifier is reasonable - maybe there’s an issue with the lambda dependence here for the shift?
A agree with the first paragraph. I now also realized that the whole dispersion correction approach anyhow assumes homogeneous particle distribution beyond the cut-off. The question is then if there is a bias in the distribution within the cut-off. This can very well be, as certain atom types prefer to interact with certain other atom types and also only certain atom types are bonded to each other.
This should be easy to check numerically for a few systems of interest.