Description of a free energy calculation output

Free Energy Perturbation

When running FEP, the alchOutFile contains electrostatic and van der Waals energy data calculated for alchLambda and alchLambda2, written every alchOutFreq 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. Upon completion of alchEquilSteps steps, the calculation of dE_avg and dG is restarted. The accumulated net free energy change is written at each lambda value and at the end of the simulation.

Whereas the FEP module of NAMD supplies free energy differences determined from equation (54), the wealth of information available in alchOutFile may be utilized profitably to explore different routes towards the estimation of $\Delta A$. Both BAR and SOS methods, which combine advantageously direct and reverse transformations to improve convergence and accuracy of the calculation, represent relevant alternatives to brute-force application of the FEP formula [38].

Within the SOS framework, the free energy difference between states $\lambda_i$ and $\lambda_{i+1}$ is expressed as:

$$\exp(-\beta \Delta A_{i \rightarrow i+1}) = \frac{\displaystyle \... ...cal H}({\bf x}, {\bf p}_x; \lambda_{i+1}) \right] \right\} \right\rangle_{i+1}}$$ (56)

and can be readily used with the statistical information provided by the forward and the backward runs.

Thermodynamic Integration

When running TI free energy calculations, the elec_dU/dl and vdW_dU/dl values reported in tiOutFile are the derivatives of the internal energy with respect to $\lambda$ -- i.e.  $\frac{\partial U}{\partial\lambda}$ for electrostatics and, van der Waals, respectively. dU/dl values are averages over the last tiOutFreq steps. Cumulative averages for each component are reported alongside in the _avg columns.

The electrostatics and vdW are separated following a partition scheme -- that is, the ``appearing'' and the ``disappearing'' atoms are accounted for separately. ``Partition 1'' contains those atoms whose interactions are switched up as $\lambda$ increases -- i.e. flagged with 1 in the alchFile. ``Partition 2'' represents those atoms whose interactions are switched down as $\lambda$ increases -- i.e. flagged with -1. $\Delta A$ values for each component are obtained by integrating from $\lambda = 0$ to $1$ using the respective ELEC / VDW LAMBDA listed for each partition after the title.

Analysis is handled by the NAMD_ti script, available from

http://www.ks.uiuc.edu/Research/namd/utilities/

Although the output format of NAMD_ti.pl may appear to lend itself easily to interpretation of the individual contributions to the free energy total (elec and vdW for each partition), this is rarely appropriate as these values are path-dependent. For example, an output such as

          |-----------------------------------------------|
          |         |    elec   |    vdW    |   Subtotal  |
          |-----------------------------------------------|
          | Part. 1 |   -0.5748 |   -6.3452 |    -6.9200  |
          | Part. 2 |    0.5391 |    4.9387 |     5.4778  |
          |-----------------------------------------------|
          | Subtotal|    0.6048 |    0.3293 |   -12.3978  |
          |-----------------------------------------------|
          Total deltaG for transition lambda 0 -> 1: -12.3978

may encourage interpretations along the lines of "the free energy for switching on the van der Waals interactions for the atoms of partition 1 was -6.35kcal/mol". This is only correct in the very narrow context of the simulation setup and parameters used in this case and is not informative in a broader sense.

Figure: Sample TI data ( $log(\left <\frac {\partial U}{\partial \lambda }\right >)$ against $\lambda$). The blue shaded area shows the integral with fine sampling close to the end point. The red area shows the difference when $\lambda$ values are more sparse. In this example, insufficient sampling before $\lambda$ $\simeq$0.1 can result in a large overestimation of the integral. Beyond $\simeq$0.2, sparser sampling is justified as dE/d$\lambda$ is not changing quickly.

The choice of $\lambda$ values will depend on the application, but in general it is important to examine the shape of the curve to ensure that sampling is adequate to give a good estimate of the integral. In particular, it will be necessary to sample more finely towards the end points in order to accurately account for the strong repulsive van der Waals forces encountered when inserting particles into a system (see Figure 10).