Free energy estimation via Thermodynamic Integration

Timing: ~30min

In this tutorial, we won't run any calculations, because required computation time would be too long. We will discuss how to use the previously learned simulation techniques to perform thermodynamic integration, and will analyze precalculated trajectories.

Theory: vibrational free energy

While the potential energy surface is shaped by the 3$N$ coordinates (degrees of freedom) that describe a molecule, the free energy is a function of thermodynamic variables (temperature, pressure, entropy, volume, etc.). This is the quantity of fundamental interest for comparison with experiments and the one that rules all equilibrium properties of the system. In the harmonic oscillator approximation, the partition function can be written as a product of multiple independent 1D harmonic oscillators, and the free energy can be defined as

$$ F^c(\beta) \approx F_{h}^c(\beta)+E_{\text{total}}, $$

where $E_{total}$ is the total energy at the potential energy surface of the system, and the vibrational contribution to $F_h$ in classical mechanics is

$$ F_{h, \text{vib}}^c(\beta) = \frac{1}{\beta} \sum_{i=1}^{3N-6} \text{ln} (\beta\hbar\omega_i), $$

where $\omega_i$ are the frequencies of the normal modes of vibration of the molecule, $\beta = 1/k_B T$, and the sum runs over all modes except the ones corresponding to translations and rotations. ¹ We learned in exercise 5 how to calculate normal modes of a molecular system.

There are limitations to the use of the harmonic approximation. It is not expected to be valid at all temperatures: for example, in simple Morse-like potentials, when atoms start to explore higher regions of the potential well, it cannot be well approximated by a parabola. Anharmonic vibrational levels are more closely spaced than what is estimated by the harmonic approximation, as can be pictorically seen in Figure 1. One way to obtain anharmonic corrections to the harmonic free-energy will be the subject of this exercise.

Morse potential

Figure 1. Schematic drawing of electronic and vibrational states of a hypothetical molecule. The Morse-like curves in black and pink represent the electronic ground-state and the first excited state, respectively. The levels drawn in each of them correspond to the possible vibrational frequencies $v_{0n}$, $v_{1n}$. The dashed parabola-like blue line represents the harmonic approximation, with equally spaced vibrational levels.

Theory: thermodynamic integration

In the canonical (NVT) ensemble, one can write

\[\label{Fbeta_differ} \frac{\partial \beta F}{\partial \beta} = -\frac{\partial \text{ln}Q}{\partial \beta} = -\frac{1}{Q}\frac{\partial Q}{\partial \beta} = \frac { \int \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta) e^{-\beta \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta)}d\textbf{p}^N d\textbf{r}^N } { \int e^{-\beta \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta)}d\textbf{p}^N d\textbf{r}^N } \]

where $Q$ is the partition function and $\mathcal{H}$ is the Hamiltonian. This leads to:

\[\label{eq:Fbeta_differ_cont} \frac{\partial \beta F}{\partial \beta} = \Big< \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta) \Big>_{\beta} \]

where $\langle..\rangle_{\beta}$ denotes ensemble average at inverse temperature $\beta$. We can then evaluate the free energy difference between inverse temperature $\beta_1$ and $\beta_2$ by integrating the equation above:

\[ \beta_2 F(\beta_2) - \beta_1 F(\beta_1) = \int_{\beta_1}^{\beta_2} \Big< \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta) \Big>_{\beta} d\beta \]

Or equivalently, $$ F(\beta_2) = \frac{\beta_1}{\beta_2}F(\beta_1) + \frac{1}{\beta_2}\int_{\beta_1}^{\beta_2} \Big< \mathcal{H}(\textbf{p}^N,\textbf{r}^N, \beta) \Big>_{\beta} d\beta $$

The averaged Hamiltonian $\langle \mathcal{H}\rangle_{\beta}$ can be written as a sum of two contributions, namely, $\langle \mathcal{H}\rangle_{\beta} = \langle U \rangle_{\beta}$ + $\langle K \rangle_{\beta}$, being $\langle U \rangle_{\beta}$ and $\langle K \rangle_{\beta}$ the potential and kinetic energy average values respectively. In order to improve the convergence of this equation, one can consider that the ensemble average of kinetic energy in an isolated non-linear polyatomic molecule is $\langle K \rangle_{\beta}=\frac{3N}{2\beta}$ if no degrees of freedom are constrained. We will constrain the center of mass motion, so we will have $\langle K \rangle_{\beta}=\frac{3N-3}{2\beta}$. Moreover, if the potential were fully harmonic, $\langle U \rangle_{\beta}$ would be $\frac{3N-6}{2\beta}$. In the anharmonic case we can write:

\[\label{eq:U_full} \langle \mathcal{H} \rangle_{\beta} = E_{\text{total}} + \underbrace{\frac{3N-6}{\beta} + \frac{3}{2\beta}}_{\langle \mathcal{H}^h \rangle_{\beta}} + \langle \mathcal{H}^{ah} \rangle_{\beta} \]

where $\langle H^{ah} \rangle_{\beta}$ is the ensemble average of the anharmonic contribution to the Hamiltonian. To be clear, we define

