GROMACS Tutorial

Step Two: Modify the Topology

The lipid bilayer we will be simulating is DPPC (dipalmitoyl phosphatidylcholine), 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?

There is no universal force field for simulating all proteins, nucleic acids, lipids, carbohydrates, and arbitrary small molecules. Since we are simulating a system that contains both protein and lipids, how do we solve the force field problem? Under GROMACS, the most widely-used parameters for the lipid component of membrane-protein simulations are commonly called "Berger 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 gromos53a6.ff/forcefield.itp. Make a new directory in your working directory called "gromos53a6_lipid.ff" and copy the following files from gromos53a6.ff into it:


Next, create a forcefield.doc file that contains a description of the force field parameters in it. Mine contains something like:

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

Next, copy and paste the entries in the [ atomtypes ], [ nonbond_params ], and [ pairtypes ] sections from lipid.itp into the appropriate 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. These nonbonded combinations are specific to the deprecated ffgmx 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. Finally, change the #include statement in your from:

#include "gromos53a6.ff/forcefield.itp"


#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, 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 green 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"

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

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

It is also possible to use the OPLS-AA force field to represent your protein. 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:

For a bit more theory on force field combination rules and theory, consult this document. Thanks to Chris Neale for providing this link.

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.

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 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

Bevan Lab Homepage

Virginia Tech Homepage

Virginia Tech Biochemistry

Site design and content copyright 2008-2015 by Justin Lemkul
Problems with the site? Send them to the Webmaster