Gromos 54A7 ff mdp neighborsearch and nonbonded parameters

User discussions
1

GROMACS version: 2019.4
GROMACS modification: No

Sorry for “trivial” question, which seems to be answered elsewhere, but I am not sure which one is correct.

Question:
Could anyone please put here neigborsearch and nonbonded parameters for Gromos 54A7 ff (for modern version of Gromacs, 2018. and newer)?*

Explanation:
I know those parameters ideally should be the same as ones with which the force field is parametrised.
After consulting literature, those should be (charge groups, twin-range cutoff-scheme, reaction field for long range electrostatics…):

cutoff-scheme = group
nstlist = 5
ns_type = grid
pbc = xyz
rlist = 0.80

coulombtype = Reaction-Field
fourierspacing = 0.18
optimize_fft = yes
rcoulomb = 1.40
epsilon-rf = 61.0

vdwtype = cut-off
rvdw = 1.40
dispcorr = no

But, modern versions of Gromacs recommends Verlet (instead of group), sigle (instead od twin-range) cutoff-scheme, and PME (instead of Reaction Field) for electrostatics, for example.

I read dosens of articles, validations of ff, tutorials… and find a lot of confusing, even contradictory statements, explanations and examples of .mdp files.

I could elaborate my confusion more, with examples, if anyone asks.

Thanks

1 Like
2

Hi,

There have been a few recent papers on this issue but I agree, it is a bit confusing at the moment. I’ve also been doing some work on it too with one bit of work looking at membranes almost finished and the other on protein simulations is getting there too. Hopefully these will add to the current literature and make things a bit clearer for people simulating recent GROMOS force fields in GROMACS.

What type of system are you simulating? You can safely use a single-range 1.4 nm cut-off with RF or PME for most simulations.

Cheers

Tom

1 Like
3

Thank you dr Piggot for your response.

I also read recent papers, which unfortunately leaves me with more questions than answers :).

I work on three simple systems: 1. peptides interaction with lipid bilayer and 2. water+TFE solutions 3. Peptide in water+TFE solution. No need to mention that your work on membranes helped me a lot, thank you!

Don’t want to take your time, but I have short question:
When you say one can use 1.4 single cut-off, does it refers to rcoulomb or rvdw, or both?

Back to my original question, this is set of parameters I came up with for now (with Gromos ff and Gromacs 2019):

cutoff-scheme = Verlet
nstlist = 5
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005 ; Default value
;rlist = (Not needed with verlet-buffer-tolerance)

coulombtype = PME
fourierspacing = 0.18
rcoulomb = 1.0 ; Probably increased automaticaly with PME

vdwtype = cut-off
rvdw = 1.40
dispcorr = no

Thanks

4

You should set both the van der Waals and Coulomic cut-offs to 1.4 nm.

I’m not sure where you’re getting the TFE parameters from but I would advise you run some tests with this cut-off setup to make sure you’re close to reported properties (e.g. for things like heats of vapourization) in either the publication of the parameters, or experimental values.

The peptides should behave fine using either RF or PME with the 1.4 nm cut-off. I’ve a whole load of simulations done to test different setups but have only had the time to do some basic analysis at the moment. That said, these setups seem to behave without problems and (with RF) are as close as you can get to the original force field settings in recent GROMACS versions.

Which exact GROMOS membrane force field/parameters are you using? Some won’t work well with the single-range Verlet 1.4 nm setup.

A final thought, and without knowing more about exactly what you are doing, is that it might be easier and computationally faster to use another force field family (e.g. CHARMM or AMBER).

Cheers

Tom

5

Hi,
I’m coming late, but I’d like to add something. I understand this is quite confusing, I’ll try to explain quickly the situation and give possible solutions.

First, the twin-range cutoff (TR), which has been classically used in the parameterization of GROMOS force field (FF), cannot be used in GROMACS starting from version 4.5 [1]. TR is a time saving approximation but one other option is to use a single-range (SR) cutoff (see below).
Second, another important aspect in GROMOS parameterization is the use of the charge-group scheme for neighbour searching which is no longer supported in recent GROMACS versions (starting from 2020).

