Alchemical Free Energy Perturbation Calculations

This feature has been contributed to NAMD by the following authors:

Surjit B. Dixit and Christophe Chipot

Equipe de chimie théorique,
Institut nancéien de chimie moléculaire,
UMR CNRS/UHP 7565,
Université Henri Poincaré,
BP 239,
54506 Vanduvre-lès-Nancy cedex, France

Introduction and theoretical background

A method to perform alchemical free energy perturbation (FEP)  [21,4,20,19,12,9,13,8] within NAMD has now been implemented. Within FEP, the difference in free energy between two states, $a$ and $b$, is expressed by:

$$\Delta A_{a \rightarrow b} = -k_B T \ \ln \left< \exp\left[-... ...}) - {\cal H}_a({\bf r}, {\bf p})} {k_B T}\right] \right>_a$$ (3)

wherein $k_B$ is the Boltzmann constant, $T$ is the kinetic temperature, and ${\cal H}_a({\bf r}, {\bf p})$ and ${\cal H}_b({\bf r}, {\bf p})$ are the Hamiltonians characteristic of states $a$ and $b$, respectively. $\left< \cdots \right>_a$ denotes an ensemble average over configurations representative of the initial state, $a$. In practice, the transformation between the two thermodynamic states is replaced by a series of transformations between non-physical, intermediate states along a pathway that connects $a$ to $b$. This pathway is characterized by a variable, referred to as ``coupling parameter'', [4,13,11] $\lambda$, that makes the free energy a continuous function of this parameter between $a$ and $b$:

$$\Delta A_{a \rightarrow b} = -k_B T \ \sum_{k = 1}^N \ln \le... ...cal H}({\bf r}, {\bf p}; \lambda_k)} {k_B T}\right] \right>_k$$ (4)

Here, $N$ stands for the number of intermediate states, or ``windows'' between the initial and the final states.

In a typical FEP setup involving the transformation of one chemical species into another one in the course of the simulation, the atoms in the molecular topology can be classified into three groups. A group of atoms that do not change during the simulation -- e.g.the environment, the atoms describing the initial state, $a$, of the system, and, last, the atoms that correspond to the final state, $b$, at the end of the alchemical transformation. The atoms representative of state $a$ do not interact with those of state $b$ throughout the entire molecular dynamics simulation. Such a setup, in which atoms of both the initial and the final states of the system are present in the molecular topology file -- i.e. the psf file -- is characterisitic of the so-called ``dual topology'' protocol. [2] The hybrid Hamiltonian of the system, which is a function of the coupling parameter $\lambda$, that smoothly connects state $a$ to state $b$, is calculated as:

$${\cal H}(\lambda) = {\cal H}_0 + \lambda {\cal H}_a + (1-\lambda) {\cal H}_b$$ (5)

where ${\cal H}_a$ is the Hamiltonian of the group of atoms representative of the initial state, $a$, and ${\cal H}_b$ characterizes the final state, $b$. ${\cal H}_0$ is the Hamiltonian for those atoms that do not undergo any transformation during the MD simulation.

In the present implementation of FEP in NAMD, we employ a hamiltonian scaling procedure as is done in the ``dual topology'' approach, [15] i.e.instead of scaling the non-bonded parameters of states $a$ and $b$, namely the net atomic charges, together with the Lennard-Jones parameters, as a function of the coupling parameter, $\lambda$:

$$q_i = \lambda_i \ q_i$$

$$\epsilon_{ij} = \sqrt{\lambda_i \ \epsilon_{ii} \ \lambda_j \ \epsilon_{jj}}$$ (6)

$$\sigma_{ij} = \frac{\lambda_i \ \sigma_{ii} + \lambda_j \ \sigma_{jj}}{2}$$

where $\lambda_i$ and $\lambda_j$ take the value of $\lambda$ or $(1-\lambda)$, depending upon state $a$ or $b$, to which $i$ or $j$ belong.

For instance, in a transformation involving the mutation of an alanine side chain into that of glycine, using the FEP, the topology of both the methyl group borne by the C$_\alpha$ in alanine, and the hydrogen of glycine co-exist throughout the simulation (see Figure 5).

Figure 5: Dual topology description for an alchemical simulation. Case example of the mutation of alanine into glycine. The lighter color denotes the non-interacting, alternate state.

The charge and Lennard-Jones parameters of the alanine and the glycine side chains are defined as a function of $\lambda$, in such a fashion that the interaction of the methyl group of alanine with the rest of the protein is effective at the beginning of the simulation, i.e.$\lambda$ = 0, while the glycine C$_\alpha$ hydrogen does not interact with the rest of the protein, and vice versa at the end of the simulation, i.e.$\lambda$ = 1. For intermediate values of $\lambda$, both the alanine and the glycine side chains participate in non-bonded interactions with the rest of the protein, scaled on the basis of the current value of $\lambda$. It should be emphasized that these side chains, however, do not interact with themselves. It is not necessary to explicitly exclude those atoms that are created from those that will be annihilated in the course of the FEP calculation.

It is also worth noting that the free energy calculation does not alter intramolecular potentials, i.e.bond stretch, valence angle deformation, torsions, etc..., during the simulation. In calculations targetted at the estimation of free energy differences between two states characterized by distinct environments -- e.g.a ligand, bound to a protein in the first simulation, and solvated in water, in the second -- as is the case for most free energy calculations that make use of a thermodynamic cycle, perturbation of intramolecular terms can be safely avoided. [5]

Implementation of free energy perturbation in NAMD

The procedure implemented in NAMD is particularly adapted for performing free energy calculations that split the $\lambda$ reaction path into a number of non-physical, intermediate states, or ``windows''. Seperate simulations can be started for each window. Alternatively, the Tcl scripting ability of NAMD can be employed advantageously to perform the complete simulation in a single run. An example making use of such script is supplied at the end of this user guide.

The following keywords can be used to control the alchemical free energy calculations.

• fep $<$ Is alchemical FEP to be performed? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Turns on parameter scaling and ensemble averaging for alchemical FEP.

• lambda $<$ Coupling parameter value $>$
  Acceptable Values: positive decimal between 0.0 and 1.0
  Description: The coupling parameter value determining the progress of the perturbation. The non-bonded parameters of the atoms vanishing in the course of the MD simulation are scaled by lambda, while the parameters of those atoms grown are scaled by (1-lambda).

• lambda2 $<$ Coupling parameter comparison value $>$
  Acceptable Values: positive decimal between 0.0 and 1.0
  Description: The lambda2 value corresponds to the coupling parameter to be used for sampling in the next window. The free energy difference between lambda2 and lambda is calculated. Through simulations at progressive values of lambda and lambda2 the total free energy difference may be determined.

• fepEquilSteps $<$ Number of equilibration steps in the window, before data collection $>$
  Acceptable Values: positive integer less than numSteps or run
  Default Value: 0
  Description: In each window fepEquilSteps steps of equilibration can be performed before ensemble averaging is initiated. The output also contains the data gathered during equilibration and is meant for analysis of convergence properties of the FEP calculation.

• fepFile $<$ pdb file with perturbation flags $>$
  Acceptable Values: filename
  Default Value: coordinates
  Description: pdb file to be used for indicating the FEP status for each of the atoms pertaining to the system. If this parameter is not declared specifically, then the pdb file containing the initial coordinates specified by coordinates is utilized for this information.

• fepCol $<$ Column in the fepFile that carries the perturbation flag $>$
  Acceptable Values: X, Y, Z, O or B
  Default Value: B
  Description: Column of the pdb file to use for retrieving the FEP status of each atom, i.e.a flag that indicates which atom will be perturbed in the course of the simulation. A value of -1 in the specified column indicates the atom will vanish during the FEP calculation, whereas a value of 1 indicates that the atom would grow.

• fepOutFreq $<$ Frequency of FEP energy output in time-steps $>$
  Acceptable Values: positive integer
  Default Value: 5
  Description: Every fepOutFreq number of MD steps, the output file fepOutFile is updated by dumping energies that are used for ensemble averaging. This variable could be set to 1 to include all the configurations for ensemble averaging. Yet, it is recommended to update fepOutFile energies with a higher frequency to avoid large correlation between consecutive configurations.

• fepOutFile $<$ FEP energy output filename $>$
  Acceptable Values: filename
  Default Value: outfilename
  Description: An output file named fepOutFile.fep, generated by NAMD, contains the FEP energies, dumped every fepOutFreq steps.

Examples of input files for running FEP alchemical calculations

The first example illustrates the use of Tcl scripting for running the alchemical FEP feature of NAMD:

fep		on  
fepfile		ion.fep
fepCol		X
fepOutfile	ion.fepout
fepOutFreq	5
fepEquilSteps	5000

set step 0.0
set dstep 0.1

while {$step <= 0.9} {
 lambda $step
 lambda2 [expr $step+$dstep]
 run  10000
 set step [expr $step+$dstep]
}

Here, the pdb file read by NAMD to retrieve the information about perturbed atoms is biotin.fep. The pertinent information is present in the X column. The output file of the free energy calculation is biotinr.fepout, in which energies are written every 5 steps. $\delta \lambda$, the width of the windows, is set to 0.1. 5000 MD steps are performed in each window to equilibrate the system. In this particular instance, the current value of $\lambda$ is controlled by the syntax set step. The FEP calculation is run until $\lambda$ reaches the value 0.9. In every window, 10000 MD steps are being performed.

In the second example, each $\lambda$-state is declared explicitly, avoiding the use of Tcl scripting:

fep		on  
fepfile		ion.fep
fepCol		X
fepOutfile	ion.fepout
fepOutFreq	5
fepEquilSteps	5000

lambda		0.0
lambda2		0.1
run		10000

lambda		0.1
lambda2		0.2
run		10000

$\vdots$

lambda		0.8
lambda2		0.9
run		10000

lambda		0.9
lambda2		1.0
run		10000

The FEP calculation is carried out from $\lambda$ = 0.0 to 0.9. In each new window, 10000 MD steps are performed.

Description of FEP simulation output

The fepOutFile contains electrostatic and van der Waals energy data calculated for $\lambda$ and $lambda2$, written every fepOutFreq steps. The column dE is the energy difference of the single configuration, dE_avg and dG are the instantaneous ensemble average of the energy and the calculated free energy at the time step specified in column 2, respectively. The temperature is specified in the penultimate column. On completion of fepEquilSteps steps, the calculation of dE_avg and dG is restarted. The accumulated net free energy change is output and the end of the simulation at each lambda value. The cummulative average energy dE_avg value may be summed using the trapezoidal rule to obtain an approximate TI estimate for the free energy change during the run.