How to prepare a system

From ChengLab
Revision as of 01:43, 29 April 2020 by Kevin (talk | contribs) (Protein-membrane systems)
Jump to: navigation, search

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

H++

ProPka

PDB2PQR

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

References: [1] [2] [3] [4]

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
}