Coverage for gpaw/test/test_jellium.py: 96%
28 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1import pytest
2import numpy as np
3from ase import Atoms
4from ase.units import Bohr
5from gpaw.jellium import JelliumSlab
6from gpaw import GPAW, Mixer
8rs = 5.0 * Bohr # Wigner-Seitz radius
9h = 0.24 # grid-spacing
10a = 8 * h # lattice constant
11v = 3 * a # vacuum
12L = 8 * a # thickness
13k = 6 # number of k-points (k*k*1)
15ne = a**2 * L / (4 * np.pi / 3 * rs**3)
18@pytest.mark.libxc
19def test_jellium(in_tmp_dir, gpaw_new):
20 x = h / 4 # make sure surfaces are between grid-points
21 bc = JelliumSlab(ne, z1=v - x, z2=v + L + x)
23 surf = Atoms(pbc=(True, True, False),
24 cell=(a, a, v + L + v))
25 params = dict(
26 mode='fd',
27 poissonsolver={'dipolelayer': 'xy'},
28 xc='LDA_X+LDA_C_WIGNER',
29 eigensolver='dav',
30 kpts=[k, k, 1],
31 h=h,
32 maxiter=300,
33 convergence={'density': 1e-5},
34 mixer=Mixer(0.3, 7, 100),
35 nbands=int(ne / 2) + 15)
36 if gpaw_new:
37 surf.calc = GPAW(**params, environment=bc)
38 else:
39 surf.calc = GPAW(**params, background_charge=bc)
41 _ = surf.get_potential_energy()
43 # Get the work function
44 v_r = surf.calc.get_electrostatic_potential()
45 efermi = surf.calc.get_fermi_level()
46 phi = v_r[:, :, -1].mean() - efermi
47 assert phi == pytest.approx(2.715, abs=1e-3)
48 # Reference value: Lang and Kohn, 1971, Theory of Metal Surfaces:
49 # Work function
50 # DOI 10.1103/PhysRevB.3.1215
51 # r_s = 5, work function = 2.73 eV
53 surf.calc.write('jellium.gpw')