Project

General

Profile

Feature #773

Implement GROMOS96 54A7

Added by Justin Lemkul over 6 years ago. Updated almost 4 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
preprocessing (pdb2gmx,grompp)
Target version:
Difficulty:
uncategorized
Close

Description

The new GROMOS96 54A7 parameter set was recently published. It involves minor changes relative to 53A6, but they are significant in that they improve upon 53A6 in meaningful ways. I have a working copy locally that I am testing. It would be nice to include this force field in an upcoming release, and I am willing to commit my files once I'm satisfied I did everything right :)

Associated revisions

Revision 2e7ace75 (diff)
Added by Justin Lemkul over 4 years ago

Added GROMOS96 54A7 force field files.

Files are taken from ATB, which is hosted by Alan Mark, who
was an author on the GROMACS96 54A7 paper, which is cited in
forcefields.doc. So these files are effectively equivalent to
the reference implementation.

Fixes #773.

Change-Id: I800da3d3535a09cd5beb7f2b4d1b42422cc4329f

History

#1 Updated by Thomas Piggot about 6 years ago

Hi,

I have a GROMACS version of the 54A7 files I can provide, if you wish to compare with your files?

The ones I made were modified from the older format GROMACS 54A7 files available from the ATB website (http://compbio.biosci.uq.edu.au/atb/forcefield/gromos54a7_for_gromacs3_and_4.0.tar.gz). Having quickly checked the provided files I think there are some problems with the entries for the new sodium and chloride ions in the [ nonbond_params ] section that I have hopefully fixed correctly (these entries had not been modified from the GROMOS 53A6 files). However, from some initial simulation tests I have been seeing a problem with the ion parameters leading to the unnatural formation of salt crystals (in both the provided and my modified files). I haven't contacted the person who made the initial converted files yet as I haven't had time to fully explore this.

It would be good to be able to ensure that my conversion is correct, if you are happy to share your files too.

Cheers

Tom

#2 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

Hi,

I have a GROMACS version of the 54A7 files I can provide, if you wish to compare with your files?

The ones I made were modified from the older format GROMACS 54A7 files available from the ATB website (http://compbio.biosci.uq.edu.au/atb/forcefield/gromos54a7_for_gromacs3_and_4.0.tar.gz). Having quickly checked the provided files I think there are some problems with the entries for the new sodium and chloride ions in the [ nonbond_params ] section that I have hopefully fixed correctly (these entries had not been modified from the GROMOS 53A6 files). However, from some initial simulation tests I have been seeing a problem with the ion parameters leading to the unnatural formation of salt crystals (in both the provided and my modified files). I haven't contacted the person who made the initial converted files yet as I haven't had time to fully explore this.

It would be good to be able to ensure that my conversion is correct, if you are happy to share your files too.

Cheers

Tom

I just ran a 10-ns simulation of 100 mM NaCl and didn't see any sign of crystallization beyond transient ion pairings. I did not get my files from ATB; I created them by modifying the existing 53A6 files according to the literature. If you'd like to try my force field files, I can send them to you privately if you like. I do not want to post them here yet since they are not official.

-Justin

#3 Updated by Thomas Piggot about 6 years ago

Yes, it would be good to have a look at your files. As I mentioned I have not explored this problem fully, I only made the files for a colleague who saw salt crystals forming after 10-20 ns in a membrane-protein system with 1 M NaCl (and in a couple of repeat simulations).

Thanks

Tom

#4 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

Yes, it would be good to have a look at your files. As I mentioned I have not explored this problem fully, I only made the files for a colleague who saw salt crystals forming after 10-20 ns in a membrane-protein system with 1 M NaCl (and in a couple of repeat simulations).

I think this result is more due to the fact that probably all biomolecular force fields will fail at such high salt concentrations. Doubtful it's limited to only 54A7. In any case, I'll send you a tarball of my files that you can try.

-Justin

#5 Updated by Thomas Piggot about 6 years ago

I had wondered about this, but multiple simulations of the exact same system using the GROMOS 53A6 force field didn't cause the same problem. Maybe at these a high concentrations the new parameters are just more likely to form salt crystals?

From the files you sent (thanks by the way!) most things are the same between my files and your files and I have actually spotted another couple of problems with the initial files I downloaded. There are, however, still differences in the nonbonded interactions. From what I could see these are (apart from things like any rounding differences):

1. I think you have forgotten to convert the NA+ and CL- interaction parameters for types CChl onwards (until NUrea) in the [ nonbond_params ] section. Currently these are just the same as the 53A6 interactions.

2. You do not have interactions between AR and NA+ or CL-. You also have interactions between NA+ and CL- with SI twice.

3. I think the C6 value of the NA+ and CU2+ interaction may not be correct in your files.

4. The final difference is that I had made a mistake with the NA+ to OM interaction (I had mistakenly read the OM - NA+ interaction from table 8 of the Oostenbrink paper GROMOS 53A6 paper as to use the OM C12 interaction 2 and not interaction 3 as it should be). You do actually mention something about this entry with the "53A6 has (NA+)--OM, paper has --OM" comment, however I don't really follow the comment.

Let me know if you agree with these points and I can also send you my files if you want to take a look? I wouldn't think that the NA+ to OM mistake will have caused the salt crystals issue, however I shall run a simulation to check.

Cheers

Tom

#6 Updated by Thomas Piggot about 6 years ago

Sorry, I just made an error in the message I sent (it has been a long week!).

It was to do with the mistake I thought I had made (point 4. in the message above). I think that this mistake is actually in your files where you have used the OM C12 interaction type 2 when it should be 3 for this interaction. I am really sorry about getting this the wrong way round. Your comment makes more sense now, but I think you have misread the table in the Oostenbrink 2004 53A6 paper, where the C12 OM interaction with NA+ is specified as type 3 and you have used type 2.

Cheers

Tom

#7 Updated by Thomas Piggot about 6 years ago

Apologies for all the messages but I have spotted what I think is another issue with your non-bonded interactions. The C12 value for the new N to O interaction should be 1.943e-06 and not 1.523e-06.

Hopefully that is everything now!

Cheers

Tom

#8 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

I had wondered about this, but multiple simulations of the exact same system using the GROMOS 53A6 force field didn't cause the same problem. Maybe at these a high concentrations the new parameters are just more likely to form salt crystals?

From the files you sent (thanks by the way!) most things are the same between my files and your files and I have actually spotted another couple of problems with the initial files I downloaded. There are, however, still differences in the nonbonded interactions. From what I could see these are (apart from things like any rounding differences):

1. I think you have forgotten to convert the NA+ and CL- interaction parameters for types CChl onwards (until NUrea) in the [ nonbond_params ] section. Currently these are just the same as the 53A6 interactions.

Oops, I'll fix these.

2. You do not have interactions between AR and NA+ or CL-. You also have interactions between NA+ and CL- with SI twice.

I'll add the AR interactions. I didn't realize before that all the C12 interactions were the same; since they were missing from Table 8 of the 53A6 paper I left them out. I've killed the duplicate SI entries.

3. I think the C6 value of the NA+ and CU2+ interaction may not be correct in your files.

Yep, screwed that up, didn't I?

4. The final difference is that I had made a mistake with the NA+ to OM interaction (I had mistakenly read the OM - NA+ interaction from table 8 of the Oostenbrink paper GROMOS 53A6 paper as to use the OM C12 interaction 2 and not interaction 3 as it should be). You do actually mention something about this entry with the "53A6 has (NA+)--OM, paper has --OM" comment, however I don't really follow the comment.

Here's where I find a problem. The Table itself is inconsistent. For Row 37, Column 2, the value is 2, but for Row 2, Column 37, it is 3. I suspect someone needs to resolve this, since the order (I,J) vs (J,I) should not matter.

Let me know if you agree with these points and I can also send you my files if you want to take a look? I wouldn't think that the NA+ to OM mistake will have caused the salt crystals issue, however I shall run a simulation to check.

Cheers

Tom

#9 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

Apologies for all the messages but I have spotted what I think is another issue with your non-bonded interactions. The C12 value for the new N to O interaction should be 1.943e-06 and not 1.523e-06.

I don't follow - why is type II C12 for N involved if the interaction has been modified to be type I in 54A7?

-Justin

Hopefully that is everything now!

Cheers

Tom

#10 Updated by Thomas Piggot about 6 years ago

Hi,

I think your confusion over both the last two points (C12 for the N to O interaction and C12 for the NA+ to OM interaction) are down to how you are interpreting table 8 from the Oostenbrink paper. It can be bit confusing at first (which is indeed why I thought I had made the mistake before). Take, for example, the interaction between NA+ and OM. For the value to use for NA+ when combining the different parameters, first you need to go to the NA+ row (row 37) and then along to the OM column (column 2) and this number is 2 (so you use the second entry for NA+ in the C12 section of table 7). For OM you also need to do the same, so go to the OM row (row 2) and then along to the NA+ column (column labelled 37) and this entry is 3 (so you need to use the third entry in the OM row in table 7). You always first go down to the correct row for the atomtype of interest and then along to the column of the other atomtype in the interaction.

Moving onto the new N to O interaction, I think the modification made for GROMOS 54A7 is simply to change the entry on row 1, column 6 of table 8 from 2 to 1. This means that the O value to use for the combination of parameters is changed from the second entry (1.130) in table 7 to the first one (1.000). In contrast, the value to use for N remains as the second entry (1.943) in table 7 (as row 6, column 1 of table 8 is unchanged).

Hopefully this makes sense. As to why some of the interactions are given as type 3 when there are no corresponding entries for type 3 in table 7, or why different C12 types are given when the entries for these two are identical in table 7, I have no idea (as a guess maybe different parameters once existed in an ancient GROMOS force field?).

Cheers

Tom

#11 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

Hi,

I think your confusion over both the last two points (C12 for the N to O interaction and C12 for the NA+ to OM interaction) are down to how you are interpreting table 8 from the Oostenbrink paper. It can be bit confusing at first (which is indeed why I thought I had made the mistake before). Take, for example, the interaction between NA+ and OM. For the value to use for NA+ when combining the different parameters, first you need to go to the NA+ row (row 37) and then along to the OM column (column 2) and this number is 2 (so you use the second entry for NA+ in the C12 section of table 7). For OM you also need to do the same, so go to the OM row (row 2) and then along to the NA+ column (column labelled 37) and this entry is 3 (so you need to use the third entry in the OM row in table 7). You always first go down to the correct row for the atomtype of interest and then along to the column of the other atomtype in the interaction.

OK, thanks for the education! This was not at all apparent to me from reading the table, and it is easy to go cross-eyed when reading it ;) I think I have everything fixed with respect to these incorrect entries.

Moving onto the new N to O interaction, I think the modification made for GROMOS 54A7 is simply to change the entry on row 1, column 6 of table 8 from 2 to 1. This means that the O value to use for the combination of parameters is changed from the second entry (1.130) in table 7 to the first one (1.000). In contrast, the value to use for N remains as the second entry (1.943) in table 7 (as row 6, column 1 of table 8 is unchanged).

Got it.

Hopefully this makes sense. As to why some of the interactions are given as type 3 when there are no corresponding entries for type 3 in table 7, or why different C12 types are given when the entries for these two are identical in table 7, I have no idea (as a guess maybe different parameters once existed in an ancient GROMOS force field?).

I wondered as well, maybe some ancient construct in the way they've done their code that is still hanging around. Oh well.

I'll send you an updated ffnonbonded.itp and we'll see if we've arrived at the same parameters.

-Justin

#12 Updated by Thomas Piggot about 6 years ago

Hi,

Two final issues and then we arrive at the same files :). I have tested this using a single point energy calculation of a solvated membrane-protein system, where after the two modifications below our files give the same energies.

