-
-
Notifications
You must be signed in to change notification settings - Fork 46
Tutorial Using ACPYPE for GROMACS
Alan Silva edited this page Feb 5, 2022
·
4 revisions
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 be aware of how ACPYPE works. If you have suggestions about how to improve this tutorial, please send a comment to [email protected].
NB: Besides acpype, antechamber and openbabel, you will need GROMACS, which comes with AMBER force fields now.
See HERE the several ways you may obtain the tools needed above.
Install GROMACS. Something like:
-
conda install -c bioconda gromacs # if you use conda, or -
sudo apt-get install gromacs # if you use Ubuntu Linux, or -
brew install gromacs # if you use Mac
Should do the trick.
For this tutorial, I'm using GROMACS 2019.1, AmberTools 21.11 and OpenBabel 3.1.0, all installed via conda. And the latest ACPYPE, version 2021.12.14.
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., obabel -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 obabel does.
In a script-like way:
mkdir acpype_tutorial
cd acpype_tutorial
# Assuming Complex.pdb (= 1BVG.pdb), split it in Protein.pdb and Ligand.pdb
wget http://www.pdb.org/pdb/files/1BVG.pdb
grep 'ATOM ' 1BVG.pdb>| Protein.pdb
grep 'HETATM' 1BVG.pdb>| Ligand.pdb
cp Protein.pdb ProteinAmber.pdb
# Process with pdb2gmx and define water
# It's a simple tutorial, you should always check for CYS bridges, HIS and other residues that
# may have different protonation states.
gmx pdb2gmx -ff amber99sb -f ProteinAmber.pdb -o Protein2.pdb -p Protein.top -water spce -ignh
# Generate Ligand topology file with acpype (GAFF2 by default now)
acpype -i Ligand.pdb
# Open Ligand.acpype/acpype.log if you want to know the execution details
# 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
# `#include "Ligand.itp"` has to be inserted right after
# `#include "amber99sb.ff/forcefield.itp"`
# line and before `Protein_*.itp` line in _Complex.top_.
# The command below will do it automatically.
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
gmx editconf -bt triclinic -f Complex.pdb -o Complex.pdb -d 1.0
gmx solvate -cp Complex.pdb -cs spc216.gro -o Complex_b4ion.pdb -p Complex.top
# Create ions.mdp file
cat << EOF >| ions.mdp
; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 200 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbours list and long range forces
cutoff-scheme = Verlet ; Buffered neighbours searching
ns_type = grid ; Method to determine neighbours list (simple, grid)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
EOF
# Create em.mdp file
cat << EOF >| em.mdp
; em.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 200 ; Maximum number of (minimization) steps to perform (should be 50000)
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbours list and long range forces
cutoff-scheme = Verlet ; Buffered neighbours searching
ns_type = grid ; Method to determine neighbours list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
EOF
# Create md.mdp file
cat << EOF >| md.mdp
;define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 1000 ; 2 * 1000 = 2 ps (should be 50000: 100 ps)
dt = 0.002 ; 2 fs
; Output control
nstxout-compressed = 2 ; save compressed coordinates every 20 fs
nstxout = 0 ; save coordinates
nstvout = 0 ; save velocities
nstenergy = 10 ; save energies every 20 fs
nstlog = 10 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbour searching
ns_type = grid ; search neighbouring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
EOF
# Setup ions
gmx grompp -f ions.mdp -c Complex_b4ion.pdb -p Complex.top -o Complex_b4ion.tpr
cp Complex.top Complex_ion.top
# SOL group
echo 15| gmx genion -s Complex_b4ion.tpr -o Complex_b4em.pdb -neutral -conc 0.15 -p Complex_ion.top
mv Complex_ion.top Complex.top
# Run minimisation
gmx grompp -f em.mdp -c Complex_b4em.pdb -p Complex.top -o em.tpr
gmx mdrun -ntmpi 1 -v -deffnm em
# Run a short simulation
gmx grompp -f md.mdp -c em.gro -p Complex.top -o md.tpr # -r em.gro
gmx mdrun -ntmpi 1 -v -deffnm md
# Create vmd.tcl file
cat << EOF >| vmd.tcl
display projection Orthographic
display rendermode GLSL
mol modselect 0 0 protein
mol modstyle 0 0 NewCartoon 0.300000 10.000000 4.100000 0
mol modcolor 0 0 Structure
mol addrep 0
mol modselect 1 0 noh and resname DMP
mol modstyle 1 0 Licorice 0.300000 12.000000 12.000000
mol modcolor 1 0 Type
color Type C white
mol addrep 0
mol modselect 2 0 noh same resid as within 8 of resname DMP
mol modstyle 2 0 Licorice 0.200000 12.000000 12.000000
mol addrep 0
mol modselect 3 0 noh same resid as within 8 of resname DMP
mol representation Licorice 0.200000 12.000000 12.000000
mol modstyle 3 0 HBonds 3.000000 20.000000 1.000000
mol modcolor 3 0 ColorID 0
mol modstyle 3 0 HBonds 3.000000 20.000000 6.000000
mol modcolor 3 0 ColorID 4
mol smoothrep 0 0 5
mol smoothrep 0 1 5
mol smoothrep 0 2 5
mol smoothrep 0 3 5
EOF
# Visualise with VMD
vmd md.gro md.xtc -e vmd.tclVoilà!
(updated on 05 Feb 2022)