The two-dimensional electron density of gallium arsenide on a plane through the Ga-As bond, computed with Quantum ESPRESSO alone -- a type 1, independent calculation with the QELK suite (no Elk involved here; the Elk-only counterpart is the GaAs_density2D_elk example, and the two codes side by side are in GaAs_density2D_elk_vs_qe).
The input file
We need only the input file
GaAs_density2D_qe.qelk
version 1
cell bohr
-5.3415 0.0000 5.3415
0.0000 5.3415 5.3415
-5.3415 5.3415 0.0000
end
atoms crystal
Ga 0.00 0.00 0.00
As 0.25 0.25 0.25
end
species Ga
elk_file Ga.in
qe_pseudo Ga.pbe-dn-kjpaw_psl.0.2.upf
qe_mass 69.723
end
species As
elk_file As.in
qe_pseudo As.pbe-n-kjpaw_psl.0.2.upf
qe_mass 74.922
end
step scf qe
from initial
calculation scf
ecutwfc 40
kpoints 6 6 6
occupations fixed
end
# pp.x: a 'density' step continues from the scf and takes only the region.
step dens2D_qe density
from scf
density 2d
npoints2d 100 100
plane
-0.45784271 -0.45784271 0.39068542
0.14215729 0.14215729 0.99068542
0.10784271 0.10784271 -0.74068542
end
end
Note: the cell is the zincblende primitive cell with a = 10.683 bohr (5.653 Å), written as 2 x 5.3415, with Ga at (0,0,0) and As at (1/4,1/4,1/4). The two species blocks point to the PBE PAW pseudopotentials (the Elk file names are not used by this run). GaAs is an insulator, so occupations fixed is enough; the 6 x 6 x 6 mesh is a density, not a DOS, calculation.
Note: QE needs two steps: scf (pw.x, the charge density) and dens2D_qe (pp.x, which reads that density and samples it on the plane). The pp.x step has from scf and takes no ecutwfc or k-points of its own -- those belong to the scf step. The line density 2d restricts the step to the 2D region; without a density line the step would do the 1D, 2D and 3D densities, and would need both a line and a plane block. The default quantity is plot_num = 0, the total valence charge density.
Note: the plane block holds three vertices in fractional cell coordinates: the origin corner, the end of the first edge and the end of the second edge. pp.x needs the two edges to be perpendicular in Cartesian space; qelk checks this before anything runs. The plane is the (110) plane that contains the Ga-As bond (a <111> direction). The bond length is √3/4 · a = 4.626 bohr. Edge 1 follows the bond, direction (-1,1,1), for 2.4 bond lengths (11.10 bohr); edge 2 is perpendicular to it inside the plane, direction (1,-1,2), for 1.6 bond lengths (7.40 bohr). The origin lies 0.7 bond lengths before Ga along edge 1 and 0.8 bond lengths below it along edge 2, so both atoms are well inside the rectangle.
Calculation
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_density2D_qe> qelk GaAs_density2D_qe.qelk
A first check with --dry-run writes pw.in and the pp.x input without launching any binary:
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_density2D_qe> qelk GaAs_density2D_qe.qelk --dry-run
The steps run in file order: run/00_scf and run/01_dens2D_qe. Important files are also copied next to the .qelk file, prefixed with the step id: scf.pw.in, scf.pw.out, dens2D_qe.density_pp_2d.in, dens2D_qe.density_pp_2d.out and the data file dens2D_qe.density_2d.dat (also in run/01_dens2D_qe/, together with rho.pp).
The pp.x input that qelk wrote (dens2D_qe.density_pp_2d.in):
&INPUTPP prefix = 'hybrid' outdir = '.../out' filplot = 'rho.pp' plot_num = 0 / &PLOT iflag = 2 x0(1) = 0.04748738, x0(2) = -0.04748738, x0(3) = -0.64748737 e1(1) = -0.84852814, e1(2) = 0.84852814, e1(3) = 0.84852814 e2(1) = 0.40000000, e2(2) = -0.40000000, e2(3) = 0.79999999 nx = 100, ny = 100 output_format = 7 fileout = 'density_2d.dat' /
Note: this is the translation of the three fractional vertices into what pp.x wants: iflag = 2 is a 2D plot, x0 the origin, e1 and e2 the two edge vectors (Cartesian, in units of alat, which for ibrav = 0 is the length of the first cell vector, 7.554 bohr), nx and ny the number of points, and output_format = 7 the gnuplot x, y, f layout. You never write this by hand: the vertices of the .qelk file are the only description of the plane.
SCF results
grep -E "convergence has been achieved|^!|number of electrons|number of k points|highest occupied" scf.pw.out grep -E "Min, Max|Plot Type" dens2D_qe.density_pp_2d.out
What this run gave:
number of electrons = 18.00
number of k points= 16
highest occupied level (ev): 7.6247
! total energy = -453.96263609 Ry
convergence has been achieved in 9 iterations
Min, Max, imaginary charge: 0.002612 1.512113 0.000000
Plot Type: 2D contour Output format: gnuplot x,y,f
The 18 electrons are 13 from Ga (the 3d shell is in the valence of this pseudopotential) and 5 from As. The 6 x 6 x 6 mesh reduces to 16 k-points, and the SCF converges in 9 iterations to -453.96263609 Ry. pp.x then wrote density_2d.dat, whose density runs from 0.0026 to 1.5121 electrons/bohr3, with no imaginary part.
Plot
density_2d.dat has three columns, x, y and rho: the two coordinates along the edges of the plane, in bohr (pp.x already gives bohr in 2D), and the density in electrons/bohr3, with a blank line between rows. The plot script reads it from the .qelk file; --code qe draws only the QE panel:
plot_qe_density.py --qelk GaAs_density2D_qe.qelk --code qe -o density2D_qe.png
Note: the figure has a filled-contour map on top and a 3D relief below (height and colour are the density). The atoms on the plane are marked and named, the axes are labelled with the crystal directions of the plane (x along [1 1 1], y along [1 1 -2] in the lattice basis), and the black lines are the cross-section of the unit cell with the plane. Useful options: --style contour or --style relief for one panel only, --log for a logarithmic colour scale (it makes the weak bonding region visible), --vmax to set the top of the scale, --cmap for another colormap.
Analysis
The map shows 100 x 100 = 10000 points. Ga is at x = 3.25 bohr and As at x = 7.85 bohr on the same row, 4.63 bohr apart. The pseudo-valence density has its maximum, 1.512 e/bohr3, on Ga: the 3d electrons are in the valence of the Ga pseudopotential, and they are compact. On As the density is only 0.12 e/bohr3 at the nucleus, and the pseudo-valence density has no cusp or core peak anywhere. This is what pp.x with plot_num = 0 gives: the valence density of the pseudo-wavefunctions, not the all-electron density, so values inside the atomic spheres should not be read as the real density of the atom.
Along the row through the two atoms the density drops from Ga to a minimum of 0.0705 e/bohr3 at x = 5.38 bohr, and is 0.0740 at the grid point nearest the bond centre (x = 5.61 bohr); it then rises again toward As. A second, shallower minimum (0.045 at x = 7.18 bohr) lies on the As side, close to the nucleus, where the pseudo-valence density is not the real density; do not read it as a bond feature.
Parameters to test next: ecutwfc and the k mesh (check that the bond region does not change), and npoints2d for a smoother map. To compare with an all-electron reference, run the same plane in Elk (GaAs_density2D_elk) or the comparison example.
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
0 Comments