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 <0,1>$ and determined $\langle V_{1}-V_{0,\mathbf {x} }\rangle$