PBE0 calculations for bulk silicon

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

PBE and PBE0 band gaps

The band structure can be calculated like this:

# Creates: si-gaps.csv
from __future__ import print_function
import numpy as np
from ase.build import bulk
from ase.parallel import paropen
from gpaw.xc.exx import EXX
from gpaw.xc.tools import vxc
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 = 'Si-{0}'.format(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 = []
    pbeeigs = []
    for kpt in [(0, 0, 0), (0.5, 0.5, 0)]:
        # Find k-point index:
        i = abs(ibzkpts - kpt).sum(1).argmin()
        kpt_indices.append(i)
        pbeeigs.append(si.calc.get_eigenvalues(i)[n1:n2])

    # DFT eigenvalues:
    pbeeigs = np.array(pbeeigs)
    
    # PBE contribution:
    dpbeeigs = vxc(si.calc, 'PBE')[0, kpt_indices, n1:n2]
        
    # Do PBE0 calculation:
    pbe0 = EXX(name + '.gpw',
               'PBE0',
               kpts=kpt_indices,
               bands=[n1, n2],
               txt=name + '.pbe0.txt')
    pbe0.calculate()
    
    dpbe0eigs = pbe0.get_eigenvalue_contributions()[0]
    pbe0eigs = pbeeigs - dpbeeigs + dpbe0eigs

    print('{0}, {1:.3f}, {2:.3f}, {3:.3f}, {4:.3f}'
          .format(k,
                  pbeeigs[0, 1] - pbeeigs[0, 0],
                  pbeeigs[1, 1] - pbeeigs[0, 0],
                  pbe0eigs[0, 1] - pbe0eigs[0, 0],
                  pbe0eigs[1, 1] - pbe0eigs[0, 0]),
          file=fd)
    fd.flush()

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.xc.exx import EXX
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 = 'si-{0:.2f}-{1}'.format(a, k)
        si.calc.write(name + '.gpw', mode='all')
        pbe0 = EXX(name + '.gpw', 'PBE0', txt=name + '.pbe0.txt')
        pbe0.calculate()
        epbe0 = pbe0.get_total_energy()
        
        con.write(si, a=a, k=k, epbe0=epbe0)
        del con[id]

Plot the results like this:

# Creates: si-a.png
from __future__ import division
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')
../../_images/si-a.png