\[\label{eq:anharm} \langle \mathcal{H}^{ah} \rangle_{\beta} = \langle \mathcal{H} \rangle_{\beta} - \langle \mathcal{H}^h \rangle_{\beta} -E_{\text{total}} \]

Note that if different constraints are applied to the system, the equations above have to be modified accordingly.

We can now make an approximation that simplifies the calculation of these quantities. that our system can be approximated as harmonic at $\beta_1$, so that $F(T_1) \approx F_h^c(T_1)$. This approximation does its best when $T_1$ is a very low temperature. By realizing that $$ \int_{\beta_1}^{\beta_2} \langle \mathcal{H}^h \rangle_{\beta} d \beta = \beta_2 F_c^h(T_2) - \beta_1 F_c^h(T_1) $$ we can finally write $$ \begin{aligned} F^c(T_2) = E_{\text{total}} + F_c^h(T_2) + k_B T_2 \int_{\beta_1}^{\beta_2} \langle U_{ah} \rangle_{\beta} d\beta = E_{\text{total}} + F_c^h(T_2) - T_2 \int_{T_1}^{T_2} \frac{\langle U_{ah} \rangle_{T}}{T^2} dT \end{aligned} $$

Can you demonstrate this equation by yourself?

where a change of variables was performed in the last line.

Note

In this exercise, we provide trajectories for you to analyze because running all of them in a laptop would take many hours. However, if you want to run them by yourself to see how it goes, we also provide inputs that you can use.

For the analysis, instead of providing you a single blackbox program that performs all the post processing, we provide you several scripts (exercise_6/scripts/) that perform specific tasks, so that you understand the procedure step by step. You can find the description of the scripts at the end of this document.

Exercise Part 1

Go to the scripts/ directory for all scripts necessary to complete this task.
Use f_harm_c.py script to calculate harmonic, classical free energy as it is defined in the equation for $F_h^c$. This script takes as the input file containing vibrational frequencies in phonons/ directory and generates $F_h^c(T)$ data file that will be used in next steps. You have learned how to perform a vibrational analysis in the last exercise, and in fact you could just use that result here too.
Calculate the average potential energy of MD simulations data using average_U_anh.py. This script uses as its input the main directory containing MD runs (like the solutions/long_runs/fhi-aims/ directory) and produces data file containing the average anharmonic part of the potential energy together with its statistical error at different temperatures.
Use the averaged anharmonic part of the potential energy calculated in previous step to perform the integration using integrate_U.py.
Finally, use the files from the second and fourth step to calculate anharmonic free energy using f_anharm_c.py script.
You can now visualise the data using plot_Fh_vs_Fa_c.py and plot_Fh_vs_Fa_c_and_q.py.
Which differences between harmonic and anharmonic free energies do you see at this point? How does the difference behave with increasing temperature? How can you understand this? Do you spot anything unphysical with the curves that are plotted? Can you guess why this is happening?

Answer

The anharmonic free energy is below the harmonic one. This is a consequence of the anharmonicity of the potential, which effectivelly becomes "softer" (less curvature) in the anharmonic case, as compared to the harmonic one. This deviation becomes more and more pronounced as the amplitude of atomic motion increases with temperature, exploring more anharmonic parts of the potential.

However, there is something very odd with this curve: The entropy, given by $-dF/dT$ would be negative in this case! In fact if you plot the free energy curve up to 1000 K in the harmonic case, you will see that there is an inflection point and the free energy starts to have a negative slope (positive entropy) at higher temperatures. This is an intrisic failure of classical mechanics for the harmonic oscillator because for this molecule, there are only few frequencies of vibrations and none of them are very low in frequency (as you have seen in the previous exercise). This means that the classical assumption in which $k_BT \approx \hbar \omega$ is only fulfilled at quite high temperatures even for the lower frequencies of vibration. We will see how in quantum mechanics everything is well in the next part of the exercise.

plot_F_c

So far, we have been treating all nuclei as classical particles. In order to include their quantum nature we will implement a fairly simple approach, where the harmonic part of the free energy is taken as the one corresponding to quantum harmonic oscillators.

$$\label{eq:F_anharm_q} F^q(T) \approx F^c(T) - F_{h,vib}^c(T) + F_{h,vib}^q(T), $$

where

$$ F_{h,vib}^q (T) = \sum_{i=1}^{3N-6} \left[ \underbrace{\frac{\hbar \omega_i}{2}}_{\text{Zero Point Energy}} + \frac{1}{\beta} \ln\left(1-\exp^{-\hbar \omega_i \beta}\right)\right]. $$

It is worth pointing out that this is an approximate value, where the anharmonic part is still treated within classical approximation. One possibility to include the quantum contribution to the anharmonic free energy would be to calculate $\langle \mathcal H \rangle$ from path integral molecular dynamics. This approach can be easily simulated with i-PI, but we are not going to learn it in this tutorial. In case you are interested, please consult the tutors for more information.

Exercise Part 2

