# Adiabatic switching and thermodynamic integration

he free-energy perturbation approach evokes a physical picture in which configurations sampled from the canonical distribution of state $$\cal A$$ are immediately switched'' to the state $$\cal B$$ by simply changing the potential from $$U_{\cal A}$$ to $$U_{\cal B}$$. Such instantaneous'' switching clearly represents an unphysical path from one state to the other, but we need not concern ourselves with this because the free energy is a state function and, therefore, independent of the path connecting the states. Nevertheless, we showed that the free-energy perturbation theory formula, eqn. (6), is only useful if the states $$\cal A$$ and $$\cal B$$ do not differ vastly from one another, thus naturally raising the question of what can be done if the states are very different.

The use of a series of intermediate states, by which eqn. (7) is derived, exploits the fact that any path between the states can be employed to obtain the free energy difference. In this section, we will discuss an alternative approach in which the system is switched slowly or adiabatically from one state to the other, allowing the system to fully relax at each point along a chosen path from state $$\cal A$$ to state $$\cal B$$, rather than instantaneously switching the system between intermediate states, as occurs in eqn. (7). In order to effect the switching from one state to the other, we will employ a common trick in the form of an external'' switching parameter, $$\lambda$$. This paramter is introduced by defining a new potential energy function

 $U({\textbf r}_1,...,{\textbf r}_N,\lambda) \equiv f(\lambda)U_{\cal A} ({\text r}_1,...,{\textbf r}_N)+ g(\lambda)U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)$ (8)

The functions $$f (\lambda)$$ and $$g (\lambda )$$ are referred to as switching functions, and they required to satisfy the conditions $$f(0) = 1, f (1) = 0$$, corresponding to the state $$\cal A$$, and $$g (0) = 0, g (1) = 1$$, corresponding to the state $$\cal B$$. Apart from these conditions, $$f (\lambda)$$ and $$g (\lambda )$$ are completely arbitrary. The mechanism embodied in eqn. (8) is one in which some imaginary external controlling influence (hand of God''), represented by the $$\lambda$$ parameter, starts the system off in state $${\cal A} (\lambda = 0 )$$ and slowly switches off the potential $$U_{\cal A}$$ while simultaneously switching on the potential $$U_{\cal B}$$. The process is complete when $$\lambda = 1$$, when the system is in state $$\cal B$$. A simple choice for the functions $$f (\lambda )$$ and $$g (\lambda )$$ is, for example, $$f (\lambda) = 1 - \lambda$$ and $$g (\lambda ) = \lambda$$.

In order to see how eqn. (8) can be used to compute the free energy difference $$A_{\cal A \cal B}$$, consider the canonical partition function of a system described by the potential of eqn. (8) for a particular choice of $$\lambda$$:

 $Q(N,V,T,\lambda) = C_N\int\;d^N{\textbf p}\;d^N{\textbf r}\; exp\left \{- \beta \left [\sum_{i=1}^N{{\textbf P}^2_i \over 2m_i}+ U({\textbf r}_1,...,{\textbf r}_N,\lambda)\right]\right\}$ (9)

This partition function leads to a free energy $$A (N, V, T, \lambda)$$ via

 $A(N,V,T,\lambda) = -kT\ln Q(N,V,T,\lambda)$ (10)

Recall, however, that the derivatives of the free energy with repsect to $$N$$, and $$V$$ and $$T$$ lead to the chemical potential, pressure and entropy, respectively. What does the derivative of the free energy $$A (N, V, T, \lambda )$$ with respect to $$\lambda$$ represent? According to eqn. (10)

 ${\partial A \over \partial \lambda} =-{kT\over Q}{\partial Q \over \partial \lambda} = -{kT \over Z}{\partial Z \over \partial \lambda }$ (11)

The reader should check that the expressions involving $$Q$$ and $$Z$$ are equivalent. Computing the derivative of $$Z$$ with respect to $$\lambda$$, we find

 $${kT \over Z}{\partial Z \over \partial \lambda}$$ $${kT \over Z}{\partial \over \partial \lambda}\int\;d^N{\textbf r}\;e^{-\beta U({\textbf r}_1,...,{\textbf r}_N,\lambda)}$$ $${kT \over Z}\int\;d^N{\textbf r}\;\left(-\beta {\partial U \over \partial \lambda}\right)e^{-\beta U({\textbf r}_1,...,{\textbf r}_N,\lambda)}$$ $$-\left<{\partial U \over \partial \lambda}\right>$$ (12)

