Calculation of the layer-by-layer magnetic moments of a Ni(001) slab with QE code

The layer-by-layer magnetic moments of a Ni(001) slab, computed with Quantum ESPRESSO alone -- a type 1, independent calculation with the QELK suite. A spin-polarised five-layer slab with vacuum on both sides shows how the moment of the surface layers differs from that of the layers below. Like the Cu_slab_workfunction_qe example, it is a non-bulk structure: the cell is periodic along z, and the slab is the crystal with a large empty region.

 

Before you start: one small update of the suite

 

The input uses mixing_beta, the charge-mixing parameter of pw.x, which the suite did not pass through. It is part of the same in-place update script as the spin-density example (apply_spin_density_patch.py; it keeps a FILE.bak_spin backup of each file and can be run twice safely):

algerien1970@linux-ml0a:~> python3 apply_spin_density_patch.py ~/abinitio/QELK

 

The input file

 

We need only the input file

Ni_slab_surface_moment_qe.qelk

version 1

cell bohr
  4.7037  0.0000   0.0000
  0.0000   4.7037  0.0000
  0.0000   0.0000   28.0000
end

atoms crystal
  Ni  0.0000  0.0000  0.26242857
  Ni  0.5000  0.5000  0.38121429
  Ni  0.0000  0.0000  0.50000000
  Ni  0.5000  0.5000  0.61878571
  Ni  0.0000  0.0000  0.73757143
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         50
  ecutrho         200
  kpoints         8 8 1
  smearing        mv
  mixing_beta     0.2
  degauss         0.02
  extra_system    nspin=2,starting_magnetization(1)=0.3
end

Note: The slab geometry. Ni is fcc with a = 6.652 bohr (3.52 Å). The (001) surface cell is a square of side a/√2 = 4.7037 bohr with one atom per layer; successive layers are a/2 = 3.326 bohr apart and shifted by (1/2, 1/2) in the surface plane. Five layers stand in the middle of a cell of c = 28 bohr, with the middle layer at z = 1/2, so atoms 1 and 5 are the two equivalent surfaces; the vacuum between the periodic images is 28 - 4 x 3.326 = 14.7 bohr thick. The layers are at their bulk positions: no relaxation.

Note: The magnetism is switched on in QE with extra_system nspin=2,starting_magnetization(1)=0.3. All five atoms are one species and start with the same moment. The mesh is 8 x 8 x 1 (no sampling along the vacuum direction) and the cutoff is 50 Ry for the wavefunctions and 200 Ry for the density, which is enough for this norm-conserving pseudopotential, so ecutrho is set by hand instead of the 8 x ecutwfc that the suite writes by default for ultrasoft and PAW files.

Note: Getting the SCF to converge. A magnetic metal slab is a harder SCF problem than the bulk. The first run (Methfessel-Paxton smearing, starting magnetization 0.5, mixing_beta 0.7, the settings of the bulk examples) did not converge: the estimated accuracy stayed between 0.1 and several Ry and the total magnetization wandered between 4 and 6 Bohr magnetons in 30 iterations. A second try with Marzari-Vanderbilt smearing and a starting magnetization of 0.3 was still not converging after 21 iterations when we stopped it. The input above, which also lowers mixing_beta to 0.2, converged in 38 iterations. We changed these parameters together and have not tested which one is responsible.

 

Calculation

 

algerien1970@linux-ml0a:~/abinitio/QELK/examples/Ni_slab_surface_moment_qe> qelk Ni_slab_surface_moment_qe.qelk

One step, in run/00_scf; scf.pw.in and scf.pw.out are copied next to the .qelk file. On one core it took about six minutes.

 

SCF results

 

grep -E "convergence has been achieved|^!|number of k points|Fermi energy" scf.pw.out
grep -E "total magnetization|absolute magnetization" scf.pw.out | tail -2
grep -A5 "Magnetic moment per site" scf.pw.out | tail -6

     number of k points=    15  Marzari-Vanderbilt smearing, width (Ry)=  0.0200
     the Fermi energy is    10.6347 ev
!    total energy              =   -1631.58228862 Ry
     convergence has been achieved in  38 iterations
     total magnetization       =     5.03 Bohr mag/cell
     absolute magnetization    =     5.36 Bohr mag/cell

The last "Magnetic moment per site" block gives the charge and the moment inside a sphere around each atom (radius 0.412 in units of alat, about 1.94 bohr):

atom	layer	charge in sphere	moment (muB)
1	surface	16.0716	1.1242
2	subsurface	16.2260	0.9612
3	centre	16.2401	0.9654
4	subsurface	16.2260	0.9612
5	surface	16.0716	1.1242

 

Plot

 

The script layer_moments_qe.py (in the example folder) reads the last "Magnetic moment per site" block of scf.pw.out and plots the moment of each layer:

./layer_moments_qe.py

Analysis

 

1. The surface moment is larger. The two surface atoms carry 1.124 Bohr magnetons in the sphere, against 0.961 and 0.965 for the subsurface and central layers: about 16% more. The two surfaces are equivalent by symmetry and come out identical, as they should. The central layer and the layer below it differ by only 0.004, so by the second layer the moment is already bulk-like in this slab; five layers cannot say more than that.

2. The charge in the sphere is lower at the surface (16.07 against 16.23 to 16.24 electrons): the surface atoms have a free side, and part of their charge extends into the vacuum, outside the sphere.

3. The moments of the spheres add to more than the total. The five site moments sum to 5.14 Bohr magnetons, the total magnetization of the cell is 5.03, so the region outside the spheres carries about -0.1: the same slightly negative spin density between the atoms that we saw in the bulk Ni example.

4. Absolute values need care. The moments in the sphere depend on its radius, and QE's moment for bulk Ni in the Ni_spin_density2D_elk_vs_qe example is 0.84 per atom against 0.69 from Elk; this norm-conserving pseudopotential has no nonlinear core correction, which is a likely cause of the overestimate (not tested). The layer-to-layer trend is therefore more reliable than the numbers. The sphere moment of the bulk example, 0.873 for the same radius, comes from a different smearing, cutoff and mesh and is not a fair reference for the central layer here (0.965): for a real comparison, run bulk Ni with the same settings.

5. Limits. One slab thickness, one vacuum thickness, one cutoff, one mesh and no relaxation of the surface layers. Each is a one-line change and one more run; a thicker slab would also show how far the surface perturbation reaches.