A python script that draws a band diagram from the VASP output file EIGENVAL
happy New Year. This is Hashimoto.
The band diagram is important information for understanding the dynamic properties of material properties. Electrons are indispensable especially in the field of semiconductor materials, such as which surface is easy to conduct electricity and which surface should be joined to make a device. (Maybe)
Recently, it has become possible to obtain band diagrams directly with ARPES, but (https://arpes.phys.tohoku.ac.jp/contents/study.html) Prediction by numerical calculation is reasonable, and the first It can be said that simulation by one-principles calculation is still active. In addition, the first-principles calculation package VASP provides a calculation method that uses Green's functions to improve accuracy, and vigorous research is underway.
In this article, I will introduce a Python script that draws a band diagram from the "EIGENVAL" file, which is the output file of VASP.
0. Premise
If you are not interested in manually organizing the contents of the results, for example, use pymatgen "https://pymatgen.org/" or p4vasp "http://www.p4vasp.at/#/". You can output a band diagram with a black box.
For pymatgen, the official reference and the article by RKS WEBSITE "http://ryokbys.web.nitech.ac.jp/pymatgen.html" will be helpful. The method for p4vasp is described on the VASP official website "https://www.vasp.at/wiki/index.php/Fcc_Si_bandstructure".
That said, honestly, it's subtle when asked if it's easy to use. Especially when there are more than 100 bands, I can't draw well. In this article, the goal is "Because it is not such a difficult format, you can extract EIGENVAL with python and draw it manually".
1. Get EIGENVAL
You need to run the calculation twice (or three times) to get the band diagram.
Step1.1 First calculation
Prepare the files required for the calculation.
・ "POSCAR": After optimizing the structure, or creating the structure you want to calculate, prepare a file with a fixed structure.
・ "KPOINTS": Prepare a file with enough k points for the accuracy of the calculation to converge.
・ "POTCAR": Prepare your favorite pseudopotential file.
・ "INCAR": Please refer to "https://qiita.com/youkihashimoto3110/items/de92172e0b5e9f3872d3" for the meaning of each tag used in the calculation setting file INCAR.
The first calculation is a new turn of the charge optimization (Electronic Relaxation) calculation. Comment out all the structural optimization parameters. Required tags are
ISTART = 0
PREC = Accurate
EDIFF = 1E-5
LREAL = Auto
ALGO = VeryFast
is. In the case of metal
ISMEAR = 1 #Or 2,-5
SIGMA = 0.5
In the case of insulators and semiconductors
ISMEAR = 0 #Or-5
SIGMA = 0.01
Use to adjust the SIGMA value so that the outcar entropy term is less than 1 meV per atom.
___ Do not use `ISMEAR> 0``` for insulators / semiconductors! ___ Also, using `ISMEAR = -5``` (tetrahedron method) improves accuracy, but is not suitable for most structures. (???)
Also, use "effective measures in case of non-convergence" as appropriate.
After completing the settings, execute the first calculation.
step2. Second calculation
When the calculation is completed, we will start the second calculation. In the second calculation, the "INCAR" file and "KPOINTS" are rewritten in order to perform charge optimization in the continuous job with the respecified k-mesh.
The changes from the first calculation file
ISTART = 1
ICHARG = 11
ISYM = 0
LORBIT = 11
is. Let's recalculate the others in the same way. The meaning of each is described in the template.
In the worst case, the "KPOINTS" file can be left as it is, but I think it is easier to see the phenomenon if you cut it with a mesh on the symmetric Bruyan zone.
You can specify it manually, but I don't feel like motivating something that is just 2D and in 3D, so I borrow the power of pymatgen.
Download and install Anaconda from "https://www.anaconda.com/products/individual" and install pymatgen with the conda command `conda install --channel conda-forge pymatgen`.
Create a python file and do the following:
generate_kpoints.py
from pymatgen.io.vasp.inputs import Kpoints
from pymatgen.core import Structure
from pymatgen.symmetry.bandstructure import HighSymmKpath
struct = Structure.from_file("POSCAR")
kpath = HighSymmKpath(struct)
kpts = Kpoints.automatic_linemode(divisions=40,ibz=kpath)
kpts.write_file("KPOINTS_nsc")
The POSCAR is read and a K mesh "KPOINT_nsc" that passes through k points with good symmetry is generated. This "KPOINT_nsc" will be replaced with a new "KPOINTS".
Start the job when you are ready.
2. Read EIGENVAL
Open the "EIGENVAL" output in the second calculation with any text editor.
The contents of EIGENVAL are as shown in the figure below.
I'm sorry it's hard to see.
Also, the example shown above does not take spin into account in the calculation. In the case of calculation considering spins, it is also output for each spin.
3. Output the result with python
That is, "There are 6 lines of fixed character strings, and there are 3 coordinates in reciprocal space and ??. The numerical value of ?? does not change even if the coordinates of the K point are different. Below that, the band number, energy, and occupancy rate are the next K. It continues to the coordinates. And there are three coordinates in the next reciprocal space and ??. The numerical value of ?? does not change even if the coordinates of the K point are different. Below that, the band number, energy, occupancy rate ... " You can see that it is a structure. In addition, it is necessary to skip blank line breaks and blanks.
If you program based on these, it will be as follows.
band.py
import csv
import matplotlib.pyplot as plt
import numpy as np
import random
eigen_name = "EIGENVAL" #EIGENVAL filename for VASP output
kpoint_name = "KPOINTS" #VASP KPOINTS filename
save_csv_name = "eigenvalue.csv" #Eigenvalue storage file
save_band = "band.csv" #Save the band diagram
#File name text_Returns a list with name
def file_to_list(text_name):
with open(text_name, newline='') as texts:
read = csv.reader(texts, delimiter=' ', skipinitialspace=True)
lists = [i for i in read]
return lists
#Save file name csv_Save the list with save and list lists
def list_to_csv(csv_name,lists):
with open(csv_name, mode='w', newline='') as csvname:
write = csv.writer(csvname)
write.writerows(lists)
#List all_Specify the start point index of list, cut out the group and store it for each array
def list_groups(all_list,index):
groups_all = []
while True:
index = index + 1
temp = []
count = 0
temp.append(all_list[index])
while True:
index = index + 1
count = count + 1
temp.append(all_list[index])
if len(all_list)==index+1:
break
if len(temp[count]) != len(all_list[index+1]):
groups_all.append(temp)
index = index + 1 #To skip blank line breaks
break
if len(all_list)==index+1:
break
return groups_all
#When kpoint is automatically mesh cut
def k_auto_band(groups_all,k_file):
spin_up = []
spin_do = []
for i in range(len(groups_all)):
temp_up = []
temp_do = []
for j in range(len(groups_all[i])):
if j == 0: #Put a K point in the first column
k_vec = [float(f) for f in groups_all[i][j][0:3]]
label = k_file_search(k_file,k_vec)
for k in groups_all[i][j]:
label = label +";"+str(k)
temp_up.append(label)
if len(groups_all[i][j])==4:
temp_do.append(label)
else:
temp_up.append(groups_all[i][j][1])
if len(groups_all[i][j])==4:
temp_do.append(groups_all[i][j][2])
spin_up.append(temp_up)
if len(groups_all[i][j])==4:
spin_do.append(temp_do)
return spin_up, spin_do
def k_file_search(k_file,lists):
if k_file[1][0]=="0":
return ""
else:
temp = []
for i in range(len(k_file)):
if len(k_file[i]) >= 4:
for j in range(3):
temp.append(float(k_file[i][j]))
if temp == lists:
return k_file[i][4]
else:
return ""
#plot. Excel can only display up to 256 y series.
def graf_show(lists):
fig = plt.figure()
list_t = np.array(lists).T
x = range(len(list_t[0]))
for i in range(len(lists[0])):
plt.plot(x,list_t[i])
fig.savefig("img.png ")
#I can't draw even with the ↑ method, so I randomly select 200 series.
def random_choice(lists):
list_t = np.array(lists).T
random_c = random.sample(range(len(list_t[1:])), 200)
new_list = []
new_list.append(list_t[0])
for i in random_c:
new_list.append(list_t[i])
return np.array(new_list).T.tolist()
#Main
eigen_list = file_to_list(eigen_name)
kpoint_list = file_to_list(kpoint_name)
list_to_csv(save_csv_name,eigen_list)
all_eigen_list = list_groups(eigen_list,6)
spin_up, spin_do = k_auto_band(all_eigen_list,kpoint_list)
list_to_csv(save_band,spin_up)
#graf_show(spin_up) #Plot with matplotlib (fault)
#random_c = random_choice(spin_up) #When there are more than 200 bands, randomly select 200 from the bands
#list_to_csv("radom.csv",random_c) #Save the band diagram randomly selected in ↑.
It's a terrible program, but I don't think it's hard to read. After that, you can plot the band diagram by plotting with Excel etc.
Also, when I tried to plot about 800 bands using the matplotlib module, the processing stopped, so I am looking for a lighter method.
Since Excel can only plot up to 250 series, we have prepared a method to randomly extract 200 bands and save them in csv format. For systems with more than 200 bands, uncomment `random_c = random_choice (spin_up)` and `list_to_csv (" radom.csv ", random_c)`.
An example using the program is shown above.
Recommended Posts