Protein-ligand with LPH atoms BFE calculations (Single Trajectory method) -- CHARMMff files¶

Sample input file for PB calculation with halogenated compounds

&general
sys_name="PB_Halogens",
PBRadii=7,
/
&pb
inp=1, saopt=1, radiopt=0, cavity_surften=0.005, cavity_offset=0.0000
/

Requirements¶

In this case, gmx_MMPBSA requires:

-> Must be defined -- -> Optional, but recommended -- -> Optional

See a detailed list of all the flags in gmx_MMPBSA command line here

In order to generate the corresponding files (The MD Structure+mass(db), index, trajectory and the topology files) without the LPH particles, it's necessary to run a few commands. Bear with me!

Let's generate the index file first:

gmx make_ndx -f com.tpr -o index_mod_gromacs.ndx

  0 System              : 70483 atoms
  1 Protein             :  5580 atoms
  2 Protein-H           :  2817 atoms
  3 C-alpha             :   334 atoms
  4 Backbone            :  1002 atoms
  5 MainChain           :  1335 atoms
  6 MainChain+Cb        :  1654 atoms
  7 MainChain+H         :  1654 atoms
  8 SideChain           :  3926 atoms
  9 SideChain-H         :  1482 atoms
 10 Prot-Masses         :  5580 atoms
 11 non-Protein         : 64903 atoms
 12 Other               : 64903 atoms
 13 3G5                 :    32 atoms
 14 CLA                 :    62 atoms
 15 SOD                 :    63 atoms
 16 TIP3                : 64746 atoms

Splitting the ligand (group 13) by atoms
>splitat 13

Grouping both LPH particles
>47|48

Excluding both LPH particles from the ligand
>13&!49

Naming ligand as lig
>name 50 lig

Grouping rec and lig
>1|50

Cleaning
>del 17-49

save and quit
>q

This is how it should look like at the end

  0 System              : 70483 atoms
  1 Protein             :  5580 atoms
  2 Protein-H           :  2817 atoms
  3 C-alpha             :   334 atoms
  4 Backbone            :  1002 atoms
  5 MainChain           :  1335 atoms
  6 MainChain+Cb        :  1654 atoms
  7 MainChain+H         :  1654 atoms
  8 SideChain           :  3926 atoms
  9 SideChain-H         :  1482 atoms
 10 Prot-Masses         :  5580 atoms
 11 non-Protein         : 64903 atoms
 12 Other               : 64903 atoms
 13 3G5                 :    32 atoms
 14 CLA                 :    62 atoms
 15 SOD                 :    63 atoms
 16 TIP3                : 64746 atoms
 17 lig                 :    30 atoms
 18 Protein_lig         :  5610 atoms

Let's generate the MD Structure+mass(db) file:

echo 18 | gmx trjconv -s com.tpr -f traj_fit.xtc -dump 0 -o str_noLP.pdb -n index_mod_gromacs.ndx

Open str_noLP.pdb in your favorite visualizer and see it doesn't contain the LPH particles. Now, let's generate the trajectory with no LPH particles:

echo 18 | gmx trjconv -s com.tpr -f traj_fit.xtc -o com_traj.xtc -n index_mod_gromacs.ndx

Finally, let's edit the topology file. Go inside the toppar folder and open the HETA.itp file. As you will see, we deleted all the information related with LPH particles (atom numbers 31, and 32 respectively). In this case, we deleted the information for LPH particles in atoms (lines 47, 48) and pairs (lines 124, 132, 150, 153, 154, 157, 158, 159). Besides, delete the whole [ virtual_sites3 ] (lines 296-299) and [ exclusions ] (lines 301-318) fields. The original .itp (HETA_original_with_LPH_info.itp) is included for comparison purposes.

Command-line¶

That being said, once you are in the folder containing all files, the command-line will be as follows:

gmx_MMPBSA -O -i mmpbsa.in -cs str_noLP.pdb -ci index_mod_gromacs.ndx -cg 1 17 -ct com_traj.xtc -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv

mpirun -np 2 gmx_MMPBSA MPI -O -i mmpbsa.in -cs str_noLP.pdb -ci index_mod_gromacs.ndx -cg 1 17 -ct com_traj.xtc -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv

where the mmpbsa.in input file, is a text file containing the following lines:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

Sample input file for PB calculation
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters 
according to what is better for your system.

&general
sys_name="Prot-Lig-ST",
startframe=5,
endframe=9,
solvated_trajectory=0,
# In gmx_MMPBSA v1.5.0 we have added a new PB radii set named charmm_radii. 
# This radii set should be used only with systems prepared with CHARMM force fields. 
# Uncomment the line below to use charmm_radii set
#PBRadii=7,
/
&pb
# radiopt=0 is recommended which means using radii from the prmtop file
# for both the PB calculation and for the NP calculation

istrng=0.15, fillratio=4, radiopt=0, inp=1,
/

See a detailed list of all the options in gmx_MMPBSA input file here as well as several examples

Considerations¶

In this case, a single trajectory (ST) approximation is followed, which means the receptor and ligand structures and trajectories will be obtained from that of the complex. To do so, an MD Structure+mass(db) file (str_noLP.pdb), an index file (index_mod_gromacs.ndx), a trajectory file (com_traj.xtc), and both the receptor and ligand group numbers in the index file (1 17) are needed. The mmpbsa.in input file will contain all the parameters needed for the MM/PB(GB)SA calculation. A topology file is also needed (mandatory) in this case to generate the topology files in amber format with all the terms for CHARMM force field.

A plain text output file with all the statistics (default: FINAL_RESULTS_MMPBSA.dat) and a CSV-format output file containing all energy terms for every frame in every calculation will be saved. The file name in '-eo' flag will be forced to end in [.csv] (FINAL_RESULTS_MMPBSA.csv in this case). This file is only written when specified on the command-line.