PBE0 calculations for bulk silicon¶

This tutorial will do non self-consistent PBE0 based on self-consistent PBE.

PBE and PBE0 band gaps¶

The band structure can be calculated like this:

# web-page: si-gaps.csv
from ase.build import bulk
from ase.parallel import paropen
from gpaw.hybrids.eigenvalues import non_self_consistent_eigenvalues
from gpaw import GPAW, PW

a = 5.43
si = bulk('Si', 'diamond', a)

fd = paropen('si-gaps.csv', 'w')

for k in range(2, 9, 2):
    name = f'Si-{k}'
    si.calc = GPAW(kpts={'size': (k, k, k), 'gamma': True},
                   mode=PW(200),
                   xc='PBE',
                   convergence={'bands': 5},
                   txt=name + '.txt')
    si.get_potential_energy()
    si.calc.write(name + '.gpw', mode='all')

    # Range of eigenvalues:
    n1 = 3
    n2 = 5

    ibzkpts = si.calc.get_ibz_k_points()
    kpt_indices = []
    for kpt in [(0, 0, 0), (0.5, 0.5, 0)]:  # Gamma and X
        # Find k-point index:
        i = abs(ibzkpts - kpt).sum(1).argmin()
        kpt_indices.append(i)

    # Do PBE0 calculation:
    epbe, vpbe, vpbe0 = non_self_consistent_eigenvalues(
        name + '.gpw',
        'PBE0',
        n1, n2,
        kpt_indices,
        snapshot=name + '.json')

    epbe0 = epbe - vpbe + vpbe0

    gg = epbe[0, 0, 1] - epbe[0, 0, 0]
    gx = epbe[0, 1, 1] - epbe[0, 0, 0]
    gg0 = epbe0[0, 0, 1] - epbe0[0, 0, 0]
    gx0 = epbe0[0, 1, 1] - epbe0[0, 0, 0]

    print(f'{k}, {gg:.3f}, {gx:.3f}, {gg0:.3f}, {gx0:.3f}', file=fd)
    fd.flush()

assert abs(gg - 2.559) < 0.01
assert abs(gx - 0.707) < 0.01
assert abs(gg0 - 3.873) < 0.01
assert abs(gx0 - 1.828) < 0.01

These are the resulting \(\Gamma\)-\(\Gamma\) and \(\Gamma\)-\(X\) gaps for PBE and PBE0 in eV:

k-points	PBE(G-G)	PBE(G-X)	PBE0(G-G)	PBE0(G-X)
2	2.445	0.598	4.187	2.153
4	2.538	0.688	3.954	1.913
6	2.553	0.701	3.892	1.849
8	2.555	0.704	3.874	1.831

Lattice constant and bulk modulus¶

Here is how to calculate the lattice constant:

import ase.db
from ase.build import bulk
import numpy as np
from gpaw.hybrids.energy import non_self_consistent_energy as nsc_energy
from gpaw import GPAW, PW

a0 = 5.43

con = ase.db.connect('si.db')

for k in range(2, 9):
    for a in np.linspace(a0 - 0.04, a0 + 0.04, 5):
        id = con.reserve(a=a, k=k)
        if id is None:
            continue
        si = bulk('Si', 'diamond', a)
        si.calc = GPAW(kpts=(k, k, k),
                       mode=PW(400),
                       xc='PBE',
                       eigensolver='rmm-diis',
                       txt=None)
        si.get_potential_energy()
        name = f'si-{a:.2f}-{k}'
        si.calc.write(name + '.gpw', mode='all')
        epbe0 = nsc_energy(name + '.gpw', 'PBE0').sum()

        con.write(si, a=a, k=k, epbe0=epbe0)
        del con[id]

Plot the results like this:

# web-page: si-a.png
import matplotlib.pyplot as plt
import ase.db
from ase.eos import EquationOfState


def lattice_constant(volumes, energies):
    eos = EquationOfState(volumes, energies)
    v, e, B = eos.fit()
    a = (v * 4)**(1 / 3)
    return a


con = ase.db.connect('si.db')
results = []
K = list(range(2, 9))
A = []
A0 = []
for k in K:
    rows = list(con.select(k=k))
    V = [row.volume for row in rows]
    E = [row.energy for row in rows]
    E0 = [row.epbe0 for row in rows]
    A.append(lattice_constant(V, E))
    A0.append(lattice_constant(V, E0))

print(K, A, A0)

plt.plot(K, A, label='PBE')
plt.plot(K, A0, label='PBE0')
plt.xlabel('number of k-points')
plt.ylabel('lattice constant [Ang]')
plt.savefig('si-a.png')