LID Phonon (lid_phonon.py)

Background

lid_phonon is a code which performs the LID method, but only at second order. The advantage of having a standalone code for second order is that one can basically hard code the group theory for this particular case, which ensures that the code will be relatively fast when generating runs for large FTGs.

lid_phonon does the necessary group theory, sets up runs using xtal_derivative.py, and then will gather all the results to form the dynamical matrix.

Here we give a basic description of how to run lid_phonon.py on a simple test case of UO2. The actual runs can be found here:

/home/cam1/dft_runs/uo2/pbe/cut450_k16/ferro/u0/lid_phonon

We start with a YAML input file (e.g. inputs.yaml) containing the necessary keywords:

fpinp: ../static/             # relative path to first-principles (DFT code) input files
pgname: Oh                    # point group
schema: generate              # tells code to generate runs
undistorted_path: ../static/  # path to undistorted results (if needed)
supa: -1,1,1;1,-1,1;1,1,-1    # 3x3 matrix defining the supercell
xtal: ../static/POSCAR        # path to file defining the crystal structure
delta: 0.01,0.04,0.01         # squence of displacement amplitudes to run: 0.01, 0.02, ...
target: force                 # take derivatives of forces (or energy)

We can now execute lid_phonon.py

lid_phonon.py -i input.yaml

This will create the following directories:

$ ls */

q1__000:
  irr_T1u.0  irr_T2g.0

q2__110:
  irr_1.Eu.0  irr_A1g  irr_A2u  irr_B1u  irr_Eg.0  irr_Eu.0

Given that the determinant of our supercell matrix is 4, there will be four q-points in the discretized first Brillouin zone. However, the Oh point symmetry dictates that only two of them will be irreducible. Therefore, we have two q-point directories as shown above. The integer next to the “q” is the denominator, and the three integers to the right of the underscore are the numerators: we have q-points (0,0,0) and (1/2,1/2,0) in lattice coordinates. Given that UO2 has 9 degrees of freedom in the primitive cell (i.e. 3 atoms with 3 degrees of freedom each), the listed irreducible representation directors should have dimension which adds to 9. This is not true for the Gamma point, which only has two three-fold irreducible representation directories. This is due to homogeneity of free space, meaning that we do not need to compute the derivatives of the shift modes, which are zero by construction.

Within each irreducible representation directory, there is a directory for each displacement amplitude, and within that directory, there are subdirectories for each multiple of that amplitude which is needed for the derivative at hand:

$ ls q1__000/irr_T1u.0/
delta_1.0000  delta_2.0000  delta_3.0000

$ ls q1__000/irr_T1u.0/delta_1.0000/
run_m1  run_p1

The notation “m1” and “p1” indicate that we are displacing -1*delta and +1*delta, respectively. These directories contain the actual DFT files that were copied from the template directory fpinp provided in the input file.

You are now free to run all the DFT calculations in each subdirectory. When the runs are finished, you then need to run lid_phonon.py again to extract all of the data, perform the central finite difference calculations of the derivtives, and construct the dynamical matrices. This is done by simply changing the schema variable from generate to evaluate. After doing so, the code will create an output file in each q-point directory:

$ ls q*
q1__000:
Dq.pkl  Dq_symm.txt  Dq.txt  irr_T1u.0  irr_T2g.0  Uk.txt

q2__110:
Dq.pkl  Dq_symm.txt  Dq.txt  irr_1.Eu.0  irr_A1g  irr_A2u  irr_B1u  irr_Eg.0  irr_Eu.0  Uk.txt

The file Dq_symm.txt contains the dynamical matrix in the symmetrized basis, while Dq.txt contains the usual standard basis; and Dq.pkl is a pickle that will be read when generating the phonon bands, etc. We will now need to interface with the interpolation code, fourier_interp_Dqn.py.

In order to interface with fourier_interp_Dqn.py, we need to creat the proper input file, called fi_info. The following is an example to create this file given a particular FTG, given by nq, and the Oh point group. The value of delta is irrelavant. This command should be run in a subdirectory as it will create a bunch of directories that are not needed: only fi_info is needed from this run.

mkdir generate_fi_info
cd generate_fi_info
fourier_interp_Dqn.py nosym=False sym=False use_symadt=False nq='-1 1 1 1 -1 1 1 1 -1' pos=POSCAR delta=0.01 order=2 pgn=Oh -gen
cd ..
ln -s fi_info .

Please note that I have copied the POSCAR file into this subdirectory. This is because the POSCAR needs to follow a deprecated convention where the label on each atom needs to be followed by a “:p”, which tells the group theory code that each atom is associated with three displacement vectors.

$ more POSCAR
UO2 Fluorite
1.0
0.00000000   2.71091838   2.71091838
2.71091838   0.00000000   2.71091838
2.71091838   2.71091838   0.00000000
1 2
d
0.00000000   0.00000000   0.00000000  U:p
0.25000000   0.25000000   0.25000000  O:p
0.75000000   0.75000000   0.75000000  O:p

Now that fi_info is present in the current directory, you can simply run fourier_interp_Dqn.py without any commands, and it will generate PhiN_WS.pkl, which are the force constants that have been repacked into the Wigner-Seitz cell such that they can be used to compute the phonons for any q-point.

If you want to create a band structure or DOS, you can now use tb_tools.py, and a typical command would be

tb_tools.py  pos=POSCAR ham=PhiN_WS.pkl -mass nk=10,10,10 -bands units=THz

Do not forget the mass key, because otherwise the code will think you are doing electrons. The path in k-space is defined in the input file called kinput, which contains a list of q-points.

$ more kinput
0 0 0
0.5 0.5 0

The output is a file called bands.dat, in addition to a file called plt_bands.gr which can be called using our interface to xmgrace called gracey.py (found in the utils directory).

gracey.py -i plt_bands.gr

[ ]:



[ ]: