Free Energy Perturbation

Table of Contents

1 Introduction

This example shows how to analyze the data generated by free energy perturbation (FEP) simulations by using UWHAM or ST-SWHAM.

1.1 System

We calculate the solvation free energy of a water molecule in pure solvent by using the slow-growth method. The simulations were performed by using the GROMACS 5.1.2 and the TIP3P water model. In this example, we ran \(11\) independent parallel simulations for the pure solvent with a tagged water molecule fixed at \(r_{6d}\). Here \(r_{6d}\) includes the coordinates of the oxygen atom of the tagged water molecule \((x,y,z)\) and the orientation of the tagged water molecule \((\alpha,\beta,\gamma)\). Each simulation in the canonical ensemble follows the Hamiltonian (potential) function

\begin{equation} H_{i}(r_{6d},\{\vec{x}\})=U_{vv}(\{\vec{x}\})+E_{v\bar{v}}(r_{6d},\{\vec{x}\},\gamma_{i})\,, \end{equation}

where \(\{\vec{x}\}\) are the coordinates of the other water molecules and \(E_{v\bar{v}}(r_{6d},\{\vec{x}\},\gamma_{i})\) is the soft-core total interaction energy between the tagged water molecule and the other water molecules at the \(i^{th}\) \(\lambda\)-state. \(E_{v\bar{v}}(r_{6d},\{\vec{x}\},\gamma_{i})\) change from \(0\) to \(U_{v\bar{v}}(r_{6d},\{\vec{x}\})\) when \(\gamma_{i}\) changes from \(0\) to \(1\), where \(U_{v\bar{v}}(r_{6d},\{\vec{x}\})\equiv E_{v\bar{v}}(r_{6d},\{\vec{x}\},1)\) is the full interaction energy between the tagged water molecule and the other water molecules. \(U_{vv}(\{\vec{x}\})\) is the total interaction energy of all the other (untagged) water molecules with each other. The Coulombic interaction and van der Waals interaction between the tagged water molecule and the other water molecules were gradually turned off together using \(11\) \(\lambda\)-states. The chosen \(\lambda\) values are \(0.0\), \(0.1\), \(0.2\), \(0.3\), \(0.4\), \(0.5\), \(0.6\), \(0.7\), \(0.8\), \(0.9\), \(1.0\). All independent simulations were run at \(300\ K\) and each lasted \(20\ ns\). The potential energies of the system were recorded every \(0.1\ ps\). Therefore, there are \(11\) states in this system, and the simulation at each state generated \(200,000\) data points.

2 Files

2.1 Data Files

  • data.tgz: Unzipping this file generates a directory called "Data" which contains \(11\) data files, which were generated by GROMACS directly. The first \(31\) lines of each data file are comments beginning with either a "#" or a "@" character. After the comments there is a \(200,00 \times 13\) matrix. Each row of the matrix corresponds to an observation at this \(\lambda\)-state. The first column is a time index, and the second column is the derivative which can be used for thermodynamic integration (TI) calculations. These two columns are not used by the UWHAM analysis. The matrix element \(p[j][k+2]\) is the potential energy of the \(j^{th}\) observation calculated by using the Hamiltonian function of the \(k^{th}\) \(\lambda\)-state. Note that in the output of GROMACS, the potential energies has been shifted so that \(p[j][k+2]\) is zero when the \(k^{th}\) \(\lambda\)-state is the state at which this observation was observed during the FEP simulations. See the "-d" option of UWHAM for more discussions about the input for the UWHAM program.

3 UWHAM Analysis

Please read the manual of the UWHAM program before practicing.

  1. The first step is to construct the biasing matrix and create the input file for UWHAM. For the data files generated by FEP simulations with GROMACS, one can simply combine them together by

    cat Data/* > input4uwham.dat
    
  2. Run UWHAM

    /path of UWHAM/uwham.o -d input4uwham.dat -g 0 -t 300 -o 1 -q 5 -m 1 -e 1.0e-10
    

    The explanation for each option is:

    -d input4uwham.dat read the file "input4uwham.dat" for the biasing matrix
    -g 0 the input data file is a combination of GROMACS output files, NVT simulation
    -t 300 all the simulations were run at \(300\ K\)
    -o 1 the unit of the output data is \(kcal/mol\)
    -q 5 -m 1 run the self-consistent iteration 5 times, then switch to the Convex Newton-Raphson iteration
    -e 1.0e-10 stop after the changes the unitless free energy estimates are smaller than 1.0e-10

    If nothing goes wrong, the output of UWHAM program will show these lines:

    The free energy differences in unit of kcal/mol:
    0    2.12802    1.32699    -1.27085    -2.56216    -2.30276    -1.21943    0.319232    2.12711    4.10132    6.17941
    

    The free energy difference between the first and the last \(\lambda\)-states, \(-6.17941\ kcal/mol\) is the estimate of the solvation free energy for growing a water molecule in pure solvent.

4 ST-SWHAM Analysis

  1. As can be seen, it is easy to prepare the input file for UWHAM based on the FEP simulation data generated by GROMACS. However, the "input.dat" file analyzed by UWHAM with the option "-g" cannot be used by SWHAM program directly. We can construct a new input file by running

    for i in {01..11}
    do
    echo "# state$i.dat" >> input4swham.dat
    perl -wnal -e ' /^#/ or /^@/ or print ' Data/state$i.xvg > tmp.dat
    perl -wnal -e ' shift(@F); shift(@F); print(join " ", @F) ' tmp.dat >> input4swham.dat
    rm tmp.dat
    done
    

    In fact, this "input4swham.dat" can be used as the input for UWHAM. Run UWHAM

    /path of UWHAM/uwham.o -d input4swham.dat -i 2 -o 1 -m 1 -q 5 -e 1.0e-10
    

    Do you get the same free energy difference between the first and the last \(\lambda\)-states this time?

  2. Run ST-SWHAM to estimate the free energy differences

    /path of SWHAM/swham.o -d input4swham.dat -u 0 -i 2 -o 1 -t 300 -q 1 -n 10 -s 1-11 -g 1 -m
    

    The explanation for each option is:

    -d input4swham.dat read the file "input4swham.dat" for the information of observations
    -u 0 use the type 0 function implemented in the SWHAM program
    -i 2 -o 1 the input data are in units of \(kJ/mol\); and the output data are in units of \(kcal/mol\)
    -t 300 all the simulations were run at \(300\ k\).
    -q 1 -n 10 -m equilibrate the system for 1,000,000 cycles first; then run ST-SWHAM for 10,000,000 cycles
    -g 1 -m every 1,000,000 cycle, print out the free energy estimates
    -s 1-11 print out the free energy estimates for all states

    Do the ST-SWHAM results agree with the previous UWHAM results?

Author: Bin Zhang

Created: 2018-09-05 Wed 15:07

Validate