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

User discussions
1

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.

2

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?

3

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.