GROMACS Tutorial

Step Two: Modify the Topology

The lipid bilayer we will be simulating is DPPC (dipalmitoylphosphatidylcholine), so we need parameters for it as well. But if pdb2gmx is designed to handle only proteins, nucleic acids, and a finite amount of cofactors, where do we get parameters for this molecule, and how do we apply them to our system?

For this tutorial, we will use a united-atom force field to describe the lipids derived by Berger, Edholm, and Jähnig (ref). These parameters can be combined with a GROMOS representation of the protein.

The so-called "Berger lipids" are somewhat of a hybrid between GROMOS atomtypes and OPLS partial charges. Since the long alkane chains are poorly represented by GROMOS bonded parameters, a Ryckaert-Bellemans dihedral potential is used, and a scaling factor of 0.125 is applied to Lennard-Jones 1-4 interactions. These lipid parameters are distributed by D. Peter Tieleman, through his website. While you are there, download the following files:

  • dppc128.pdb - the structure of a 128-lipid DPPC bilayer
  • dppc.itp - the [moleculetype] definition for DPPC
  • lipid.itp - Berger lipid parameters

So what exactly is lipid.itp, and how do you use it? Think of this analogy: gromos53a6.ff/forcefield.itp is to your protein what lipid.itp is to the lipids. Essentially, lipid.itp contains all the atom types, nonbonded parameters, and bonded parameters for a large class of lipids, much like how forcefield.itp, ffnonbonded.itp, and ffbonded.itp function in relation to a protein. That said, we cannot simply #include "lipid.itp" within our topology, since it is at the same level (precedence) as forcefield.itp.

To use the parameters in lipid.itp, we will have to make some changes to our pre-packaged GROMOS96 53A6 force field files (in $GMXLIB/gromos53a6.ff). Make a copy of this directory in your working directory called "gromos53a6_lipid.ff" (assuming you have GROMACS installed in /usr/local/gromacs):

cp -r /usr/local/gromacs/share/gromacs/top/gromos53a6.ff/ ./gromos53a6_lipid.ff

Inside this directory, you will find the following files:

aminoacids.c.tdb
aminoacids.hdb
aminoacids.n.tdb
aminoacids.r2b
aminoacids.rtp
aminoacids.vsd
atomtypes.atp
ff_dum.itp
ffbonded.itp
ffnonbonded.itp
forcefield.doc
forcefield.itp 
ions.itp
spc.itp
spce.itp
tip3p.itp
tip4p.itp
watermodels.dat

Next, modify forcefield.doc to update the description of the force field parameters in it based on the modifications we will make in this tutorial. Mine contains something like:

GROMOS96 53A6 force field, extended to include Berger lipid parameters

