pymatgen自动切割表面-高通量计算初探

  1. 练习一:切割金属表面Au(111)面
  2. 练习二:切割复杂化合物表面,LiFePO4-(001)-Fe-terminal

本文章为原创,版权归作者刘锦程所有,文章转载请先取得作者的同意,非常欢迎转发文章链接!严禁以任何方式挪用本文内容,用于以盈利为目的各种活动。

相关内容

pymatgen API的使用-高通量计算初探

pymatgen自动切割表面-高通量计算初探

pymatgen自动生成表面吸附模型-高通量计算初探

pymatgen画能带图方法-高通量计算初探

进行表面slab模型的搭建需要在晶体conventional cell的基础上,切相应miller指数的晶面,手动切面需要用到Materials studio、p4vasp、VESTA等程序,这些程序难以实现大规模批处理。如果需要批处理表面模型,实现自动化表面切割就很重要了。

切割表面要考虑以下几个问题:

  1. 表面slab的厚度,
  2. 真空层的厚度,
  3. 表面supercell的大小,
  4. 表面切断的位置(暴露的表面原子),
  5. 表面的重构。

使用到的模块是pymatgen.core.surface,该模块非常智能,这些内容全都可以通过pymatgen程序实现。

练习一:切割金属表面Au(111)面

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
from pymatgen.analysis.adsorption import *
from pymatgen.core.surface import Slab, SlabGenerator, generate_all_slabs, Structure, Lattice, ReconstructionGenerator
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.core.structure import Structure
from pymatgen.ext.matproj import MPRester
from matplotlib import pyplot as plt
from pymatgen.io.vasp.inputs import Poscar
mpr = MPRester('yourPrivateKey')

mp_id = "mp-81"
struct = mpr.get_structure_by_material_id(mp_id)
struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()
slab = SlabGenerator(struct, miller_index=[1, 1, 1], min_slab_size=8.0,\
min_vacuum_size=15.0, center_slab=True)

for n, slabs in enumerate(slab.get_slabs()):
slabs_bak = slabs.copy()
slabs.make_supercell([[2,0,0],
[0,2,0],
[0,0,1]])
fig = plt.figure()
ax = fig.add_subplot(111)
plot_slab(slabs, ax, adsorption_sites=True)
plt.savefig(str(n) + '-Au-111.png', img_format='png')
open('POSCAR' + mp_id + '-' + str(n), 'w').write(str(Poscar(slabs)))

mpr = MPRester('yourPrivateKey')位置输入你的material project的私钥,就可以用MPRester导入任意

material project上有的结构。

切表面之前要用 SpacegroupAnalyzer(struct).get_conventional_standard_structure() 确定已经将晶胞转化成了惯用晶胞。

然后用SlabGenerator产生所有可能截断的表面,先读入一个structure类的结构,确定晶面指数:miller_index=[1, 1, 1],确定slab层的最小厚度:min_slab_size=8.0 (unit: Angstrom),确定真空层的厚度:min_vacuum_size=15.0 (unit: Angstrom)。

然后在循环所有可能的截断的表面,如果是Au(111)的话只有一个可能的暴露表面,如果是化合物,可能有多种表面。

strucutre类里的make_supercell函数可以方便的给表面晶胞进行扩包,只需要读入一个 3*3 的数组:

1
2
3
slabs.make_supercell([[2,0,0],
[0,2,0],
[0,0,1]])

对角元表示a,b,c 每个晶格矢量方向上扩大的倍数,此时保持非对角元为0,如果想要对晶格矢量进行旋转可以对非对角元进行改动比如:

1
2
3
slabs.make_supercell([[1,1,0],
[2,1,0],
[0,0,1]])

代表新生成的晶格矢量:A = a + b; B = 2a - b; C = c

这两个建出来的模型如下($\times$代表可能的吸附位点):

这是 $\operatorname{Au}(111)-(2 \times 2)$ 表面:

这是$\operatorname{Au}(111)-(\sqrt{3} \times \sqrt{3}) R 30^{\circ} $表面:

把结构文件以POSCAR的格式导出:

1
open('POSCAR' + mp_id + '-' + str(n), 'w').write(str(Poscar(slabs)))

在VESTA里显示的$\operatorname{Au}(111)-(\sqrt{3} \times \sqrt{3}) R 30^{\circ} $表面:

如果想要从cif文件里导入结构,用IStructure.from_file也可以导入structure类的结构。例如:

1
2
from pymatgen.core.structure import Structure, IStructure
struct = IStructure.from_file("LaSiRu_mp-4579_primitive.cif")

到这里,我们不借助于MS、p4vasp就建立好了表面模型。

练习二:切割复杂化合物表面,LiFePO4-(001)-Fe-terminal

对于复杂化合物,切表面时候往往不只一种表面暴露情况,比如LiFePO4,就可以是Li暴露,Fe暴露,P暴露,O暴露的,也可能表面有一些分子钝化表面,比如OH,H集团来饱和表面的原子。

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
from pymatgen.analysis.adsorption import *
from pymatgen.core.surface import Slab, SlabGenerator, generate_all_slabs, Structure, Lattice, ReconstructionGenerator
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.core.structure import Structure
from pymatgen.ext.matproj import MPRester
from matplotlib import pyplot as plt
from pymatgen.io.vasp.inputs import Poscar
mpr = MPRester('yourPrivateKey')


mp_id = "mp-19017"
struct = mpr.get_structure_by_material_id(mp_id)
struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()
slab = SlabGenerator(struct, miller_index=[0, 0, 1], min_slab_size=8.0,\
min_vacuum_size=15.0, center_slab=True)


for n, slabs in enumerate(slab.get_slabs(bonds = {('P', 'O'): 3})):
slabs_bak = slabs.copy()
slabs.make_supercell([[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
fig = plt.figure()
ax = fig.add_subplot(111)
plot_slab(slabs, ax, adsorption_sites=True)
plt.savefig(str(n) + '-LiFePO4-001.png', img_format='png')
open('POSCAR' + mp_id + '-' + str(n), 'w').write(str(Poscar(slabs)))

其中在get_slabs()中,添加了bonds = {('P', 'O'): 3},对于P-O键在3 A以内的键,切表面的时候不允许截断,这在多酸化合物中是个准则。最后得到了两种表面:

第一种是上表面暴露O,下表面暴露Fe;第二种是上下表面都暴露Fe。两种情况体相中PO4都得以保留。


欢迎在评论区讨论问题、补充内容、指出错误。请勿发送占楼、沙发、点赞等无意义的回复。
目录