7.4. SrTiO₃ : A self-consistent phonon example

This page explains how to run an self-consistent phonon (SCP) calculation using ALAMODE. The target materials is cubic SrTiO₃, which is a well-known perovskite solid exhibiting soft phonons.

The example input files are provided in example/SrTiO3/reference.

Let’s move to the example directory

$ cd ${ALAMODE_ROOT}/example/SrTiO3

1. Prepare force constants and BORNINFO

In this tutorial, we assume that the harmonic and anharmonic force constants (up to 4th order) are already calculated and only focus on the SCP calculation part. The tutorial for calculating anharmonic force constants is given in another page.

To start an SCPH calculation, please copy the force constant files and BORNINFO as follows:

$ cp reference/STO_anharm.xml.bz2 .
$ bunzip2 STO_anharm.xml.bz2
$ cp reference/BORN .

Note

The anharmonic force constants were computed using the following settings:

VASP, PBEsol, ENCUT = 550 eV
2x2x2 supercell of cubic SrTiO₃ (40 atoms)
Harmonic force constants were computed using ordinary least squares
Anharmonic force constants were estimated using LASSO from 40 training structures

For more details, please see the original paper.

2. Phonons within the harmonic approximation

Before moving onto an SCP calculation, let’s first confirm that the cubic STO is dynamically unstable within the harmonic approximation.

Please prepare an input file for band structure calculation (phband.in) as follows:

&general
 PREFIX = STO222_NA3
 MODE = phonons
 NKD = 3; KD = Sr Ti O
 FCSXML = STO_anharm.xml 
 NONANALYTIC = 3; BORNINFO = BORN
/

&cell
7.363
1.0 0.0 0.0 
0.0 1.0 0.0
0.0 0.0 1.0
/

&kpoint
1
G 0.0 0.0 0.0 X 0.5 0.0 0.0 51
X 0.5 0.0 0.0 M 0.5 0.5 0.0 51
M 0.5 0.5 0.0 G 0.0 0.0 0.0 51
G 0.0 0.0 0.0 R 0.5 0.5 0.5 51
R 0.5 0.5 0.5 M 0.5 0.5 0.0 51
/
 

Please run anphon and plot the band structure as

$ ${ALAMODE_ROOT}/anphon/anphon phband.in > phband.log
$ python ${ALAMODE_ROOT}/tools/plotband.py STO222_NA3.bands

It is clear that phonons are dynamically unstable at the \(\Gamma\) and R points.

3. Phonons within the self-consistent phonon method

Now, let’s include the finite-temperature effect using the SCP method.

Please copy phband.in to scph.in and edit scph.in as follows:

&general
 PREFIX = STO_scph2-2
 MODE = SCPH
 NKD = 3; KD = Sr Ti O
 FCSXML = STO_anharm.xml 
 NONANALYTIC = 3; BORNINFO = BORN
 TMIN =0; TMAX = 1000; DT = 50
/

&cell
7.363
1.0 0.0 0.0 
0.0 1.0 0.0
0.0 0.0 1.0
/

&kpoint
1
G 0.0 0.0 0.0 X 0.5 0.0 0.0 51
X 0.5 0.0 0.0 M 0.5 0.5 0.0 51
M 0.5 0.5 0.0 G 0.0 0.0 0.0 51
G 0.0 0.0 0.0 R 0.5 0.5 0.5 51
R 0.5 0.5 0.5 M 0.5 0.5 0.0 51
/

&scph
 SELF_OFFDIAG = 0
 MAXITER = 500
 MIXALPHA = 0.2
 KMESH_INTERPOLATE = 2 2 2
 KMESH_SCPH = 2 2 2
/

Here, the &scph field is added to the input. The important parameters for the SCP calculation are KMESH_INTERPOLATE and KMESH_SCPH. The former defines the \(\boldsymbol{q}\)-point mesh used to solve the SCP calculation and perform (inverse) Fourier transformation of dynamical matrices, while the latter defines the \(\boldsymbol{q}\)-point mesh used to compute the renormalization term associated with the loop self-energy diagram. Also, SELF_OFFDIAG option controls whether the off-diagonal components of the loop self-energy are considered or not. In the above example, we neglect the off-diagonal components to make the calculation faster. We will see the effect of the off-diagonal components later.

