如何优雅地画态密度和能带结构?

基于Pymatgen和plotly的DOS和Band作图

Requirements

  • Pymatgen
  • plotly
  • matplotlib

DOS作图

这篇paper将通过两个简单的例子,展示如何通过pymatgen画图

  • pymatgen: (Python Materials Genomics:Python材料基因模块),是一个强大的,适用于材料分析研究的开源的库。
  • plotly: 是一个可交互的画图工具

首先导入必要的package

1
2
3
4
5
6
7
8
9
10
11
%matplotlib inline
from pymatgen import Element
from pymatgen.io.vasp import Vasprun # 读取vasp作业(vasprun.xml)
from pymatgen.io.vasp import Vasprun, BSVasprun
from pymatgen.electronic_structure.plotter import DosPlotter
from pymatgen.electronic_structure.core import Spin, OrbitalType
import plotly.plotly as pltly
import plotly.tools as tls
import plotly.graph_objs as go

简单体系 $Al_1$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
dosrun = Vasprun('dos_and_band_jobs/dos/Al/vasprun.xml')
spd_dos = dosrun.complete_dos.get_spd_dos()
trace_tdos = go.Scatter(
x=dosrun.tdos.densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode='lines',
name='Total DOS',
line=go.Line(color='#444444'),
fill='tozeroy')
# 3s contribution to the total DOS
trace_3s = go.Scatter(
x=spd_dos[OrbitalType.s].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="3s",
line=go.Line(color="red"))
# 3p contribution to the total DOS
trace_3p = go.Scatter(
x=spd_dos[OrbitalType.p].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="3p",
line=go.Line(color="green"))
dosdata = go.Data([trace_tdos, trace_3s, trace_3p])
dosxaxis = go.XAxis(
title="Density of states",
showgrid=True,
showline=True,
range=[.01, 3],
mirror="ticks",
ticks="inside",
linewidth=2,
tickwidth=2
)
dosyaxis = go.YAxis(
title="$E - E_f \quad / \quad \\text{eV}$",
showgrid=True,
showline=True,
ticks="inside",
mirror='ticks',
linewidth=2,
tickwidth=2,
zerolinewidth=2
)
doslayout = go.Layout(
title="Density of states of Silicon",
xaxis=dosxaxis,
yaxis=dosyaxis
)
dosfig = go.Figure(data=dosdata, layout=doslayout)
plot_url = pltly.plot(dosfig, filename="DOS_Al", auto_open=False)
tls.embed(plot_url)

Eaxmple Example (稍微复杂的体系) $H{18}C{70}N_{67}F$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
dosrun = Vasprun('dos_and_band_jobs/dos/HCNF/vasprun.xml')
spd_dos = dosrun.complete_dos.get_spd_dos()
trace_tdos = go.Scatter(
x=dosrun.tdos.densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode='lines',
name='Total DOS',
line=go.Line(color='#444444'),
fill='tozeroy')
# 3s contribution to the total DOS
trace_3s = go.Scatter(
x=spd_dos[OrbitalType.s].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="3s",
line=go.Line(color="red"))
# 3p contribution to the total DOS
trace_3p = go.Scatter(
x=spd_dos[OrbitalType.p].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="3p",
line=go.Line(color="green"))
dosdata = go.Data([trace_tdos, trace_3s, trace_3p])
dosxaxis = go.XAxis(
title="Density of states",
showgrid=True,
showline=True,
range=[.01, 50],
mirror="ticks",
ticks="inside",
linewidth=2,
tickwidth=2
)
dosyaxis = go.YAxis(
title="$E - E_f \quad / \quad \\text{eV}$",
showgrid=True,
showline=True,
ticks="inside",
mirror='ticks',
linewidth=2,
tickwidth=2,
zerolinewidth=2
)
doslayout = go.Layout(
title="Density of states",
xaxis=dosxaxis,
yaxis=dosyaxis
)
dosfig = go.Figure(data=dosdata, layout=doslayout)
plot_url = pltly.plot(dosfig, filename="$H_{18}C_{70}N_{67}F_1$", auto_open=False)
tls.embed(plot_url)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
# DOS 投影到不同元素
cdos = dosrun.complete_dos
element_dos = cdos.get_element_dos()
trace_tdos = go.Scatter(
x=dosrun.tdos.densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode='lines',
name='Total DOS',
line=go.Line(color='#444444'),
fill='tozeroy')
# H contribution to the total DOS
trace_H = go.Scatter(
x=element_dos[Element('H')].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="H",
line=go.Line())
# C contribution to the total DOS
trace_C = go.Scatter(
x=element_dos[Element('C')].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="C",
line=go.Line())
# N contribution to the total DOS
trace_N = go.Scatter(
x=element_dos[Element('N')].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="N",
line=go.Line())
# F contribution to the total DOS
trace_F = go.Scatter(
x=element_dos[Element('F')].densities[Spin.up],
y=dosrun.tdos.energies - dosrun.efermi,
mode="lines",
name="F",
line=go.Line())
dosdata = go.Data([trace_tdos, trace_H, trace_C, trace_N, trace_F])
dosxaxis = go.XAxis(
title="Density of states",
showgrid=True,
showline=True,
range=[.01, 50],
mirror="ticks",
ticks="inside",
linewidth=2,
tickwidth=2
)
dosyaxis = go.YAxis(
title="$E - E_f \quad / \quad \\text{eV}$",
showgrid=True,
showline=True,
ticks="inside",
mirror='ticks',
linewidth=2,
tickwidth=2,
zerolinewidth=2
)
doslayout = go.Layout(
title="Density of states",
xaxis=dosxaxis,
yaxis=dosyaxis
)
dosfig = go.Figure(data=dosdata, layout=doslayout)
plot_url = pltly.plot(dosfig, filename="$H_{18}C_{70}N_{67}F_1$", auto_open=False)
tls.embed(plot_url)

