Thermodynamic integration with harmonic reference

The Helmholtz free energy ($A$) of a fully interacting system (1) can be expressed in terms of that of system harmonic in Cartesian coordinates (0,$\mathbf{x}$) as follows

$A_{1} = A_{0,\mathbf{x}} + \Delta A_{0,\mathbf{x}\rightarrow 1}$

where $\Delta A_{0,\mathbf{x}\rightarrow 1}$ is anharmonic free energy. The latter term can be determined by means of thermodynamic integration (TI)

$\Delta A_{0,\mathbf{x}\rightarrow 1} = \int_0^1 d\lambda \langle V_1 -V_{0,\mathbf{x}} \rangle_\lambda$

with $V_i$ being the potential energy of system $i$, $\lambda$ is a coupling constant and $\langle\cdots\rangle_\lambda$ is the NVT ensemble average of the system driven by the Hamiltonian

$\mathcal{H}_\lambda = \lambda \mathcal{H}_1 + (1-\lambda)\mathcal{H}_{0,\mathbf{x}}$

Free energy of harmonic reference system within the quasi-classical theory writes

$A_{0,\mathbf{x}} = A_\mathrm{el}(\mathbf{x}_0) - k_\mathrm{B} T \sum_{i = 1}^{N_\mathrm{vib}} \ln \frac{k_\mathrm{B} T}{\hbar \omega_i}$

with the electronic free energy $A_\mathrm{el}(\mathbf{x}_0)$ for the configuration corresponding to the potential energy minimum with the atomic position vector $\mathbf{x}_0$, the number of vibrational degrees of freedom $N_\mathrm{vib}$, and the angular frequency $\omega_i$ of vibrational mode $i$ obtained using the Hesse matrix $\underline{\mathbf{H}}^\mathbf{x}$. Finally, the harmonic potential energy is expressed as

$V_{0,\mathbf{x}}(\mathbf{x}) = V_{0,\mathbf{x}}(\mathbf{x}_0) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_0)^T \underline{\mathbf{H}}^\mathbf{x} (\mathbf{x} - \mathbf{x}_0)$

Thus, a conventional TI calculation consists of the following steps:

1. determine $\mathbf{x}_0$ and $V_{0,\mathbf{x}}(\mathbf{x}_0)$ in structural relaxation
2. compute $\omega_i$ in vibrational analysis
3. use the data obtained in the point 2 to determine $\underline{\mathbf{H}}^\mathbf{x}$ that defines the harmonic forcefield
4. perform NVT MD simulations for several values of $\lambda \in \langle0,1\rangle$ and determine $\langle V_1 -V_{0,\mathbf{x}} \rangle$
5. integrate $\langle V_1 -V_{0,\mathbf{x}} \rangle$ over the $\lambda$ grid and compute $\Delta A_{0,\mathbf{x}\rightarrow 1}$