GROMACS version: 2020.3
GROMACS modification: Yes/No
Here post your question
Dear all,
I’ve done a MD simulation of a protein-ligand complex with gormacs 2020.3. I wonder whether I can calculate the binding energy between the ligand and the protein with MM/PBSA method. I find the MMPBSA 2.1 written by A. Spitaleri, but it can work with gromacs 4, not with gromacs 2020 versions. Do you know whether there is any MMPBSA program that can work with gromcas 2020?
Best regards
Hi,
GMXPBSA can work with the last versions of Gromacs if you edit gmxpbsa0.sh a little bit (for me it worked with Gromacs2019.4)
.
*line 134 becomes: set_variable_default2 "gmx_prefix" "gmx_AVX2"; (gmx_AVX2 is the name of my executable, change with wathever you want)
*line 140 becomes: pdb2gmx="${gmx_prefix} pdb2gmx${gmx_suffix}"; trjconv="${gmx_prefix} trjconv${gmx_suffix}"; mdrun="${gmx_prefix} mdrun${gmx_suffix} -ntmpi 1"; grompp="${gmx_prefix} grompp${gmx_suffix}"; editconf="${gmx_prefix} editconf${gmx_suffix}"; tpbconv="${gmx_prefix} tpbconv${gmx_suffix}"; gmx="${gmx_prefix}" (I added the spaces after gmx_prefix)
*line 149, 151, 155, 168, 172 : change $editconf with $gmx in the conditions ; e.g. for line 149, which becomes: nf=$(which $gmx 2>/dev/null | awk -F / '{print NF-1 }')
*line 198: the mdrun command was changed to mdrun="${gmx} mdrun -ntmpi 1"
Hope that helps.
Nicolas
Hi Jason.
I came across the same issue a few months ago. I tried 3 approaches: the first two (gmxpbs2.1 and g_mmpbsa) failed miserably (I sucessfuly edited the scripts of one and compiled the other, however with unacceptable results after the calculation), while the last (ambertools) while still failing in the Pbsa was marvelous within Gbsa.
If you or someone else wishes, I can later provide my “recipe” in here.
Also, if someone succeeded on the other approaches, I’m also interested in listening. :)
Best wishes
After editing the gmxpbsa0.sh as you have written, I still get an error:
Fatal error:
The group cutoff scheme has been removed since GROMACS 2020. Please use the
Verlet cutoff scheme.
I guess there is still something in the script that needs to edit. Do you have any idea how to solve this matter?
Best regards
If the script is providing you with .mdp files, they are designed for older versions of GROMACS, so you will either need to use something outdated or adapt the .mdp files to be functional with version 2020.
You are right, I had forgotten some modifications I made:
In function_gmx.dat, line 270: I have added a ‘-maxwarn 1’ at the end (i.e. it becomes … -o $name.tpr -maxwarn 1 >>STD_ERR0 2>&1)
In print_files.dat:
Nicolas
PS : @jdandrade If you have a protocol that works with AmberTools, I think it could be of interest to some people that you post it.
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:
load the environment variables from both gromacs and ambertools:
source …/bin/GMXRC
test -f …/ambertools/amber.csh && source …/ambertools/amber.csh
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
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
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:
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.
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
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
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:
“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.
“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).
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
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).
Hi jdandrade.
I’ve faced similar situation. Does the script you wrote can apply the gromacs54a7 united-atom model?
Thank you :)
Hi @Jinyoung
I don’t know. As far as I can imagine, it would only depend on the ambertools: if he can or cannot deal with a UA force-field. I have no idea.
You can use the commands described above: there is not much more to make a working (csh) script.
(a bit more work may be needed to make it work on bash: variable definitions and “echos” new lines are different)
Please tell me if I can be of assistance! :)
Best regards
Thank you jdandrade very much.
I try to run the amber MMPBSA.py with gromacs UA model and No matter there is problem or not, I tell you the result that it can be run or not.
Thank you:)
Hollo ncheron,
I modified the two .dat files and the problem about Mn.mdp file has gone, but a new error appears:
FatalFatal error:
Topology include file “./charmm36-mar2019.ff/forcefield.itp” not found
The input trajectory file was generated with the charmm36-mar2019.ff and I put the charmm36-mar2019.ff directory under the GMXPBSA root directory (the same directory where the receptor.itp and the ligand.itp files are placed). Is this error caused but misplacing of the ff directory? If so where should it be places?
Best regards.
./ means the force field needs to be in the working directory.
I have also tried to put the “./charmm36-mar2019.ff” directory in the working directory (in which the “$GMXPBSA/gmxpbsa0.sh” command is run, but the same error comes out.
If I remember, GMXPBSA creates subfolders and perform the calculations from these subfolders. Thus, placing the .ff directory in the same directory as gmxpbsa0.sh or from your working directory will not work. You could create in advance the directories that will be created by GMXPBSA and place the ff folder there ; or edit gmxpbsa0.sh so that a symbolic link is created from each folder ; or place the .ff folder in the gromacs folder (in {Gromacs folder}/share/gromacs/top/) ; or probably something else.
Nicolas
Nicolas
Hi ncheron,
During my GMXPBSA running, it creates a directory named “RUN1_2myN2-ZMR”. Creating this directory in advance and place the ff folder there doesn’t work, and it will give the error:
The directory RUN1_2myN2-ZMR is already present in /path/to/ZMR-md_2. Please remove it or choose an appropriate value for the variable “run” in the INPUT.dat file.
Putting the ff directory in {Gromacs folder}/share/gromacs/top/ also doesn’t work. The ff directory can still not be found by the scripts:
Fatal error:
Topology include file “./charmm36-mar2019.ff/forcefield.itp” not found
I cannot figure out where the symbolic link should be created in the gmxpbsa0.sh. Could you help with this?
Best regards
Did you try to comment out the function check_run on line 47 of the file function_base.dat? Regarding symbolic link, searching “mkdir” in gmxpbsa0.sh seems to point towards line 409 and 412. I would add a line somewhere there.
It turns out that there are mistakes in gmxpbsa0.sh:
Line 412 “mkdir $FoLdeR” should be “mkdir $RUN/$FoLdeR”
Line 436 “cp -r ./.ff …/$RUN/$FoLdeR" should be "cp -r …/.ff …/$RUN/$FoLdeR” (the ff directory is placed under the root directory [the same directory where the input.dat is]
parmed gromber doesn’t support gromacs UA model.
@Jinyoung: Well, now we know: the limit is the need for AA simulations. Sorry.
@Jason and @ncheron: Despite me having issues with that perl script, I don’t remember needing so much editing to make it run (my issues were with non-sense final results). The major issue was the “old way” of several executables instead of the actual “gmx something”, and pretty small modifications at the “mdp model” files they used in order to “make” the mdp.
I tested it on 2018 gmx version, but I wouldn’t expect too many differences.
Resurrecting this old conversation as the advice on updating GMXPBSA has been useful but I’ve hit a new stumbling block.
After editing the gmxpbsa0.sh script as ncheron/Nicolas suggested to bring it in line with GROMACS 2020+ versions, the script succeeds in running the single point energy calculations for my 100 frames. So far, so good.
Unfortunately when it comes to calculating the LJ parameters, it throws the following ERROR:
ERROR: Variable L14CC is empty
The log files, eg. comp0.log look normal and contain the LJ-14 block:
Step Time
0 0.00000Energies (kJ/mol)
Bond Angle Proper Dih. Per. Imp. Dih. LJ-14
4.85399e+03 1.19592e+04 1.48530e+04 7.07028e+02 5.34866e+03
Coulomb-14 LJ (SR) Coulomb (SR) Potential Kinetic En.
6.03493e+04 -1.03224e+04 -1.05587e+05 -1.78383e+04 4.07081e+02
Total Energy Temperature Pressure (bar)
-1.74312e+04 5.82864e+00 -8.45444e+01
I think the relevant section is lines 1081-1243 - is there something I’m missing/needs updating here? Any help will be most welcome given the difficulty of finding MMPBSA programs compatible with recent GROMACS versions.
Running GMXPBSA2.1 from github; GROMACS 2023.
Thanks,
Natalie