solvation free energy with couple-intramol set to no
proposed fix for redmine bug #3403
Fixes #3403 by preventing the free energy kernel from skipping pairs
that are beyond the cutoff but also in the exclusion list generated by
couple-intramol=no. Forces calculation of the reciprocal-space Ewald
component for those pairs in the free energy kernel. Also extends the
size of the relevant tables to rcoulomb+tabext to avoid getting nonsense
energies/forces beyond the cutoff. Molecules being decoupled in this way
should be smaller than rcoulomb+tabext.
#1 Updated by Vytautas Gapsys 6 months ago
A possible issue when calculating solvation free energies and using the option: couple-intramol = no
The observation is that the overall electrostatic energy (and, in turn, the final free energy) depend on the rcoulomb cutoff, in spite of using PME.
Attached are the examples illustrating the issue. The system contains a solvated ligand. All, except for two charges in the system are set to 0.0. The two atoms on the molecule which is being decoupled via the .mdp settings have charges of +1 and -1. The distance between these atoms is 0.74 nm. Several scenarious to illustrate the issue are shown where overall electrostatic energy is calculated at the first MD step:
noFE_rcoul0.7: no free energy activated, rcoulomb set to 0.7 nm, energy -188.755 kJ/mol
noFE_rcoul0.8: no free energy activated, rcoulomb set to 0.8 nm, energy -188.745 kJ/mol
withFE_rcoul0.7: free energy activated, couple-intramol=no, rcoulomb 0.7 nm, energy -375.565 kJ/mol
withFE_rcoul0.8: free energy activated, couple-intramol=no, rcoulomb 0.8 nm, energy -188.745 kJ/mol
A possible explanation to what might be happening: for this simulation, intramolecular interactions are replaced by explicit pair interactions, however, this apparently is done for both states - coupled and decoupled. This way, in the coupled state those intramolecular electrostatic interactions that are outside of the rcoulomb cutoff are counted twice (when using PME): first, they are counted explicitly ignoring the rcoulomb cutoff, second time they appear in the reciprocal part of PME. This is verified by extending rcoulomb up to the point where all the intramolecular interactions are within the rcoulomb cutoff. When extending the cutoff further from that point the overall electrostatic energy becomes rcoulomb independent.
Such a dependence on rcoulomb is problematic, as tunepme can change rcoulomb value during runtime.
- Status changed from New to Accepted
This indeed seems to be an issue. The radius of the molecule should be smaller than the cut-off radius. We should document this and ideally also add a check for this.
Note that PME tuning is not an issue. This can only increase the cut-off, so never causes issues.
#4 Updated by Yuriy Khalak 5 months ago
I managed to trace the problem to the free energy kernel not applying the long range part of the potential for pairs that are in the exclusion list (as part of the molecule being decoupled) and beyond the Coulombic cutoff.Attached is a diff against gromacs 2019.4 that fixes this behavior by:
- adding long range interactions for those specific pairs
- making the relevant potential table extend to rlist + table-extension (like in the case of the group cutoff scheme) so that interactions beyond the cutoff can still be evaluated.
Do note that this moves the problem from applying to molecules larger than rcoulomb to those larger than rlist + table-extension. Ideally grompp should give a warning if this is about to occur. However, table-extension can be tuned from the mdp file without affecting performance by much.
For a system of 6593 atoms, 40 of which are being decoupled, this fix leads to a performance loss of ~5%, so there is room for optimization.
Hope this helps.
#6 Updated by Yuriy Khalak 5 months ago
I meant the need to subtract (not add) the "reciprocal-space Ewald component" to interactions within the decoupled molecule that are beyond the cutoff. This is already done for all interactions within cutoff here: https://github.com/gromacs/gromacs/blob/328a18d71dda42ca67edf76b2f93781dab6fdf9d/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp#L783-L820
As those interactions are on the exclusion list due to couple-intramol = no, they ignore cutoff in https://github.com/gromacs/gromacs/blob/328a18d71dda42ca67edf76b2f93781dab6fdf9d/src/gromacs/listed-forces/pairs.cpp#L96-L143 and get a 1/r potential there.
So what ends up happening to short range (PME LR not considered) interactions is:
#7 Updated by Vytautas Gapsys 5 months ago
Berk Hess wrote:
Note that table-extension does nothing in release-2020 which only supports the Verlet cutoff-scheme.
Could you explain a little more why table-extension does nothing?
In the attached examples I used gmx release-2020. In one case in the .mdp I set table-extension=1 (table1.tar.gz), in another case table-extension=5 (table5.tar.gz). Verlet cutoff scheme is used in both cases. For the shorter table extension I get warnings like this one:
"WARNING: Listed nonbonded interaction between particles 99 and 361
at distance 2.517 which is larger than the table limit 2.441 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions."
For the longer table extension, the warnings are not present. Would this mean that table extension makes a difference?