GROMACS Tutorial

Step Three: The Workflow

The free energy change of transforming a system from state A to state B, ΔGAB, is a function of a coupling parameter, λ, which indicates the level of change that has taken place between states A and B, that is, the extent to which the Hamiltonian has been perturbed and the system has been transformed. Simulations conducted at different values of λ allow us to plot a ∂H/∂λ curve, from which ΔGAB is derived. The first step in planning free energy calculations is how many λ points will be used to describe the transformation from state A (λ = 0) to state B (λ = 1), with the goal of collecting adequate data to exhaustively sample phase space and produce a reliable ∂H/∂λ curve.

λ = 0 λ = 0.5 λ = 1

For decoupling Coulombic interactions, which depend linearly upon λ, an equidistant λ spacing from 0 to 1 should usually suffice, with λ spacing values of of 0.05-0.1 being common. Note that this is a broad generalization, and in fact many molecules will require a great deal more finesse, such as those that interact very strongly with the surrounding environment through hydrogen bonding. In this case, λ spacing may need to be decreased, such that more points are used, perhaps even in an asymmetric fashion.

For decoupling van der Waals interactions, λ spacing can be much trickier. For reasons described by Shirts et al. and elsewhere, a great deal of λ points may be necessary to properly describe the transformation. Clustering λ points in regions where the slope of the ∂H/∂λ curve is steep may be necessary. For the purposes of this tutorial, we will use equidistant λ spacing of 0.05, but in many cases, one may need to use different spacing, especially clustering values in the range of 0.6 ≤ λ ≤ 0.8.

For each value of λ, a complete workflow (energy minimization, equilibration, and data collection) must be conducted. I generally find it most efficient to run these jobs as batches, passing the output of one step directly into the next. The workflow utilized here will be:

  1. Steepest descents minimization
  2. L-BFGS minimization
  3. NVT equilibration
  4. NPT equilibration
  5. Data collection under an NPT ensemble

Two energy minimization steps are used in order to better minimize the starting structure. I have found that L-BFGS converges very prematurely, resulting in unstable systems, but when used in conjunction with steepest descents (which has its own problems finding "true" energy minima), the results are quite satisfactory.

The .mdp files for the different steps of the tutorial (for λ = 0 only) are provided at the following links:

I am also making available two accessory Perl scripts, write_mdp.pl and write_sh.pl, that you can use to very quickly produce all of the input files needed for running these simulations.

You can also download a shell script (job.sh) for running these jobs according to the workflow I will describe.

If you issue:

perl write_mdp.pl em_steep.mdp

You will receive files named em_steep_0.mdp, em_steep_1.mdp, ... em_steep_20.mdp, with the relevant values of init_lambda_state inserted into the filename. Analogously, you can produce the rest of your .mdp files in the same way (em_l-bfgs.mdp, nvt.mdp, npt.mdp, and md.mdp).

The write_sh.pl script works in the same fashion, generating job_*.sh files for each value of init_lambda_state to execute the commands in the workflow. You will need to change the values of $FREE_ENERGY and $MDP, otherwise they likely won't work.


Notes on .mdp settings

Before proceeding, it is important to understand the free energy-related settings in the .mdp files (using em_steep_0.mdp as an example). They are the same for all steps in the workflow. I will assume that you are familiar with the other settings found in the .mdp file. If this is not the case, please refer to some more basic tutorial material before proceeding.


free_energy = yes Indicates that we are doing a free energy calculation, and that interpolation between the A and B states of the chosen molecule (defined below) will occur.

init_lambda_state = 0 In previous GROMACS versions, the "init_lambda" keyword specified a single value of λ directly. As of version 5.0, λ is now a vector that allows for transformation of bonded and nonbonded interactions. With the init_lambda_state keyword, we specify an index (starting from zero) of the vector to be utilized in the simulation (more on this later).

delta_lambda = 0 The value of λ can be incremented by some amount per timestep (i.e., δλ/δt) in a technique called "slow growth." This method can have significant errors associated with it, and thus we will make no time-dependent changes to our λ values.

fep_lambdas = (nothing) You will note that this keyword is not specified. In previous GROMACS versions, the use of init_lambda and foreign_lambda controlled the value of λi and the additional values of λ for which energy differences would be evaluated for configurations at λi. This is no longer the case. One can explicitly set values of λ in the fep_lambdas keyword, but instead we allow the calc_lambda_neighbors keyword (see below) to automatically determine these additional values.

calc_lambda_neighbors = 1 The number of neighboring windows for which energy differences are computed with respect to λi. For instance, if init_lambda_state is set to 10, then energy differences with respect to λ states 9 and 11 are computed during the run with calc_lambda_neighbors = 1.

vdw_lambdas = ... An array of λ values for the transformation of van der Waals interactions.

coul_lambdas = ... An array of λ values for the transformation of Coulombic (electrostatic) interactions.

