Anharmonicity and thermodynamic integration

If the displacements are not small, for example because the temperature is high, then the harmonic expansion of the potential energy function  may not be accurate enough and higher order terms would be needed (see Fig. ).

Figure: Potential energy $U(x)$ (black solid line) and its harmonic approximation (red dashed line).

In general it is difficult to take into account these terms, and in any case eventually every term included in the expansion becomes inaccurate for large enough displacements. In addition, for large displacements the assumption of constant underlying potential also breaks down, and so its Taylor expansion becomes meaningless.

An alternative procedure to compute the full free energy of the system goes under the name of thermodynamic integration. To explain how this works, let us consider the total potential energy of a crystal, and write it as:

$$U(\{{\bf r}\}) = U^h (\{{\bf r}\}) + U^\prime (\{{\bf r}\}),$$ (7.72)

where $U^h$ is the harmonic potential, defined as the sum of the first two terms in , and $U^\prime$ is the reminder. Note that here $U$ is the total potential energy function, which is not necessarily the one determined by the atoms in their equilibrium positions $\{{\bf r}^0\}$. As a consequence, $U^\prime$ is not necessarily equal to the $o(u^3)$ term in .

In general we do not have an analytic expression of $U^\prime$, but we might be able to sample it at some values of $\{{\bf r}\} = {\bf r}_1,\dots,{\bf r}_N$. Now define:

$$U_\lambda(\{{\bf r}\}) = U^h (\{{\bf r}\}) + \lambda U^\prime (\{{\bf r}\}),$$ (7.73)

where $\lambda$ is a parameter that we let vary between zero and one. For $\lambda = 0$ we have $U_{\lambda = 0} = U^h$ and for $\lambda = 1$ we have $U_{\lambda = 1} = U$. By switching the parameter $\lambda$ between zero and one we can therefore switch the potential energy function from $U^h$, which we know analytically, to $U$, which we may not fully know. The full free energy of the system can be written as:

$$F = F^h + (F - F^h),$$ (7.74)

where $F^h$ is the free energy of the harmonic system. The difference $F - F^h$ can be written as:

$$F - F^h = \int_0^1 d\lambda \frac{dF_\lambda}{d\lambda},$$ (7.75)

where $F_\lambda = -k_{\rm B}T \ln Z_\lambda$ is the free energy of the system with potential energy $U_\lambda$ and $Z_\lambda$ is the partition function, given by the usual sum of Boltzmann factors over the available microstates $i$:

$$Z_\lambda = \sum_i e^{-E^i_\lambda / k_{\rm B}T},$$ (7.76)

where $E^i_\lambda$ are the corresponding energies of the microstates. We therefore obtain:

$$\frac{dF_\lambda}{d\lambda} = -k_{\rm B}T \frac{1}{Z_\lambda}\fra... ...eft \langle \frac{\partial E_\lambda}{\partial \lambda} \right \rangle_\lambda,$$ (7.77)

where $\langle \cdot \rangle_\lambda$ denotes the (thermal) average over the canonical ensemble generated by the system with potential energy $U_\lambda$.

If the system is classical the partition function can be separated into a kinetic and a potential energy part, $Z_\lambda = Z_P Q_\lambda$, with

$$Q_\lambda = \frac{1}{V^N}\int_V d^3\{{\bf r}\} e^{-\beta U_\lambda(\{{\bf r}\}) },$$ (7.78)

where we have written $d^3\{{\bf r}\}$ as a shorthand for $d^3{\bf r}_1 \dots d^3{\bf r}_N$. Since only the potential energy term depends on $\lambda$, we have:

$$\frac{dF_\lambda}{d\lambda} = -k_{\rm B}T \frac{d \ln Q_\lambda}{... ...eft \langle \frac{\partial U_\lambda}{\partial \lambda} \right \rangle_\lambda.$$ (7.79)

We see that the canonical average in a classical system only depends on the space of configurations, as the momenta are integrated out and do not contribute. We can therefore write our final expression for  and :

$$F = F^h + \int_0^1 d\lambda \left \langle \frac{\partial E_\lambda}{\partial \lambda} \right \rangle_\lambda,$$ (7.80)

which for the classical case simplifies to:

$$F = F^h + \int_0^1 d\lambda \left \langle \frac{\partial U_\lambda}{\partial \lambda} \right \rangle_\lambda.$$ (7.81)

Because of the simple form  that we chose for $U_\lambda$ this simply amounts to calculating the canonical average of $U^\prime$. If we do not have an analytic expression for $U^\prime$ we cannot compute its average exactly, however, if we have some means to compute $U^\prime$ for some values of $\{{\bf r}\}$ distributed according to the canonical ensemble, then we can can obtain its average with a sampling error. As we mentioned in Sec. , for large systems it does not matter much if the average is taken in the canonical or the microcanonical ensemble, and so for simplicity we will discuss the latter first.

Note that although we described the thermodynamics integration procedure in the context of the anharmonic contribution to the free energy of a crystal, the procedure is completely general, and can be used to compute free energy differences between any two systems. The only requirement is that we must be able to integrate the quantity $dF_\lambda/d\lambda$ in the range $0\le \lambda \le 1$.