I want to calculate energy binding with the LIE method in Gromacs, but Gromacs 2021 cannot read the energy file correctly (ener.edr). For this purpose, I will explain the simulation step by step so that if there is a problem, help me solve this problem. CNT and inhibitor molecules are simulated in water. An additional simulation (Inhibitor in water) is required to calculate energy binding with the LIE method. According to the Gromacs Manual, I went through all the steps as follows:
Add in .mdp file “energygrps = CNT inhibitor” (in CNT-inhibitor-Water simulation)
Add in .mdp file “energygrps = inhibitor SOL” (in inhibitor-Water simulation)
The following commands were executed to simulate the inhibitor-Water:
At this point, Gromacs cannot extract the relevant data from the ener.edr file and the numbers it reports are unrealistic and incorrect. Thank you for helping me solve this problem.
If ener.edr is from the original simulation without different energygrps specified, then certainly the output of gmx lie will be nonsensical. Do you get an error message, or just wild values?
Potentially related question:
Why are you overriding up to 13 warnings? What is grompp complaining about? These are meant to stop you from doing something unphysical and there’s almost never a reason to use -maxwarn at all, let alone to override this many issues.
I considered maxwarn 13 because of the overwrite of some atoms. Of course, I considered the number too much.
The calculated values are not logical. When I do Sand-inhibitor-water simulation, the values are the same. It does not matter if it is Sand or CNT, or any other molecule; the answer is the same.
My .mdp files are as follows:
binding_energy1.mdp
integrator = md
nsteps = 5000000
dt = 0.002
nstlist = 25
rlist = 1.2
cutoff-scheme = Verlet
vdw-type = cut-off
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
tcoupl = v-rescale
tc-grps = system
tau-t = 0.1
ref-t = 298
gen-vel = yes
gen-temp = 298
pcoupl = Parrinello-Rahman
tau-p = 1.0
compressibility = 1e-5 1e-5 1e-5 0 0 0
ref-p = 1.0
refcoord-scaling = all
pbc = xyz
nstxtcout = 50
nstxout = 50
nstenergy = 50
energygrps = CNT inhibitor
constraints = all-bonds
My question comes from your nomenclature. You’re calling ener.edr in the gmx lie command, but the interaction energies should be stored in ener2.edr - is this just a typo, or are you actually using a different file? Can you please define what you mean by “values are not logical?” What does that mean? Can you extract reasonable interaction energy values from ener2.edr using gmx energy?
After simulating CNT-inhibitor-water in a separate folder, I rerun the final simulation file run.xtc. ener.edr file was created (energygrps = CNT inhibitor in mdp).
After simulating inhibitor-water in a separate folder, I rerun the final file run.xtc. After rerun in this step, the ener2.edr file was created (energygrps = inhibitor SOL in mdp).
I had nothing to do with the energy output file of both simulations and left it unchanged (name of the simulation energy file run.edr). For each simulation I rerun again and get the edr files.The rerun CNT-inhibitor-SOL output energy file is called ener .edr, and the rerun inhibitor-SOL output energy file is called ener2.edr.
The energy binding in the case of CNT or Sand or any other molecule is 204 kj / mol.The binding energy remains constant with the change of molecule (204 kj / mol). any molecule-inhibitor-Water energy binding is equal (204 kj / mol). Do you think the problem is due to incorrect input of edr files?
I got the parameters Elj = 29.338, Eqq = -419.44 from the ener2.edr file.Are these numbers wrong (inhibitor-water)?
-Elj : Lennard-Jones interaction between ligand and solvent
-Eqq : Coulomb interaction between ligand and solvent
What are the values of interaction energy between the CNT and the inhibitor? The formula is straightforward in terms of reproducing what the code is doing, but the values of the inhibitor interaction energies sound reasonable, so what do the interactions with the CNT look like? I suspect there’s essentially no Coulombic contribution, and all LJ. That may mean the LIE assumption breaks down.
The low interaction energy between the species is due to ions in the solution. Because10 Mg, Ca, and 34 Na are present in the 5x5x5 nm3 box reduce interaction energy between species. Does Interaction energy between species correct (204 KJ/mol)?
Those interaction energies suggest the small molecules have essentially no interaction with the CNT or Sand. So yes, the high repulsive value generated by LIE seems reasonable here. The attraction is weak and the desolvation penalty is high (based on the strength of the inhibitor-water interactions).
(Though the value can be different for different system)
If this is the case, the results obtained are incorrect.
If it is necessary for you to perform LIE analysis at this time, you can use the recently released 2022 beta 2. You do not need to perform the simulation again, just use the fixed gmx lie program contained in this version.
Thank you very much. Gromacs 2021 had bugs. When I checked with Gromacs 2022, I got the answer. But there is a question for me that you might help with. The binding energy of two of my simulations is negative. Is the negative binding energy correct? The other two simulations have positive binding energies. Does negative or positive connection energy have a specific meaning?
There is a question for me that you might help with. The binding energy of two of my simulations is negative. Is the negative binding energy correct? The other two simulations have positive binding energies. Does negative or positive connection energy have a specific meaning?
What you refer to as “energy of binding” is the free energy \Delta G for the reaction Protein + Ligand -> Protein-Ligand complex, and it has the same properties with regards to the meaning of its sign as any other chemical reaction.
Dear sir @jalemkul ,
I want to calculate binding free energies for docked complexes. However, protein is a membrane protein. Can I run the simulations for complexes in water to calculate binding free energies with gmx lie using another method in GROMACS?
Thanks