bonded_lambdas = ... An array of λ values for the transformation of bonded interactions.

restraint_lambdas = ... An array of λ values for the transformation of restraints.

mass_lambdas = ... An array of λ values for the transformation of atomic masses; used if tranforming the chemical nature of the molecule.

temperature_lambdas = ... An array of λ values for the transformation of temperatures; used only for simulated tempering.

sc-alpha = 0.5 The value of the α scaling factor used in the "soft-core" Lennard-Jones equation. See equations 4.124 - 4.126 in the GROMACS manual, section 4.5.1, for a complete description of this term, as well as the next three.

sc-coul = no Transform Coulombic interactions linearly. This is the default behavior, but is written out for clarity.

sc-power = 1.0 The power for λ in the soft-core equation.

sc-sigma = 0.3 The value of σ assigned to any atom types that have C6 or C12 parameters equal to zero or σ < sc-sigma (typically H atoms). This value is used in the soft-core Lennard-Jones equation.

couple-moltype = Methane The name of the [moleculetype] in that will have its topology interpolated from state A to state B. Note that the name given here must match a [moleculetype] name, and not the residue name. In some cases these may be the same, but I have maintained different names for the [moleculetype] and component residue for instructive purposes.

couple-lambda0 = vdw The types of nonbonded interactions that are present in state A between the interpolated [moleculetype] and the remainder of the system. The value "vdw" indicates that only van der Waals terms are active between methane and water; there are no solute-solvent Coulombic interactions.

couple-lambda1 = none The types of nonbonded interactions that are present in state B between the interpolated [moleculetype] and the remainder of the system. The value "none" indicates that both van der Waals and Coulombic interactions are off in state B. Relative to couple-lambda0, this indicates that only van der Waals terms have been turned off.

couple-intramol = no Do not decouple intramolecular interactions. That is, the λ factor is applied to only solute-solvent nonbonded interactions and not solute-solute nonbonded interactions.

Setting "couple-intramol = yes" is useful for larger molecules that may have intramolecular interactions occuring at longer distance (e.g. peptides or other polymers). Otherwise, distal parts of the molecule will interact via explicit pair interactions in an unnaturally strong manner (since solute-solvent interactions are weakened as a function of λ, but the intramolecular terms would not be), giving rise to configurations that will incorrectly bias the resulting free energy.

nstdhdl = 10 The frequency with which ∂H/∂λ and ΔH are written to dhdl.xvg. A value of 100 would probably suffice, since the resulting values will be highly correlated and the files will get very large. You may wish to increase this value to 100 for your own work.

There are also a number of differences in the .mdp settings that will be used here relative to the settings used by Shirts et al. and those in the Mobley tutorial. Among these:

  1. rlist=rcoulomb=rvdw=1.2. In order to use PME, rlist must be equal to rcoulomb, a check that was enforced in grompp sometime after the release of version 3.3.1. The Verlet cutoff scheme introduced in version 4.6 that allows for buffered neighbor lists also requires rvdw=rcoulomb. The value of rlist will be tuned during the run for the purposes of energy conservation, thus providing the necessary buffer for the switched van der Waals interactions. The switching range (from 1.0-1.2 nm) agrees with the methods of Shirts et al. for the treatment of solute-solvent interactions.
  2. Temperature coupling is handled through the use of the Langevin piston method (via the "sd" integrator), rather than an Andersen thermostat. There is no Andersen thermostat yet implemented in GROMACS.
  3. The temperature is set to 300 K, as in the .mdp files provided in the Mobley tutorial. The original study conducted by Shirts et al. set the reference temperature to 298 K, the thermodynamic standard state. In the interest of reproducing the results of the original tutorial, I leave the temperature set to 300 K. The difference in the resulting free energy at 298 K should be relatively small and will be discussed later.
  4. tau_t = 1.0. In the Mobley tutorial, the inverse friction coefficient was set to 0.1, but upon recommendation of one of the GROMACS developers, this term has been increased to 1.0 to avoid over-damping the dynamics of water.

Let us take a moment to dissect the λ vectors a bit more closely. As an example, for init_lambda_state = 0, that means we are specifying that the state in the transformation corresponds to the vector with index 0 in the *_lambdas keywords, like so:

; init_lambda_state        0    1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18   19   20
vdw_lambdas              = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

Setting initial_lambda_state = 1 would correspond to the next column to the right (λ for van der Waals = 0.05), while initial_lambda_state = 20 would specify the final column (λ vector), in which the value of λ for van der Waals = 1.0. For the purposes of this tutorial, we are only transforming van der Waals interactions, leaving everything related to charges, bonded parameters, restraints, etc. alone. Please note that, in this example, where charges are zero in the topology and the .mdp settings never indicate that charges should be present (neither couple-lambda0 nor couple-lambda1 include 'q'), the choice of values for coul_lambdas is irrelevant, but this is not always the case!

Back: Topology Next: Energy Minimization

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