The two-dimensional electron density of gallium arsenide on a plane through the Ga-As bond, computed with Elk alone -- a type 1, independent calculation with the QELK suite (no Quantum ESPRESSO involved here; the QE-only counterpart is the GaAs_density2D_qe 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_elk.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
# One Elk step: task 0 is the SCF, task 32 the 2D density.
step dens2D_elk elk
from initial
task 0 32
npoints2d 100 100
plane
-0.45784271 -0.45784271 0.39068542
0.14215729 0.14215729 0.99068542
0.10784271 0.10784271 -0.74068542
end
rgkmax 7.0
ngridk 6 6 6
xctype 20
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 species blocks name Elk's Ga.in / As.in; the qe_pseudo and qe_mass lines are not used by this run, but they keep the species block identical to the QE and comparison examples, so the three files differ only in their steps. xctype 20 is Elk's PBE, the functional of the PBE pseudopotentials used in the QE example.
Note: a single Elk step does both jobs: the task list 0 32 runs the ground state (task 0) and then the 2D density (task 32). The dimension is chosen by the task, not by a keyword: task 31 is a line (1D), 32 a plane (2D), 33 the whole cell (3D), and task 0 31 32 33 would give all three from one SCF. qelk also writes tshift .false. in elk.in, so that Elk keeps the atoms where the .qelk file puts them; otherwise Elk shifts the crystal and the plane would no longer pass through the atoms.
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. The two edges must be perpendicular in Cartesian space; qelk checks this before anything runs. Here 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 goes along 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 sit well inside the rectangle, with room around them.
Note: npoints2d 100 100 is what you ask for; Elk samples the plane at i/n, so qelk writes plot2d ... 99 99 in elk.in and RHO2D.OUT holds 99 x 99 = 9801 points (the same spacing that pp.x gives with 100 x 100, which is what makes the comparison example point-by-point).
Calculation
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_density2D_elk> qelk GaAs_density2D_elk.qelk
A first check with --dry-run writes elk.in without launching Elk:
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_density2D_elk> qelk GaAs_density2D_elk.qelk --dry-run
There is one step, so everything is in run/00_dens2D_elk. Important files are also copied next to the .qelk file, prefixed with the step id (dens2D_elk.elk.in, dens2D_elk.INFO.OUT, dens2D_elk.TOTENERGY.OUT, dens2D_elk.GEOMETRY.OUT, dens2D_elk.elk.stdout). The density itself, RHO2D.OUT, stays in run/00_dens2D_elk/.
SCF results
cat dens2D_elk.elk.stdout grep -E "Total number of k-points|Muffin-tin radius times|Convergence targets" dens2D_elk.INFO.OUT wc -l dens2D_elk.TOTENERGY.OUT tail -1 dens2D_elk.TOTENERGY.OUT tail -1 run/00_dens2D_elk/GAP.OUT
What this run gave (abridged):
Info(checkmt): reduced muffin-tin radius of species 1 (Ga) from 2.4000 to 2.2879
Info(checkmt): reduced muffin-tin radius of species 2 (As) from 2.4000 to 2.2879
Info(elk): current task : 32
Info(rhoplot): 2D density plot written to RHO2D.OUT
Total number of k-points : 22
Muffin-tin radius times maximum |G+k| : 7.000000000
Convergence targets achieved
20 dens2D_elk.TOTENERGY.OUT
-4205.24022721
0.187999762681E-01
Elk reduced the muffin-tin radii of both atoms from the 2.4 bohr of the species files to 2.288 bohr, because the two spheres would otherwise overlap. The 6 x 6 x 6 mesh reduces to 22 k-points with the 24 crystal symmetries; the SCF converged in 20 iterations to a total energy of -4205.2402 Ha. This is an all-electron energy (core states included), so it cannot be compared with the pseudopotential total energy of QE. The last line of GAP.OUT is the band gap of the converged SCF: 0.0188 Ha, about 0.51 eV -- the usual PBE underestimate of the GaAs gap, which is of no consequence for a density map.
Plot
RHO2D.OUT has the columns x, y, rho: the two coordinates along the edges of the plane, in bohr, and the density in electrons/bohr3. The plot script finds the file from the .qelk file; with a single Elk step, use plot_elk_density.py:
plot_elk_density.py --qelk GaAs_density2D_elk.qelk -o density2D_elk_full_scale.png
On the full scale the figure is nearly empty: Elk's all-electron density reaches 1967.6 electrons/bohr3 at the nuclei, while the bonding region is around 0.07, so the colour bar is spent on the two nuclear spikes. The fix is to cut the colour scale with --vmax:
plot_elk_density.py --qelk GaAs_density2D_elk.qelk --vmax 1.5 -o density2D_elk.png
Note: the figure has a filled-contour map on top and a 3D relief below (height and colour are the density). The atoms that lie 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. --style contour or --style relief keeps only one of the two panels, --log uses a logarithmic colour scale (an alternative to clipping), and --cmap changes the colormap.
Analysis
RHO2D.OUT has 99 x 99 points; the density runs from 0.0027 e/bohr3 at the lowest point of the plane to 1967.6 at the nucleus of As (about 1867 at the grid point nearest the nucleus of Ga). The two atoms are on the same row, 4.63 bohr apart, with Ga at x = 3.25 bohr and As at x = 7.85 bohr; the other two blobs at the edges of the map are periodic images of the atoms of the neighbouring cells.
Along that row, the density falls from the Ga nucleus to a minimum of 0.0702 e/bohr3 at x = 5.38 bohr, close to where the two muffin-tin spheres meet (the radii are 2.288 bohr each, 2 x 2.288 = 4.576 bohr against a bond length of 4.626 bohr), and then rises again into As. The minimum lies on the Ga side of the bond centre (x = 5.55 bohr); the density at the grid point nearest the centre (x = 5.61 bohr) is 0.0748, higher than at x = 5.38 bohr. This points to charge displaced toward As, which is what one expects for the polar Ga-As bond; a single row through the nuclei is only a hint, so look at the whole map before drawing conclusions.
To go further: change the three vertices to look at another plane (keep the two edges perpendicular), or use a miller h k l line instead of the plane block when the plane is a standard one -- the GaAs_density_planes_elk_vs_qe example does this for the (110) and (111) planes. For the whole cell, task 0 33 gives RHO3D.OUT, which the suite can convert to a .cube file for VESTA.
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
0 Comments