布里渊区及能带路径的可视化

先前本公众号已经介绍"五种方法生成能带结构计算的高对称点"，使用里面的方法即可获取相应晶体对应的布里渊区图。在VASPKIT 1.20新版本中即将实现布里渊区及能带路径的可视化。相应的脚本文件(见文末代码)已经开放内测，感兴趣的可以到VASPKIT FAQs QQ 群331895604群文件下载。相比其它软件，VASPKIT可以更加方便实现自动标记对应能带计算的k点路径。

晶体结构的第一布里渊区可通过调用Voronoi图算法实现，特别感谢中科大郑奇靖老师提供了关于实现算法的建议。接下来我们演示如何调用VASPKIT结合python脚本实现布里渊区及能带路径的可视化：

第一步：我们以Cu和graphene为例，准备好POSCAR文件，注意一定是原胞结构(primtive cell)。然后分别运行vaspkit -303 (Cu 体相)和vaspkit -302 (二维体系graphene)生成包含推荐能带路径的KPATH.in文件和高对称点HIGH_SYMMETRY_POINTS文件;

===================== Structural Options ======================== 01) VASP Input Files Generator 02) Elastic-Properties 03) K-Path Generator 04) Structure Editor 05) Catalysis-ElectroChem Kit 06) Symmetry Search ===================== Electronic Options ======================== 11) Density-of-States 21) DFT Band-Structure 23) 3D Band-Structure 25) Hybrid-DFT Band-Structure 26) Fermi-Surface 28) Band-Structure Unfolding =========== Charge & Potential & Wavefunction Options =========== 31) Charge & Spin Density 42) Potential-Related 51) Wave-Function Analysis ====================== Misc Utilities =========================== 71) Optical-Properties 72) Molecular-Dynamics Kit 73) VASP2other Interface 91) Semiconductor Calculator 92) 2D-Materials Kit 0) Quit------------>>302 +-------------------------- Warm Tips --------------------------+ See An Example in vaspkit/examples/seek_kpath/graphene_2D. More details are given in arXiv:1806.04285 (2018) and refs. This feature is still experimental & check PRIMCELL.vasp file. +---------------------------------------------------------------+-->> (01) Reading Structural Parameters from POSCAR File... +-------------------------- Summary ----------------------------+ The vacuum slab is supposed to be along z axis Prototype: AB2 Total Atoms in Input Cell: 3 Lattice Constants in Input Cell: 3.190 3.190 12.000 Lattice Angles in Input Cell: 90.000 90.000 120.000 Total Atoms in Primitive Cell: 3 Lattice Constants in Primitive Cell: 3.190 3.190 12.000 Lattice Angles in Primitive Cell: 90.000 90.000 120.0002D Bravais Lattice: Hexagonal Point Group: 26 [ D3h ] International: P-6m2 Symmetry Operations: 12 Suggested K-Path: (shown in the next line) [ GAMMA-M-K-GAMMA ] +---------------------------------------------------------------+-->> (02) Written PRIMCELL.vasp file.-->> (03) Written HIGH_SYMMETRY_POINTS File for Reference.-->> (04) Written KPATH.in File for Band-Structure Calculation. +----------------------------WARNING----------------------------+ | Do NOT forget to copy PRIMCELL.vasp to POSCAR unless you know | | what you are doing. Otherwise you might get wrong results! | +---------------------------------------------------------------+`

第二步，终端运行python3 visualize_BZ.py (代码见下面)得到下图：

import numpy as np, numpy.linalg as nplfrom scipy.spatial import Voronoifrom itertools import productfrom matplotlib.patches import FancyArrowPatchfrom mpl_toolkits.mplot3d import proj3dimport matplotlib as mplmpl.rcParams['font.size'] = 12.

def read_poscar(poscar, species=None): poscar = open(poscar,'r') title = poscar.readline().strip() scale = float(poscar.readline().strip()) s = float(scale) lattice_vectors = [[ float(v) for v in poscar.readline().split() ], [ float(v) for v in poscar.readline().split() ], [ float(v) for v in poscar.readline().split() ]] lattice_vectors = np.array(lattice_vectors) reciprocal_lattice_vectors= np.linalg.inv(lattice_vectors).T reciprocal_lattice_vectors=reciprocal_lattice_vectors*np.pi*2return reciprocal_lattice_vectors

def read_high_symmetry_points(): f = open("HIGH_SYMMETRY_POINTS",'r') fa=f.readline() n=0 hpts=[] while True: fa=f.readline() hpts.append(fa) n=n+1if fa.find('use') > 0: break f.close() kpts=[] klabels=[] for i in range(0,n-3): kpts.append(hpts[i].split()[0:3]) klabels.append(hpts[i].split()[3]) kpts=np.array(kpts,dtype=np.float64)return kpts, klabels

def is_greek_alphabets(klabels): Greek_alphabets=['Alpha','Beta','Gamma','Delta','Epsilon','Zeta','Eta','Theta', 'Iota','Kappa','Lambda','Mu','Nu','Xi','Omicron','Pi','Rho','Sigma','Tau','Upsilon','Phi','Chi','Psi','Pega'] group_labels=[]for i in range(len(klabels)): klabel=klabels[i] for j in range(len(Greek_alphabets)):if (klabel.find(Greek_alphabets[j].upper())>=0): latex_exp=r''+'$\\'+str(Greek_alphabets[j])+'$' klabel=klabel.replace(str(Greek_alphabets[j].upper()),str(latex_exp))if (klabel.find('_')>0): n=klabel.find('_') klabel=klabel[:n]+'$'+klabel[n:n+2]+'$'+klabel[n+2:] group_labels.append(klabel) klabels=group_labelsreturn klabels

def read_kpath(): kpath=np.loadtxt("KPATH.in", dtype=np.string_,skiprows=4)#print(kpath) kpath_labels = kpath[:,3].tolist() kpath_labels = [i.decode('utf-8','ignore') for i in kpath_labels]for i in range(len(kpath_labels)):if kpath_labels[i]=="Gamma": kpath_labels[i]=u"Γ" kpaths=np.zeros((len(kpath_labels),3),dtype=float)for i in range(len(kpath_labels)): kpaths[i,:]=\ [float(x) for x in kpath[i][0:3]]return kpath_labels, kpaths

def get_Wigner_Seitz_BZ(lattice_vectors):# Inspired by http://www.thp.uni-koeln.de/trebst/Lectures/SolidState-2016/wigner_seitz_3d.py # Inspired by https://github.com/QijingZheng/VASP_FermiSurface/blob/master/fs.py latt = [] prefactors = [0., -1., 1.]for p in prefactors:for u in lattice_vectors: latt.append(p * u) lattice = []for vs in product(latt, latt, latt): a = vs[0] + vs[1] + vs[2]if not any((a == x).all() for x in lattice): lattice.append(a) voronoi = Voronoi(lattice) bz_facets = [] bz_ridges = [] bz_vertices = []for pid, rid in zip(voronoi.ridge_points, voronoi.ridge_vertices):if(pid[0] == 0 or pid[1] == 0): bz_ridges.append(voronoi.vertices[np.r_[rid, [rid[0]]]]) bz_facets.append(voronoi.vertices[rid]) bz_vertices += rid bz_vertices = list(set(bz_vertices))return voronoi.vertices[bz_vertices], bz_ridges, bz_facets

class Arrow3D(FancyArrowPatch):def __init__(self, xs, ys, zs, *args, **kwargs): FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs) self._verts3d = xs, ys, zsdef draw(self, renderer): xs3d, ys3d, zs3d = self._verts3d xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M) self.set_positions((xs[0], ys[0]), (xs[1], ys[1])) FancyArrowPatch.draw(self, renderer)

def visualize_BZ_matplotlib(points,ridges,facets,reciprocal_lattice_vectors,kpts,klabels,kpaths):import matplotlib as mplimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection='3d') basis_vector_clrs = ['r', 'g', 'b'] basis_vector_labs = ['x', 'y', 'z']for ii in range(3): arrow = Arrow3D([0,reciprocal_lattice_vectors[ii, 0]], [0,reciprocal_lattice_vectors[ii, 1]], [0,reciprocal_lattice_vectors[ii, 2]], color=basis_vector_clrs[ii], mutation_scale=20,lw=1,arrowstyle="->") ax.add_artist(arrow) ax.text(reciprocal_lattice_vectors[ii, 0], reciprocal_lattice_vectors[ii, 1],reciprocal_lattice_vectors[ii, 2], basis_vector_labs[ii])for ir in ridges: ax.plot(ir[:, 0], ir[:, 1], ir[:, 2], color='k', lw=1.5,alpha=0.5)for i in range(len(klabels)): kpt=np.dot(kpts[i,:],reciprocal_lattice_vectors) ax.scatter(kpt[0], kpt[1], kpt[2],c='b', marker='o',s=20,alpha=0.8) ax.text(kpt[0], kpt[1], kpt[2],klabels[i],c='b')for i in range(kpaths.shape[0]): kpaths[i,:]=np.dot(kpaths[i,:],reciprocal_lattice_vectors)for i in range(0,kpaths.shape[0],2): arrow = Arrow3D([kpaths[i,0],kpaths[i+1,0]],[kpaths[i,1],kpaths[i+1,1]],[kpaths[i,2],kpaths[i+1,2]],mutation_scale=20,lw=1.5,arrowstyle="->", color="magenta") ax.add_artist(arrow) ax.set_axis_off() ax.view_init(elev=12, azim=23) plt.savefig('Brillouin_Zone.pdf',dpi=100) plt.show()

def welcome(): print('') print('+---------------------------------------------------------------+') print('| A VASPKIT Plugin to Visualize Brillouin Zone Using Matplotlib |') print('| Written by Vei WANG (wangvei@icloud.com) |') print('+---------------------------------------------------------------+') print('')

if __name__ == "__main__": welcome() reciprocal_lattice_vectors=read_poscar('POSCAR') lattice_vectors=[np.array(reciprocal_lattice_vectors[0,:]),np.array(reciprocal_lattice_vectors[1,:]),np.array(reciprocal_lattice_vectors[2,:])] kpts,klabels=read_high_symmetry_points() klabels=is_greek_alphabets(klabels) kpath_labels,kpaths=read_kpath() points, ridges, facets = get_Wigner_Seitz_BZ(lattice_vectors)visualize_BZ_matplotlib(points, ridges, facets, reciprocal_lattice_vectors,kpts,klabels,kpaths)