GROMACS version:
GROMACS modification: Yes/No
Hi @jalemkul@hess
I want to calculated the energy interaction between two protein groups upon MD simulation is finished. My question is: which .gro and .cpt file should I use for the energy calculation? Should I use production files (pro1.part0001.gro and pro1.part0001.cpt) or NPT files (NPT.gro and NPT.cpt)?
Also for -rerun is 200 ps sufficient or I have to run at least for 10 ns?
When creating your new TPR file with grompp you should use the NPT files, as you are setting up a rerun of your production. Remember to define the groups of interest in your system in your new parameter file, e.g. energygrps = protein_1 protein_2. Afterward, execute mdrun on this new TPR file with the flag -rerun pro1.xtc (assuming based on your GRO file name).
Your question about the time is not straightforward to answer, if you want to have thermodynamically sound results you should calculate the interaction energies from the moment your simulation’s potential energy enters equilibrium (stops decreasing) until the end of your simulation.
Note that the interaction energies between parts of the system often do not provide useful information, as the behavior of groups is mediated by many other interactions which are all coupled.
Are you suggesting that I should obtain the potential energy for the entire simulation and then select the time span where the energy is least volatile or at its lowest?
Currently, I have been using RMSD (with the reference being frame 0) to choose the time span to calculate the interaction energy. I select the continuous time frames where the RMSD is more stable.
Any suggestions on this approach would be greatly appreciated.
Thanks …
Exactly. Using RMSD is a common practice for deciding when a system is stable. Still, it is not a thermodynamically meaningful quantity, so I’d say it is more advisable to monitor the system’s potential energy (which you can easily do with gmx energy).
As Hess said, depending on what you want to do and your chosen force field, the interaction energies might not be useful. In the best-case scenario, you have a good approximation of the binding enthalpy while neglecting the entropy cost of binding.
If you want to assess something like binding free energy (dG) or the difference in binding free energy upon mutation of a residue (ddG) you are better off using free energy calculation methods, which albeit more complicated will be much more useful.
As an extreme example: the free-energy of an Na+ Cl- ion pair at contact in water is between -10 and -7 kJ/mol. The electrostatic interaction energy between the two ions is about -500 kJ/mol, the LJ energy is around 20 kJ/mol. So looking at the potential energy alone overestimates the interaction by a factor 50 in this case. This is because interactions in water have high enthalpy-entropy compensation.
If I want to calculate only protein-protein interactions, should I specify “Protein Protein” instead of “Protein Non-Protein” in the “comm_grps” section of the “.mdp” file?
The comm (center of mass motion removal) groups have nothing to do with energy output.
I assume you have also “Protein” in your energy groups. Note that the grid part of the PME electrostatics will not be decomposition. And the usual caveat I mentioned above is also present.
Does this mean that choosing group names (Protein Non-protein or Protein Protein) in the “comm_grps” setting doesn’t affect energy calculations, whether for a Protein-Ligand complex or Protein-Protein complex, as long as the groups are present in the index.ndx file?
Again, the comm groups have nothing to do with the energy caculations (apart from that the kinetic energy of the COMs will be zero). I don’t understand why you write “as long as the groups are present”. To select groups there either need to be present in the index file or they need to be part of the auto-generated groups when not using an index file.
PS: We do not recommend the use of multiple comm groups, except for special cases.
I mentioned this because when I executed gmx grompp with Protein Non-protein specified in the comm_grps, I encountered an error after deleting the SOL entries from the index file for convenience as I was looking only for Protein Protein interaction.