3 PHONON TRANSITIONS CODE

**********************************************************************

NOTE: 

This program is a work in process!! It is slow, annoying etc.
Please bear with me.

It now supports frequencies in cm-1, THz, eV and meV. 
It cannot handle symmetry.

Future plans: Parallelise the code to speed up operation.
Rewrite in python, add symmetry.
Any suggestions? :)

J. Buckeridge (j.buckeridge@lsbu.ac.uk) Nov. 2019

**********************************************************************

DESCRIPTION:

A program to determine the allowed phonon decays and collisions, given
a set of phonon frequencies calculated on a grid of q-points. The
format for input is a mesh.yaml file created by the Phonopy code.

The point of this code is to ascertain which phonon modes are involved
in transitions in 3-phonon processes that determine thermal
conductivity. Whether the aao, aaa, aoo etc processes (a = acoustic, o
= optical) dominate in a 3-phonon process could be inferred from the
results of this program, but of course the strengths of the
transitions are unknown. A full 3-phonon calculation needs to be
carried out to determine these strengths.

See the paper PRB 91 094306 (2015), where A. Togo et al describe the
theory behind the 3-phonon interaction calculations encoded in
Phono3py. What this code determines is the relation \Delta(q + q' +
q") in Eq. 10, and the delta function delta(omega - omega' - omega")
in Eq. 11. In this way, conditions for which phonon decays are allowed
are determined. Of course, when a phonon decay is allowed, so is a
phonon collision (with e.g. omega and omega' swapped).

This code reads in the q vectors and frequencies omega at each q, so
that it finds a set of (q,omega) points. It also reads in the
reciprocal lattice vectors. It then finds the combinations of q that
satisfy the quasimomentum conservation criterion:

q - q' - q" - G = 0,

where G is a reciprocal lattice vector. The max G checked: G = +/- (3
b1 + 3 b2 + 3 b3), where bi are the reciprocal lattice vectors.

When quasimomentum conservation is obeyed, the code then checks for
energy conservation omega = omega' + omega". If this is obeyed, a
transition has been found.

As symmetry is not used, redundancies occur in the output. These are
checked for and removed.

--------------------
UNITS:

The default unit for frequency is cm-1, but the units can be specified
in a file 'unit.dat', which the code will search for. It will 
recognise cm-1, THz, eV and meV.

--------------------
BAND STRUCTURE PLOTTING:

The code can also read in a band.yaml file, which it assumes is
consistent with the mesh.yaml (i.e. is for the same system, but just
has the frequencies along a high-symmetry line in the BZ. The band and
mesh files must use the same frequency units). For the allowed 
transitions that has been determined, it checks if the all q points 
are included in the BZ path in the band.yaml file. If so, it produces
a file 'bz-transitions.dat', which contains data to plot lines on a
phonon band structure diagram. The x-axis, which corresponds to 
q-points in the BZ, is given in units of "distance" that are included
in the band.yaml file. 

In the package, I have also included a simple bash script (again which
is a bit slow), that converts the band.yaml file to a file 'disp.dat'
which can be plotted in xmgrace. By running 'xmgrace -nxy disp.dat',
then importing the 'bz-transitions.dat' file (as single source), you 
can plot the band structure with the lines indicating the allowed
transitions. The script is called 'lattice-vib-disp.sh'.

---------------------
JOINT DOS:

The code can calculate the joint density of states (JDOS) for phonons,
eqns 19 and 22-24 of Togo's paper mentioned above, PRB 91 094306 (2015). 
It will do so if it finds a file 'temperature.dat', which contains a 
temperature value. The JDOS N2 depends on the occupations of the modes, 
hence requires a temperature. The output is written to the file 
'jdos.dat'. 

**********************************************************************

INPUT:

The code reads in a mesh.yaml file produced by Phonopy. It cannot
handle symmetry, so to produce the appropriate mesh.yaml, you need to
have MESH_SYMMETRY = .FALSE. in your mesh.conf file.

The default frequency unit is cm-1, but the units can be specified
by including the unit.dat file, in which you can write THz, eV, meV
or cm-1.

A band.yaml file can be read in, and the allowed transitions are 
checked to see if the (q,omega) points are included in the list in 
that file.

To calculate the JDOS, the file temperature.dat must be present. In
it the temperature is given, which is necessary to compute the N2
JDOS (see Eq. 22-24 of Togo's paper PRB 91 094306 (2015)).

**********************************************************************

OUTPUT:

The code produces a file, 'allowed-omega.dat' which contains the
values of omega, omega' and omega" for each allowed transition, as 
well as the degeneracy of the transition (which can be large for a 
highly symmetric system, as there may be many equivalent q-points in
the mesh). One can then plot omega' and omega" against omega in a 
2D plot by extracting the first three columns and running e.g. 
xmgrace -nxy allowed-omega-extract.dat (use only points, not lines to
make the output clear). You should then be able to see if collisions
are of aaa, aao, aoo, ooo nature. The degeneracies could be added as
a heat map using other plotting software.

If band.yaml was included, and some transitions were found that involved
(q,omega) points in the range of the band structure, a file containing
lines to be plotted on top of the band structure is outputted, called
'bz-transitions.dat'. The x-axis units are distance along the BZ path
(as given in the band.yaml file), while the y-axis is phonon frequency. 
The bash script 'lattice-vib-disp.sh' produces a file 'disp.dat' from
the band.yaml file, which is compatible with the 'bz-transitions.dat'
file, so both can be plotted together. You will need to put in the 
high symmetry points in the BZ path yourself.

If 'temperature.dat' is present, then the JDOS will be calculated and 
output to the file 'jdos.dat'. The JDOS are given for each value of
q and omega.

**********************************************************************

INSTALL:

It is a FORTRAN 90 code and should work with the gnu compiler. There
is a Makefile, so put in your preferred fortran compiler there and
type 'make'.

**********************************************************************

FINAL REMARKS:

Imaginary modes are ignored, i.e. set to zero and not included in the
checks. There is a bit of ambiguity in determining the minimum
frequency at Gamma, as the three acoustic modes there should be
zero. I'm still working on that too, so if anything weird happens with
imaginary modes and/or the three lowest modes at Gamma, let me know!

This code will hopefully be useful if you would like to get an idea of
how many 3-phonon interactions will be allowed in your system and what
there nature will be, without having to do the very expensive full
3-phonon calculation (using phono3py).

Using an MIT licence.

Email me please with all criticisms issues bugs complaints etc!

j.buckeridge@lsbu.ac.uk