Thermodynamic integration calculations

Thermodynamic integration (TI) is an approach where two different calculations (e.g., different k-point meshes, machine learning force fields vs. ab-initio) are run in parallel with linear combination of the forces and stresses. This creates allows controlled integration between the two, from a reference system to and interacting system. The free energy of a fully interacting system can be written as the sum of the free energy a non-interacting reference system and the difference in the free energy of the fully interacting system and the non-interacting system

[math]\displaystyle{ F_{1} = F_{0} + \Delta F_{0\rightarrow 1} }[/math].

Using thermodynamic integration the free energy difference between the two systems is written as

[math]\displaystyle{ \Delta F_{0\rightarrow 1} = \int\limits_{0}^{1} d\lambda \langle U_{1}(\lambda) - U_{0}(\lambda) \rangle_{\lambda} }[/math].

Here [math]\displaystyle{ U_{1}(\lambda) }[/math] and [math]\displaystyle{ U_{0}(\lambda) }[/math] describe the potential energies of a fully-interacting and a non-interacting reference system, respectively. The coupling strength of the systems is controlled via the coupling parameter [math]\displaystyle{ \lambda }[/math]. It is neccessary that the connection of the two systems via the coupling constant is reversible. The notation [math]\displaystyle{ \langle \ldots \rangle_{\lambda} }[/math] denotes an ensemble average of a system driven by the following classical Hamiltonian

[math]\displaystyle{ H_{\lambda}= \lambda H_{1} + (1-\lambda) H_{0} }[/math].

VASP supports three approaches to thermodynamic integration (TI):

• TI between any two states using VCAIMAGES.
• TI with a harmonic solid or ideal gas as a reference state using SCALEE.
• TI with a harmonic solid as a reference state using TILAMBDA.

Details on choosing the ensemble size and how to perform the integration are described in the main text and especially the supplementary information of reference . Caution: the tag ISPECIAL=0 used in that reference is not valid anymore, instead the tag PHON_NSTRUCT=-1 is used.

TI using VCAIMAGES

Main tags:

• VCAIMAGES
• NCORE_IN_IMAGE1.

Directory structure:

• A parent directory with two subdirectories: 01 and 02 (these represent two images).
• All three directories must include the VCAIMAGES and NCORE_IN_IMAGE1 tags in their INCAR files, with identical values.
• Subdirectories must contain identical POSCAR and POTCAR files. Other files (e.g., KPOINTS) can differ, enabling TI between distinct calculation setups (e.g., different k-point meshes, machine learning force fields vs. ab-initio).
• The parent directory's files should match those in 01, except for the INCAR, which only needs the VCAIMAGES and VCAIMAGES entries.
• Parameters controlling atomic motion (e.g., IBRION, ISIF, POTIM, MDALGO) must be identical in both subdirectories.
• Force calculation parameters may differ.

How it works:

• Similar to nudged elastic band calculations, multiple calculations are run (subdirectories 01 and 02).
• The VCAIMAGES tag linearly combines the forces and stresses:
  • Weight for 01: VCAIMAGES
  • Weight for 02: 1 - VCAIMAGES
  • Thus, VCAIMAGES must be between 0 and 1.
• NCORE_IN_IMAGE1 sets the number of cores used for 01; remaining cores are used for 02.

Warning: Both directories must contain the same positions after applying the forces and the thermostat contributions. So thermostats containing random numbers cannot be used (MDALGO must be 2).

TI using SCALEE

Main tags:

• SCALEE: Sets the coupling parameter [math]\displaystyle{ \lambda }[/math] and determines the Hamiltonian used in the calculation.

Directory structure:

• Thermodynamic integration (TI) with SCALEE is performed in a single directory.
• Optionally, the calculation can read the file DYNMATFULL if present (see details below).

How it works:

• The main control tag is SCALEE, which sets the coupling parameter [math]\displaystyle{ \lambda }[/math] and determines the Hamiltonian used in the calculation.

Available options for reference system

Ideal gas: By default the thermodynamic integration is carried out from the ideal gas to the fully interacting case (in the case when no DYNMATFULL is present in the calculation folder). Usually the Stirling approximation is used for the free energy of the ideal gas written as

[math]\displaystyle{ F = -\frac{1}{\beta} \mathrm{ln} \left[ \frac{V^{N}}{\Alpha^{3N} N!} \right] }[/math]

where [math]\displaystyle{ V }[/math] is the volume of the system, [math]\displaystyle{ N }[/math] is the number of particles in the system and [math]\displaystyle{ \Alpha }[/math] is the de Broglie wavelength. The Stirling approximation applies in principle only in the limes of infinitely many particles. In reference the exact ideal gas equation was used since it helped to speed up the convergence of the final free energy of liquid Si with respect to the system size.

Harmonic solid: If the file DYNMATFULL exists in the calculation directory and SCALEE[math]\displaystyle{ \ne }[/math]1, the second order Hessian matrix is added to the force and thermodynamic integration from a harmonic model to a fully interacting system is carried out. The DYNMATFULL file stores the eigenmodes and eigenvalues from diagonalizing the dynamic matrix. This file is written by a previous calculation using the INCAR tags IBRION=6 and PHON_NSTRUCT=-1. This calculaion runs in a single folder. It optionally reads in a DYNMATFULL file in the calculation directory (for more details see below). The tag SCALEE sets the coupling parameter [math]\displaystyle{ \lambda }[/math] and hence controls the Hamiltonian of the calculation. By default SCALEE=1 and the scaling of the energies and forces via the coupling constant is internally skipped in the code. To enable the scaling SCALEE[math]\displaystyle{ \lt }[/math]1 has to be specified.

TI using TILAMBDA

Main tags:

• TILAMBDA: Sets the value for the coupling parameter [math]\displaystyle{ \lambda }[/math].

Directory structure:

• Thermodynamic integration (TI) with TILAMBDA is performed in a single directory.
• Internal coordinates for the TI calculation are defined in the ICONST file, with the status set to 3.
• The Hessian matrix in Cartesian coordinates, [math]\displaystyle{ \mathbf{\underline{H}}^\mathbf{x} }[/math], must be provided in the HESSEMAT file. The calculation automatically transforms this matrix to internal coordinates, [math]\displaystyle{ \mathbf{\underline{H}}^\mathbf{q} }[/math].

How it works:

• The potential energies for both systems (system 1 and system 0) are computed in the internal coordinate representation, [math]\displaystyle{ \mathbf{q} }[/math].
• These values are used to evaluate the integrand [math]\displaystyle{ \langle V_1 - V_{0,\mathbf{q}} \rangle }[/math] in the TI expression [math]\displaystyle{ \Delta F_{0,\mathbf{q} \rightarrow 1} }[/math].
• The TI calculations are performed in the NVT ensemble using any available thermostat.

Output:

• The required energies are written in the REPORT file. Look for lines beginning with the string e_ti>.