In
the following tutorial it is shown how to calculate electronic
structure of graphene with PWSCF code included in the QUANTUM ESPRESSO
package installed on CRESCO.
The details on keywords are found in the QUANTUM ESPRESSO web-page and pw.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.
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).
You can also add the option to limit time, with -W TIME, where TIME is the time in minutes.
Variable Cell Relax and Relax
1)
Build up the graphene structure with your preferred atomistic model
editor. Make note of the unit cell parameters and atomic coordinates.
2)
Make a relax of the unit cell and atomic coordinates at the same time
with a variable cell relax (vc-relax). Create your pwscf input file graphene.pbe.vc-relax.in as follows:
&CONTROL
calculation='vc-relax'
title='graphene'
prefix='graphene'
verbosity='high'
restart_mode='from_scratch'
nstep=1000
iprint=1
tprnfor=.true.
outdir='./tmp'
disk_io='default'
pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'
tstress=.true.
forc_conv_thr=1.0d-4
etot_conv_thr=1.0d-5
/
&SYSTEM
ibrav = 12,
a = 2.460,
b = 2.460,
c= 20.000,
cosab=-0.500000,
nat = 2,
ntyp = 1,
ecutwfc = 40.0 ,
ecutrho = 400.0 ,
input_DFT = 'PBE' ,
occupations = 'smearing' ,
degauss = 1.0d-4 ,
smearing = 'marzari-vanderbilt' ,
/
&ELECTRONS
electron_maxstep = 100,
conv_thr = 1.0d-10 ,
mixing_mode = 'plain' ,
mixing_beta = 0.3d0 ,
/
&IONS
ion_dynamics='bfgs'
upscale=20.0
/
&CELL
press_conv_thr = 0.5D0
press = 0.D0
cell_dynamics = bfgs,
cell_dofree = '2Dxy'
cell_factor = 1.5D0
/
ATOMIC_SPECIES
C 12.011 C.pbe-van_ak.UPF
ATOMIC_POSITIONS crystal
C 0.000000000 0.000000000 0.250000000
C 0.333333333 0.666666667 0.250000000
K_POINTS automatic
12 12 1 0 0 0
Notes:
i) The atomic coordinates have been inserted in
crystal units. However, you can also use coordinates in Angstrom. You
can easily switch between the units used for atomic coordinates by using
XCrySDen.
2b) Launch the calculation with EXE=pw .
3) Open the output file graphene.pbe.vc-relax.out and look for the final coordinates as follows:
End of BFGS Geometry Optimization
Final enthalpy = -22.7899403251 Ry
Begin final coordinates
new unit-cell volume = 711.58442 a.u.^3 ( 105.44593 Ang^3 )
CELL_PARAMETERS (alat= 4.64872629)
1.002996671 -0.000000000 0.000000000
-0.501498335 0.868620597 0.000000000
0.000000000 0.000000000 8.130081301
ATOMIC_POSITIONS (crystal)
C 0.000000000 0.000000000 0.250000000
C 0.333333333 0.666666667 0.250000000
End final coordinates
Extract
the new unit cell parameters a, b and the atomic coordinates from the
output file. You can also use Xcrysden to have them.
4) Make a further relax calculation. Create your pwscf input file graphene.pbe.relax.in as follows:
calculation='relax'
title='graphene'
prefix='graphene'
verbosity='high'
restart_mode='from_scratch'
nstep=1000
iprint=1
tprnfor=.true.
outdir='./tmp'
disk_io='default'
pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'
tstress=.true.
forc_conv_thr=1.0d-4
etot_conv_thr=1.0d-5
/
&SYSTEM
ibrav = 0,
celldm(1) = 4.64872629
nat = 2,
ntyp = 1,
ecutwfc = 40.0 ,
ecutrho = 400.0 ,
input_DFT = 'PBE' ,
occupations = 'smearing' ,
degauss = 1.0d-4 ,
smearing = 'marzari-vanderbilt' ,
/
&ELECTRONS
electron_maxstep = 100,
conv_thr = 1.0d-10 ,
mixing_mode = 'plain' ,
mixing_beta = 0.3d0 ,
/
&IONS
ion_dynamics='bfgs'
upscale=20.0
/
ATOMIC_SPECIES
C 12.011 C.pbe-van_ak.UPF
CELL_PARAMETERS alat
1.002996671 -0.000000000 0.000000000
-0.501498335 0.868620597 0.000000000
0.000000000 0.000000000 8.130081301
ATOMIC_POSITIONS crystal
C 0.000000000 0.000000000 0.250000000
C 0.333333333 0.666666667 0.250000000
K_POINTS automatic
12 12 1 0 0 0
4b) Launch the calculation with EXE= pw .
5) Open the output file graphene.pbe.relax.out and look for the final coordinates as follows:
End of BFGS Geometry Optimization
Final energy = -22.7899971025 Ry
Begin final coordinates
ATOMIC_POSITIONS (crystal)
C 0.000000000 0.000000000 0.250000000
C 0.333333333 0.666666667 0.250000000
End final coordinates
Self-Consistent Calculation
1) Prepare the input file graphene.pbe.scf.in for the self-consistent scf calculation as follows:
&CONTROL
calculation='scf'
title='graphene'
prefix='graphene'
verbosity='high'
restart_mode='from_scratch'
nstep=1000
iprint=1
tprnfor=.true.
outdir='./tmp'
wf_collect=.true.
disk_io='default'
pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'
tstress=.true.
/
&SYSTEM
ibrav = 12,
a = 2.4674,
b = 2.4674,
c= 20.000,
cosab=-0.500000,
nat = 2,
ntyp = 1,
ecutwfc = 40.0 ,
ecutrho = 400.0 ,
input_DFT = 'PBE' ,
occupations = 'smearing' ,
degauss = 1.0d-4 ,
smearing = 'marzari-vanderbilt' ,
/
&ELECTRONS
electron_maxstep = 100,
conv_thr = 1.0d-10 ,
mixing_mode = 'plain' ,
mixing_beta = 0.3d0 ,
/
ATOMIC_SPECIES
C 12.011 C.pbe-van_ak.UPF
ATOMIC_POSITIONS angstrom
C 0.000000000 0.000000000 5.000000000
C 0.000000000 1.424500000 5.000000000
K_POINTS automatic
12 12 1 0 0 0
2) Launch the calculation with EXE= pw .
When the calculation is done, you can extract the total energy with the command:
tail -200 graphene.pbe.scf.out | grep !
You will get:
! total energy = -22.78999768 Ry
Post-Processing: Potential and Charge Density
1)
The wave functions and the charge density are now available in
outdir='./tmp' and can be used for the post-processing. First, we
collect the charge density map in the file graphene.pbe.rho. Create the
post-processing input file graphene.pbe.rho.in as follows:
&inputpp
prefix='graphene'
outdir='./tmp'
filplot='graphene.pbe.rho'
plot_num=0
/
&plot
/
2) Launch the post-processing with EXE= pp .
3)
We collect the electrostatic potential (V_bare + V_H) grid in the file
graphene.pbe.pot. Create the post-processing input file graphene.pbe.pot.in as follows:
&inputpp
prefix='graphene'
outdir='./tmp'
filplot='graphene.pbe.pot'
plot_num=11
/
&plot
/
4) Launch the post-processing with EXE= pp .
5) Make readable in a suitable format the electrostatic potential and the charge density.
5a) Create the post processing file graphene.pbe.rho.3D.in:
&inputpp
/
&plot
nfile=1
filepp(1)='graphene.pbe.rho'
weight(1)=1.0
iflag=3
output_format=3
e1(1)=2.0, e1(2)=0.0, e1(3)=0.0
e2(1)=0.0, e2(2)=2.0, e2(3)=0.0
e3(1)=0.0, e3(2)=0.0, e3(3)=8.130081301
x0(1)=0.0, x0(2)=0.0, x0(3)=0.0
nx=50, ny=50, nz=50
fileout='graphene.pbe.rho.xsf'
/
5b) Launch with EXE= pp (1 processor is recommended).
5c) Create the post processing file graphene.pbe.rho.3D.in:
&inputpp
/
&plot
nfile=1
filepp(1)='graphene.pbe.pot'
weight(1)=1.0
iflag=3
output_format=3
e1(1)=2.0, e1(2)=0.0, e1(3)=0.0
e2(1)=0.0, e2(2)=2.0, e2(3)=0.0
e3(1)=0.0, e3(2)=0.0, e3(3)=8.130081301
x0(1)=0.0, x0(2)=0.0, x0(3)=0.0
nx=50, ny=50, nz=50
fileout='graphene.pbe.pot.xsf'
/
Notes: The xsf is a 3D data file that can be read with xcrysden.
5d) Launch with EXE= pp (1 processor is recommended).
6) The charge density and the electrostatic potential can be visualized wit XCrySDen by opening with it the two files graphene.pbe.rho.xsf and graphene.pbe.pot.xsf.
Density of States
1)
Now we calculate the density of states. First we perform a non
self-consistent calculation (nscf) based on the results of the
self-consistent calculation in the previous step. In the following input
file graphene.pbe.dos.nscf.in the number of bands (nbnd) and the k-point mesh are changed:
calculation='nscf'
title='graphene'
prefix='graphene'
verbosity='high'
restart_mode='from_scratch'
nstep=1000
iprint=1
tprnfor=.true.
outdir='./tmp'
wf_collect=.true.
disk_io='default'
pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'
tstress=.true.
/
&SYSTEM
ibrav = 12,
a = 2.4674,
b = 2.4674,
c= 20.000,
cosab=-0.500000,
nat = 2,
ntyp = 1,
ecutwfc = 40.0 ,
ecutrho = 400.0 ,
nbnd = 8 ,
input_DFT = 'PBE' ,
occupations = 'tetrahedra' ,
/
&ELECTRONS
electron_maxstep = 100,
conv_thr = 1.0d-10 ,
mixing_mode = 'plain' ,
mixing_beta = 0.3d0 ,
/
ATOMIC_SPECIES
C 12.011 C.pbe-van_ak.UPF
ATOMIC_POSITIONS angstrom
C 0.000000000 0.000000000 5.000000000
C 0.000000000 1.424500000 5.000000000
K_POINTS automatic
24 24 1 0 0 0
Notes:
I) nbnd is set to 1.2 times the number of occupied bands (the half of the number of electrons).
II) A dense mesh of k-points is recommended.
III) occupations = 'tetrahedra' is best suited for DOS calculations.
2) Launch the nscf calculation with EXE= pw .
3)
To calculate the DOS here we use projwfc.x in the QUANTUM ESPRESSO
package that will allow, in general, to calculate also the projected
density of states. The input file graphene.pbe.proj.in is:
&projwfc
prefix='graphene'
outdir='./tmp'
ngauss=0
degauss=0.01
Emin=-10.0, Emax=10.0, DeltaE=0.1
filproj='graphene.pbe.proj'
filpdos='graphene.pbe.pdos'
/
4) Launch with EXE= projwfc .
Notes:
i)
We can launch the job on a different number of processor (from 32 to
16) because we set the keyword wf_collect=.true. . Otherwise, the same
number of processors would be mandatory. However, it can be useful to
decrease the number of processors when the "too many processors" error
is found.
The output files will be the total DOS graphene.pbe.pdos.pdos_tot and the projected DOS on the atomic orbitals graphene.pbe.pdos.pdos_atm#*(C)_wfc#*(*).
You can plot them with xmgrace.
For
example, by setting nbnd=16, the mesh of k-points to 48 48 1 for scf
calculation and 96 96 1 for nscf calculation, you get a more accurate
DOS as in the figure below (we set degauss=0.02 in the projwfc input).
(Then, try to increase the mesh of k-points further and check convergence)
Band Structure
1)
Now we calculate the band structure. We perform again a non
self-consistent calculation (nscf) where k-point mesh is substituted
with paths along the lines connecting the high symmetry points in the
Brillouin zone. The input file graphene.pbe.bands.nscf.in follows:
calculation='bands'
title='graphene'
prefix='graphene'
verbosity='high'
restart_mode='from_scratch'
nstep=1000
iprint=1
tprnfor=.true.
outdir='./tmp'
wf_collect=.true.
disk_io='default'
pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'
tstress=.true.
/
&SYSTEM
ibrav = 12,
a = 2.4674,
b = 2.4674,
c= 20.000,
cosab=-0.500000,
nat = 2,
ntyp = 1,
ecutwfc = 40.0 ,
ecutrho = 400.0 ,
nbnd = 16 ,
input_DFT = 'PBE' ,
occupations = 'smearing' ,
degauss = 0.001 ,
smearing = 'marzari-vanderbilt' ,
/
&ELECTRONS
electron_maxstep = 100,
conv_thr = 1.0d-10 ,
mixing_mode = 'plain' ,
mixing_beta = 0.3d0 ,
/
ATOMIC_SPECIES
C 12.011 C.pbe-van_ak.UPF
ATOMIC_POSITIONS angstrom
C 0.000000000 0.000000000 5.000000000
C 0.000000000 1.424500000 5.000000000
K_POINTS crystal_b
4
-0.333333333 0.666666667 0.000000000 50 ! K
0.0000000000 0.000000000 0.000000000 50 ! G
0.000000000 0.500000000 0.000000000 50 ! M
-0.333333333 0.666666667 0.000000000 50 ! K
2) Launch the bands nscf calculation with EXE= pw .
3) Now, the bands are to be organized for plotting with the bands.x program. Create the following input file graphene.pbe.bands.in:
&bands
outdir = './tmp'
prefix = 'graphene'
filband = 'graphene.pbe.bands'
/
4) Launch the bands calculation with EXE= bands .
5) Plot the bands from the graphene.gnu output file with gnuplot or xmgrace (energy units are eV).
For
PWSCF versions older than 6.0 you have to generate a plot of the bands.
You can use plotband.x inserting graphene.pbe.bands as input file.
/afs/enea.it/software/qu_esp/qe600/cresco4-ser/plotband.x
You can also create plotb.in input file:
graphene.pbe.bands
-25 10
graphene.pbe.bands.dat
graphene.pbe.bands.ps
-2.3577
5 0
and run:
/afs/.enea.it/software/qu_esp/qe600/cresco4-ser/plotband.x < plotb.in > plotb.out
Hello, my name is Dr Abderrahmane Reggad. I'm a 53 year old and I am a researcher in materials' sciences. 
0 Comments