How to extract the Potential Energy of Protein and Ligand from the Entire System after MD Simulation

GROMACS version: 2022.2
GROMACS modification: No

Hello, dear GROMACS users!

I would appreciate any assistance for a GROMACS beginner.
I have constructed the entire system with a cell membrane + water + ligand + protein, and then proceeded with the MD simulation using the following mdp.

integrator = md
dt = 0.002
nsteps = 500000
nstxout = 50000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
tcoupl = Nose-Hoover
tc_grps = SOLU MEMB SOLV
tau_t = 1.0 1.0 1.0
ref_t = 303.15 303.15 303.15
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
nstcomm = 100
comm_mode = linear
comm_grps = SOLU_MEMB SOLV

Afterwards, to extract the energy specific to the protein and ligand rather than the entire system, I performed the steps as follows.

  1. gmx make_ndx -f step8_601.tpr -o index_2.ndx
    (Here, select atoms for protein and ligand to create a new index.)
  2. gmx trjcat -f step8*.trr -o combined.trr
  3. gmx convert-tpr -s step8_601.tpr -n index_2.ndx -o energy_tpr.tpr
  4. gmx trjconv -s energy_tpr.tpr -f combined.trr -n index_2.ndx -o group_trajectory.trr
  5. gmx mdrun -s energy_tpr.tpr -rerun group_trajectory.trr -deffnm rerun_specific -nb gpu
  6. gmx eneconv -f rerun_specific.edr -o final_energy.edr
  7. gmx energy -f final_energy.edr -o energy_data.xvg

Could there be a problem with my method?

I would be grateful for your response.

It’s difficult to tell from just the commands alone, was there any specific issue you ran into when running the commands that you needed help with?

The only problem I will mention is not a technical one, but a theoretical one. Extracting individual energies of various components has no physical meaning. No force field is parametrized in such a way that the absolute potential energies are targeted. You can calculate the quantity as shown, but the results will not tell you anything.