How to prepare a system

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

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
}