A Practical Guide to Frozen Phonon Calculations
Density functional theory (DFT) provides a way to get 0 Kelvin quantum mechanical energies and forces for atoms in crystals. However, if you want first principles thermodynamic data for crystals at finite temperatures, then you need more than just static or geometrical relaxation calculations. One important contribution to finite temperature thermodynamics of crystals is the energies and entropies of vibrations. These can be calculated by integrating a function called the phonon density of states. Phonons are quantizations of the normal modes of atomic oscillations in a crystal. We can use DFT and some post-processing software to calculate the vibrational modes of a crystal, calculate the phonon density of states, and integrate that to get the vibrational entropy, enthalpy, and heat capacity of the crystal. This guide will detail the steps required to successfully perform a set of phonon calculations with DFT. As the Vienna Ab-initio Simulation Package (VASP) is the DFT computational work-horse of the Wolverton Research Group, examples will take the form of VASP input and output files.
Theory of Phonons
Frozen Phonons
One way to implement the theory of crystal vibrations in a practical way is to explicitly calculate the forces between every atom in the crystal and construct the force constant matrix of the crystal. This force constant matrix then allows us to calculate the normal modes of at any particular wavevector, q. To calculate the forces caused by an atom i, we displace atom i, and then use DFT to calculate the forces on every atom using the Hellman-Feynman theorem. This method of calculating the force constant matrix by explicitly displacing atoms is called the Frozen-Phonon method. It is generally quicker and computationally cheaper than the Linear Response method, which utilizes density functional perturbation theory to calculate forces.
One drawback of the frozen-phonon method is that large supercells are needed to accurately calculate the force constant matrix. The periodic boundary conditions used in DFT calculations cause problems for the frozen phonon method. Displacing one atom in a small unit cell creates forces on all the atoms in the same unit cell, but also on the periodic images of these atoms.
Setting Up the Calculations
The theory of phonons assumes that the crystal is in equilibrium; the energy is minimized, and there are no forces on atoms and no stresses in the unit cell. So, the first step to implementing this theory in DFT is to relax the crystal and bring it to its equilibrium geometry. This should first be done for the primitive cell of the crystal structure of interest. However, accurate frozen phonon calculations require that (relatively) large supercells be used for calculations to ensure that the force constants between distant atoms go to zero. This means that relaxations should also be performed for any supercell that will be used in the frozen-phonon calculation. A good utility for manipulating VASP structure files, such as creating supercells from primitive cells, is convasp.
Relaxation
Relaxations should be carried out with respect to forces, not energies. A good initial guess (such as from a previous energy relaxation) will speed up the calculation, however. The VASP manual suggests some input parameters which will aid in finding an accurate geometry. These include:
EDIFFG = -1E-2 (-1E-3 if you really want those forces zeroed)
IBRION = 1
EDIFF = 1E-5
NELMIN = 6
ADDGRID = .TRUE.
ENAUG = 1500
MAXMIX = 40
The force relaxation calculations should be performed multiple times (say, at least twice), with fixed cutoff energy each time, as changes in volume result in changes in the basis set resulting in artificial Pulay stresses. The VASP manual has more details on this phenomenon. In addition to the relaxation calculations, a static calculation (NSW = 0, ISIF = 2) should be performed after the forces have been sufficiently converged. To speed up the frozen phonon calculations which will follow, this static calculation should write out the charge density to the CHGCAR file (LCHARG = .TRUE.). This will provide a good starting guess for the charge densities of the (many) calculations which still need to be done.
GoBaby Inputs
There are many codes which can be used to perform and automate the task of setting up and calculating frozen phonons. Some examples include, ATAT, Phonopy, Phonons (not free), and GoBaby. This guide will focus on GoBaby, which is widely used with the Wolverton and Ozolins research groups, but is not available for public use. Nevertheless, some of the advice and discussion here may carry over to using other phonon software packages.
The input file for GoBaby is apos.dat. It contains information on the structure to be calculated as well as the specific q-points at which phonons should be calculated. A bash script make_apos.sh has been attached to this document, which creates an apos.dat file from existing CONTCAR and POTCAR files. Note that this script was written for VASP 4.x files. In order to make it work for VASP 5.x, the 6th line of the CONTCAR file (the line with the symbols of the elements) must be deleted. An example apos.dat file for a 3x3x3 supercell of the primitive cell of PbTe follows:
0.0000000000 3.2814333620 3.2814333620 (First three lines are unit cell vectors)
3.2814333620 0.0000000000 3.2814333620
3.2814333620 3.2814333620 0.0000000000
3 3 3 (Fourth line has supercell dimensions)
2 1 (Fifth line has number of atom types {n} and number of q-points {m})
207.200 127.600 (Sixth line has atomic masses)
1 1 (Seventh line has number of atoms of each type)
0.0000000000 0.0000000000 0.0000000000 Pb (Next n lines have cartesian atomic coordinates.)
3.2814333659 3.2814333659 3.2814333659 Te (Atomic symbols not necessary.)
0.0 0.0 0.0 1.0 (Next m lines have q-points and weights)
In addition to the apos.dat file, GoBaby requires POTCAR, KPOINT, and INCAR files, all with the suffix _phon appended to them. In addition, CHGCAR_phon and WAVCAR_phon files can be used to create CHCCAR and WAVCAR files for every phonon calculation, speeding up each calculation.
The INCAR_phon file needs to be modified slightly for use with GoBaby. GoBaby uses VASP's MD mode to send each atom on a trajectory and calculate the forces at several points along this trajectory. This improves the fit of the force constant matrix relative to a fit using only one displacement. Correspondingly, the INCAR_phon file needs to be setup to perform an MD calculation. The tags to add are:
ENAUG = 1500
IBRION = 0
SMASS = -2
POTIM = 1.0
ISIF = 0
In addition, the NSW line should be removed, as GoBaby will add that back in later.
To construct the apos.dat file for an arbitrary supercell, I recommend the following procedure: i) Relax the primitive cell of the structure. ii) Create a supercell of the structure and relax it as well, creating a CHGCAR file from the relaxed structure. iii) Using convasp, convert the supercell back to a primivite cell. iv) Create an apos.dat file from the new primitive cell from convasp, putting the supercell size in the apos.dat file.
The benefits of this procedure are several. First, the forces will be converged for any supercell size. If the k-point mesh used for the primitive cell is not a multiple of the supercell size, then there will be residual forces on the supercell of the relaxed primitive cell. However, by re-relaxing the supercell, this potential problem is dealt with. Second, a CHGCAR file is created for the supercell. This will speed up the next step considerably, especially if the system has low symmetry and many calculations will have to be performed. Third, the symmetry finding of GoBaby is not the best, and more calculations will be set up for a supercell that is explicitly put into the apos.dat file than for a supercell that has been converted back to a primitive cell and made a supercell again by GoBaby. This procedure avoids this problem by converting the supercell back to a primitive cell before inputing it into the apos.dat file.
Calculating Force Constants
The frozen phonon method calculates the force constant matrix by displacing each atom in a unit cell and calculating the resulting forces on every other atom. Ideally, this is done for an infinite crystal, and the induced forces between atoms go to zero as the distance between the involved atoms increases. In practice however, we must make do with a finite sized cell, impose periodic boundary conditions, and deal with the consequences. The consequence for phonons is that the force constant matrix, instead of being a 3N x 3N matrix with N ~ 1023 whose off diagonal elements all go to 0 very quickly, is a 3n x 3n matrix with n ~ 100 whose off-diagonal elements do not necessarily disappear.
Running GoBaby
There are three components to the GoBaby code, a set of Fortran programs which do all of the heavy lifting, and two bash scripts GoBaby-Pt1 and GoBaby-Pt2. As of this writing, the Fortran code has only been compiled on one of the Wolverton Group computer clusters, Victoria.
To run GoBaby, a directory should be created containing the input files mentioned above. Running GoBaby-Pt1 from this directory will create a number of subdirectories, each containing the files needed to run a VASP calculation. The calculations will be MD runs, with only one atom moving along one direction per calculation. The symmetry of the unit cell will determine how many calculations need to be carried out (see caveat above on GoBaby's symmetry finding). Next, run VASP on each of the subdirectories as you normally do. A script phonon_job.sh that will copy and submit queue files from a directory to all subdirectories containing the same prefix is attached to this document. Once all of the phonon MD calculations have finished successfully, the file kpgen.dat should be created. This file is used as input to make phonon pDOS files. An example kpgen.dat file is shown below:
1 1 1 0 0 0 (The first line contains the q-point mesh phonons were calculated on, and any shift to that mesh)
T (The second line is T if interpolation of q-points should be done, otherwise F)
1 1 1 (The third line gives the q-point mesh for DOS interpolation)
The first line should always be 1 1 1 0 0 0, unless you used a larger set of q-points in the apos.dat file (unlikely, and not recommended). The third line, is an interpolation mesh for calculating the DOS at q-points not explicitly calculated. This set of numbers can be varied and converged with respect to the shape of the phonon DOS.
Imaginary Modes
The results from frozen phonon calculations can sometimes give normal modes with negative eigenvalues. This corresponds to having imaginary frequencies, a normal mode which decreases energy along its displacement vector, and an unstable structure. That is, if the mode is actually unstable and not just a numerical instability. To explicitly check if a mode is imaginary, one can perform static calculations along a trajectory of the phonon displacement vector. If the mode is actually unstable, then the energy of this displacement should decrease along the trajectory. If the mode is only numerically unstable, then the energy should increase along the trajectory. GoBaby has some tools and output files that make calculating phonon mode displacements straightforward.
The output file sc_pos.dat contains a list of the atoms on their initial sites. The file disp.dat contains a list of phonon mode displacement vectors. These are 3Nx1 complex vectors (though they should all be real at the gamma point) that are listed in a format:
x1 x2 x3 y1 y2 y3
...
where for each