Update/ Please use the binder only for input generation, to run the simulation please use a third-party OpenMM platform. This is kept here for my convenience but may be not so consistent to this repo. The simulation package is designed for OpenMM.
conda create -n openmm-env python=3.7
conda activate openmm-env
conda install -c omnia openmm-setup
################################################################## Recommended workflow.
-
Docking with pl3 as per https://github.com/quantaosun/pl3
-
MD input files generation with Charmm gui input generator (solution builder), with files generated by pl3 (By converting Docked1.pdb to Docked1.sdf with obabel)
-
MD simulation with local Gromacs or gmx inside a gmx_mmpbsa environment
Modify REAMD for loop to only run the first 2, with step5.production.mdp nsteps extended 10 times (from 1 to 10 ns) of the original value
An example of README on normal local Linux without GPU could be
#!/bin/csh
set init = step3_input
set mini_prefix = step4.0_minimization
set equi_prefix = step4.1_equilibration
set prod_prefix = step5_production
set prod_step = step5
# Minimization
# In the case that there is a problem during minimization using a single precision of GROMACS, please try to use
# a double precision of GROMACS only for the minimization step.
gmx grompp -f ${mini_prefix}.mdp -o ${mini_prefix}.tpr -c ${init}.gro -r ${init}.gro -p topol.top -n index.ndx -maxwarn -1
gmx mdrun -v -deffnm ${mini_prefix}
# Equilibration
gmx grompp -f ${equi_prefix}.mdp -o ${equi_prefix}.tpr -c ${mini_prefix}.gro -r ${init}.gro -p topol.top -n index.ndx
gmx mdrun -v -deffnm ${equi_prefix}
# Production
set cnt = 1
set cntmax = 2
while ( ${cnt} <= ${cntmax} )
@ pcnt = ${cnt} - 1
set istep = ${prod_step}_${cnt}
set pstep = ${prod_step}_${pcnt}
if ( ${cnt} == 1 ) then
set pstep = ${equi_prefix}
gmx grompp -f ${prod_prefix}.mdp -o ${istep}.tpr -c ${pstep}.gro -p topol.top -n index.ndx
else
gmx grompp -f ${prod_prefix}.mdp -o ${istep}.tpr -c ${pstep}.gro -t ${pstep}.cpt -p topol.top -n index.ndx
endif
gmx mdrun -v -deffnm ${istep}
@ cnt += 1
end
An example on normal local Linux with GPU could be (-ntmpi added)
#!/bin/csh
set init = step3_input
set mini_prefix = step4.0_minimization
set equi_prefix = step4.1_equilibration
set prod_prefix = step5_production
set prod_step = step5
# Minimization
# In the case that there is a problem during minimization using a single precision of GROMACS, please try to use
# a double precision of GROMACS only for the minimization step.
gmx grompp -f ${mini_prefix}.mdp -o ${mini_prefix}.tpr -c ${init}.gro -r ${init}.gro -p topol.top -n index.ndx -maxwarn -1
gmx mdrun -v -deffnm ${mini_prefix} -ntmpi 1
# Equilibration
gmx grompp -f ${equi_prefix}.mdp -o ${equi_prefix}.tpr -c ${mini_prefix}.gro -r ${init}.gro -p topol.top -n index.ndx
gmx mdrun -v -deffnm ${equi_prefix} -ntmpi 1
# Production
set cnt = 1
set cntmax = 2
while ( ${cnt} <= ${cntmax} )
@ pcnt = ${cnt} - 1
set istep = ${prod_step}_${cnt}
set pstep = ${prod_step}_${pcnt}
if ( ${cnt} == 1 ) then
set pstep = ${equi_prefix}
gmx grompp -f ${prod_prefix}.mdp -o ${istep}.tpr -c ${pstep}.gro -p topol.top -n index.ndx
else
gmx grompp -f ${prod_prefix}.mdp -o ${istep}.tpr -c ${pstep}.gro -t ${pstep}.cpt -p topol.top -n index.ndx
endif
gmx mdrun -v -deffnm ${istep} -ntmpi 1
@ cnt += 1
end
#!/bin/bash
#
# Generated by CHARMM-GUI (http://www.charmm-gui.org) v3.7
#
# This folder contains GROMACS formatted CHARMM36 force fields, a pre-optimized PDB structure, and GROMACS inputs.
# All input files were optimized for GROMACS 2019.2 or above, so lower version of GROMACS can cause some errors.
# We adopted the Verlet cut-off scheme for all minimization, equilibration, and production steps because it is
# faster and more accurate than the group scheme. If you have a trouble with a performance of Verlet scheme while
# running parallelized simulation, you should check if you are using appropriate command line.
# For MPI parallelizing, we recommand following command:
# mpirun -np $NUM_CPU gmx mdrun -ntomp 1
init="step3_input"
mini_prefix="step4.0_minimization"
equi_prefix="step4.1_equilibration"
prod_prefix="step5_production"
prod_step="step5"
# Minimization
# In the case that there is a problem during minimization using a single precision of GROMACS, please try to use
# a double precision of GROMACS only for the minimization step.
#gmx grompp -f "${mini_prefix}.mdp" -o "${mini_prefix}.tpr" -c "${init}.gro" -r "${init}.gro" -p topol.top -n index.ndx -maxwarn -7
#gmx mdrun -v -deffnm "${mini_prefix}" -ntmpi 1
# Equilibration
#gmx grompp -f "${equi_prefix}.mdp" -o "${equi_prefix}.tpr" -c "${mini_prefix}.gro" -r "${init}.gro" -p topol.top -n index.ndx -maxwarn -7
#gmx mdrun -v -deffnm "${equi_prefix}" -ntmpi 1
# Production
cnt=1
cntmax=2
while [ $cnt -le $cntmax ]
do
pcnt=$((cnt-1))
istep="${prod_step}_${cnt}"
pstep="${prod_step}_${pcnt}"
if [ $cnt -eq 1 ]
then
pstep="${equi_prefix}"
gmx grompp -f "${prod_prefix}.mdp" -o "${istep}.tpr" -c "${pstep}.gro" -p topol.top -n index.ndx -maxwarn -7
else
gmx grompp -f "${prod_prefix}.mdp" -o "${istep}.tpr" -c "${pstep}.gro" -t "${pstep}.cpt" -p topol.top -n index.ndx -maxwarn -7
fi
gmx mdrun -v -deffnm "${istep}" -ntmpi 1
cnt=$((cnt+1))
done
mdp file after 10 time logner than default (nsteps has been modified with one more 0)
integrator = md
dt = 0.002
nsteps = 5000000
nstxtcout = 50000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
rlist = 1.2
rcoulomb = 1.2
coulombtype = PME
;
tcoupl = Nose-Hoover
tc_grps = SOLU SOLV
tau_t = 1.0 1.0
ref_t = 303.15 303.15
;
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-5
ref_p = 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = SOLU SOLV
;
- MMPBSA calcultion of protein-ligand based on MD trajectories (and per-residue decomposition)
- Visualisation and/or chart/imge making with gmx_MMPBSA_ana
Elaboration on step 3, modify README to only run 1 nano second instead of 10, so the final trajectory only named as step5_1.xtc, and all mmpbsa calculation will be on the basis of step5_1
Elaboration on step 4 (When there is no halogen atoms, for halogen contained situation, see next below). To know the index for protein and ligand, run
gmx make_ndx -f step5_1.gro -o index2.ndx
1 | 13
To run mmpbsa calculation with decompositon
gmx_MMPBSA -O -i mmpbsa.in -cs step5_1.tpr -ci index2.ndx -cg 1 13 -ct step5_1.xtc -cp topol.top -o FINAL_RESULTS.dat -eo FINAL_RESULTS.csv
where mmpbsa.in is
&general
sys_name="Prot-Lig-ST",
startframe=1,
endframe=12,
/
&gb
igb=5, saltcon=0.150,
/
&decomp
idecomp=2, dec_verbose=3,
print_res="within 4"
/
To determine how many frames in your trajectory file, use pymol
load step5_1.gro
load step5_1.xtc
In case the small molecules contains F Cl Br, it is necessary to delete the virtual site not supported by gmx_MMPBSA.
There is a much easier work around of halogen atom issue, whcih is just using OpenFF instead of Charmm, all the below discussion is only necessary if you used Charmm for the small molecule.
See details from https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/examples/Protein_ligand_LPH_atoms_CHARMMff/, note the modification of topol.top is extremely easy to go wrong, be careful.
- Similar to above, generate a new index file
Command line:
gmx make_ndx -f step5_1.tpr -o index2.ndx
Reading structure file
Reading file step5_1.tpr, VERSION 2022.2 (single precision)
Reading file step5_1.tpr, VERSION 2022.2 (single precision)
Going to read 0 old index file(s)
Analysing residue names:
There are: 361 Protein residues
There are: 32371 Other residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
0 System : 102688 atoms
1 Protein : 5924 atoms
2 Protein-H : 2952 atoms
3 C-alpha : 361 atoms
4 Backbone : 1083 atoms
5 MainChain : 1443 atoms
6 MainChain+Cb : 1781 atoms
7 MainChain+H : 1788 atoms
8 SideChain : 4136 atoms
9 SideChain-H : 1509 atoms
10 Prot-Masses : 5924 atoms
11 non-Protein : 96764 atoms
12 Other : 96764 atoms
13 UNL : 36 atoms
14 POT : 91 atoms
15 CLA : 100 atoms
16 TIP3 : 96537 atoms
nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups
'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues
't': atom type '|': or 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
'ri': residue index
> splitat 13
0 System : 102688 atoms
1 Protein : 5924 atoms
2 Protein-H : 2952 atoms
3 C-alpha : 361 atoms
4 Backbone : 1083 atoms
5 MainChain : 1443 atoms
6 MainChain+Cb : 1781 atoms
7 MainChain+H : 1788 atoms
8 SideChain : 4136 atoms
9 SideChain-H : 1509 atoms
10 Prot-Masses : 5924 atoms
11 non-Protein : 96764 atoms
12 Other : 96764 atoms
13 UNL : 36 atoms
14 POT : 91 atoms
15 CLA : 100 atoms
16 TIP3 : 96537 atoms
17 UNL_O_5925 : 1 atoms
...................................
49 UNL_H11_5957 : 1 atoms
50 UNL_H12_5958 : 1 atoms
51 UNL_H13_5959 : 1 atoms
52 UNL_LP1_5960 : 1 atoms
nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups
'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues
't': atom type '|': or 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
'ri': residue index
> 13&!52
Copied index group 13 'UNL'
Copied index group 52 'UNL_LP1_5960'
Complemented group: 102687 atoms
Merged two groups with AND: 36 102687 -> 35
53 UNL_&_!UNL_LP1_5960 : 35 atoms
> name 53 lig
0 System : 102688 atoms
1 Protein : 5924 atoms
2 Protein-H : 2952 atoms
3 C-alpha : 361 atoms
4 Backbone : 1083 atoms
5 MainChain : 1443 atoms
6 MainChain+Cb : 1781 atoms
7 MainChain+H : 1788 atoms
8 SideChain : 4136 atoms
9 SideChain-H : 1509 atoms
10 Prot-Masses : 5924 atoms
11 non-Protein : 96764 atoms
12 Other : 96764 atoms
13 UNL : 36 atoms
14 POT : 91 atoms
15 CLA : 100 atoms
16 TIP3 : 96537 atoms
17 UNL_O_5925 : 1 atoms
18 UNL_C1_5926 : 1 atoms
....................................
50 UNL_H12_5958 : 1 atoms
51 UNL_H13_5959 : 1 atoms
52 UNL_LP1_5960 : 1 atoms
53 lig : 35 atoms
> 1|53
Copied index group 1 'Protein'
Copied index group 53 'lig'
Merged two groups with OR: 5924 35 -> 5959
54 Protein_lig : 5959 atoms
0 System : 102688 atoms
1 Protein : 5924 atoms
2 Protein-H : 2952 atoms
3 C-alpha : 361 atoms
4 Backbone : 1083 atoms
5 MainChain : 1443 atoms
6 MainChain+Cb : 1781 atoms
7 MainChain+H : 1788 atoms
8 SideChain : 4136 atoms
9 SideChain-H : 1509 atoms
10 Prot-Masses : 5924 atoms
11 non-Protein : 96764 atoms
12 Other : 96764 atoms
13 UNL : 36 atoms
14 POT : 91 atoms
15 CLA : 100 atoms
16 TIP3 : 96537 atoms
17 UNL_O_5925 : 1 atoms
...................................
50 UNL_H12_5958 : 1 atoms
51 UNL_H13_5959 : 1 atoms
52 UNL_LP1_5960 : 1 atoms
53 lig : 35 atoms
54 Protein_lig : 5959 atoms
> del 17-52
>
0 System : 102688 atoms
1 Protein : 5924 atoms
2 Protein-H : 2952 atoms
3 C-alpha : 361 atoms
4 Backbone : 1083 atoms
5 MainChain : 1443 atoms
6 MainChain+Cb : 1781 atoms
7 MainChain+H : 1788 atoms
8 SideChain : 4136 atoms
9 SideChain-H : 1509 atoms
10 Prot-Masses : 5924 atoms
11 non-Protein : 96764 atoms
12 Other : 96764 atoms
13 UNL : 36 atoms
14 POT : 91 atoms
15 CLA : 100 atoms
16 TIP3 : 96537 atoms
17 lig : 35 atoms
18 Protein_lig : 5959 atoms
So we have generated group 17 which is a new version of 13 but without any lone pair halogen atoms, and a protein-ligand complex without LPH, next we also need to generate a pdb file without LPH as the -cs input, and generate a new trajectory file without LPH, at last update the topol file which is refering to /toppar/UNL.itp
PDB file
echo 18 | gmx trjconv -s step5_1.tpr -f step5_1.xtc -dump 0 -o step5_1_noLP.pdb -n index2.ndx
Trajectory file
echo 18 | gmx trjconv -s step5_1.tpr -f step5_1.xtc -dump 0 -o step5_1_noLP.xtc -n index2.ndx
Next,deleted all the information related with LPH particles (atom numbers 52). Dleted the information for LPH particles in atoms ( atom 52 ) and pairs (4 lines need to be deleted, in other case if there are 2 LPH you need delete 8 lines). Besides, delete the whole [ virtual_sites ] and [ exclusions ]
Finally run
gmx_MMPBSA -O -i mmpbsa.in -cs step5_1.tpr -ci index2.ndx -cg 1 17 -ct step5_1.xtc -cp topol.top -o FINAL_RESULTS.dat -eo FINAL_RESULTS.cs
where mmpbsa.in is exactly the same as above (doubt check the frame number for the newly generated trajectory file though)
&general
sys_name="Prot-Lig-ST",
startframe=1,
endframe=2,
/
&gb
igb=5, saltcon=0.150,
/
&decomp
idecomp=2, dec_verbose=3,
print_res="within 4"
/
An example of the final result of such a calculation (per-residue decomposition of binding free energy)
