# How to find out the total functional form of the forcefield that is being used for a current simulation in GROMACS?

This is a doubt regarding the usage of various forcefields together in MD simulation of hydrates. For Eg., we use TIP4P/ice for water molecules, OPLS-UA for methane and EPM2 for carbon dioxide and so on. How exactly does GROMACS work in this scenario? Can it work with different functional forms or is it just different parameters for the same functional form? Because, in the GROMACS topology or .itp file, we sometimes use different numbers for different functional forms for different molecules under one same directive.(Eg: func 1 for OPLS for one molecule and func 2 for GROMOS for some other molecule under the same 'bonds' directive)

Majority of the research work in Molecular dynamics simulations of hydrates is done with such combinations. For example, this paper uses such different forcefields for different molecules in one simulation.

Also in this link, Justin clearly dismisses the idea of using different forcefields together.

It will be helpful if you can clear this confusion. Is there a way to check the functional form of the forcefields that is being used to calculate, in this software?

Thanks and Regards

• Thanks for coming back to post this after the Beta got started, and sorry for those downvoters on Area51! I have messaged Meilani on Facebook so that she knows to copy and paste her answer from here: area51.meta.stackexchange.com/a/31131/190792 – Nike Dattani Apr 29 '20 at 20:49
• @NikeDattani, thankyou for being so patient and kind. – Kavya Mrudula Apr 29 '20 at 21:19
• Maybe you need a bit of clarification about functional. There is the general functionnal of the force field and there may be different functionals for bond, angle, dihedrals in the same force field. From experience, you can use different parameters or different functionals of bond, angle, dihedral in the same general force field. But I am curious about the different general force fields. – BND Apr 30 '20 at 8:52
• @BND That's interesting. Thank you. Can you please provide an example/reference of different functionals being used for the same general functional of a forcefield, so that I can understand on how to search for compatibilities. – Kavya Mrudula Apr 30 '20 at 13:25

I am a regular GROMACS user. I do not think you can just easily have the total functional form "printed out" from GROMACS. However, you can obtain it from the topology file (top and itps) and the mdp file. Let me elaborate a little below.

The functional forms you are referring to in the GROMACS top/itp file are the ones used in bonded interactions - bonds, angles, and (proper and improper) dihedrals - and they indeed may vary across force fields, mostly for historical reasons. For example, GROMOS uses cosine based angle potentials (see GROMACS manual - bonded interactions for the functional form), while other force fields use harmonic angle potentials. Another example is OPLS using Fourier functions (GROMACS manual - bonded interactions for the functional form) for proper dihedrals, while other force fields use other functions. However, most of this functional forms are very similar, or in some cases equivalent: for example, at least in the case of dihedrals, you can actually have the same dihedral potential expressed via a Fourier function or by a Ryckaert-Bellemans function. As you recognized, the functional form of a particular bonded interaction in GROMACS is defined in the itp file by a functional form "ID" - called function type: i.e., the func 1, func 2 you mentioned. This Table from the GROMACS manual contains a complete list of the different functions available in GROMACS for the bonded terms, their "itp ID" (function type), and the parameters each of them requires.

Besides bonded interactions, the total functional form is made up of a second part: the nonbonded interactions (GROMACS manual - nonbonded interactions; see also, e.g., van Gusteren et al). This part is possibly more relevant in the context of your sub-question about force field compatibility. In classical MD, nonbonded interactions most commonly means describing van der Waals and Coulomb interactions. However, it is in how such nonbonded terms are specifically treated which can make a big difference and renders force field mixing not recommendable (i.e., Justin's replies on Researchgate here and here) unless a proper validation is carried out (see below). This is the case because, the specific force fields have been all parametrized with different nonbonded settings. What I refer to here is to the exact treatment of Lennard-Jones and Coulomb interactions. For example,  while the use of Lennard-Jones potential itself is pretty much ubiquitous, the scaling of the 1-4 interactions, i.e., the interactions between atom i and atom j which is 3 bonds away from i, varies considerably across force fields, see, e.g., Table 2 of S. Riniker. This difference (and others) impacts the overall functional form of the force field, and in turn the parameters which make up the force field. So, different force fields end up having different balances of Lennard-Jones and Coulomb interactions: e.g., some force fields have on average partial charges with larger magnitudes than others, and this usually means that the Lennard-Jones parameters will be somewhat less cohesive to compensate for that. The paper of S. Riniker gives a detailed overview of the four major force fields historically developed for biomolecules - AMBER, CHARMM, GROMOS, OPLS -, their refinement over time, and all the different settings which come with that.

Also note that, despite the fact that water models, such as TIP3P or SPC, can be considered a "separate" force field, they are compatible with certain force fields because such force fields have been developed using them. In other words, each force field chose to use a certain water model - possibly because they were the developers of that model - for their force field. The specific force field was then parametrized using that chosen water model, e.g., TIP3P for OPLS and SPC for GROMOS. This then means that you should use TIP3P with OPLS and SPC with GROMOS and not vice versa, etc.

A last note on the validation: force fields have historically been validated in various ways, that is, by using different parametrization targets. For example, GROMOS 53A6 relies primarily on solvation free energies, the coarse-grain Martini force field on free energies of transfer, etc. The chosen parametrization targets reflect the "philosophy" behind the force field, and tell you also about the domain of application of the force field. Now, this means that if you create a new molecule (a new itp file) within that force field by following the guidelines (e.g., use a particular method to compute the partial charges, use the force field default atom type for the Lennard-Jones parameter of the aromatic atoms, etc.), the idea is that you get a molecule which is compatible with the force field and which, in principle, should give you, for example in the case of GROMOS, solvation free energies in agreement with experimental data. Now, you see that if you mix in topologies from different force fields, there's no guarantee that this is the case. Potentially, for a specific application, if you really need to mix in different force field, you should first carefully validate that properties which are relevant for your application and for that system are reproduced well (solvation free energy, mass density, etc.). But this is very application specific, and should be checked very carefully. Now, I think I should leave it here and not go much further off topic. :)

• +100 !! What a great first answer! Welcome to the site and thank you for your contributions! I hope to see more of you! – Nike Dattani May 17 '20 at 20:32
• @NikeDattani Thank you! – Riccardo Alessandri May 17 '20 at 20:35
• @RiccardoAlessandri long time no see! I should have invited you much earlier, my bad... Anyway, glad that Kavya got more detailed information from the right person. – Meilani Wibowo May 17 '20 at 21:36
• @MeilaniWibowo long time indeed :) thanks! No problem - and I should actually spread the word about this SE a bit more myself. – Riccardo Alessandri May 18 '20 at 7:46
• @RiccardoAlessandri Thank you so much for such a clear and detailed explanation. – Kavya Mrudula May 18 '20 at 11:21

I am not a GROMACS user, but have used it once or twice during my course. As far as I can remember and understand, it is common to use different force field in MD if one simulate something solvated in something, e.g. formaldehyde (OPLSAA force field) in water (SPC,TIP3P, etc.). I guess we can also do that in GROMACS, but I do not really know how to specify this in the topology file (I remember that there must be directory which lists available force fields and how to call them in the topology file). There must be some discussions already in the mailing list of GROMACS (http://www.gromacs.org/Support/Mailing_Lists).

For other processes, such as the interaction between organic molecule and DNA, I do not think that you can mix the force fields, e.g. using MM-type for organic molecule and AMBER for DNA. One reason that I can think of is because those force fields are parameterised differently and you might need to do validation too, which I am not really sure how to do that. Also, you cannot mix up force fields for MD and coarse grain.