In this example, we will calculate the phonon dispersion relations of Cu2O at the DFT-PBE0/TZVP level of theory.
Table of contents
To calculate the phonon dispersion relations of Cu2O at the DFT-PBE0/TZVP level of theory, we need the optimized crystal structure of Cu2O at the same level of theory. For details, see the page Geometry optimization of Cu2O.
Phonopy reads in the optimized geometry from a CRYSTAL output file. There are several possible choices for the output file:
- Output from IR and Raman spectra calculation (recommended). The output includes dielectric tensor and Born effective charges, which are required to calculate the non-analytical term correction to phonon dispersion relations. If you don't need the Raman spectrum, you can save time by leaving out keyword
- Output from harmonic frequency calculation at the Γ-point (recommended only if you are sure that you will not need the non-analytical term correction).
- Output from geometry optimization (not recommended). It is always a good practice to check that the structure is a true local minimum based on the Γ-point frequencies before carrying out much more expensive phonon dispersion relation calculations.
Copy the CRYSTAL output file to a new directory and rename the file to
crystal.o. This is the default filename that Phonopy expects for a CRYSTAL calculation.
2. Phonon supercell dimensions
In the directory with the file
crystal.o, create a file called
DIM that contains the dimensions of the phonon supercell. Here, we use a rather small 2x2x2 phonon supercell and the contents of the
DIM file would be the following line:
2 2 2
3. TEMPLATE file for phonon supercell calculations
By default, Phonopy creates just the phonon supercell geometries and you need to manually complete the input files. If you provide a file called
TEMPLATE in the directorywhere you are running the phonon calculations, Phonopy will prepare complete input files automatically. The file must contain the basis set block for CRYSTAL input and the computational parameters block without final
END. Here is an example for Cu2O:
Pay attention on SHRINK
Note how the reciprocal space k-sampling has been halved in the TEMPLATE file from the original 8x8x8 to 4x4x4:
4 4 4
The reason is that the real-space unit cell is now a 2x2x2 phonon supercell, and therefore the reciprocal space k-mesh can be halved.
4. Creating phonon supercells
Phonopy is operated from the command line with the command
phonopy and appropriate flags. In the directory with the files
phonopy -v --crystal --dim="$(<DIM)" -d
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. 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.
5. Completing the phonon supercell input files
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 DFT-PBE0/TZVP level are provided in the Optional files above if you do not wish to calculate them yourself. Please note that if you use the output files provided here, they may 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.
6. Submitting the phonon supercell input files
The phonon supercells are typically fairly large and may require several gigabytes of memory per CPU. Here we use 3 GB per core:
jsub -np 12 -smp -mem 3G crystal supercell-001.d12
jsub -np 12 -smp -mem 3G crystal supercell-002.d12
7. Processing the phonon supercell outputs
When all the force calculations have succesfully finished, we feed them to Phonopy with the command
phonopy -v --crystal -f supercell-*o
Phonopy reads the CRYSTAL output files and creates the file FORCE_SETS.
8. 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.
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 manual.
9. Calculating phonon dispersion relations
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 phonon 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.
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 --readfc -p -s --factor=521.47083 band.conf
will tell Phonopy to calculate and plot the data (
-p), save it into a PDF (
-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 keywords from band.conf. Phonopy again prints out the structure and the path in the reciprocal space, as shown below.
10. 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 2. 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
11. Useful tips:
11.1. 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.
11.2. 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.
11.3. 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.
12. Summary of the commands
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