Recently, two refs [2, 3] compared single-range (SR) and twin-range (TR) cutoff schemes on various properties (such as amino acid free energy of solvation, etc.). The authors observed minor differences between TR and SR. All their simulations were performed using the GROMOS software with the GROMOS force field (54A8). Some past tests using GROMACS (version 4.5) showed that using SR was also fine to reproduce some results on lipids (with GROMOS 53A6L) [1].
In refs [2, 3], the atom-based cutoff scheme (such as Verlet) was also compared to the group-based cutoff scheme for neighbour searching. Slightly larger differences were observed compared to TR/SR. Citing the conclusions “Thus, in codes where a GROMOS-type twin-range is not implemented users should use a group-based, single-range cutoff to obtain equivalent results”.
A recent review [4] explains in detail how these assumptions / approximations (truncation schemes, multiple time-step algorithms, etc.) work and their outcome on simulation results (see sections 2.5, 2.6, 2.7 and 2.8). There also is recommended to use a group-based cutoff scheme when using a GROMOS force field in another simulation software.

Regarding Tom’s comment, I’m not sure to which GROMOS membrane force field/parameters he is referring to, but in principle using SR vs TR should have only minor effects if the charge-group is used.

Now with that in mind, what versions of GROMACS are usable? The group-based cutoff scheme (cutoff-scheme = group) is deprecated since version 5.1 but is present up to version 2019. In version 2020 and 2021, it is no longer supported.

So here are your possible choices :

  • You can use SR cutoff with the group-based cutoff (cutoff-scheme = group) with rlist = rcoulomb = rvdw = 1.4 (plus coulombtype = Reaction-Field for electrostatics as indicated in your parameters above). This is possible up to GROMACS version 2019. It may be slow but will guarantee you to be as close as possible to the original FF as demonstrated in [1, 2].
  • If you want to use a more recent GROMACS version (2020 or 2021), you’ll have to use SR cutoff with the Verlet neighbour-searching scheme (cutoff-scheme = Verlet). But in this case, you may have to show that the properties you’re interested in are not affected compared to the charge-group cut-off scheme.

Since you were talking about using GROMACS 2019, the first solution would be probably the best solution.

I also advice you to use the thermostat / barostat settings of the original parameterization paper if you can.

Best,

Patrick

[1] Bug #1400: Can't reproduce results on DPPC with reaction field under version 4.5.3 compared to 4.0.7 - GROMACS - GROMACS development
[2] Diem and Oostenbrink, JCC, 2020, https://www.doi.org/10.1002/jcc.26426
[3] Diem and Oostenbrink, JCTC, 2020, https://www.doi.org/10.1021/acs.jctc.0c00509
[4] van Gunsteren et al., CPC, 2021, https://www.doi.org/10.1002/cphc.202000968

1 Like
6

Hi Patrick,

Thank you for adding to this discussion as I think this is important for the future of the GROMOS force fields and their inclusion in GROMACS.

What you say is correct when you look at the current literature that has been published. However, I’ve also been working on further simulations over the past few years to test some of the different GROMOS force field cut-off settings in GROMACS for both lipid and protein simulations. I’d hoped to have the first of two works (the one on lipids) published by now but with COVID, home-schooling kids, etc. I’ve unfortunately not managed to get the last 10% of it done yet.

