How to prepare a system
Contents
Fetch PDB
You can fetch pdb using
wget https://files.rcsb.org/download/6I1K.pdb
ParmEd
ParmEd is handful for converting common file formats such as psf/pdb to top/gro (Gromacs) or top/inpcrd (Amber).
import parmed as pmd gmx = pmd.load_file('topol.top', xyz='system.gro') gmx.save('system.top', format='amber') gmx.save('system.inpcrd', format='rst7') gmx.save('ionized.psf')
Please note that the Gromacs output comes WITHOUT forcefields.
Protein-membrane systems
CHARMM to Gromacs
File conversion from CHARMM to Gromacs is often encountered, because CHARMM-GUI is very convenient for creating a membrane system. There is a script psf2itp provided by CHARMM-GUI (modified by Kev 2019) to convert a psf to topology file for Gromacs with CHARMM36m ff. You may replace the toppar to change the force-field (but still have to be CHARMM). CHARMM36m ff for Gromacs could be downloaded here.
python psf2itp.py toppar.str PSFfile
Insert your protein into membrane
A simple VMD script deletes "bad lipids" and place your favourite protein inside the membrane:
set POPC "resname POPC" set all [atomselect top all] $all set beta 0 set seltext1 "$POPC and same residue as (within 0.6 of protein)" set sel1 [atomselect top $seltext1] $sel1 set beta 1 set badlipid [atomselect top "name P and beta > 0"] $badlipid get resid set sel [atomselect top "beta<1 "] $sel writepdb prot_memb_hole.pdb
Remove water inside membrane
This VMD script solvates the system and removes unnecessary water inside the membrane
package require solvate solvate 5kxi_popc.psf 5kxi_popc.pdb -o 5kxi_popc_water_TEMP -b 1.5 -minmax {{-38 -38 -39} {39 39 50}} set all [atomselect top all] $all set beta 0 set seltext "segid WT1 to WT99 and same residue as abs(z) <25" set sel [atomselect top $seltext] $sel set beta 1 set badwater [atomselect top "name OH2 and beta > 0"] set seglist [$badwater get segid] set reslist [$badwater get resid] mol delete all package require psfgen resetpsf readpsf 5kxi_popc_water_TEMP.psf coordpdb 5kxi_popc_water_TEMP.pdb foreach segid $seglist resid $reslist { delatom $segid $resid } writepdb 5kxi_popcw.pdb writepsf 5kxi_popcw.psf file delete kcsa_popc_water_TEMP.psf file delete kcsa_popc_water_TEMP.pdb
Protonation states of amino acids
PlayMolecule (Using ProPka)
Improper dihedral problem in Gromacs' Amber force-field
In ffbonded.itp:
CT CV CC NA 4 180.00 4.60240 2
However pdb2gmx will make
CT CC CV NA
for HID, causing error message of "No default Proper Dih. types". Therefore you have to manually find the improper dihedral terms in the topology (usually topol_Protein.itp made by pdb2gmx) and exchange the second and the third indices. You could find the line number in the error message.
Ligand topology - using Amber force-fields in Gromacs
antechamber -i lig.mol2 -fi mol2 -o LIG.mol2 -fo mol2 -j 4 -at gaff -c bcc -nc 0 parmchk -i LIG.mol2 -f mol2 -o LIG.frcmod tleap -s -f tleap.all acpype -p com_solvated.top -x com_solvated.crd -b complex
You can find the tleap.all here.
Ligand topology - using CHARMM force-fields in Gromacs
Just use SwissParam
RNA
Rosetta: rna_helix.py stepwise
ModeRNA: [5]
Some useful VMD scripts
Re-number residues
foreach ii {A B C D} { mol new $ii.pdb set sel [atomselect top all] set oldNumbering [$sel get residue] set newNumbering {} foreach resid $oldNumbering { lappend newNumbering [expr {$resid + 1}] } $sel set resid $newNumbering $sel writepdb $ii.pdb }
Make rotations of coordinates
set sel [atomselect top "segname P1 and resid 1 to 13"] set selcenter [atomselect top "segname P1 and resid 14 and name N"] set centerX [list [$selcenter get x] [$selcenter get y] [$selcenter get z]] set matrixX [trans center $centerX axis y -10] $sel move $matrixX
set sel [atomselect top "protein"] set selcenter [atomselect top "protein"] set centerX [measure center $selcenter] for {set ii 0} {$ii < 360} {incr ii 30} { set matrixX [trans center $centerX axis x -10] $sel move $matrixX }