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.