To keep things simple, I will just focus on single-range simulations (although I do have results for some TR simulations too). As you mentioned, some of the recent work from the Oostenbrink lab, plus work from Silva et al. (https://pubs.acs.org/doi/abs/10.1021/acs.jctc.8b00758) and Gonçalves et al. (https://pubs.acs.org/doi/10.1021/acs.jctc.8b00425 ; Patrick, I know you’re an author on this one but I’ve included the link for others), have looked at a comparison of multiple-atom (MA) versus single-atom (SA) charge group cut-offs with GROMOS force fields. For properties used in the force field parameterisation, MA and SA cut-offs give fairly similar results. However, simulations of membranes (e.g. see the DMPC results in Silva et al. that weren’t really discussed in much detail) give more substantial differences.

To examine this further, I’ve performed a whole series of predominantly DPPC simulations with a few different GROMOS lipid parameter sets to try and work out what is going on. To summarise quite a lot of work in a short space is difficult, but the main result is that the SA/MA differences (seen by Silva for DMPC and reproduced in my simulations) really only arise with the GROMOS 53A6L force field parameters of Poger and Mark (i.e. the default GROMOS 54A7 lipid parameters). In particular, a higher APL close to experiment only arises with MA charge groups and RF, when using a SR cut-off (so there is no possible impact of the multiple-timestepping algorithm that was blamed in the work of Reißer et al. https://pubs.acs.org/doi/abs/10.1021%2Facs.jctc.7b00178). PME or RF with SA charge groups (using either the Verlet scheme, or done manually) give a lower APL. The same issue doesn’t arise with the GROMOS-CKP parameters and I therefore believe that this is an issue with this particular force field and the necessity of RF with MA charge groups being something that has been accidentally parameterised into the force field.

I’ve also been doing a whole series of long-timescale protein tests too. I’ve not done full and extensive analysis of these simulations yet but what I’m 99% sure is happening is a similar issue with the GROMOS 54A7 protein force field, although reversed in terms of matching experiment. When performing simulations (multiple repeats across several proteins), those with RF and MA charge groups eventually result in a loss of protein structure. PME, RF and SA cut-offs (either made manually or with the Verlet scheme) are all fine. In some regards this is a comforting result, with the more accurate methods (see the Hess et al. preprint On The Importance of Accurate Algorithms for Reliable Molecular Dynamics Simulations) maintaining protein structure.

To me, I think this makes it pretty clear that if you are wanting to simulate either a protein or a membrane (or indeed both) using the GROMOS force fields in GROMACS you should not use your first suggested setup. Rather you should use the Verlet scheme with either RF or PME (i.e. the second setup you mention). This is also nice in that these are currently maintained methods in GROMACS and is why I’d strongly argue to keep the GROMOS force fields as part of GROMACS. I’d also advise caution when using the GROMOS 53A6L parameters too.

I should stress what I said above, I’ve only tested all this for lipid and protein systems. Things will quite possibly change for other molecules. For example, in this particular case, TFE should definitely be validated as it has a low dielectric constant and the use of SA charge groups with RF will quite likely have a substantial impact on the RDF.

Cheers

Tom

7

Hi Tom,

Thanks for your message, it is indeed useful we have this discussion. Hopefully we’ll understand what is going on and find a solution so that GROMOS FFs can be used in recent GROMACS versions. Thanks also for sharing the conclusions of your tests, I’m looking forward to having a detailed read of them. What I find difficult in this story is that there are two things interlaced (and maybe more!), single-range vs twin-range on one hand and charge-group vs atomistic group cutoff (MA and SA in your notation) on the other hand. Not easy to tease apart the effect of both of them.

I comment some of your points below.

In the work from Diem and Oostenbrink it was also noted that the use of an atomistic cutoff scheme (SA in your notation) leads to an artificial structure in the solvent at the cutoff distance (though small). I’m not surprised that switching from MA to SA has more dramatic consequences on lipids.

This is interesting, but I have a hard time explaining it. The GROMOS-CKP lipids also use charge-group in their parameterization, right? If that is the case, I don’t really understand why the behaviour would be different compared to 53A6L.

I’m suprised by this result. It looks like it contradicts the paper from Diem and Oostenbrink (DOI: 10.1002/jcc.26426). Did you run much longer simulations?

Given what you said above, I understand your statement. But to agree with this conclusion, I would need to see more details in the different parameters you use. I’m notably thinking here about the thermostat and barostat. In general, the GROMOS FFs were parameterised with the weak coupling algorithm with a very tight coupling (0.1 ps) (also with a tight pressure coupling of 0.5 ps). Changing the tightness of the coupling or going to another thermostat might have important effects.

I definitely agree about the need to validate the TFE model in this case.

Best,

Patrick

8

Hi,

Yes, there are some other differences with twin and single-range that do vary with changes in GROMACS version (as has been discussed in the literature). However, as I said in my email, all of what I discussed here was for single-range simulations. In fact, to ensure nothing else was impacting on things, the neighbour list in the group cut-off simulations (both SA and MA) were updated every step. As you can image, manual single-atom charge group simulations with a single-range 1.4 nm cut-off using the group scheme and the neighbour list updated every step aren’t the fastest to run! Particularly for the lengths I’ll mention below.

With regards to your points, these are all things that I’ve either looked at or are looking at:

  1. On the differences between GROMOS-CKP and GROMOS53A6L lipids, the differences between the force fields are really small (the DPPC from GROMOS-CKP is that originally proposed by Chandrasekhar et al. https://link.springer.com/article/10.1007/s00249-002-0269-4 ; the CKP name comes from other modifications/additions to this). The main differences between the two are the CH0 atomtype for the carbonyl carbons (versus C in 53A6L) and the modified choline methyl to phosphate non-ester oxygen interactions in the 53A6L head group. Because the differences are so small, I’ve ran simulations to work out exactly what causes the charge group dependence. It turns out that it is the CH3-O LJ interactions in the head group. If you look at the size of the choline and phosphate charge groups (which are big) you can see why differences can occur as in MA charge group simulations all of these modified interactions (designed to increase the APL) are either applied or not depending on whether the charge group is included or not, while in the SA charge group simulations this is not the case.

  2. For the proteins and simulation lengths, yes. I’ve tested things on three different systems: villin head piece, lysozyme and ubiquitin. Simulation lengths respectively for all simulations of each protein are 1, 5, and 10 microseconds. I should also say that this is all in GROMACS. We (the royal we, as I didn’t run these couple of simulations) have done a few simulations of DPPC using the GROMOS 11 software too and the results seem to match. However, given how long the simulations take to run in GROMOS, I’m less confident with these results as they are shorter with less repeats. We’ve not done any protein simulations with GROMOS 11 and so this could be another source of difference. Another potential source of differences might be the fact that slightly different protein force fields (54A7 and 54A8) have been used. However, that said, on the hundreds of ns timescale you cannot see the differences between the MA and SA protein simulations and you do require these long simulations (which is maybe why it’s not been seen before).

  3. All simulations use Nose-hoover thermostat and the Parinello-Rahman barostat (coupling constants of 2 and 5 ps respectively). That said, I’ve also ran some simulations using the Berendsen weak coupling thermostat with tight coupling (0.1 ps) as I thought this might impact things. For the lipids this didn’t change anything. For the proteins (well just villin in fact), I still need to analyse these simulations but it may well be that the tight temperature coupling can stop the proteins losing structure. That said, I’d be worried if to allow the use of RF with multiple-atom charge groups you’d require such specific temperature coupling.

So I’m still pretty confident in my recommendation, that for protein and lipid simulations in GROMACS with the force fields mentioned above, you should use the Verlet scheme.

Cheers

Tom

9

Hi Tom,

thanks for your detailed answer. Sure I well noted that all you describe above is about the SR case. My comment was just to say there are multiple things coming into play.

About the three points, here are my comments :

  1. Thanks for all the details, this is very interesting. First, I guess you call these lipids CKP for Chandrasekar-Kukol-Poger, right? I re-read your paper https://doi.org/10.1021/ct3003157 which was useful for reminding me the details.
    I understand this story of big charge groups in the polar head between choline CH3 and phosphate O. Given their big size, it makes sense they behave differently between charge-group (MA) and atomistic (SA) cutoff schemes (I guess it’s an all-or-nothing behavior). By having a look at the topologies, I see the charge-groups are defined equally in both FFs. But, as you said, I also see that the CH3-OM LJ parameters are different for the C12. This leads to more repulsive interactions for 53A6L which should increase the APL more in the case of a charge-group (MA) cutoff scheme than for the atomistic (SA) one. But then, what I’m not sure of is how the CKP reaches the right APL if the choline-phosphate is more attractive, is this due to the difference in the carbonyl that you mentionned?

  2. Regarding your simulation of proteins it is the first time I hear about so long simulations with the GROMOS FF. I’m a bit surprised that RF / charge-group cutoff (MA) simulations eventually lead to a loss of protein structure, but not the atomistic (SA) ones. That would be nice to understand the origin of that. Could it be that the increase of ordering with the atomistic cutoff (SA), as observed in https://doi.org/10.1002/jcc.26426, plays a role?
    Indeed that would be very interesting to test similar long simulations with the GROMOS software.

  3. As for the thermostat, the tighter coupling that is needed is in line with the higher noise observed when using the charge-group scheme (MA) compared to the atomistic (SA) one (https://doi.org/10.1002/jcc.26426). In fact, I guess it is not the thermostat itself which matters but the tightness of the coupling.

Overall, regarding your recommendation on the cutoff scheme, after our discussion, I’d be less categorical. Of course, I understand your point regarding proteins, one doesn’t want them to unfold in very long simulations. On the other hand, the atomistic cutoff (SA) scheme comes with some artificial ordering at the cutoff distance. This has been observed for water, maybe it also occurs on lipids.
And according to what you said for 53A6L, one should use the charge-group (MA) cutoff as this is the only one that gives a higher APL.
So for now, I’d see this more as a compromise. Some other investigations, including yours, will probably shed light on what to do while we understand better what is going on.

Anyway, thanks again for this enlightening discussion.

Best,

Patrick

10

Hi Patrick,

No problem. As I said before, I’d hoped to have published this stuff by now but unfortunately haven’t managed to do so. And looking back on my last message, it came across a bit harsh regarding my point about it being only SR I was talking about! I completely agree there is lots of other stuff going on with TR setups too that makes it all quite complex. I do have a few TR simulations too that will go in the membrane paper but I won’t go into detail on them now. As to your comments:

  1. It’s actually Chandrasekar-Kukol-Piggot (excuse the vanity here!). Although Kukol did more things wrong than right, it felt appropriate to include the name as the DPPC parameters are the same as he used. The “C” also includes Chiu, given the use of the near ubiquitous united-atom lipid charges originally published 35+ years ago. As to the differences, you’re quite right, it is the CH0 atomtype that compensates. And this modification isn’t charge-group dependent, which is why the MA and SA results are the same with CKP. If I use both the CH0 atomtype and the Poger LJ parameters, I can get an even higher APL for DPPC but that is now dependent on the use of MA charge groups.

  2. For the proteins, I also agree that there may be some artefacts using the SA charge-groups with RF due to the ordering at the cut-off. As suggested to me by Chris Oostenbrink, I will be looking at this, however, I’ve not done this analysis yet. All I can say at the moment is that if this does cause artefacts, it doesn’t result in destabilisation of the protein structure. I’m also try to look into what is causing this issue but I am unsure as yet. It may well be that PME will end up being the best solution but until I’ve done the analysis I can’t say more than what I’ve seen from some basic analysis of the protein simulations. The good thing about ubiquitin and lysozyme in particular is that there is plenty of experimental data (e.g. NMR) to compare to (as done in typical GROMOS protein FF parameterisation/testing).

  3. Agreed here too. But to go back to my point, I’d be hesitant to use a force field (or setup) that required a specific thermostat/setting to be applied to dampen the errors that otherwise occur.

And you’re right with the 53A6L lipids. If you do really want to use these parameters in GROMACS, you should use MA charge groups and RF with the group cut-off scheme. However, along the lines of the thermostat, I’d be hesitant to use a force field that required such a specific cut-off scheme. Particularly given the results seen with the long protein simulations using the same setup and the fact that they eventually become unstable. Another potential issue with the Poger 53A6L lipids is that because the correction/fudge to increase the area is in the head group, you can’t simulate non-PC phospholipids either. So, in my mind, this limits this force field to simulations of PC membranes, using and older version of GROMACS, and without being able to confidently simulate proteins. Well, that’s unless you do something else to increase the APL such as playing around with the cut-off. This gets you into another whole world of pain but is something I’ve also been working on and will include in these couple of papers (but won’t go into now otherwise I won’t have much to write about that I haven’t talked about in this thread!).

And, as per above with the proteins, I am also looking at potential artificial ordering with the SA RF membrane simulations at the cut-off (this analysis is one of the few bits left to do and, again, Chris Oostenbrink gets the credit for the suggestion here). You mentioned the artefacts with an atomistic cut-off and water. If possible can you give some references here? I know that for liquids with low dielectric permittivity you definitely get a substantial error but from some of the papers that I’ve read, I thought that with water/SPC this artefact was relatively small (e.g. smaller than the other MA/SA differences when comparing to something like PME as a reference). I might well have missed something or be misinterpreting something in the literature here.

I’m glad the discussion has been useful for you, and hopefully also for others out there too.

Cheers

Tom

11

Hi Tom,
thanks for your complete answer. Sure, you already revealed many things and I won’t ask you any other question ;-). But the discussion was very enriching, so thanks again. I will read your future papers describing those things with interest! And I think you deserve the letter P in the CKP lipids, no vanity here ;-).
Regarding your question about some artefacts with an atomistic cut-off and water, I was just thinking to the work of Chris Oostenbrink (https://www.doi.org/10.1002/jcc.26426). There, a radial distribution function on a simple water box shows some (slight) peaks around the cut-off for the atomistic cut-off (Figure 4). It doesn’t depend on the water model. These are some results I could also reproduce with a recent GROMACS version (2019.6).
Indeed, Chris’ idea to check whether such ordering occurs in proteins will be interesting to look at. Maybe also on lipids.
Best,

Patrick