Worked example: Phonon dispersion
Table of Contents
*Necessary files: crystal.o
To begin, you need to have the output of a CRYSTAL calculation in the folder. If you do not have that ready and don't know how to operate CRYSTAL, see the wiki's CRYSTAL page for more information. All the used example files are also available for download without needing to reproduce them yourself.
Reading the structure and creating supercells
Phonopy is operated from the command line with the command 'phonopy' and adding appropriate flags to tell Phonopy what to do. To get started, have the CRYSTAL output named crystal.o in the folder and enter the command
phonopy --crystal -v -d --dim="2 2 2"
Here we tell phonopy to invoke CRYSTAL mode (
--crystal), tell it to show extended output (
-v, stands for verbose), create the necessary displacements of atoms (
-d, displacements) and do all this in a supercell that is of size 2x2x2 the primitive cell of Cu2O. You can also feed Phonopy an arbitrarily named CRYSTAL output with
-c your_output.o The choice of 2x2x2 is chosen to be computationally effective, see Useful tips 3.
The output printed in the terminal (see below) will show the lattice parameters and atom positions of the supercell as well as say "disp.yaml and supercells have been created." The folder now contains the input files for all the supercells that need to be calculated.
Figure 1. Output of Phonopy for Cu2O when creating 2x2x2 sized supercells. Figure: Jarno Linnera.
Performing the force calculations
By default, the prepared .d12-files contain only the keywords "GRADCAL", "EXTERNAL" and "TOLDEE 10", telling CRYSTAL to perform force calculations on a geometry that is specified in an .ext-file and to use a strict SCF energy convergence threshold. Information about the basis set and DFT parameters should then be added to the CRYSTAL input files. This can be automated, as the interface also looks for a file called TEMPLATE in the folder. By having the basis set and DFT parameters in the file TEMPLATE, they are added in all the supercell-***.d12-files which are then ready to be calculated right off the bat. The TEMPLATE file used in this example can be found here.
The supercell calculation outputs at the PBE0/TZVP level are provided in the Optional files if you do not wish to calculate them yourself. Please note that if you use the output files provided here, they differ from the ones that would be produced by the TEMPLATE file given here. The example outputs were produced and calculated with older versions of the CRYSTAL-Phonopy interface and CRYSTAL.
When all the force calculations have succesfully finished (the supercell without a number does not need to be calculated), we feed them to Phonopy with the line
phonopy -v --crystal -f supercell-*o
Phonopy reads the CRYSTAL output files and creates the file FORCE_SETS.
[jarno@wihuri wiki_test]$ phonopy -v --crystal -f supercell-*o _ _ __ | |__ ___ _ __ ___ _ __ _ _ | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | | | |_) | | | | (_) | | | | (_) || |_) | |_| | | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, | |_| |_| |___/ 1.13.2 1. Drift force of "supercell-001.o" to be subtracted 0.00000000 0.00000000 0.00000000 2. Drift force of "supercell-002.o" to be subtracted 0.00000000 0.00000000 -0.00000000 FORCE_SETS has been created. _ ___ _ __ __| | / _ \ '_ \ / _` | | __/ | | | (_| | \___|_| |_|\__,_|
Calculating force constants
You are now ready to calculate and plot the phonon dispersion. It is, however, a good idea to calculate and write the force constants to a file, so they can be read from the file whenever needed instead of computing again. This is accomplished with the command
phonopy -v --crystal --dim="2 2 2" --writefc
The output should contain the structure and the following block in the end.
---------------------------------------------------------------------------- Computing force constants... Max drift of force constants: 0.001669 (zz) -0.000000 (zx) Force constants are written into FORCE_CONSTANTS. ---------------------------------------------------------------------------- One of the following run modes may be specified for phonon calculations. - Mesh sampling (MP, --mesh) - Q-points (QPOINTS, --qpoints) - Band structure (BAND, --band) - Animation (ANIME, --anime) - Modulation (MODULATION, --modulation) - Characters of Irreps (IRREPS, --irreps) - Create displacements (-d) ---------------------------------------------------------------------------- _ ___ _ __ __| | / _ \ '_ \ / _` | | __/ | | | (_| | \___|_| |_|\__,_|
Phonopy also prints the maximum drift of force constants for given directions. It can be described as a measure of accuracy for the calculated force constants. If the drift is large, the results may not be reliable. The drift shown here for Cu2O is very small and warrants no action. If the system shows large drift, that may be overcome with symmetrization of the force constants. For more information, see the Phonopy home page.
Calculating phonon dispersion
As Phonopy suggests, we may use one of the listed run modes. We will be using the flag
--band for the phonon dispersion (phonon band structure). We also need to provide Phonopy with information on the reciprocal space paths along which the dispersion is calculated and plotted. This can be given in the command when running the program, but it is much more convenient to write the information in a file, such as band.conf used here. An example for the contents of band.conf for Cu2O is given below and a matching file is found here. The nomenclature of the band points is based on the paper of Setyawan and Curtarolo.
DIM = 2 2 2
BAND = 0 0 0 0 1/2 0 1/2 1/2 0 0 0 0 1/2 1/2 1/2 0 1/2 0 1/2 1/2 0 1/2 1/2 1/2
BAND_LABELS = $\Gamma$ X M $\Gamma$ R X M R
BAND_POINTS = 101
Here we specify the supercell dimensions with
DIM, the path in the reciprocal space with
BAND, names of the points where band segments meet with
BAND_LABELS and the number of points sampled between each band segment with
BAND_POINTS. Now the command
phonopy -v --crystal -p -s --factor=521.47083 --readfc band.conf
will tell Phonopy to calculate and plot the data (
-p), save it into a pdf-file (
-s), use a factor of 521.47083 for converting energy units to cm-1 (default is THz), read the force constants from the file rather than re-calculate them (
--readfc) and finally to read all the other information from band.conf. Phonopy again prints out the structure and the path in the reciprocal space, as shown below.
Max drift of force constants: 0.001669 (xx) -0.000000 (zx) Reciprocal space paths in reduced coordinates: [ 0.00 0.00 0.00] --> [ 0.00 0.50 0.00] [ 0.00 0.50 0.00] --> [ 0.50 0.50 0.00] [ 0.50 0.50 0.00] --> [-0.00 -0.00 0.00] [ 0.00 0.00 0.00] --> [ 0.50 0.50 0.50] [ 0.50 0.50 0.50] --> [-0.00 0.50 -0.00] [ 0.00 0.50 0.00] --> [ 0.50 0.50 0.00] [ 0.50 0.50 0.00] --> [ 0.50 0.50 0.50] _ ___ _ __ __| | / _ \ '_ \ / _` | | __/ | | | (_| | \___|_| |_|\__,_|
Plotting tools and useful tips
Now the folder should contain the files band.yaml, which contains the information about the phonon dispersion, and the plot in band.pdf. The .pdf should look like this.
Figure 52. Phonon dispersion of Cu2O. Figure: Jarno Linnera.
If you wish to modify the plot, you can write your own scripts to read the data from the files or use the plotting script bandplot. For usage and functionalities, enter
1. SCF convergence issues with displaced supercells
Even though the displacements are small, sometimes the supercells can suffer from SCF convergence issues. This can be, in most cases, salvaged by giving CRYSTAL a better initial guess for the wavefunction. A good candidate for it is obtained by doing a single-point calculation for the regular supercell, where no atoms are displaced (remove the keyword GRADCAL, forces are not needed). Copying the wavefunction file supercell.w as supercell-***.w for the needed cell(s) and inserting the keyword GUESSP in the SCF-section of supercell-***.d12 will result in CRYSTAL using the pre-calculated wavefunction as the initial guess.
2. Reducing the CPU time of the force calculations
Following the previous advice, one can use the wavefunction from the undisplaced supercell as an initial guess for all of the displaced supercells. Whether or not this is useful, depends on the situation. If the number of cells is large, shaving some cycles from the SCF procedure can result in a decent reduction of CPU time.
3. Ensuring convergence with respect to interatomic distances
The supercell size of 2x2x2 has been tested to produce the correct dispersion for Cu2O although the lattice constant is only roughly 8.6 Å. Generally, in order to ensure the force constants between a pair of atoms vanish within the supercell, distances of 10 - 15 Å are recommended. This is to make sure all interactions of relevant magnitude are taken into account. To make sure the results are converged, up to a 4x4x4 supercell would be preferable to check for Cu2O.
Expect something to go wrong if you don't know what you're doing
1. phonopy -v --crystal --dim="2 2 2" -d
2. Supercell calculations
3. phonopy -v --crystal -f supercell-*o
4. phonopy -v --crystal --dim="2 2 2" -p -s band.conf