The MMPBSA.py program distributed Open Source in the AmberTools21 package is a powerful tool for end-point free energy calculations on molecular dynamics simulations. In its most simple application, MMPBSA.py is used to calculate the free energy difference between the bound and unbound states of a protein-ligand complex. In order to use it, however, you need to have an Amber-compliant trajectory file, which means you need to setup and run your simulation fairly carefully.

While the Amber Manual and the MMPBSA tutorial provide lots of helpful information, putting everything together into a full pipeline taking you from structure to a free energy is another story. The goal for this guide is to provide a schematic you can follow to get started. This guide assumes you are familiar with molecular dynamics simulations and the theory of MMPBSA.

The easiest way I have found to do this, using only Open Source software, is:

(1) Download your raw PDB file. If you are lucky and it contains a complete set of heavy atoms (excepting perhaps a terminal OXT here and there, which tleap will add for you in step 3) you are good to go.

(2) Use the H++ webserver to determine the protonation states of each residue and add hydrogens as needed. This webserver is particularly convenient because it will allow you to directly download a PQR file that you can use to generate your starting topology and coordinates. Note that you have various options to choose the pH and internal/external dielectric constants for the calculation.

(3) Use tleap to generate your topology (prmtop) and coordinate (mdcor) files for your simulations. Do not forget that you will need not only the prmtop for the solvated complex, but also a dry prmtop for each of the complex, receptor, and ligand. Load the PQR file from H++ and do not forget to set PBRadii *to the same value for all prmtops*. A typical tleap script for setting up your solvated complex would look something like:

source leaprc.protein.ff14SB
source leaprc.water.tip3p
set default PBRadii mbondi2
pdb = loadpdb amber_pdbs/7bbg.pdb
center pdb
alignAxes pdb
solvatebox pdb TIP3PBOX 14
addIons pdb Cl- 98
addIons pdb Na+ 96
saveAmberParm pdb prmtop_files/7bbg.prmtop initial/7bbg.mdcor
quit

Note that 7bbg.pdb is just the renamed output PQR from H++.

(4) Minimize, heat, and equilibrate your structure as desired in OpenMM. You can, of course, use Amber/sander for this part, but I like OpenMM because it keeps things Open Source.

(5) Process your raw trajectory files for MMPBSA.py analysis. While it’s convenient to write in binary for (4), MMPBSA.py will expect Amber format .mdcor files in ASCII for analysis. You can use cpptraj, also included with AmberTools21, to preprocess your DCDs. A cpptraj script for this purpose will look something like:

parm prmtop_files/7bbg.prmtop
parm prmtop_files/7bbg.prmtop [top1]
trajin output/7bbg/1_7bbg_prod.dcd 501 last 10
trajout output/7bbg/1_7bbg_prod_501-2500_skip10.mdcrd [top1]

In which I am loading the binary DCD frames 501 to end (in this case that means 501-2500) with a skip of 10 and then writing the downsampled trajectory to the new .mdcrd file. It’s important that we remove all frames that we do not want to analyze, remembering that adjacent frames/time points are correlated.

(6) Run the MMPBSA.py analysis. This step requires another input file that, for a GB run, will look something like this:

&general
startframe=1,
endframe=6000,
interval=1,
verbose=2,
keep_files=2,
/
&gb
igb=2,
saltcon=0.15,
/

Note well: we need to make sure that our igb option is compatible with the PBRadii used in model building; checking the Amber21 manual (p. 71) PBRadii=mbondi2 works with igb = 2 or 5.

And… that’s pretty much it. You can now run your MMPBSA.py calculation by submitting a command like

mpirun -np 20 --use-hwthread-cpus MMPBSA.py.MPI -O \
-i /data/dnissley/mmpbsa/scripts/7bbg_gb.in \
-o /data/dnissley/mmpbsa/analysis/7bbg/gb_6000frame/FINAL_MMPBSA.dat \
-cp /data/dnissley/mmpbsa/prmtop_files/7bbg_noWAT.prmtop \
-lp /data/dnissley/mmpbsa/prmtop_files/7bbg_TCRm_noWAT.prmtop \
-rp /data/dnissley/mmpbsa/prmtop_files/7bbg_pMHC_noWAT.prmtop \
-sp /data/dnissley/mmpbsa/prmtop_files/7bbg.prmtop \
-y /data/dnissley/mmpbsa/output/7bbg/*_7bbg_prod_501-2500_skip10.mdcrd

Good luck and happy free energy calculations.