The first issue is a simple one and is not in the ffnonbonded.itp but rather something in your aminoacids.c.tdb I missed before. You have forgotten to update the N-C-CA-O2 dihedral in this file to the new dihedrals.

The second is regarding the ffnonbonded.itp and the new CH3p atomtype. You have only explicitly included the interaction of CH3p to OM as this is the only interaction which uses a different C12 parameter to the standard one (i.e. type 1 in table 7 of the Oostenbrink paper). The idea to let GROMACS calculate the other interactions from the entries in the [ atomtypes ] section does not, however, work for all of the atomtypes. This is because some of the atomtpyes (OA, NT and NR) don't use type 1 from table 7 in their interactions with themselves but rather type 2 (and so their entries in the [ atomtypes ] section are for the C12 type 2). This means that the interactions GROMACS calculates for these atomypes with CH3p uses type 2 for these atomtypes and not the desired type 1. To get solve this problem the CH3p to OA, NT and NR interactions need to be explicitly given. To make things less confusing I think it might just be easier to explicitly give all of the CH3p parameters in the [ nonbond_params ] and [ pairtypes ] sections.

If you want I can send you the corrected aminoacids.c.tdb file and also the CH3p nonbonded parameters. You can also get these parameters in the initial GROMACS conversion files available from the ATB. Now that I am happy I know what the problems are in these files I will contact the person who provided them to let him know of the problems.

