1import pyscf.pbc.tools.pyscf_ase as pyscf_ase
2import pyscf.pbc.gto as pbcgto
3import pyscf.pbc.dft as pbcdft
4
5import matplotlib.pyplot as plt
6
7from ase.lattice import bulk
8from ase.dft.kpoints import sc_special_points as special_points, get_bandpath
9
10c = bulk('C', 'diamond', a=3.5668)
11print(c.get_volume())
12
13cell = pbcgto.Cell()
14cell.atom = pyscf_ase.ase_atoms_to_pyscf(c)
15cell.a = c.cell
16
17cell.basis = 'gth-szv'
18cell.pseudo = 'gth-pade'
19cell.verbose = 5
20cell.build(None,None)
21
22points = special_points['fcc']
23G = points['G']
24X = points['X']
25W = points['W']
26K = points['K']
27L = points['L']
28band_kpts, kpath, sp_points = get_bandpath([L, G, X, W, K, G], c.cell, npoints=50)
29band_kpts = cell.get_abs_kpts(band_kpts)
30
31#
32# band structure from Gamma point sampling
33#
34
35mf = pbcdft.RKS(cell)
36print(mf.kernel())
37
38e_kn = mf.get_bands(band_kpts)[0]
39vbmax = -99
40for en in e_kn:
41    vb_k = en[cell.nelectron//2-1]
42    if vb_k > vbmax:
43        vbmax = vb_k
44e_kn = [en - vbmax for en in e_kn]
45
46#
47# band structure from 222 k-point sampling
48#
49
50kmf = pbcdft.KRKS(cell, cell.make_kpts([2,2,2]))
51print(kmf.kernel())
52
53e_kn_2 = kmf.get_bands(band_kpts)[0]
54vbmax = -99
55for en in e_kn_2:
56    vb_k = en[cell.nelectron//2-1]
57    if vb_k > vbmax:
58        vbmax = vb_k
59e_kn_2 = [en - vbmax for en in e_kn_2]
60
61au2ev = 27.21139
62
63emin = -1*au2ev
64emax = 1*au2ev
65
66plt.figure(figsize=(5, 6))
67nbands = cell.nao_nr()
68for n in range(nbands):
69    plt.plot(kpath, [e[n]*au2ev for e in e_kn], color='#87CEEB')
70    plt.plot(kpath, [e[n]*au2ev for e in e_kn_2], color='#4169E1')
71for p in sp_points:
72    plt.plot([p, p], [emin, emax], 'k-')
73plt.plot([0, sp_points[-1]], [0, 0], 'k-')
74plt.xticks(sp_points, ['$%s$' % n for n in ['L', r'\Gamma', 'X', 'W', 'K', r'\Gamma']])
75plt.axis(xmin=0, xmax=sp_points[-1], ymin=emin, ymax=emax)
76plt.xlabel('k-vector')
77
78plt.show()
79
80