Introduction

This tutorial is to show how to prepare a system to run on GROMACS, starting with a PDB file for a complex protein/ligand.

It is a mere proof of concept. One should see how ACPYPE is trying to do it in

HowAcpypeWorks.

NB: Besides acpype, antechamber and babel, you will need GROMACS with ffAMBER. GROMACS 4.5.x comes with AMBER force fields now.

Getting GROMACS

Install GROMACS. Something like:

• sudo apt-get install gromacs # if you use Ubuntu Linux, or

• fink install gromacs # if you use Mac

should do the trick.

Getting ffAMBER

If you are installing GMX 4.5.x, you don't need ffAMBER.

For GMX 4.0.x fetch the appropriate version of ffAMBER. Got it with PDFs. You definitely need to read those papers.

To summarise in a script-like way:

# Do this block only if not using GMX 4.5.x
# assuming you are at acpype installation folder
gmx_top_dir=/sw/share/gromacs/top
wget -c http://ffamber.cnsm./ffamber_v4.0-doc.tar.gz
tar xvfz ffamber_v4.0-doc.tar.gz
cd ffamber_v4.0/
sudo cp -bv ffamber_* ${gmx_top_dir}
sudo cp -bv ffamber*/* ${gmx_top_dir}
sudo cp -bv vdwradii.dat ${gmx_top_dir}
sudo cp -bv aminoacids-NA.dat ${gmx_top_dir}/aminoacids.dat
# Get ACPYPE additions:
# go to acpype folder
sudo cp -vb ffamber_additions/* ${gmx_top_dir}

NB: /sw/share/gromacs/top/ is pertinent to Mac. Find the equivalent GROMACS top folder within your platform.

Now you should have an operative GROMACS with ffAMBER.

Running an Example

This is for protein 1BVG.pdb (get it at PDB), a homodimer (HIV protease) with a ligand called DMP. We will use force field Amber99SB.

Luckily, this pdb file has all hydrogens for the ligand, which is necessary for antechamber. One can use either, e.g., babel -h _mol_w/o_H_.pdb _mol_with_H.pdb or YASARA View to automatically add missing hydrogens to your compound. The former just puts 'H' for atom names while the latter puts more meaningful atom name, e.g., 'HCA' for a H bonded to a CA and not a simply 'H' as babel does.

In a script-like way:

# Assuming Complex.pdb (= 1BVG.pdb), split it in Protein.pdb and Ligand.pdb
wget -c "http://www./pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1BVG" -O 1BVG.pdb
grep 'ATOM  ' 1BVG.pdb>| Protein.pdb
grep 'HETATM' 1BVG.pdb>| Ligand.pdb

# If using GMX 4.0.x
# Edit Protein.pdb according to ffAMBER (http://ffamber.cnsm./#usage)
sed s/PRO\ A\ \ \ 1/NPROA\ \ \ 1/g Protein.pdb | sed s/PRO\ B\ \ \ 1/NPROB\ \ \ 1/g | sed s/PHE\ A\ \ 99/CPHEA\ \ 99/g | sed s/PHE\ B\ \ 99/CPHEB\ \ 99/g | sed s/O\ \ \ CPHE/OC1\ CPHE/g | sed s/OXT\ CPHE/OC2\ CPHE/g | sed s/HIS\ /HID\ /g | sed s/LYS\ /LYP\ /g | sed s/CYS\ /CYN\ /g >| ProteinAmber.pdb

# If using GMX 4.5.x
\cp Protein.pdb ProteinAmber.pdb

# Process with pdb2gmx and define water
pdb2gmx -ff amber99sb -f ProteinAmber.pdb -o Protein2.pdb -p Protein.top -water spce -ignh

# Generate Ligand topology file with acpype (GAFF)
acpype -i Ligand.pdb

# Merge Protein2.pdb + updated Ligand_NEW.pdb -> Complex.pdb
grep -h ATOM Protein2.pdb Ligand.acpype/Ligand_NEW.pdb >| Complex.pdb

# Edit Protein.top -> Complex.top
\cp Ligand.acpype/Ligand_GMX.itp Ligand.itp
\cp Protein.top Complex.top

# See NB(1) below

# If using GMX 4.0.x
cat Complex.top | sed '/\#include\ \"ffamber99sb\.itp\"/a#include "Ligand.itp"
' >| Complex2.top

# If using GMX 4.5.x
cat Complex.top | sed '/forcefield\.itp\"/a#include "Ligand.itp"
' >| Complex2.top

echo "Ligand   1" >> Complex2.top
\mv Complex2.top Complex.top

# Setup the box and add water
editconf -bt triclinic -f Complex.pdb -o Complex.pdb -d 1.0
genbox -cp Complex.pdb -cs spc216.gro -o Complex_b4ion.pdb -p Complex.top

# Create em.mdp file
cat << EOF >| em.mdp
define                   = -DFLEXIBLE
integrator               = cg ; steep
nsteps                   = 200
constraints              = none
emtol                    = 1000.0
nstcgsteep               = 10 ; do a steep every 10 steps of cg
emstep                   = 0.01 ; used with steep
nstcomm                  = 1
coulombtype              = PME
ns_type                  = grid
rlist                    = 1.0
rcoulomb                 = 1.0
rvdw                     = 1.4
Tcoupl                   = no
Pcoupl                   = no
gen_vel                  = no
nstxout                  = 0 ; write coords every # step
optimize_fft             = yes
EOF

# Create md.mdp file
cat << EOF >| md.mdp
integrator               = md
nsteps                   = 1000
dt                       = 0.002
constraints              = all-bonds
nstcomm                  = 1
ns_type                  = grid
rlist                    = 1.2
rcoulomb                 = 1.1
rvdw                     = 1.0
vdwtype                  = shift
rvdw-switch              = 0.9
coulombtype              = PME-Switch
Tcoupl                   = v-rescale
tau_t                    = 0.1 0.1
tc-grps                  = protein non-protein
ref_t                    = 300 300
Pcoupl                   = parrinello-rahman
Pcoupltype               = isotropic
tau_p                    = 0.5
compressibility          = 4.5e-5
ref_p                    = 1.0
gen_vel                  = yes
nstxout                  = 2 ; write coords every # step
lincs-iter               = 2
DispCorr                 = EnerPres
optimize_fft             = yes
EOF

# Setup ions
grompp -f em.mdp -c Complex_b4ion.pdb -p Complex.top -o Complex_b4ion.tpr
\cp Complex.top Complex_ion.top

# If using GMX 4.0.x
echo 13| genion -s Complex_b4ion.tpr -o Complex_b4em.pdb -neutral -conc 0.15 -p Complex_ion.top -norandom

# If using GMX 4.5.x
echo 15| genion -s Complex_b4ion.tpr -o Complex_b4em.pdb -neutral -conc 0.15 -p Complex_ion.top -norandom

\mv Complex_ion.top Complex.top

# Run minimisaton
grompp -f em.mdp -c Complex_b4em.pdb -p Complex.top -o em.tpr
mdrun -v -deffnm em

# Run a short simulation
grompp -f md.mdp -c em.gro -p Complex.top -o md.tpr
mdrun -v -deffnm md

# or with openmpi, for a dual core
grompp -f em.mdp -c Complex_b4em.pdb -p Complex.top -o em.tpr
om-mpirun -n 2 mdrun_mpi -v -deffnm em
grompp -f md.mdp -c em.gro -p Complex.top -o md.tpr
om-mpirun -n 2 mdrun_mpi -v -deffnm md

# Visualise with VMD
vmd md.gro md.trr

NB(1): #include "Ligand.itp" has to be inserted right after ffamber**.itp line and before Protein_*.itp line in Complex.top.

Voila!