How to Calculate the linear interaction energy

GROMACS version: 2021
GROMACS modification: Yes/No

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:

gmx grompp -f binding_energy2.mdp -c run.gro -p topol.top -o binding.tpr -maxwarn 13

gmx mdrun -s binding.tpr -rerun run.xtc –nb cpu –deffnm ener2

gmx energy -f ener2.edr -o bindingenergy.xvg -b 9000 > binding.log

Select the inhibitor-SOL option

Elj=29.338 ,Eqq=-419.446

The following commands were executed to simulate the CNT-inhibitor-Water:

gmx grompp -f binding_energy1.mdp -c run.gro -p topol.top -o binding.tpr -maxwarn 13

gmx mdrun -s binding.tpr -rerun run.xtc –nb cpu –deffnm ener

In the last step, the “gmx lie” command was used to calculate energy binding.

gmx lie -f ener.edr -b 9000 -Elj 29.338 -Eqq -419.446 -ligand inhibitor -o lie.xvg

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

binding_energy2.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 = inhibitor SOL
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

Yes, either ener.edr or ener2.edr are logical. The main problem is in output numbers from gmx lie.

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.

Coul-SR:CNT-Inhibitor -0.0842631 Kj/mol
LJ-SR:CNT-Inhibitor -0.0527407 Kj/mol

Coul-SR:SAND-Inhibitor 0.000447952 Kj/mol
LJ-SR:SAND-Inhibitor -0.00142884 Kj/mol

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).

1 Like

Thank you so much for taking the time to answer me patiently.

There is a bug in gmx lie starting from gromacs 2018 or 2019. It will be fixed in 2022.

You are affected by the bug if, in the ouput xvg file (lie.xvg), you see exactly the same value on each line. For instance:

         0     48.5812
        50     48.5812
       100     48.5812
       150     48.5812
       200     48.5812

(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.

1 Like

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.