IntroductionThis 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.If you have suggestions about how to improve this tutorial, please send a comment (at the bottom of the page). NB: Besides acpype, antechamber and babel, you will need GROMACS with ffAMBER. GROMACS 4.5.x comes with AMBER force fields now. Getting GROMACSInstall GROMACS. Something like:
should do the trick. Getting ffAMBERIf 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 foldergmx_top_dir=/sw/share/gromacs/topwget -c http://ffamber.cnsm.csulb.edu/ffamber_v4.0-doc.tar.gztar xvfz ffamber_v4.0-doc.tar.gzcd 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 foldersudo 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 ExampleThis 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.pdbwget -c "http://www.pdb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1BVG" -O 1BVG.pdbgrep 'ATOM ' 1BVG.pdb>| Protein.pdbgrep 'HETATM' 1BVG.pdb>| Ligand.pdb# If using GMX 4.0.x# Edit Protein.pdb according to ffAMBER (http://ffamber.cnsm.csulb.edu/#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 waterpdb2gmx -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.pdbgrep -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.xcat Complex.top | sed '/\#include\ \"ffamber99sb\.itp\"/a#include "Ligand.itp"' >| Complex2.top# If using GMX 4.5.xcat Complex.top | sed '/forcefield\.itp\"/a#include "Ligand.itp"' >| Complex2.topecho "Ligand 1" >> Complex2.top\mv Complex2.top Complex.top# Setup the box and add watereditconf -bt triclinic -f Complex.pdb -o Complex.pdb -d 1.0genbox -cp Complex.pdb -cs spc216.gro -o Complex_b4ion.pdb -p Complex.top# Create em.mdp filecat << EOF >| em.mdpdefine = -DFLEXIBLEintegrator = cg ; steepnsteps = 200constraints = noneemtol = 1000.0nstcgsteep = 10 ; do a steep every 10 steps of cgemstep = 0.01 ; used with steepnstcomm = 1coulombtype = PMEns_type = gridrlist = 1.0rcoulomb = 1.0rvdw = 1.4Tcoupl = noPcoupl = nogen_vel = nonstxout = 0 ; write coords every # stepoptimize_fft = yesEOF# Create md.mdp filecat << EOF >| md.mdpintegrator = mdnsteps = 1000dt = 0.002constraints = all-bondsnstcomm = 1ns_type = gridrlist = 1.2rcoulomb = 1.1rvdw = 1.0vdwtype = shiftrvdw-switch = 0.9coulombtype = PME-SwitchTcoupl = v-rescaletau_t = 0.1 0.1tc-grps = protein non-proteinref_t = 300 300Pcoupl = parrinello-rahmanPcoupltype = isotropictau_p = 0.5compressibility = 4.5e-5ref_p = 1.0gen_vel = yesnstxout = 2 ; write coords every # steplincs-iter = 2DispCorr = EnerPresoptimize_fft = yesEOF# Setup ionsgrompp -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.xecho 13| genion -s Complex_b4ion.tpr -o Complex_b4em.pdb -neutral -conc 0.15 -p Complex_ion.top -norandom# If using GMX 4.5.xecho 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 minimisatongrompp -f em.mdp -c Complex_b4em.pdb -p Complex.top -o em.tprmdrun -v -deffnm em# Run a short simulationgrompp -f md.mdp -c em.gro -p Complex.top -o md.tprmdrun -v -deffnm md# or with openmpi, for a dual coregrompp -f em.mdp -c Complex_b4em.pdb -p Complex.top -o em.tprom-mpirun -n 2 mdrun_mpi -v -deffnm emgrompp -f md.mdp -c em.gro -p Complex.top -o md.tprom-mpirun -n 2 mdrun_mpi -v -deffnm md# Visualise with VMDvmd 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! |
联系客服
