Free energy calculations, soft-core potential
Current Gromacs implementation of the soft-core potential for free energy calculations implements a version described by Beutler et al, 1994. Potential of this form has been shown to suffer from undesired minima along the alchemical path when van der Waals and electrostatic interactions are switched simultaneously. For certain free energy calculation protocols, e.g. discrete TI, BAR, a workaround of separate switching of the electrostatic and van der Waals interactions has become a standard practice. For the single topology alchemical methods, however, such an approach introduces the requirement of generating several separate topologies. Also, for different free energy calculation protocols, e.g. non-equilibrium methods, performing non-bonded interaction switching simultaneously is more practical.
To alleviate the problems of the current soft-core potential and allow for a simultaneous Coulombic and van der Waals switching, an alternative formulation of the potential was proposed https://dx.doi.org/10.1021/ct300220p
The new soft-core version sets the force to a finite value between the alchemically transformed particles at short distances. In contrast to the Beutler et al, 1994 soft-core, this finite, but small force allows avoiding the unwanted minima along the transition path. The new soft-core version has been implemented in Gromacs 4.6.2 and extensively tested by calculating free energy changes for the amino acid, nucleic acid and ligand modifications.
The current aim is to incorporate the new soft-core version into the Gromacs master branch. The new function will not replace the current soft-core implementation, but allow the user to choose between the soft-core versions.
People involved in implementation: Luka Stanisic, Vytautas Gapsys, Carsten Kutzner, Markus Rampp, Bert de Groot
#4 Updated by Mark Abraham 12 months ago
Sounds like a good thing to do. Currently we have a few other things on the roadmap before we'd be keen on new code here (e.g. add verlet table support, kill the group scheme, port the perturbed-interaction kernel to use verlet-scheme neighbor list data structures, remove vestiges of the group scheme), but I'm sure we can have a good discussion e.g. in Gottingen in October.
Could you please provide the formulas for this new softcore potentiaL? We need to discuss those first before we even start looking at code. From Vytas I understood there were several parameters that needed to be optimized, but I can not find any in the code (and the are no comments and there is no documentation).