We will begin by generating a topology for ethanol and solvating it in a box of water. For this tutorial, we will be using the CHARMM force field, specifically the CHARMM General Force Field to represent ethanol. Download a PDB file of ethanol here. You will also need the CHARMM force field in GROMACS format, which you can download from Prof. Alex MacKerell's website.

Generate the topology for ethanol and center it in a cubic box of water with the following commands:

gmx pdb2gmx -f etoh.pdb -o etoh.gro

gmx editconf -f etoh.gro -o etoh_box.gro -c -box 4.0

gmx solvate -cp etoh_box.gro -cs spc216.gro -o etoh_solv.gro -p topol.top

The free energy change of transforming a system from state A to state B, ΔG_AB, is a function of a coupling parameter, λ, which indicates the level of change that has taken place between states A (λ = 0) and B (λ = 1). λ therefore takes on values within {0..1} and represents the extent to which the Hamiltonian has been perturbed and the system has been transformed. The Hamiltonian at a given value of λ is therefore a combination of the A- and B-state energy functions:

Simulations conducted at different values of λ allow us to plot a ∂H/∂λ curve, from which ΔG_AB is obtained. 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.

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. A larger number 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. NVT equilibration
3. NPT equilibration
4. Data collection under an NPT ensemble

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

• Steepest descents EM
• NVT equilibration
• NPT equilibration
• Production MD

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

You can also download a Bash 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_30.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 (nvt.mdp, npt.mdp, and md.mdp).

In job.sh, you will need to change the values of $FREE_ENERGY and $MDP, otherwise the script certainly won't work.

Notes on .mdp settings

Before proceeding, it is important to understand the free energy-related settings in the .mdp files; this tutorial assumes that you are familiar with the basic MD settings found in the .mdp file, as well as those related to free energy calculations. If this is not the case, please refer to some more basic tutorial material and the simpler free energy tutorial before proceeding.

Given that we are transforming both Coulombic and van der Waals terms here, we will take a moment to look at 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   21   22   23   24   25   27   28   29   30
vdw_lambdas              = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 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.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; We are not transforming any bonded interactions
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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; There are restraints in this system
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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated tempering here
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 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 electrostatics = 0.10), while initial_lambda_state = 30 would specify the final column (λ vector), in which the value of λ for van der Waals and electrostatics = 1.0. The conventional approach when transforming both electrostatic and van der Waals terms is to remove the charges first, followed by van der Waals terms. If done in the opposite order, atoms would not repel from one another and charges would collapse into themselves, resulting in Coulombic singularities and the simulations will crash. Conversely, if one is transforming from a non-interacting state, the van der Waals terms should be switched on first, followed by electrostatic terms.