Now, the free energy difference $$A_{\cal A \cal B}$$ can be obtained trivially from the relation

 $A_{\cal A\cal B}= \int_0^1 {\partial A \over \partial \lambda}d\lambda$ (13)

Substituting eqns. (11) and (12) into eqn. (13), we obtain the free energy difference as

 $A_{\cal A\cal B}= \int_0^1 \left<{\partial U \over \partial \lambda}\right>_{\lambda}d\lambda$ (14)

where $$\langle\cdots\rangle_{\lambda}$$ denotes an average over the canonical ensemble described by the distribution $$\exp[-\beta U({\textbf r}_1,...,{\textbf r}_N,\lambda)]$$ with $$\lambda$$ fixed at a particular value. The special choice of $$f (\lambda ) = 1 - \lambda$$ and $$g (\lambda ) = \lambda$$ has a simple interpretation. For this choice, eqn. (14) becomes

 $A_{\cal A\cal B}= \int_0^1\left_{\lambda}d\lambda$ (15)

The content of eqn. (15) can be understood by recalling the relationship between work and free energy from the second law of thermodynamics. If, in transforming the system from state $$\cal A$$ to state $$\cal B$$, an amount of work $$W$$ is performed on the system, then

 $W \ge A_{\cal A \cal B}$ (16)

where equality holds only if the transformation is carried out along a reversible path. Since reversible work is related to a change in potential energy, eqn. (15) is actually a statistical version of eqn. (16) for the special case of equality. Eqn. (15) tells us that the free energy difference is the ensemble average of the microscopic reversible work needed to change the potential energy of each configuration from $$U_{\cal A}$$ to $$U_{\cal B}$$ along the chosen $$\lambda$$-path. Note, however, that eqn. (14), which is known as the thermodynamic integration formula, is true independent of the choice of $$f (\lambda)$$ and $$g (\lambda )$$, which means that eqn. (14) always yields the reversible work via the free energy difference. The flexibility in the choice of the $$\lambda$$-path, however, can be exploited to design adiabatic switching algorithms of greater efficiency that can be achieved with the simple choice $$f (\lambda) = 1 - \lambda, g (\lambda ) = \lambda$$.

In practice, the thermodynamic integration formula is implemented as follows: A set of $$M$$ values of $$\lambda$$ is chosen from the interval $$[0, 1]$$, and at each chosen value $$\lambda _k$$, a full molecular dynamics or Monte Carlo calculation is carried out in order to generate the average $$\left<{\partial U \over \partial \lambda_k}\right>_{\lambda_k}$$. The resulting values of $$\left<{\partial U \over \partial \lambda_k } \right>_{\lambda_k}$$, $$k = 1, \cdots , M$$ are then substituted into eqn. (14), and the resulted is integrated numerically to produce the free energy difference $$A_{\cal A \cal B}$$. Thus, we see that the selected values $$\left \{ \lambda _k \right \}$$ can be evenly spaced, for example, or they could be a set of Gaussian quadrature nodes, depending on how $$A (N, V, T, \lambda )$$ is expected to vary with $$\lambda$$ for the chosen $$f (\lambda )$$ and $$g (\lambda )$$.

As with free-energy perturbation theory, the thermodynamic integration approach can be implemented very easily. An immediately obvious disadvantage of the method, however, is the same one that applies to eqn. (7): In order to perform the numerical integration, it is necessary to perform many simulations of a system at physically uninteresting intermediate values of $$\lambda$$ where the potential $$U({\textbf r}_1,...,{\textbf r}_N,\lambda)$$ is, itself, unphysical. Only $$\lambda = 0, 1$$ correspond to actual physical states and ultimately, we can only attach physical meaning to the free energy difference $$A_{\cal A \cal B}= A(N,V,T,1)-A(N,V,T,0)$$. Nevertheless, the intermediate averages must be accurately calculated in order for the integration to yield a correct result. The approach to be presented in the next section attempts to reduce the time spent in such unphysical intermediate states and focuses the sampling in the important regions $$\lambda = 0, 1$$.