GROMACS version: 2020.3
GROMACS modification: No
Here post your question
Greetings! I’m a neophyte trying to port a new force field into gromacs. I have a quick question about atom names in the [ atomtypes ] for a ffnonbonded.itp in a force field.
In the user manual, the atom name at the very beginning of an atom definition contains a short string for identifying the atom before the interger number used to designate the element number. I’ve noticed that certain general force fields present with gromacs actually include two names here, particularly OPLS-aa. This field has an atom name, then a bond type name, and then the element number. As far as I can tell, the second name ties the entry in the ffnonbonded.itp to sister entries in the ffbonded.itp file and seems to act as a sort of wild card to allow a single bonded entry cooperate with a number of different nonbonded types.
I can’t find anything in the gromacs manual about how to instruct gromacs as to whether or not the second name should be included, so I must assume that either gromacs looks for multiple names automatically, or that these names must be defined somewhere that I haven’t found. What are the limits of how this can be used? Can I add ten names before the integer? Does it only look for the first name under certain circumstances and the second name in others, like using the first name to identify VdW interactions and the second name for occurrences in the bonded interactions? If I add a [ nonbond_params ] directive for including special non-diagonal interactions, bypassing combination rules, must it refer to the first name and not the second or does it not care?
If anyone can give me some guidance on this, I would be appreciative. Thank you very much!
Long story short, each force field is a little bit different and there is some magic within the code that determines how to interpret the contents of [atomtypes] based on the number of fields read. Most force fields do not use the alternate bonded atom type mapping; it’s basically an OPLS thing.
One nonbonded type must be mapped to one bonded type. There are far fewer of the latter so they can be reused, but the nonbonded types are unique.
This directive will use the nonbonded type (first column) since it is purely a LJ concept. It doesn’t care about bonded atom types.
I would suggest looking at basically any of the other force fields as “typical” implementations to get comfortable with these concepts before trying to work off of OPLS, which is a bit of an oddball.
The [ atomtypes ] section contains optional bonded atom type and atomic number entries which are unfortunately not listed and described in the reference manual. I now pushed up a change to add those to the description and topology format table for 2021.4.
Thank you very much, both of you. This is helpful. I will take a closer look at some of the other force fields in the package because it may give me some other ideas on how to deal with my particular issue.
I have one follow-up question. Based upon what you’ve said, I think I can guess the answer is probably ‘no’, but I’ll ask anyway.
Within a molecular topology, if I have an molecule specified in [ moleculetype ] with an [ atoms ] directive after it, can I identify the ‘unique’ atoms by the bond-type name rather than by the first identifier in the line of an [ atomtypes ] atom?
I know this seems like a bit of a crazy question, but I’ll add some background for context. We’re physicists working to simulate a liquid crystal. Right now, we have a custom force field with Buckingham potentials which is specifically designed to simulate the liquid crystal molecule we’re interested in. This force field uses a topology where bond-types, VdW types and partial charges are all fairly independently specified. In this case, there are more unique bonded interaction types than there are unique VdW types. In gromacs, if and only if the [ atoms ] specification doesn’t care which moniker I choose, if I identify the atoms by the bonded-type names and also specify the partial charge for each atomic center directly in that topology, the non-bonded [ atomtypes ] specification becomes really quite compact and more similar to the original form of the field. If I didn’t have to name by the most unique name in the [ atomtypes ], I wouldn’t have quite so many terms to specify in [ nonbond_params ] since I’m not really working with an available combination rule. If you look closely at OPLS, this is also true; while there are some 900 atom types, the actual VdW parameters being used to specify unique L-J potentials is a small fraction of that, with most of the variability coming from partial charge specification and bond-type naming combinations.
I think I can get around this by simply using only one identifier in the [ atomtypes ] which matches the identifier in the bonded interactions and then computer generating the extra non-unique [ nonbond_params ] types that result. It’s not nearly as compact a table, but it would get me to where I’m going. I think the only table that probably actually matters is the one that Gromacs produces in the background from the inputs that I never actually see;-)
If you have any other thoughts for something I’m overlooking, I’m all ears! Thank you for all your help!