Jellium

In this tutorial, we try to reproduce some old jellium calculations by Lang and Kohn [Lang70].

Bulk

Let’s do a calculation for \(r_s=5\) Bohr. We use a cubic cell with a lattice constant of \(a=1.6\) Å, 8*8*8 grid points and and a k-point sampling of 12*12*12 points.

import numpy as np
from ase import Atoms
from ase.units import Bohr
from gpaw.jellium import Jellium
from gpaw import GPAW, PW

rs = 5.0 * Bohr  # Wigner-Seitz radius
h = 0.2          # grid-spacing
a = 8 * h        # lattice constant
k = 12           # number of k-points (k*k*k)

ne = a**3 / (4 * np.pi / 3 * rs**3)
jellium = Jellium(ne)

bulk = Atoms(pbc=True, cell=(a, a, a))
bulk.calc = GPAW(mode=PW(400.0),
                 background_charge=jellium,
                 xc='LDA_X+LDA_C_WIGNER',
                 nbands=5,
                 kpts=[k, k, k],
                 h=h,
                 txt='bulk.txt')
e0 = bulk.get_potential_energy()

In the text output from the calculation, one can see that the cell contains 0.0527907 electrons.

Surfaces

Now, we will do a surface calculation. We put a slab of thickness 16.0 Å inside a box, periodically repeated in the x- and y-directions only, of size 1.6*1.6*25.6 Å and use 12*12 k-points to sample the surface BZ:

import json
from pathlib import Path
import numpy as np
from ase import Atoms
from ase.units import Bohr
from gpaw import GPAW, PW
from gpaw.core import UGArray
from gpaw.new.extensions import Jellium

rs = 5.0 * Bohr  # Wigner-Seitz radius
h = 0.2          # grid-spacing
a = 8 * h        # lattice constant
v = 3 * a        # vacuum
L = 10 * a       # thickness
k = 12           # number of k-points (k*k*1)

ne = a**2 * L / (4 * np.pi / 3 * rs**3)


class JelliumSlab(Jellium):
    def update_mask(self, mask_r: UGArray) -> None:
        z = mask_r.desc.xyz()[:, :, :, 2] * Bohr
        eps = 0.001  # displace surfaces away from grid-points
        mask_r.data[np.logical_and(z > v - eps, z < v + L + eps)] = 1.0


surf = Atoms(pbc=(True, True, False),
             cell=(a, a, v + L + v))
surf.calc = GPAW(mode=PW(400.0),
                 extensions=[JelliumSlab(charge=ne)],
                 xc='LDA_X+LDA_C_WIGNER',
                 kpts=[k, k, 1],
                 h=h,
                 convergence={'density': 0.001},
                 nbands=int(ne / 2) + 15,
                 txt='surface.txt')
e = surf.get_potential_energy()
density = surf.calc.get_pseudo_density()[0, 0]
Path('surface.json').write_text(json.dumps(density.tolist()))

The surface energy is:

>>> from ase.io import read
>>> e0 = read('bulk.txt').get_potential_energy()
>>> e = read('surface.txt').get_potential_energy()
>>> a = 1.6
>>> L = 10 * a
>>> sigma = (e - L / a * e0) / 2 / a**2
>>> print('%.2f mev/Ang^2' % (1000 * sigma))
5.39 mev/Ang^2
>>> print('%.1f erg/cm^2' % (sigma / 6.24150974e-5))
86.4 erg/cm^2

which is reasonably close to the value from Lang and Kohn: 100 ergs/cm\(^2\).

Here is the electron density profile:

# web-page: fig2.png
import json
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from ase.units import Bohr

rs = 5.0 * Bohr
density = np.array(json.loads(Path('surface.json').read_text()))
h = 0.2
a = 8 * h
v = 3 * a
L = 10 * a
z = np.linspace(0, v + L + v, len(density), endpoint=False)
# Position of surface is between two grid points:
z0 = (v + L - h / 2)
n = 1 / (4 * np.pi / 3 * rs**3)  # electron density
kF = (3 * np.pi**2 * n)**(1.0 / 3)
lambdaF = 2 * np.pi / kF  # Fermi wavelength
plt.figure(figsize=(6, 6 / 2**0.5))
plt.plot([-L / lambdaF, -L / lambdaF, 0, 0], [0, 1, 1, 0], 'k')
plt.plot((z - z0) / lambdaF, density / n)
# plt.xlim(xmin=-1.2, xmax=1)
plt.ylim(ymin=0)
plt.title(rf'$r_s={rs / Bohr:.1f}\ \mathrm{{Bohr}},\ '
          rf'\lambda_F={lambdaF / Bohr:.1f}\ \mathrm{{Bohr}}$')
plt.xlabel('DISTANCE (FERMI WAVELENGTHS)')
plt.ylabel('ELECTRON DENSITY')
plt.savefig('fig2.png')

Compare with Fig. 2 in [Lang70]:

Other jellium geometries

For other geometries, one will have to subclass Jellium, and implement the get_mask() method:

class gpaw.jellium.Jellium(charge)

The Jellium object

Initialize the Jellium object

Input: charge, the total Jellium background charge.

add_charge_to(rhot_g)

Add Jellium background charge to pseudo charge density rhot_g

get_mask()

Choose which grid points are inside the jellium.

gd: grid descriptor

Return ndarray of ones and zeros indicating where the jellium is. This implementation will put the positive background in the whole cell. Overwrite this method in subclasses.

set_grid_descriptor(gd)

Set the grid descriptor for the Jellium background charge

The JelliumSlab is one such example of a subclass of gpaw.jellium.Jellium:

class gpaw.jellium.JelliumSlab(charge, z1, z2)

The Jellium slab object

Put the positive background charge where z1 < z < z2.

z1: float

Position of lower surface in Angstrom units.

z2: float

Position of upper surface in Angstrom units.

get_mask()

Choose which grid points are inside the jellium.

gd: grid descriptor

Return ndarray of ones and zeros indicating where the jellium is. This implementation will put the positive background in the whole cell. Overwrite this method in subclasses.

---

[Lang70] (1,2)

N. D. Lang and W. Kohn, Phys. Rev. B 1, 4555-4568 (1970), Theory of Metal Surfaces: Charge Density and Surface Energy, doi: 10.1103/PhysRevB.1.4555