Anharmonicity

We have described in Sec.  the method of thermodynamic integration to compute free energy differences between two systems. We will now use this method to compute the free energy difference between the full DFT system and the harmonic system described in the previous section. Since anharmonicity is only going to be important in the high temperature limit, we can compute it assuming the classical approximation for the partition function.

If we write the potential energy function as $U = U^h + U'$, where $U^h$ is the harmonic potential that we have used in the previous section:

$$U^h(\{{\bf r}\}) = U_0 + \frac{1}{2}\sum_{i,j}{\bf u}_i \cdot {\bf\Phi}({\bf r}_i^0 - {\bf r}_j^0)\cdot {\bf u}_j,$$ (9.5)

where ${\bf u}_i = {\bf r}_i - {\bf r}_i^0$, with ${\bf r}_i^0$ the equilibrium position of particle $i$ and ${\bf r}_i$ its displaced position. The anharmonic contribution to the free energy is given by (see Eq. )

$$F' = \int_0^1 d\lambda \langle U' \rangle_\lambda,$$ (9.6)

with

$$\langle U' \rangle_\lambda = \frac{\int d^3 \{{\bf r}\} U'(\{{\bf... ...U_\lambda(\{{\bf r}\})}}{\int d^3 \{{\bf r}\}e^{-\beta U_\lambda(\{{\bf r}\})}}$$ (9.7)

and $U_\lambda = U^h + \lambda U'$.

We see that the quantities that we will need to compute are the canonical averages of $U'$, evaluated in the ensembles corresponding to $U_\lambda$ with $0\le \lambda \le 1$. We will therefore compute $\langle U' \rangle_\lambda$ at a number of chosen $\lambda$ values, and compute the integral  numerically. Canonical averages will be computed using molecular dynamics, with the temperature controlled using a Andersen thermostat. This will be done with an appropriately modified version of vasp.

To generate the harmonic potential $U^h$ we need to provide the equilibrium positions, via the SPOSCAR file, and the force constant matrix. The latter is written out by phon in the file HARMONIC by setting the variable:

LFORCEOUT =.TRUE.

To tell vasp to use the potential energy $U_\lambda$ and to write out $U'$ we need to set the following parameters in the INCAR file:

TIPO=inteharm         Switches on thermodynamic integration mode
R1=4.05   Radius of the force constant matrix. Must match value written in HARMONIC
XMIXREF=0.0          Ignore
XMIXHARM=1.0       Ignore
NCELL=8                 Number of primitive cells (a 2x2x2 supercell in this case)
XLAMBDA=0.5         Value of $\lambda$

And the following ones to perform a molecular dynamics run:

IBRION=0                     Switches on MD
NSW=1000                    Number of MD steps
LANDERSON=.TRUE.    Andersen thermostat
NANDERSON=50           Average number of steps in between thermalisations
TEBEG=300                  Temperature
POTIM=2.0                         Time step in fs
NBLOCK=10; KBLOCK=1     Output positions every NBLOCK steps
ISIF=0                                 Does not calculate pressure (can be expensive)
MAXMIX=40                        Ignore
NPAR = 2                            For parallel efficiency, test
LCHARG=.FALSE.                Ignore
LWAVE=.FALSE.                    Ignore
NWRITE=0                          Reduces output
IALGO = 48                          Faster diagonalisation algorithm

The number of MD steps determines the statistical accuracy of the run, and usually a few thousand steps are needed. To choose the time step the rule of thumb is that it should be shorter than about $1/30^{th}$ of the shortest period in the system. For the Al crystal the highest frequency is about 10 THz, which corresponds to a period of 100 fs, so a time step of 2 fs is a safe choice and could also be increased. The canonical ensemble is simulated by coupling the system with a Anderson thermostat, whereby the velocities of the atoms are periodically drawn from a Maxwellian distribution. These thermalisation events represent discontinuities in the Newtonian dynamics, but generate a set of configurations distributed according to the canonical ensemble. The intervals of time between thermalisations do not affect the the eventual generation of the canonical ensemble, but they may affect the time it take to have a sufficient representation of the ensemble. If one thermalises too often, then the particles don't move very much and sample the phase space slowly. If one thermalises too seldom, then energy transfers slowly between the system and the thermostat and the system behaves more like as in the microcanonical ensemble. Moreover, if the system is very non-ergodic, such as a closely harmonic system, energy is transferred between different modes only when the system is thermalised.

Note that to be consistent with the harmonic calculations the MD runs have to be performed with the same technical parameters, e.g. same k-point sampling etc.. Now perform the simulations at a number of $\lambda$ points and compute the anharmonic contribution to the free energy.