Thermodynamic integration calculations
VASP supports three approaches to thermodynamic integration (TI):
- TI between any two states using VCAIMAGES[1].
- TI with a harmonic solid or ideal gas as a reference state using SCALEE[1].
- 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 [1]. Caution: the tag ISPECIAL=0 used in that reference is not valid anymore, instead the tag PHON_NSTRUCT=-1 is used).
TI using VCAIMAGES
The main tags for this type of calculation are VCAIMAGES and NCORE_IN_IMAGE1. This calculation requires a parent directory and two subdirectories named 01 and 02 (two images). All thre directories must contain VCAIMAGES and NCORE_IN_IMAGE1 in the INCAR file with identical values in each of them. The subdirectories must contain identical POSCAR and POTCAR files. The rest of the files can be different. For example one can do a TI from a smaller number of k-points to a larger number (different KPOINTS file) or a TI from a machine learning force field (MLFF) to ab-initio. The files in the parent directory must be identical to the files of the 01 directory, with exception of the INCAR file. The parent directory only requires to have the VCAIMAGES and NCORE_IN_IMAGE1 in it. The INCAR tags concerning the motion atoms (e.g. IBRION, ISIF, POTIM, MDALGO, etc.) need to be identical in both subdirectories. The parameters concerning the calculation of the forces can be different (different parameters for electronic minimization, MLFF vs. ab-initio, etc.).
Similar to nudged elastic band calculations VASP will carry out two distinct calculations in both directories according to the parameters set in both directories. VCAIMAGES is used to linearly combine the forces from both directories and the combined force is used for the update of the atomic positions. The value of VCAIMAGES is the weight for the forces from directory 01 and (1-VCAIMAGES) is the weight for 02. So the value of VCAIMAGES must be between 0 and 1.
NCORE_IN_IMAGE1 sets the number of cores used for the calculation in subdirectory 01. The remaining available cores are used for the calculation ins subdirectory 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 (MD_ALGO must be 2). |
TI using SCALEE
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.
Two possible options are available for the 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 [1] 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.
TI using TILAMBDA
The TI calculations in internal coordinates are performed in NVT ensemble using any thermostat available in VASP. The coupling parameter [math]\displaystyle{ \lambda }[/math] is defined by setting the parameter TILAMBDA in the INCAR file. The set of internal coordinates used in the TI calculation are defined via the ICONST file by setting the status to 3. The Hesse matrix [math]\displaystyle{ \mathbf{\underline{H}}^\mathbf{x} }[/math] is provided in the file HESSEMAT and its transformation into [math]\displaystyle{ \mathbf{\underline{H}}^\mathbf{q} }[/math] is performed by VASP. The potential energies of the system 1 and 0,[math]\displaystyle{ \mathbf{q} }[/math], needed to compute [math]\displaystyle{ \langle V_1 -V_{0,\mathbf{q}} \rangle }[/math] used as integrant in the TI expression for [math]\displaystyle{ \Delta A_{0,\mathbf{q} \rightarrow 1} }[/math], are written in the file REPORT in lines introduce by a string "e_ti>"