Finally back to one of my initial questions, I haven't been able to find anything obvious in my conversion (which I am happy is correct now) to provide an explanation as to why the new ion parameters formed salt crystals (when the same simulation using the 53A6 parameters didn't). I assume the modifications to the LJ parameters must just make this more likely at this high concentration of 1M.

Cheers

Tom

#13 Updated by Justin Lemkul about 6 years ago

Thomas Piggot wrote:

Hi,

Two final issues and then we arrive at the same files :). I have tested this using a single point energy calculation of a solvated membrane-protein system, where after the two modifications below our files give the same energies.

The first issue is a simple one and is not in the ffnonbonded.itp but rather something in your aminoacids.c.tdb I missed before. You have forgotten to update the N-C-CA-O2 dihedral in this file to the new dihedrals.

The second is regarding the ffnonbonded.itp and the new CH3p atomtype. You have only explicitly included the interaction of CH3p to OM as this is the only interaction which uses a different C12 parameter to the standard one (i.e. type 1 in table 7 of the Oostenbrink paper). The idea to let GROMACS calculate the other interactions from the entries in the [ atomtypes ] section does not, however, work for all of the atomtypes. This is because some of the atomtpyes (OA, NT and NR) don't use type 1 from table 7 in their interactions with themselves but rather type 2 (and so their entries in the [ atomtypes ] section are for the C12 type 2). This means that the interactions GROMACS calculates for these atomypes with CH3p uses type 2 for these atomtypes and not the desired type 1. To get solve this problem the CH3p to OA, NT and NR interactions need to be explicitly given. To make things less confusing I think it might just be easier to explicitly give all of the CH3p parameters in the [ nonbond_params ] and [ pairtypes ] sections.

