# 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(restart=name + '.json')

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')