The spin density of ferromagnetic nickel on a plane, computed twice and independently, with Quantum ESPRESSO (pw.x, then pp.x for the spin-up and spin-down densities) and with Elk (task 72, the 2D magnetisation density), from one input of the QELK suite -- a type 3, comparison calculation. It extends the density step of the suite (the GaAs 2D density examples) to spin-polarised systems.
Before you start: one small update of the suite
This example needs two small additions to QELK: a spin_component keyword for the QE density step, and the Elk tasks 72 and 73 (magnetisation density) in the density machinery. They are delivered as a script that edits a few lines in three files in place, with a backup of each file:
algerien1970@linux-ml0a:~> python3 apply_spin_density_patch.py ~/abinitio/QELK
Each file is saved as FILE.bak_spin before it is changed, and running the script a second time changes nothing. If the Interfaces package is installed in editable mode (pip install -e) there is nothing more to do; otherwise reinstall it as you did for the first installation.
The input file
We need only the input file
Ni_spin_density2D_elk_vs_qe.qelk
version 1
# fcc Ni, a = 6.652 bohr (3.52 A): primitive vectors with a/2 = 3.326
cell bohr
-3.326 0.000 3.326
0.000 3.326 3.326
-3.326 3.326 0.000
end
atoms crystal
Ni 0.00 0.00 0.00
end
species Ni
elk_file Ni.in
qe_pseudo Ni_PBE_TM_2pj.UPF
qe_mass 58.693
end
step scf qe
from initial
calculation scf
ecutwfc 60
kpoints 12 12 12
degauss 0.02
extra_system nspin=2,starting_magnetization(1)=0.5
end
# pp.x, spin up and spin down on the same plane; the spin density is up - down
step spin_up_qe density
from scf
density 2d
spin_component 1
npoints2d 100 100
plane
-0.5 -0.5 0.25
0.5 0.5 1.25
0.5 0.5 -1.75
end
end
step spin_dn_qe density
from scf
density 2d
spin_component 2
npoints2d 100 100
plane
-0.5 -0.5 0.25
0.5 0.5 1.25
0.5 0.5 -1.75
end
end
# Elk task 72: the magnetisation density on the plane (MAG2D.OUT)
step spin_elk elk
from initial
task 0 72
npoints2d 100 100
plane
-0.5 -0.5 0.25
0.5 0.5 1.25
0.5 0.5 -1.75
end
rgkmax 8.0
ngridk 12 12 12
xctype 20
spinpol 1
bfieldc 0.0 0.0 0.01
end
Note: The magnetism is seeded as in the Ni_magnetism_elk_vs_qe and Fe examples: nspin=2,starting_magnetization(1)=0.5 in QE, spinpol 1 with the small bfieldc 0 0 0.01 in Elk (spinpol alone would stay non-magnetic). The magnetisation points along z in both codes.
Note: QE has no "spin density" quantity in pp.x; it plots the charge density of one spin channel: spin_component 1 is the spin-up density, 2 the spin-down density (with plot_num 0, the default). The spin density is their difference, so the example runs two density steps on the same plane and the plot script subtracts them. Elk has a task for it: task 72 writes the magnetisation density on the plane, task 73 on the 3D grid (the plane block is the same one used for task 32).
Note: The plane is the (110) plane of the fcc cell (it contains the z axis, the direction of the magnetisation). In fractional coordinates its three vertices are the origin corner (-0.5, -0.5, 0.25) and the ends of two edges: edge 1 is the lattice vector (1,1,1), which is [-1 1 1] in Cartesian axes and 11.52 bohr long, edge 2 is (1,1,-2), which is [1 -1 2] and 8.15 bohr long. The two edges are perpendicular, as the density step requires. The plane passes through the atom at the origin of the cell, which lies a quarter of the way along both edges, so the rectangle shows several Ni atoms.
Calculation
algerien1970@linux-ml0a:~/abinitio/QELK/examples/Ni_spin_density2D_elk_vs_qe> qelk Ni_spin_density2D_elk_vs_qe.qelk
The steps run in file order: run/00_scf, run/01_spin_up_qe, run/02_spin_dn_qe, run/03_spin_elk. The two QE densities are in density_2d.dat of their step folders; Elk's is run/03_spin_elk/MAG2D.OUT. On one core the whole run takes about three minutes.
SCF results
grep -E "number of k points|convergence has been achieved|^!|total magnetization|absolute magnetization" scf.pw.out | tail -5 grep -A1 "Magnetic moment per site" scf.pw.out | tail -1 grep -E "Min, Max" run/01_spin_up_qe/density_pp_2d.out run/02_spin_dn_qe/density_pp_2d.out grep -E "Total number of k-points|Convergence targets" run/03_spin_elk/INFO.OUT tail -1 run/03_spin_elk/MOMENTM.OUT
number of k points= 72 Methfessel-Paxton smearing, width (Ry)= 0.0200
! total energy = -329.69603739 Ry
total magnetization = 0.84 Bohr mag/cell
absolute magnetization = 0.91 Bohr mag/cell
convergence has been achieved in 10 iterations
atom 1 (R=0.413) charge= 16.3251 magn= 0.8727
run/01_spin_up_qe/density_pp_2d.out: Min, Max, imaginary charge: 0.014076 2.252026 0.000000
run/02_spin_dn_qe/density_pp_2d.out: Min, Max, imaginary charge: 0.014930 2.143502 0.000000
Total number of k-points : 72
Convergence targets achieved
0.6929144419
QE converges in 10 iterations (about 21 seconds on one core) with 72 k-points; Elk uses the same 72. QE's total magnetization is 0.84 Bohr magnetons per cell, Elk's total moment is 0.693 (the last line of MOMENTM.OUT). The two are discussed in the analysis.
Plot
The script plot_spin_density.py (in the example folder) reads the two QE files, subtracts them, reads Elk's MAG2D.OUT, and draws the two spin densities and their difference on the same plane:
./plot_spin_density.py # default colour scale: +-0.05 e/bohr^3 ./plot_spin_density.py --vmax 0.01 -o spin_density_fine.png
plane edges: 11.5216 x 8.1470 bohr; Ni atoms in the plane: 3 QE : min -0.0067 max 0.2223 plane mean 0.01634 Elk: min -0.0074 max 0.3923 plane mean 0.01756 (sign flipped) QE - Elk on the common 99 x 99 points: rms 0.03640, mean -0.00124 wrote spin_density.png
Note: How the script reads Elk's file. MAG2D.OUT holds, at each point, three numbers: the magnetisation vector projected on the plane, as its components along edge 1, along edge 2 and along the normal. For a magnet along z the spin density is the z component of that vector, and the script rebuilds it from the plane edges (read from the .qelk file). The other two components of the magnetisation vector (x, y) come out zero, as they should for a collinear magnet. Sign: with bfieldc along +z, Elk's z component came out negative; the overall sign of a spin density is a convention (reversing every spin changes nothing physical), so the script multiplies the Elk map by the sign that makes its plane average positive, as QE's majority-spin-positive map already is, and says so in the title.
Analysis
1. The spin density sits on the atoms. Both maps show the polarisation concentrated around each Ni atom, from the 3d shell, as positive blobs separated by regions where it is close to zero. The finer scale shows that between the atoms the spin density is slightly negative in both codes: the minimum is -0.0067 e/bohr3 in QE and -0.0074 in Elk, about thirty times smaller than the maximum.
2. Inside the atoms the codes do not agree, and should not. The largest values of the spin density are 0.222 e/bohr3 in QE and 0.392 in Elk. Elk is all-electron; the QE pseudopotential describes only the valence electrons, with a smooth pseudo-wavefunction inside the core region, so the spin density near the nucleus is not the same quantity. The difference map shows exactly that: a ring and a disc around each atom, and little elsewhere. The rms of the difference over the 99 x 99 points, 0.036, comes almost entirely from those regions. The plane averages are close: 0.0163 e/bohr3 (QE) and 0.0176 (Elk).
3. The total moments differ by 0.15. QE gives a total magnetization of 0.84 Bohr magnetons per atom, Elk 0.69. The measured moment of nickel is about 0.6. The QE value is high, and so was the QE value for iron in the Fe_magnetism_elk_vs_qe example. For this pseudopotential file, Ni_PBE_TM_2pj, the header says that no nonlinear core correction was used, which is known to overestimate magnetic moments, so it is a likely contributor. We have not tested it: a run with a pseudopotential that includes the core correction, or with a different smearing, would show how much of the 0.15 it explains.
4. Next. The same input with task 0 73 in Elk and a 3D density step in QE gives the spin density in the whole cell, for VESTA (RHO3D-style conversion for MAG3D.OUT is not in the suite yet).
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
0 Comments