Use f_harm_q.py script to calculate harmonic, quantum free energy $F_{h,vib}^q$ as it is defined in the equation above. This script needs as the input file containing phonon frequencies and generates $F_h^q(T)$ datafile that will be used in the next steps. Use the results from the exercise 5, or find the eigenvalues in solutions/phonons/ directory.
Using previously calculated classical harmonic and anharmonic free energy, as well as the quantum harmonic free energy, calculate the anharmonic quantum free energy $F_{anh}^q$ using f_anharm_q.py script.
Finally you can plot: harmonic, anharmonic, classical and quantum free energies with plot_free_energies_c_and_q.py script.
Inspect the plotted results. What changes in the free energy between the classical and quantum descriptions of the nuclei, in this approximation? Can you tell where it comes from? Can you guess what would happen if we would include nuclear quantum effects also in the evaluation of the anharmonic component?

Answer

The main differences are: (1) The energy is no longer zero at $T=0$, which makes perfect sense for the quantum system because of the zero-point energy contributions; (2) The change in free energy with temperature always has a negative slope, which comes from applying the right physics to the system at hand.

If nuclear quantum effects were also accounted for in the anharmonic corrections, one could expect that the anharmonic corrections would be larger. See a discussion in Ref.².

plot_F_q

All the scripts are stored in exercise_6/scripts/. A more detailed description of each one is below:

f_harm_c.py
- Syntax: python f_harm_c.py input_eigenvalues output_file
- Input: The input_eigenvalues is the name of the file containing eigenvalues of the Hessian.
- Output: Returns a file containing two columns: temperature and the vibrational harmonic free energy within the classical approximation calculated according to the equation for $F_h^c$.
f_harm_q.py
- Syntax: python f_harm_q.py input_eigenvalues output_file
- Input: The input_eigenvalues is the name of the file containing eigenvalues of the Hessian.
- Output: Returns a file containing two columns: temperature and the vibrational harmonic free energy within the quantum approximation calculated according to the equation for $F_{h,vib}^q$.
average_U_anh.py
- Syntax: python average_U_anh.py input_md_folder output_file
- Input: The input_md_folder is the name of the directory with subdirectories containing results of molecular dynamic simulations. Function automatically reads name of those subdirectories as the temperature at which the MD simulation is held.
- Output: Returns the file with three columns: temperature, the averaged anharmonic part of the potential energy and its statistical error.
integrate_U.py
- Syntax: python integrate_U.py input_average_U output_file
- Input: input_average_U - the file with potential energy at given temperatures.
- Output: Returns the file with three columns: temperature, integrated anharmonic part of the potential energy and its statistical error.
f_anharm_c.py
- Syntax: python f_anharm_c.py input_f_harm_c input_U_integrated output_file
- Input: 2 files as input, 1 for output:
  input_f_harm_c containing classical harmonic free energy,
  input_U_integrated containing integrated potential energy,
  output_file.
- Output: Returns the file with 3 columns: temperature, anharmonic classical free energy and its statistical error. The anharmonic classical free energy is calculated according to the equation for $F^c_{anh}$.
f_anharm_q.py
- Syntax: python f_harm_q.py input_f_harm_c input_f_harm_q input_f_anhharm_c output_file
- Input: 3 input files and 1 for output:
  input_f_harm_c containing classical harmonic free energy,
  input_f_harm_q containing quantum harmonic free energy,
  input_f_anhharm_c containing anharmonic classical free energy,
  output_file.
- Output: Returns a file with three columns: temperature, anharmonic quantum free energy and its statistical error. The anharmonic quantum free energy is calculated according to the equation for $F_{anh}^q$.
plot_Fh_vs_Fa_c.py
- Syntax: python plot_Fh_vs_Fa_c.py input_f_harm_c input_f_anharm_c
- Input: 2 files:
  input_f_harm_c containing harmonic classical free energy,
  input_f_anharm_c containing anharmonic classical free energy.
- Output: Returns the plot of the harmonic and anharmonic free energy, both in classical approximation.
plot_Fh_vs_Fa_c_and_q.py
- Syntax: python plot_Fh_vs_Fa_c_and_q.py input_f_harm_c input_f_anharm_c input_f_harm_q input_f_anharm_q
- Input: 4 filenames for input:
  input_f_harm_c containing harmonic classical free energy,
  input_f_anharm_c containing anharmonic classical free energy,
  input_f_harm_q containing harmonic quantum free energy,
  input_f_anharm_q containing anharmonic quantum free energy.
- Output: Plots the harmonic and anharmonic free energy in both classical and quantum approximations.

One can also include the contributions to the partition function and the free energy arising from rotations and translations (you can find possible approximations for these quantities in textbooks of statistical mechanics, e.g. see D. McQuarrie, Statistical Mechanics, University Science Books, 1st. ed., 2000), and in real systems these may even be non-separable. In this exercise we focus on vibrational contributions just for the sake of simplicity and because they represent the largest contribution to the free energy. ↩
M. Rossi, "Progress and challenges in ab initio simulations of quantum nuclei in weakly bonded systems", J. Chem. Phys. 154, 170902 (2021) ↩