The files contained in the gromos53a6_lipid.ff directory constitute the complete description of the force field. The files serve the following purposes:

  • aminoacids.c.tdb and aminoacids.n.tdb - These files are read by pdb2gmx and list the available patching actions that can be applied to protein chains
  • aminoacids.hdb - This file is read by pdb2gmx and contains a database for constructing hydrogen atoms
  • aminoacids.r2b - This file is read by pdb2gmx and contains translations between conventional residue names and any force field-specific building block names
  • aminoacids.rtp - This file is read by pdb2gmx and contains all the available protein residues (and some cofactors, water, and ions - typical contents of a coordinate file passed to pdb2gmx)
  • atomtypes.atp - This file is read by pdb2gmx and contains all the available atom types in the force field
  • ff_dum.itp - This file contains constants that are used in constructing virtual sites (more on this in a different tutorial, and not relevant for our purposes here
  • ffbonded.itp - All the parameters for bonded interactions in the force field
  • ffnonbonded.itp - All the nonbonded (Lennard-Jones) parameters in the force field
  • ions.itp - Contains [moleculetype] definitions for all monoatomic ions in the force field
  • spc.itp, spce.itp, tip3p.itp, and tip4p.itp - Water model topologies (the standard water model for GROMOS96 simulations is SPC)
  • watermodels.dat - A text file read by pdb2gmx for interactively selecting the water model to be used in the system topology

Now, to add the lipid parameters into the parent force field, we will need to copy and paste the entries in the [atomtypes], [nonbond_params], and [pairtypes] sections from lipid.itp into the corresponding headings within ffnonbonded.itp. You will find that the lipid.itp [atomtypes] section lacks atomic numbers (the at.num column), so add these in. The newly-modified lines should be:

   LO    8    15.9994      0.000     A  2.36400e-03 1.59000e-06 ;carbonyl O, OPLS
  LOM    8    15.9994      0.000     A  2.36400e-03 1.59000e-06 ;carboxyl O, OPLS
  LNL    7    14.0067      0.000     A  3.35300e-03 3.95100e-06 ;Nitrogen, OPLS
   LC    6    12.0110      0.000     A  4.88800e-03 1.35900e-05 ;Carbonyl C, OPLS
  LH1    6    13.0190      0.000     A  4.03100e-03 1.21400e-05 ;CH1, OPLS
  LH2    6    14.0270      0.000     A  7.00200e-03 2.48300e-05 ;CH2, OPLS
   LP   15    30.9738      0.000     A  9.16000e-03 2.50700e-05 ;phosphor, OPLS
  LOS    8    15.9994      0.000     A  2.56300e-03 1.86800e-06 ;ester oxygen, OPLS
  LP2    6    14.0270      0.000     A  5.87400e-03 2.26500e-05 ;RB CH2, Bergers LJ
  LP3    6    15.0350      0.000     A  8.77700e-03 3.38500e-05 ;RB CH3, Bergers LJ
  LC3    6    15.0350      0.000     A  9.35700e-03 3.60900e-05 ;CH3, OPLS
  LC2    6    14.0270      0.000     A  5.94700e-03 1.79000e-05 ;CH2, OPLS

In the [nonbond_params] section, you will find the line ";; parameters for lipid-GROMOS interactions." Delete this line and all of the subsequent lines in this section (preserving the section that starts with ";; lipid-SPC/SPCE interactions"). The protein-specific nonbonded combinations are specific to the deprecated "ffgmx" (a modified GROMOS87) force field. Removing these lines allows the interactions between the protein and the lipids to be generated by the standard combination rules of GROMOS96 53A6. Non-bonded interactions involving atom type HW are also present; since these are all zero you can delete these lines as well, or otherwise rename HW as H to be consistent with the GROMOS96 53A6 naming convention. If you do not rename these lines or remove them, grompp will later fail with a fatal error.

Append the contents of the [ dihedraltypes ] to the corresponding section of ffbonded.itp. Do not be concerned that these lines look a bit different. They are Ryckaert-Bellemans dihedrals, which differ in form from the standard periodic dihedrals used in the default GROMOS96 53A6 force field.

To recap, you must make all of the following changes for the force field to be functional:

  1. Copy [atomtypes] from lipid.itp to ffnonbonded.itp and add a column for atomic number
  2. Copy [nonbond_params] from lipid.itp to ffnonbonded.itp
  3. Remove the ;; parameters for lipid-GROMOS interactions and all subsequent lines in the [nonbond_params] section in ffnonbonded.itp
  4. Remove all lines containing "HW" in [nonbond_params] or otherwise rename them to "H"
  5. Copy [pairtypes] from lipid.itp to ffnonbonded.itp
  6. Copy [dihedraltypes] from lipid.itp to ffbonded.itp

Finally, change the #include statement in your topol.top from:

#include "gromos53a6.ff/forcefield.itp"

to:

#include "gromos53a6_lipid.ff/forcefield.itp"

Lastly, we need to include the specific parameters for our DPPC molecules. The process for doing so is quite simple. Just add the line #include "dppc.itp" in your topol.top, somewhere after the position restraints section for the protein, which defines the end of the protein [moleculetype] definition. Doing so is analogous to adding any other small molecule or solvent into the topology. In this section and throughout this tutorial, text in orange will denote lines that you should add, while other text (in black) are lines that should already be present in the topology prior to the modification indicated.

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include DPPC chain topology
#include "dppc.itp"

; Include water topology
#include "gromos53a6_lipid.ff/spc.itp"

Notes on Force Fields

Since the Berger parameters were drawn from aspects of both the GROMOS and OPLS-UA force fields, it is also possible to use the OPLS-AA force field to represent your protein and have a compatible model. You will have to make some modifications to lipid.itp so that it is consistent with OPLS conventions. For information on how to do this, read the following links:

https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2006-May/021416.html
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2006-August/023587.html
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2006-September/023761.html

Regardless of which setup you choose, you must demonstrate that your model is reasonable. The particular parameter combinations described here have been shown to work in-house and according to reports by other users. It is up to you to convince your audience (i.e., reviewers) that you know what you are doing and that your model is valid. For instance, if you intend to do long simulations, it is also known that GROMOS96 53A6 under-stabilizes α-helices and you may end up with a random coil rather than a correct helix. In this case, you may want to use the newer GROMOS96 54A7 parameter set, or another force field entirely.

If you have correctly followed all of the steps above, you will have a fully functional force field that can be used to process other membrane proteins with pdb2gmx. Doing so removes the need to manually hack topol.top after the fact. Placing the new gromos53a6_lipid.ff directory in $GMXLIB will allow you to use this force field system-wide.

Back: pdb2gmx

Next: Building the Unit Cell

Site design and content copyright Justin Lemkul

Problems with the site? Send them to the Webmaster