Modify a Topology

In this tutorial we will check out different ways to process GROMACS topology files with the command-line interface kimmdy-modify-top. Use kimmdy-modify-top -h to see the arguments you can pass to it for your kimmdy version. The python API can be found in the Reference. You need a GROMACS installation and an environment with kimmdy for the tutorial.

Parameterize a Molecule with Grappa

To parameterize a molecule with grappa using kimmdy-modify-top, you first need to have a GROMACS topology file for your system. Here, we will use a Ace-Ala-Nme peptide. First, download the peptide pdb and create a topology file using a grappa-compatible force field.

# set up your directory
mkdir <tutorial-directory>
cd  <tutorial-directory>
# use an environment with kimmdy
# download files and create topology
wget https://graeter-group.github.io/kimmdy/guide/tutorials/modify-topology-files/pep.pdb
gmx pdb2gmx -f pep.pdb -ignh  # select an appropriate force field from which partial charges and LJ parameters will be taken, e.g. amber99sbildn

Next, re-parameterize the peptide with grappa and check out the resulting topology file. Make sure to change the values for -grappa_tag and grappa_charge_model if the defaults are not applicable.

# parameterize peptide with grappa
kimmdy-modify-top topol.top topol_grappa.top -p
# check out the output file; use a text editor of your choice
gedit topol_grappa.top

For this task you can also use the command-line interface of grappa grappa-gmx. The same results as with the kimmdy-modify-top command above can be achieved with the command below.

grappa_gmx -f topol.top -o topol_grappa.top -t 'latest' -c 'amber99'

You will notice many changes compared to the topology file generated by gmx pdb2gmx. First, all includes were resolved to create a standalone topology file. That’s the reason for having the [ atomtypes ], [ bondtypes ] and other parameter definitions as well as ion and water topology definitions with [ moleculetype ], [atoms] and [bonds] sections.

The change through parameterization will become more apparent by comparing with a non-parameterized, converted topology file.

# parameterize peptide
kimmdy-modify-top topol.top topol_vanilla.top 
# compare results
vimdiff topol_grappa.top topol_vanilla.top

You should see explicit parameters for bonds, angles dihedrals (e.g 2 1 5 6 9 180.0 1.5141575 1) in the grappa parameterized file while the non-parameterized file has implicit parameters which are defined in the interaction type section (2 1 5 6 9). If that’s the case, you have successfully parameterized a molecule with grappa and are finished with this section of the tutorial!

Remove Hydrogens and parameterize for radical starting structures

An additional capability of the tool is to remove hydrogens to create radicals. Which hydrogen will be removed is defined by its atom nr in the topology file (first column). Additionally, the structure file has to be modified to get a compatible pair of structure and topology file. Look for the atom nr of the hydrogen connected to the Ala C-alpha and remove it. Keep the -p flag to reparameterize the molecule using grappa. Use a text editor or molecule visualization program to verify the successful removal of the hydrogen of choice. The nonbonded parameters all stay the same, except for the partial charge of the heavy atom from which the hydrogen is detached: The partial charge of the hydrogen will be added to the heavy atom upon removal of the hydrogen.

# remove hydrogen and re-parameterize peptide with grappa
kimmdy-modify-top topol.top topol_grappa_CA.top -p -c conf.gro -r <hydrogen nr>
# look at the result

To see how removal of a hydrogen changes the minimum energy structure, run a minimization on both the native and radical peptide. Most strikingly, you should see the radical C-alpha become planar, as expected from QM calculations. Use the mdp file mentioned below. Both structures are also shown in molstar panels.

GROMACS mdp file (minim.mdp)

minim.mdp

; minim.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       = 10.0        ; Stop minimization when the maximum force < 10.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; 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 neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor 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

# minimize radical peptide
gmx editconf -f topol_grappa_CA.gro -o topol_grappa_CA_box.gro -d 1 -bt cubic
gmx grompp -f minim.mdp -c topol_grappa_CA_box.gro -p topol_grappa_CA.top -o minim.tpr
gmx mdrun -deffnm minim
# do the same for the native peptide
# look at the result

Minimized structure of native peptide

Minimized structure with deleted hydrogen

This has only been validated for proteins but the implementation is designed to provide a general solution. Hence, non-verified and possibly wrong parameters can easily be generated. For further information, refer to appendix of the original grappa paper.

When using the option to remove hydrogens, make sure to keep a distance of at least 16-20 bonds (corresponding to 2 times the field-of-view of the grappa model) between radicals because grappa is not trained on systems with multiple radicals.

Complex Parameterization

While parsing a topology file, molecules that are neither water nor ions (as identified by the molecule name) are grouped in the Reactive molecule. This Reactive molecule is passed to grappa for parameterization. For topologies with multiple solute components such as protein:DNA complexes or membrane proteins as well as files where the default recognition of water and ions does not work, you can explicitly define included and excluded moleculetypes in the topology. The includes and excludes are supplied as .csv files.