GROMACS version:2020.3
GROMACS modification: No
Here post your question
Hi,
I am looking at interaction energies of one molecule with all other molecules by using the gmx energy and looking at (Coulomb SR + LJ-SR)/num of molecules.
My system is homogenous and only contains one type of molecule - no waters.
I have tried this with 4 forcefields: GAFF, GAFF2, Charmm36 carb and glycam06.
For Glycam06, the Coulomb potential energy is positive whereas the 3 other forcefields give negative Coulomb potential energies. This does not make sense to me, why is it different? The 4 systems have similar numbers of hydrogen bonds and charges on the oxygen atoms are similar too in the topology files.
Why do I see a difference in sign? I do not understand why this model would have a repulsive nature when the three others do not.
If I was to change the sign, the interaction energy is similar to the others for Glycam06.
Should I run md on a single molecule to calculate the energy difference between bound and unbound?
Your help in understanding this coulomb potential energy short range term and how to use it in the context of an interaction energy is greatly appreciated,
Thanks in advance
The net Coulomb energy will depend a lot on how the bonded interactions are handed. I don’t know the details of all those force-fields, but are for instance the exclusions settings the same, i.e. over how many bonds atom pairs are excluded for non-bonded interactions?
Another important effect could be PME tuning which could be different for different runs. One should always analyze the sum of SR+LR energy, not SR alone.
Thank you for answering, I have had a look in the topology files and the number of exclusions is the same but the gen-pairs comment is different: yes for gaff and no for glycam06.
In the log file, the LJ1-4 and Coulomb1-4 terms are also not present in the glycam06. Does it have something to do with this.
Another thing worth noting is that to convert glycam06 amber topology and coordinate files, I used ACPYPE because of the refitting needed.
Is it that the Coulomb SR for glycam06 has the Coulomb1-4 term included in it?
Can I ask for these terms to be calculated separately in gromacs?
gen-pairs only sets how the 1-4 parameters are generated. IIRC glycam has different scaling factors for different atom type pairs and all 1-4 parameters are looked up based on atom type. But I would still expect 1-4 energy terms.
With exclusions along the same number of bonds I would expect similar Coulomb energies.
The topology file for the glycam parameters has explicit sigma and epsilon terms for the pairs so I guess that’s why gen-pairs is off.
I have attached the list of terms I see when I use gmx energy with my enregy file from my simulation with Glycam06.
The LJ-14 and Coulomb-14 are not available. I do not know if this means I will have to manually calculate the 1-4 energy terms if I want to get the interaction energy from this simulation.
Bear in mind that ACPYPE only deals with GLYCAM06 when converting the original topology files built in AMBER to GROMACS.
ACPYPE DOES NOT generate GLYCAM06 topology directly for GMX.
I still do not understand how [pairs_nb] and [pairs] differ when I am not doing a free energy calculation and therefore do not have a A-state and B-state topology…
If all the 1-4 pairs parameters are explicitly given in the [pairtypes] section, are [pairs] and [pairs_nb] not the same, given the same function?
Ok I see, I need to add the fudgeQQ factor if I want to use [pairs] and therefore need to use function 2 because the fudge factor needs to overwrite the default fudgeQQ factor given in my topology file that matches my amber ff for my protein. Otherwise the fudgeQQ factor is applied to the charges given for each atom in the pair, is this correct?
This is actually what doglycans provides as an output gromacs topology file for the [pairs] section.
Yes. The whole idea of fudgeQQ (and fudgeLJ) is that they apply to all 1-4 interactions. That’s how most force fields were set up. For some newer force fields it would be convenient to be able to set fudgeQQ per moleculetype, but that’s not possible in GROMACS. But type 2 is a solution there.
