GROMACS can not reproduce properties with the GROMOS force fields
The GROMOS force fields have been parameterized with a physically incorrect multiple-time-stepping scheme using a twin-range cut-off. When using a single range cut-off (or a physically correct Trotter multiple-time-stepping scheme as introduced in version 4.5), physical properties, such as densities of liquids, can differ from the intended values. Thus GROMACS could no longer reproduce all properties of the GROMOS force fields since the GROMOS scheme was removed in version 4.5 released in 2010.
Note that certain properties or certain molecules might be affected less by this, but we will perform such checks. Therefore a warning will be added in grompp in version 2019.2 when a GROMOS force field from the GROMACS distribution is used. We will consider removing the GROMOS force fields in version 2020 to avoid generation of incorrect results.
#4 Updated by Mark Abraham 4 months ago
Berk Hess wrote:
We will consider removing the GROMOS force fields in version 2020 to avoid generation of incorrect results.
As one of several aspects of our long-running general move towards "correct, useful, and performant by default" I am generally on board with this. We should note a general policy in our documentation for force fields (or force field families) that we a) bundle, and b) intend to support. If e.g. CHARMM36 was a very large amount better than CHARMM22 in its GROMACS implementation, then we'd be considering deprecating our support/bundling of it, also.
As with other things we plan to remove, it would be advisable to have a period of deprecation first, to give the community a chance to observe that there are useful subsets of functionality in the GROMACS+GROMOS combination. Are there studies showing the efficacy of GROMACS implementations of GROMOS force fields in a) reproducing the GROMOS implementations (whether in older GROMACS or GROMOS itself), and b) reproducing experiment? I am aware of a few papers observing that the GROMACS did change its implementation (which was noted in the GROMACS release notes and for several years in the GROMACS manual). It's a pity some of the authors of such papers chose to speculate that this was done for performance reasons without consulting the GROMACS developers.
Note that if there is interest in continuing to have a version of the GROMOS scheme implemented using the GROMACS infrastructure, it is likely that the necessary pieces will emerge from the PRACE-funded NB-LIB project. What we do is all free and open-source, so people are welcome to step up and implement the things they way they want for their science!
#6 Updated by Thomas Piggot about 2 months ago
I have been working on testing and validating different GROMOS-based lipid (CKP and 54A7/53A6L) and protein (54A7) force fields in GROMACS for a while now. In particular I have been looking at the impact of a few different simulation settings on the properties of PC membranes and several different protein systems.
These simulations are on-going but I thought it important to comment on this bug based upon what I've done so far. My current results with several different setups are showing that for the GROMOS 54A7 protein force field, you can perform accurate protein simulations in GROMACS using recommended settings (e.g. with the Verlet cut-off scheme). For the two related lipid force fields, which it is worth noting were both originally tested using GROMACS and not the GROMOS software, the matter is slightly more complicated. The GROMOS-CKP lipids perform fine using the recommended setups in more recent GROMACS versions, as per the proteins. However, and as noted in some of recent published work out there, the 53A6L/54A7 lipid parameters do not. However, I believe that this is not a (current) GROMACS problem, rather an issue with these parameters and the use of sub-optimal simulation settings in their original testing and validation.
However, I should stress here that these above statements are based upon comparing back to some reference simulations performed in GROMACS with the group cut-off scheme, a single-range 1.4 nm cut-off, and atom-based charge groups (with both RF and PME). One of the reasons for using this as a reference setup was to remove any potential issues surrounding the different multiple-time-stepping algorithms. Please correct me if I am wrong, but wasn't the original argument for a twin-range cut-off simply one of speed rather than anything else? In other words, shouldn't the force field behave the same with a single-range cut-off as with the twin-range setup, if things are done appropriately?
As I mentioned, the results I've briefly discussed here are preliminary at the moment but I thought worth sharing. I definitely believe it is important to keep these force fields available within GROMACS.