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