Thermodynamic integration with harmonic reference

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

[math]\displaystyle{ A_{1} = A_{0,\mathbf{x}} + \Delta A_{0,\mathbf{x}\rightarrow 1} }[/math]

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

[math]\displaystyle{ \Delta A_{0,\mathbf{x}\rightarrow 1} = \int_0^1 d\lambda \langle V_1 -V_{0,\mathbf{x}} \rangle_\lambda }[/math]

with [math]\displaystyle{ V_i }[/math] being the potential energy of system [math]\displaystyle{ i }[/math], [math]\displaystyle{ \lambda }[/math] is a coupling constant and [math]\displaystyle{ \langle\cdots\rangle_\lambda }[/math] is the NVT ensemble average of the system driven by the Hamiltonian

[math]\displaystyle{ \mathcal{H}_\lambda = \lambda \mathcal{H}_1 + (1-\lambda)\mathcal{H}_{0,\mathbf{x}} }[/math]

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

[math]\displaystyle{ 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} }[/math]

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

[math]\displaystyle{ 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) }[/math]

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

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