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 (67), 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 [58].

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}}$$ (69)

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, vdW_dU/dl, and bond_dU/dl values reported in alchOutFile are the derivatives of the internal energy with respect to the scaling factors for each interaction type (i.e. electrostatics, etc.). dU/dl values are locally averaged over the last alchOutFreq steps. Cumulative averages for each component are reported alongside in the _avg columns.

The electrostatic, vdW, and bond values 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 / BOND LAMBDA listed for each partition after the title.

New as of version 2.12: The output in alchOutFile has been extensively revised and now more closely matches the NAMD standard output. Additional accounting for bonded term scaling is now also included.

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 9).