用pymatgen自带的plotter比较

1
2
3
4
5
6
from pymatgen.electronic_structure.plotter import DosPlotter
vasprun = Vasprun('dos_and_band_jobs/dos/HCNF/vasprun.xml')
tdos = vasprun.tdos
plotter = DosPlotter()
plotter.add_dos("Total DOS", tdos)
plotter.show(xlim=[-4, 4], ylim=[-80, 80])

1
2
3
4
5
cdos = vasprun.complete_dos
element_dos = cdos.get_element_dos()
plotter = DosPlotter()
plotter.add_dos_dict(element_dos)
plotter.show(xlim=[-5, 5], ylim=[-40, 40])


能带结构作图

读取vasp作业下的vasprun.xml文件

Eaxmple $Al_1$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
run = BSVasprun("dos_and_band_jobs/band/vasprun.xml", parse_projected_eigen = True) # 读取vasprun.xml
bands = run.get_band_structure(kpoints_filename="band/KPOINTS", line_mode=True, efermi=dosrun.efermi)
emin = 1e100
emax = -1e100
for spin in bands.bands.keys():
for band in range(bands.nb_bands):
emin = min(emin, min(bands.bands[spin][band]))
emax = max(emax, max(bands.bands[spin][band]))
emin = emin - bands.efermi - 1
emax = emax - bands.efermi + 1
kptslist = [k for k in range(len(bands.kpoints))]
bandTraces = list()
for band in range(bands.nb_bands):
bandTraces.append(
go.Scatter(
x=kptslist,
y=[e - bands.efermi for e in bands.bands[Spin.up][band]],
mode="lines",
line=go.Line(color="#666666"),
showlegend=False
)
)
labels = [r"$L$", r"$\Gamma$", r"$X$", r"$U,K$", r"$\Gamma$"]
step = len(bands.kpoints) / (len(labels) - 1)
# vertical lines
vlines = list()
for i, label in enumerate(labels):
vlines.append(
go.Scatter(
x=[i * step, i * step],
y=[emin, emax],
mode="lines",
line=go.Line(color="#111111", width=1),
showlegend=False
)
)
# Labels of highsymetry k-points are added as Annotation object
annotations = list()
for i, label in enumerate(labels):
annotations.append(
go.Annotation(
x=i * step, y=emin,
xref="x1", yref="y1",
text=label,
xanchor="center", yanchor="top",
showarrow=False
)
)
bandxaxis = go.XAxis(
title="k-points",
range=[0, len(bands.kpoints)],
showgrid=True,
showline=True,
ticks="",
showticklabels=False,
mirror=True,
linewidth=2
)
bandyaxis = go.YAxis(
title="$E - E_f \quad / \quad \\text{eV}$", # 费米能级对其
range=[emin, emax],
showgrid=True,
showline=True,
zeroline=True,
mirror="ticks",
ticks="inside",
linewidth=2,
tickwidth=2,
zerolinewidth=2
)
bandlayout = go.Layout(
title="Bands diagram of Al",
xaxis=bandxaxis,
yaxis=bandyaxis,
annotations=go.Annotations(annotations)
)
bandfig = go.Figure(data=bandTraces + vlines, layout=bandlayout)
plot_url = pltly.plot(bandfig, filename="Bands_Al", auto_open=False)
tls.embed(plot_url)

比较 [PyMatgen]

1
2
3
4
5
6
from pymatgen.io.vasp import Vasprun, BSVasprun
from pymatgen.electronic_structure.plotter import BSPlotter
v = BSVasprun("dos_and_band_jobs/band/vasprun.xml")
bs = v.get_band_structure(kpoints_filename="band/KPOINTS",line_mode=True)
plt = BSPlotter(bs) plt.get_plot(vbm_cbm_marker=True,ylim=(-3,3)) plt.show()

总结

  • plotly 可交互

  • Pymatgen自带,可应用于paper展示

参考

人艰不拆,生活不易