Important

All \(\boldsymbol{q}\) points generated by KMESH_INTERPOLATE should be commensurate with the supercell size (of the harmonic force constant calculation) because \(\Delta D(\boldsymbol{q}) = D_{\mathrm{SCPH}}(\boldsymbol{q}) - D_{\mathrm{HA}}(\boldsymbol{q})\) is used to interpolate dynamical matrix, for which \(D_{\mathrm{HA}}(\boldsymbol{q})\) should be exact within the harmonic approximation (HA). In the above example, we employ a 2x2x2 supercell, so we should use KMESH_INTERPOLATE = 2 2 2 (or KMESH_INTERPOLATE = 1 1 1, which is less accurate).

Important

The \(\boldsymbol{q}_1\) points generated by KMESH_SCPH can be denser than KMESH_INTERPOLATE, but the latter coarser grid must be a subset of the former denser grid. For example, when we use KMESH_INTERPOLATE = 2 2 2, we can set KMESH_SCPH = 2 2 2, KMESH_SCPH = 4 4 4, and so forth, but KMESH_SCPH = 3 3 3, KMESH_SCPH = 5 5 5, … are not allowed.

Now, let’s run anphon using MPI parallelization.

$ export OMP_NUM_THREADS=1
$ mpirun -np 4 ${ALAMODE_ROOT}/anphon/anphon scph.in > scph.log

The calculation finishes in ~2 minutes. When it is done, please check if the SCPH iteration reached convergence.

$ grep "conv" scph.log

We can see that the SCPH loop converged except for \(T = 0\) K. In the working directory, the following files are created:

STO_scph2-2.scph_dymat : Anharmonic corrections to the dynamical matrix, \(\Delta D(\boldsymbol{q}, T) = D_{\mathrm{SCPH}}(\boldsymbol{q}, T) - D_{\mathrm{HA}}(\boldsymbol{q})\), are stored. This file can be used to restart an SCP calculation with the same KMESH_* and SELF_OFFDIAG parameters.
STO_scph2-2.scph_dfc2 : Anharmonic corrrections to the second-order force constants, \(\Delta \Phi(\boldsymbol{r}(\ell), T) = M^{\frac{1}{2}} \Delta D(\boldsymbol{r}(\ell), T) M^{\frac{1}{2}}\), are stored, where \(M\) is the diagonal matrix whose elements are atomic masses, and \(\Delta D(\boldsymbol{r}(\ell), T) = N_{q}^{-1} \sum_{\boldsymbol{q}} \Delta D(\boldsymbol{q}, T) e^{-i\boldsymbol{q}\cdot\boldsymbol{r}(\ell)}\). This file is used to create effective second-order force constants via dfc2 (for ALAMODE FCSXML format) or scph_to_qefc.py (Quantum-ESPRESSO fc format).
STO_scph2-2.scph_bands : Self-consistent phonon band structures at various temperatures.

Let’s plot the finite-temperature band structures.

$ gnuplot
gnuplot> set terminal qt font "Helvetica,20"
gnuplot> set ylabel "Frequency (cm^{-1})"
gnuplot> unset key
gnuplot> plot for [col=3:17] "< awk '{if ($1!=0.0) print $0}' STO_scph2-2.scph_bands" u 2:col w l lt 1
gnuplot> replot for [col=2:16] "STO222_NA3.bands" u 1:col w l lt 2

The soft modes are stabilized by the SCP calculation and exhibit significant temperature dependences.

Note

If you are interested, please check the dependencies of the results on SELF_OFFDIAG and KMESH_SCPH values.

4. Thermodynamics from the self-consistent phonon method

We next compute various thermodynamic functions, including the vibrational free energy, mean square displacement (MSD), using the SCP calculation. Please copy scph.in to scph2.in and edit scph2.in as follows:

&general
 PREFIX = STO_scph2-2
 MODE = SCPH
 NKD = 3; KD = Sr Ti O
 FCSXML = STO_anharm.xml 
 NONANALYTIC = 3; BORNINFO = BORN
 TMIN =0; TMAX = 1000; DT = 50
 DELTA_E = 1.0
