MM/PBSA for the latest gromacs versions

User discussions
8

Hi @ncheron .

As I mentioned, my protocol only works with MMGBSA (All my three different attempts failed for MMPBSA). Also, it seems to need a fully compiled ambertools (no errors).

Given that:

A: Gromacs Files Preparation:

  1. load the environment variables from both gromacs and ambertools:
    source …/bin/GMXRC
    test -f …/ambertools/amber.csh && source …/ambertools/amber.csh

  2. Set a few variables for work:
    set JUMP = 10
    set TOP = full_system_topology_file.top
    set CTOP = dessolvated_system_topology_file.top
    set RTOP = receptor_only_topology_file.top
    set REC = receptor name in the files
    set LTOP = ligand_only_topology_file.top
    set LIG = ligand name in the files
    set TPR = tpr_file_from_simulation.tpr
    set NDX = index_file.ndx
    set TRAJ = full_trajectory_file.xtc or trr

  3. Center trajectory file on the complex and ensure non-broken molecules:
    echo “$REC”’\n’“System” | gmx trjconv -f $TRAJ -n $NDX -s $TPR -pbc whole -center -skip $JUMP -o tmp.xtc
    echo “Protein”’\n’“Protein”’\n’“System” | gmx trjconv -f tmp.xtc -n $NDX -s $TPR -pbc cluster -center -o $TRAJ

  4. Extract initial configuration non-broken configuration files:
    echo “[frames]”’\n’“1” >& first_frame.ndx
    echo “$REC”’\n’“System” | gmx trjconv -f $TRAJ -n $NDX -s $TPR -fr first_frame.ndx -pbc whole -center -o tmp.gro
    echo “Protein”’\n’“Protein”’\n’“System” | gmx trjconv -f tmp.gro -n $NDX -s $TPR -pbc cluster -center -o system.gro
    echo “$REC”’\n’“Protein” | gmx trjconv -f system.gro -n $NDX -s $TPR -pbc whole -center -o tmp.gro
    echo “Protein”’\n’“Protein”’\n’“Protein” | gmx trjconv -f tmp.gro -n $NDX -s $TPR -pbc cluster -center -o complex.gro
    echo “$REC”’\n’"$REC" | gmx trjconv -f system.gro -n $NDX -s $TPR -pbc whole -center -o receptor.gro
    echo “$LIG”’\n’"$LIG" | gmx trjconv -f system.gro -n $NDX -s $TPR -pbc whole -center -o ligand.gro

B: Gromacs to Amber Files Format Conversion:

  1. Make sure your potential does not include in the periodic dihedrals which happen to have a periodicity number equal to zero (often appears when the energy constant of that diherals is parametrized to zero): as explained in a very deeply lost email from the CCL list if I’m not mistaken, it does not make sense and leads to erros (clue, there are some in the "…/share/gromacs/top/amber*.ff/ffbonded.itp files). They have the format (one example line):
    Br CT CT Br 9 0.0 0.00000 0
    but in order do work they must be modified to:
    Br CT CT Br 9 0.0 0.00000 1
    Again, without this kind of correction in all used diheral types from your simulation, you will get an error.

  2. Set a few variables to work:
    set PRMTOP = $TOP:r:t.prmtop
    set PRMCTOP = $CTOP:r:t.prmtop
    set PRMRTOP = $RTOP:r:t.prmtop
    set PRMLTOP = $LTOP:r:t.prmtop
    set TRAJCRD = $TRAJ:r:t.mdcrd

  3. Convert configuration files:
    echo “gromber $TOP system.gro”’\n’“parmout $PRMTOP”’\n’“scnb 2”’\n’“EOF” | $AMBERHOME/bin/parmed
    echo “gromber $CTOP complex.gro”’\n’“parmout $PRMCTOP”’\n’“scnb 2”’\n’“EOF” | $AMBERHOME/bin/parmed
    echo “gromber $RTOP receptor.gro”’\n’“parmout $PRMRTOP”’\n’“scnb 2”’\n’“EOF” | $AMBERHOME/bin/parmed
    echo “gromber $LTOP ligand.gro”’\n’“parmout $PRMLTOP”’\n’“scnb 2”’\n’“EOF” | $AMBERHOME/bin/parmed

  4. Convert trajectory file:
    echo “parm $PRMTOP”’\n’“trajin $TRAJ”’\n’“trajout $TRAJCRD”’\n’“run”’\n’“quit” | $AMBERHOME/bin/cpptraj

C: Prepare the MMPBSA.py input file “mmpbsa.in” and run it:

  1. “General” (first) section of the “mmpbsa.in” file:
    &general
    endframe=50000, verbose=1, interval=2
    use_sander=1, netcdf=1,
    strip_mask=":SOL:NA"
    /
    The “endframe” does not need to match precisely the end, so make sure it is at least “after” the end or it will stop there. Set interval to your choice (it is equivalent to the “-skip” flag in gromacs programs), I usually don’t use it since I already skipped how many I want during the stage A3. The main source of error here is the “strip_mask”: to the best of my understanding, you need to include all “solution” molecules in here (I’m not an amber expert, I just needed to “make it work” in some software and this was the only one I succeced). Seems to be specially important when you have ions. Keep the rest as shown.

  2. “BD” (second) section of the “mmpbsa.in” file:
    &gb
    igb=x, saltcon=y
    /
    You can either select the model (“igb”, “x” is an integer, look at the ambertools manual to choose the appropriate model) and other variables (again, look at the ambertools manual: in “saltcon” for example “y” is a real) or keep it clean to go for the standard settings. I suggest reading and choosing the appropriate settings (I’m definitelly not an expert, I just found a way to make it work with some older models and parameters).

  3. Run MMPBSA.py:
    $AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp $PRMTOP -cp $PRMCTOP -rp $PRMRTOP -lp $PRMLTOP -y $TRAJCRD

  4. Wait: it will make some preparations, then in order will make the calculations for the complex, them the receptor and finally, the ligand. It does not have any checkpoint system in place (despite the generated files naming, disregard that), so be certain your calculation is not interrupted. Also, there seems to be some sorte of parallelization, but I was not able to make it work (I gladly accept suggestions in here).

I think I haven’t fogot anything. Hope it helps (there is not much more than it in my script dealing with it, and I preferred to explain what have to be done).

1 Like