可视化散点图怎么分析？布里渊区及能带路径的可视化
可视化散点图怎么分析？布里渊区及能带路径的可视化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(Gree
先前本公众号已经介绍"五种方法生成能带结构计算的高对称点"，使用里面的方法即可获取相应晶体对应的布里渊区图。在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)