/

&cell
7.363
1.0 0.0 0.0 
0.0 1.0 0.0
0.0 0.0 1.0
/

&kpoint
2
10 10 10
/

&scph
 SELF_OFFDIAG = 0
 MAXITER = 500
 MIXALPHA = 0.2
 KMESH_INTERPOLATE = 2 2 2
 KMESH_SCPH = 2 2 2
/

&analysis
 PRINTMSD = 1
/

We use the same PREFIX as before so that the code can restart the calculation from STO_scph2-2.scph_dymat. The code can skip the most expensive part of the SCP calculation, so we run anphon without MPI:

$ ${ALAMODE_ROOT}/anphon/anphon scph2.in > scph2.log

When the calculation finishes, the following files are created in the working directory:

STO_scph2-2.scph_dos : Phonon density of states at various temperatures
STO_scph2-2.scph_msd : Temperature-dependent MSDs for each atom
STO_scph2-2.scph_thermo : Heat capacity and vibrational free energies at various temperatures

These files can be plotted and processed easily. For example, let’s plot the MSD of Sr, Ti, and O atoms using gnuplot.

$ gnuplot
gnuplot> set terminal qt font "Helvetica,20"
gnuplot> set xlabel "Temperature (K)"
gnuplot> set ylabel "MSD (Angstrom^{2})"
gnuplot> plot "STO_scph2-2.scph_msd" u 1:2 ti "Sr, x=y=z"
gnuplot> replot "STO_scph2-2.scph_msd" u 1:5 ti "Ti, x=y=z"
gnuplot> replot "STO_scph2-2.scph_msd" u 1:8 ti "O_1, x"
gnuplot> replot "STO_scph2-2.scph_msd" u 1:9 ti "O_1, y=z"

We can see that the MSD of an oxygen is large in the plane perpendicular to the nearest Ti-O bond direction.

In the STO_scph2-2.scph_thermo file, vibrational free energies within the SCP theory are stored, which can be used to predict phase stability and thermal expansitivy of anharmonic materials. In the file, 3-5 columns give vibrational free energies.

$ head -n 10 STO_scph2-2.scph_thermo

# Temperature [K], Cv [in kB unit], F_{vib} (QHA term) [Ry], F_{vib} (SCPH correction) [Ry], F_{total} [Ry]
        0.000000      0.000000e+00      2.163071e-02     -3.390871e-04      2.129162e-02
       50.000000      1.986285e+00      2.160337e-02     -3.628028e-04      2.124057e-02
000000      5.492636e+00      2.109534e-02     -4.684166e-04      2.062692e-02
000000      8.069936e+00      1.980037e-02     -6.283993e-04      1.917197e-02
000000      9.879608e+00      1.771384e-02     -8.226380e-04      1.689120e-02
000000      1.113299e+01      1.489935e-02     -1.043142e-03      1.385621e-02
000000      1.200912e+01      1.143021e-02     -1.285467e-03      1.014474e-02
000000      1.263357e+01      7.374577e-03     -1.546631e-03      5.827946e-03
000000      1.308878e+01      2.792002e-03     -1.824440e-03      9.675618e-04

The 3rd column is the QHA-like term and the 4th column is the correction term necessary for the vibrational free energy to satisfy the thermodynamic condition of \(S_{\mathrm{vib}}=-\frac{dF_{\mathrm{vib}}}{dT}\) . The 5th column is the sum of the 3rd and 4th columns, which gives the correct vibrational free energy within the SCP scheme. So, the values in the 5th column should be used for phase stability and thermal expansivity calculations.

Important

Since the SCP calculation did not converge at 0 K, the thermodynamic quantities at 0 K are incorrect and should not be used for any purposes.

7.4. SrTiO3 : A self-consistent phonon example

1. Prepare force constants and BORNINFO

2. Phonons within the harmonic approximation

3. Phonons within the self-consistent phonon method

4. Thermodynamics from the self-consistent phonon method

7.4. SrTiO₃ : A self-consistent phonon example