## Graphene with PHONON

Tutorials

**Graphene phonons calculation using phonon (ph.x in the QUANTUM ESPRESSO package)**

In the following tutorial it is shown how to calculate phonon of graphene with Pwscf & Phonon codes included in the QUANTUM ESPRESSO package installed on CRESCO.

The details on keywords are found in the QUANTUM ESPRESSO web-page and ph.x input file description.

You can check the best suited queue and scratch directory from the Cresco website. E.g., there the information for cresco4 and cresco6 clusters.

More information about how to use the QUANTUM ESPRESSO package on CRESCO can be found at this page.

1) Build up the graphene structure with your preferred atomistic model editor.

2) Make a relax of graphene following the example of the relax calculation in the previous tutorials. Here we provide the relaxed atomic coordinates and unit cell parameters.

To run a job on CRESCO6 the command line is:

**bsub -n NP -q QUEUE -o lsf_file.out -e file.err qe 6.3.0 par EXE input.inp output.out**

where NP is the number of processors (NP=48 is one node on CRESCO6), QUEUE=cresco6_48h24 (24h queue), cresco6_48h48 (48h queue, NP >= 48), cresco6_h144 (144h queue, NP can be less than 48), EXE = executable (e.g., pw, pp, bands, projwfc, ph).

3) Create your pwscf input file

**graphane.scf.in**as follows: &CONTROL

calculation = 'scf',

restart_mode = 'from_scratch',

verbosity='high'

nstep=100

tprnfor=.true.

disk_io='default'

pseudo_dir = './',

outdir = './tmp',

prefix = 'graphene',

tstress=.true.

/

&SYSTEM

ibrav = 4,

a = 2.4629,

c = 10

nat = 2,

ntyp = 1,

occupations = 'smearing',

smearing = 'methfessel-paxton',

degauss = 0.02,

ecutwfc = 120,

nbnd = 8,

/

&ELECTRONS

mixing_mode = 'plain' ,

conv_thr = 1.0d-10,

mixing_beta = 0.7,

/

ATOMIC_SPECIES

C 12.011 C_ONCV_PBE-1.0.upf

ATOMIC_POSITIONS (crystal)

C 0.333333333 0.666666666 0.250000000

C 0.666666666 0.333333333 0.250000000

K_POINTS (automatic)

12 12 1 0 0 0

Notes:

i) The convergence threshold of SCF calculation conv_thr must be very low (1.0D-10), smaller than for total energy calculation.

ii) GGA norm-conserving pseudopotentials have been used (C_ONCV_PBE-1.0.upf generated using ONCVPSP code by D. R. Hamann).

3) Launch the scf calculation using the input file above with EXE= pw .

4) Create your phonon input file

**graphene.ph.in**as follows:phonons of graphene on a grid ! title line

&inputph

outdir="/gporq2/scratch_0/usr/buonocor/tmp/graphene-ph4",

prefix="graphene", ! must be the same than scf calculation

tr2_ph = 1d-20, ! the threshold

ldisp = .true., ! when true the calculation is on a grid

nq1=9, nq2=9, nq3=1, ! uniform q-point grid

amass(1)=12.011,

fildyn='graphene.dyn', ! name of output file

/

5) Launch the phonon calculation using the input file above with EXE= ph .

6) Check the ouput files. In the same directory than input file you will find the files graphene.dyn#, with # ranging from 0 to 12. The files from 1 to 12 are the dynamical matrices on the grid of q-points. graphane.dyn0 is the list of inequivalent q-point (12 in this case).

7) Create the following

**graphene.q2r.in**input for q2r.x:&input

fildyn='graphene.dyn', ! the dynamical matrix from phonon

zasr='simple',

flfrc='graphene.fc', ! the interatomic force constants

/

8) Launch the conversion of the dynamical matrices in a supercell 9x9x1 of the reciprocal space from the phonon calculation (q) to the interatomic force constants in a supercell 9x9x1 of real space (r), using the input file above with EXE= q2r .

9) Create the

**graphene.matdyn.disp.in**input file for matdyn.x: &input

asr='simple',

amass(1)=12.011,

flfrc='graphene.fc', ! the interatomic force constants

flfrq='graphene.freq', ! output file with frequencies

q_in_band_form=.true. ! band format of k-points

/

4

0 0 0 40 ! Gamma point

0 0.577350269 0 20 ! M point

-0.333333333 0.577350269 0 40 ! K point

0 0 0 20 ! Gamma point

i) Note: setting q_in_band_form=.true. the list of k-points must be given as follows:

# of high-symmetry k-points

k1(A) k2(A) k3(A) np(A-B)

k1(B) k2(B) k3(B) np(B-C)

.....

where np(A-B) is the number of k-points along the path from A to B, and so on, A, B, ... being the high-symmetry k-points.

10) Calculate the phonon dispersion using the interatomic force constants by calculating phonons at a generic q point, using the input file above with the EXE= matdyn .

11) Generate a plot of the phonon bands. You can use plotband.x inserting

**graphene.disp.freq**as input file. Frequency units are cm-1. Use plotband.x in /afs/enea.it/software/qu_esp/<your_qe_version>.You can plot the graphene.disp.freq.dat data file with xmgrace.

However, latest versions of QE gives the .gp file for direct plotting.

12) Create also the

**graphene.matdyn.dos.in**input file for matdyn.x:&input

asr='simple',

dos=.true., ! option for dos phonon calculation

amass(1)=12.011,

flfrc='graphene.fc',

fldos='graphene.phdos', ! output file

nk1=50,nk2=50,nk3=1, grid of q-points for DOS calculation

deltaE=1.0d0 ! interval of energy

/

13) Calculate the vibrational density of states, using the input file above with EXE= matdyn .

14) Plot the density of states from graphene.phdos output file. Frequency units are cm-1.