Agreed.

If you want I can send you the corrected aminoacids.c.tdb file and also the CH3p nonbonded parameters. You can also get these parameters in the initial GROMACS conversion files available from the ATB. Now that I am happy I know what the problems are in these files I will contact the person who provided them to let him know of the problems.

If you could send me your files, that'd be great, thanks. I had avoided asking for them before to make sure that we could arrive at the same parameters independently. Now, with just the two small tweaks, I'm happy to save a bit of typing and take yours :)

Finally back to one of my initial questions, I haven't been able to find anything obvious in my conversion (which I am happy is correct now) to provide an explanation as to why the new ion parameters formed salt crystals (when the same simulation using the 53A6 parameters didn't). I assume the modifications to the LJ parameters must just make this more likely at this high concentration of 1M.

That seems likely. For CL-, the C6 attractive coefficient was decreased slightly, but the C12 repulsive term was decreased by about 40%. For NA+, C6 became slightly more attractive, but C12 became somewhat more repulsive. So many changes make it hard to pinpoint, but since these parameters were significantly modified, I'd suspect they're the culprit. Is the crystallization with RF or PME for electrostatics?

-Justin

#14 Updated by Thomas Piggot about 6 years ago

Great, I will email you my files.

Regarding the crystallisation it is with PME. I haven't test with reaction field.

Cheers

Tom

#15 Updated by Thomas Piggot over 5 years ago

Hi,

Just so you know, I have spotted one more 'problem' with the force field files. Nothing too major, but the charges for the updated DPPC in aminoacids.rtp are still the standard GROMOS charges and not the Chiu charges, so this needs changing.

By they way, do you know if the force field is going to be incorporated into GROMACS at some point soon? If not, it might be worth adding it to the user contribution section on the website.

Cheers

Tom

#16 Updated by Justin Lemkul over 5 years ago

Thomas Piggot wrote:

Hi,

Just so you know, I have spotted one more 'problem' with the force field files. Nothing too major, but the charges for the updated DPPC in aminoacids.rtp are still the standard GROMOS charges and not the Chiu charges, so this needs changing.

By they way, do you know if the force field is going to be incorporated into GROMACS at some point soon? If not, it might be worth adding it to the user contribution section on the website.

I guess this is sort of a moot point now - you can get the 54A7 files in GROMACS format from ATB. I suppose I'll check my own files over to make sure they agree, but at this point the task is simply for pedantry. Since people can get the files from ATB, I see no real need to post them on the GROMACS site as well. I've listed the target version as 5.0 since there is a feature freeze in place for 4.6. This means it will likely be some time before the official implementation of these files, though their inclusion is largely trivial.

#17 Updated by Manuel Melo over 5 years ago

Hi,

I spotted a couple of small omissions in the .rtp file from the ATB. I write them down here just as a reminder for when the files are finally added to the tree:
- an angle potential in the ACE residue is lacking (something that had also affected the 53a6 .rtp file but that has since been corrected);
- the AIB residue, which is now present in the GROMOS files, is absent from the GROMACS version.

Cheers,
Manuel

#18 Updated by Mark Abraham over 4 years ago

What remains to do here?

#19 Updated by Justin Lemkul over 4 years ago

I'll upload the files now.

#20 Updated by Justin Lemkul over 4 years ago

  • Status changed from New to Fix uploaded

#21 Updated by Mark Abraham over 4 years ago

Looks good.

For the record, supporting new force fields is something the Stockholm crew does not want to spend any time on. Unless there's a publication, we're not interested at all (so this force field is fine). If it involves format conversion, then the person who wants it in GROMACS has the burden of showing that the conversion is accurate compared to the reference implementation (wherever that is). A file already in GROMACS format from the website of an author of the paper is good enough :-)

Otherwise, it can go on the contrib part of the website :-)

#22 Updated by Justin Lemkul over 4 years ago

Sounds like a good policy. I put this one up in particular because (1) people asked about it a lot, (2) it is very important since 53A6 has some definite deficiencies, and (3) ATB posted the files soon after I started making them, so life is easy in that sense.

#23 Updated by Justin Lemkul over 4 years ago

  • Status changed from Fix uploaded to Resolved
  • % Done changed from 0 to 100

#24 Updated by Rossen Apostolov almost 4 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF