Feature #2601

Free energy calculations, soft-core potential

Added by Vytautas Gapsys over 2 years ago. Updated 11 months ago.

Target version:


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
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

equations.pdf (57.6 KB) equations.pdf Christian Blau, 02/10/2020 11:29 AM

Related issues

Related to GROMACS - Feature #742: Enhancing the performance of the free energy codeNew


#1 Updated by Carsten Kutzner over 2 years ago

  • Target version changed from 2019 to future

#2 Updated by Mark Abraham over 2 years ago

  • Description updated (diff)

#3 Updated by Mark Abraham over 2 years ago

  • Related to Feature #742: Enhancing the performance of the free energy code added

#4 Updated by Mark Abraham over 2 years 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.

#5 Updated by David van der Spoel almost 2 years ago

The free energy kernel is also left in src/gromacs/gmxlib/nonbonded/
Has it been decided whether this is going to change? That is, does it make sense to upload patches to nb_free_energy.cpp or is someone working on this already?

#6 Updated by Berk Hess almost 2 years ago

I will soon add SIMD acceleration to this kernel. This is planned as "trivial" vectorization over j-atoms, so the data structures will remain the same. To avoid conditionals I plan to use templating.
What do you want to change?

#7 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related DRAFT patchset '1' for Issue #2601.
Uploader: Luka Stanisic ()
Change-Id: gromacs~master~I703ff4ddc57cb9f22958367d7d56bc44c2297ed0
Gerrit URL:

#8 Updated by Luka Stanisic almost 2 years ago

An initial version of the changes we intend to propose have been pushed in the patch above. Although this code is not ready to be merged in the current shape, we hope that it can provide more details, and improve discussion and code development planning.

#9 Updated by Berk Hess almost 2 years ago

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).

#10 Updated by Berk Hess almost 2 years ago

It seems like the code uses the exact formulas from the linked manuscript. Then the main question is how the parameters should be chose and how sensitive they are.

#11 Updated by Michael Shirts almost 2 years ago

If we are talking about soft core changes. I'd also mention our essentially optimal pairwise soft core form; which could have some advantages:

Part 1:
Part 2:

Also available in: Atom PDF