A Practical Tutorial



This is a tutorial for the calculation of thermal ellipsoids using VASP and Phonopy.

If you have already calculated the density of phonon states with Phonopy, please check for imaginary modes. If there are only a few imaginary modes, you can directly go to "Calculate the Thermal Displacement Parameters".

1. Pre-Optimization of the Crystal Structure


We assume you're familiar with VASP. A typical INCAR for the optimization of the crystal structure could look as follows:
[code]PREC = Accurate
ADDGRID = .TRUE.
IBRION = 2
ENCUT = 500 ## should be adapted to your POTCAR
EDIFF = 1.0e-08
EDIFFG = 1.0e-06
ISIF=3
NSW=500
ISMEAR = 0
SIGMA = 0.1
ALGO=Normal
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
[/code]
It is very important to check the convergence of the structure and energy with a growing k-point mesh. Moreover, VASP optimizes the structure with a constant basis set, so that several restarts of structure optimizations are needed (see http://cms.mpi.univie.ac.at/vasp/guide/node161.html). For further details on VASP calculations, please visit http://www.vasp.at/.

2. Calculating the Forces with the Finite Displacement Method Implemented in Phonopy and VASP as a Calculator


First, you have to create supercells with and without displacements from your optimized structure with the help of Phonopy. Please visit http://phonopy.sourceforge.net for more information about Phonopy. It is reasonable to create supercells which have three nearly identical cell lengths in all three directions of at least 15 ??.
You can easily create 2??2??3-supercells from your [kbd]POSCAR[/kbd] like this:

[code]phonopy -d --dim="2 2 3"
[/code]

2a. Run Single Point Calculations for All Supercells



A typical INCAR for such a calculation could look like this:

[code]PREC = Accurate
IBRION = 2
ENCUT = 500
EDIFF = 1.0e-08
EDIFFG = 1.0e-06
ISMEAR = 0
NSW=0
SIGMA = 0.1
ALGO=Normal
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
[/code]
Don't forget to adapt your k-point mesh to the supercell size!

2b. Calculate the Force Constants


To do so, you need the [kbd]vasprun.xml[/kbd] files from all your single-point calculations! The force calculation is done as follows:

[code]phonopy -f disp-001/vasprun.xml disp-002/vasprun.xml disp-003/vasprun.xml[/code]or
[code]phonopy -f disp-{001..003}/vasprun.xml[/code]

3. Calculate the Density of Phonon States and Check For Imaginary Modes


First, create a file (e.g. [kbd]mesh.conf[/kbd]). This includes information about the supercell size ([kbd]DIM=2 2 3[/kbd], if 2??2??3 supercells were created) and about the q-point mesh ([kbd]MP= 3 3 2[/kbd]). A first guess for this q-point mesh would be the k-point mesh from the pre-optimization.

Thus, the file looks as follows:

[code]DIM=2 2 3
MP=3 3 2
[/code]

The calculation is run with [kbd]phonopy --dos --factor=521.47 mesh.conf[/kbd]. Due to the tag [kbd]factor[/kbd] the energy is converted to a wave vector in cm-1. Plot [kbd]total_dos.dat[/kbd] with a program of your choice. To control the smearing applied to this density of phonon states include the tag [kbd]--sigma=2[/kbd] (or a similar value).
If there are only a few imaginary modes you can continue calculating the thermal displacement parameters.

4. Calculate the Thermal Displacement Parameters



Create a file tdisp.conf. It looks as follows:

[code]TDISPMAT = .TRUE.
TMIN =0
TMAX = 300
TSTEP = 10
DIM = 2 2 3
MP = 3 3 2
[/code]

DIM and MP are the same as in [kbd]mesh.conf[/kbd] which means that one uses [kbd]DIM=2 2 3[/kbd] if 2??2??3 supercells were created and the q-point mesh is given by [kbd]MP= 3 3 2[/kbd] . [kbd]TMIN[/kbd] defines the temperature at which the thermal displacement parameter calculation begins, and [kbd]TMAX[/kbd] sets the temperature where the calculation ends. [kbd]TSTEP[/kbd] defines how many temperature steps are in the total range.
Run the file with [kbd]phonopy tdisp.conf[/kbd].
Test whether your thermal displacements converge with the q-points (MP-mesh). This is very important! Numerical problems can occur and they often do occur around the Gamma-point. The frequencies at the Gamma-point are used for the ADP calculation if one uses an odd q-point mesh (e.g. 3??3??3). The displacement parameters are saved in the file [kbd]thermal_displacement_matrices.yaml[/kbd].

If there are problems with the convergence of your thermal displacement parameters, you can include the following in the [kbd]tdisp.conf[/kbd] file:

[code]CUTOFF_FREQUENCY=0.1
[/code]

[kbd]CUTOFF_FREQUENCY[/kbd] sets a cutoff below which the frequencies are suppressed in the calculation of the thermal displacement parameters. This value is given in THz. A typical cutoff frequency lies between 0.1 and 0.15 THz. You can also include the tag [kbd]FC_SYMMETRY=1[/kbd]. [kbd]FC_SYMMETRY[/kbd] resymmetrizes the force constants calculated with VASP. This may help.
Again, check the convergence of the thermal displacement parameters with a growing q-point mesh size.

5. Convert the Thermal Displacement Parameters from the Cartesian to the Crystal Coordinate System



The thermal displacement parameters are calculated in the Cartesian coordinate system by Phonopy. For crystal structures with non-orthogonal cell vectors, one needs to convert the thermal ellipsoids from the Cartesian coordinate system to the ADPs used in ".cif"-files. We provide a custom-made program for this conversion that digests the [kbd]POSCAR[/kbd] used for the building of the super-cells in Phonopy and the [kbd]thermal_displacement_matrices.yaml[/kbd]. As a result you get a [kbd].cif[/kbd]-file with space group P1, hence, independent of the space group of your crystal structure, the cell parameters, the fractional coordinates and the ADPs. This program is now available in the Download section.

last modified: 2016-02-05