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

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

相关内容

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

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

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

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

本文参考文献: A high-throughput framework for determining adsorption energies on solid surfaces

https://www.nature.com/articles/s41524-017-0017-z

Gerbrand Ceder, Kristin A. Persson等人领导团队多年来一直致力于高通量计算方法,和materials project数据库的发展,https://materialsproject.org

2017年Kristin A. Persson 和 Joseph H. Montoya发展了表面吸附模型的高通量计算模块集成到pymatgen程序包里,使得表面吸附分子的数据库构建成为可能。这篇文章就介绍一下用pymatgen实现表面高通量计算的方法。

其中有几个关键点:

  1. 通过API自动得到数据库里的晶胞结构。(以后会另外写一篇文章专门介绍Materials project的API)
  2. 自动切割表面,要考虑的问题非常多:比如,表面的厚度,最表层暴露的原子,晶面指数,真空层厚度,slab模型的位置等等。
  3. 自动识别不同的吸附位点,比如:top,bridge,hollow。识别要基于表面的对称性,这就需要程序有可靠的对称性识别算法,对于slab模型的整体对称性和不等价吸附位点能准确判断。
  4. 自动添加吸附原子/分子,并保存好POSCAR文件。
  5. 批量提交任务。

练习一:手动建立一个Ni晶胞,切(111)面,并自动识别不等价的吸附位点。

首先加载我们本文有可能用到的python包:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from pymatgen import Structure, Lattice, MPRester, Molecule
from pymatgen.analysis.adsorption import *
from pymatgen.core.surface import generate_all_slabs
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from matplotlib import pyplot as plt
from pymatgen.ext.matproj import MPRester
from pymatgen.io.vasp.inputs import Poscar

fcc_ni = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3.5), ["Ni", "Ni"], [[0, 0, 0], [0.5, 0.5, 0.5]])
slabs = generate_all_slabs(fcc_ni, max_index=1, min_slab_size=8.0, min_vacuum_size=15.0)
ni_111 = [slab for slab in slabs if slab.miller_index==(1,1,1)][0]

asf_ni_111 = AdsorbateSiteFinder(ni_111)
ads_sites = asf_ni_111.find_adsorption_sites()
# print(ads_sites)
assert len(ads_sites) == 4


fig = plt.figure()
ax = fig.add_subplot(111)
plot_slab(ni_111, ax, adsorption_sites=True)
plt.savefig('111.png', img_format='png')

运行以上代码:

Structure.from_spacegroup建立一个类,通过空间群和坐标参数生成一个Structure对象。

然后对这个晶胞结构进行自动切面,max_index=1得到最大的不等价表面,比如(1, 1, 1) (1, 1, 0) (1, 0, 0),min_slab_size=8.0, min_vacuum_size=15.0控制晶胞的厚度和真空层。

然后用AdsorbateSiteFinder找到有可能的吸附位点。把print(ads_sites)前的注释拿掉,就能打印出所有的吸附位置。

1
2
3
4
{'ontop': [array([-9.57921524e-16,  2.14330352e+00,  9.57772228e+00])], \
'bridge': [array([-6.71997040e-16, 1.69152089e-15, 9.57772228e+00])], \
'hollow': [array([1.23743687, 1.42886902, 9.57772228]), array([-4.66596435e-16, 7.14434508e-01, 9.57772228e+00])], \
'all': [array([-9.57921524e-16, 2.14330352e+00, 9.57772228e+00]), array([1.8561553 , 1.07165176, 9.57772228]), array([1.23743687, 1.42886902, 9.57772228]), array([3.57702105e-16, 7.14434508e-01, 9.57772228e+00])]}

用matplotlib就能直接画出所有的可能吸附位点了:

练习二:添加吸附分子或原子。

1
2
3
4
5
6
fig = plt.figure()
ax = fig.add_subplot(111)
adsorbate = Molecule("H", [[0, 0, 0]])
ads_structs = asf_ni_111.generate_adsorption_structures(adsorbate, repeat=[1, 1, 1])
plot_slab(ads_structs[0], ax, adsorption_sites=False, decay=0.09)
plt.savefig('111-H.png', img_format='png')

想要吸附原子或者分子,点定义一个变量:

做HER反应吸附H原子adsorbate = Molecule("H", [[0, 0, 0]])

做ORR,OER反应吸附O,OH,OOH

1
2
3
OH = Molecule("OH", [[0, 0, 0], [-0.793, 0.384, 0.422]])
O = Molecule("O", [[0, 0, 0]])
OOH = Molecule("OOH", [[0, 0, 0], [-1.067, -0.403, 0.796], [-0.696, -0.272, 1.706]])

用generate_adsorption_structures函数产生吸附结构,用matplotlib画图:

练习三:使用API,从Materials project得到晶体结构。

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
# Import statements
from pymatgen import Structure, Lattice, MPRester, Molecule
from pymatgen.analysis.adsorption import *
from pymatgen.core.surface import generate_all_slabs
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from matplotlib import pyplot as plt
from pymatgen.ext.matproj import MPRester
from pymatgen.io.vasp.inputs import Poscar
# Note that you must provide your own API Key, which can
# be accessed via the Dashboard at materialsproject.org
mpr = MPRester('yourPrivateKey')

fig = plt.figure()
axes = [fig.add_subplot(2, 3, i) for i in range(1, 7)]
mats = {"mp-23":(1, 0, 0), # FCC Ni
"mp-2":(1, 1, 0), # FCC Au
"mp-13":(1, 1, 0), # BCC Fe
"mp-33":(0, 0, 1), # HCP Ru
"mp-5229":(1, 0, 0), # Cubic SrTiO3
"mp-2133":(0, 0, 1)} # Wurtzite ZnO
for n, (mp_id, m_index) in enumerate(mats.items()):
struct = mpr.get_structure_by_material_id(mp_id)
struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()
slabs = generate_all_slabs(struct, 1, 5.0, 2.0, center_slab=True)
slab_dict = {slab.miller_index:slab for slab in slabs}
asf = AdsorbateSiteFinder.from_bulk_and_miller(struct, m_index)
plot_slab(asf.slab, axes[n])
ads_sites = asf.find_adsorption_sites()
sop = get_rot(asf.slab)
ads_sites = [sop.operate(ads_site)[:2].tolist()
for ads_site in ads_sites]
axes[n].plot(*zip(*ads_sites), color=’k’, marker=’x’,markersize=10, mew=1,\
linestyle=’’, zorder=10000)
mi_string = "".join([str(i) for i in m_index])
axes[n].set_title("{}({})".format(struct.composition.reduced_formula, mi_string))
axes[n].set_xticks([])
axes[n].set_yticks([])

这样,就从Materials Project数据库里读取到了我们想要的晶体结构,mp-23 mp-5229等都是晶体数据库里的编号,通过构建一个字典mats指定我们想要的切面,然后循环所有我们指定的表面,做模型识别吸附位点添加原子即可。

想要导出POSCAR文件用open().write()即可,这是我批量到处POSCAR的方式:

1
open('POSCAR' + mp_id + '-' + "{}({})".format(struct.composition.reduced_formula, mi_string) + '-' + str(count), 'w').write(str(Poscar(reorient_z(slabmodel))))

实战演练,金属和合金表面HER反应高通量计算

笔者读了Kristin A. Persson文章之后一直在思考一个问题,既然这个功能2017年已经开发完善了,但是这都2019年了,迟迟不见表面吸附数据库出现。经过笔者实战尝试时候终于懂了为什么。如果只考虑Ni和Pt两种元素,他们之间可能形成 Ni,Pt,NiPt,Ni3Pt,NiPt3这5中组成状况,每种晶体就算只考虑Miller最高到(1, 1, 1),这样每种晶体也有三种表面,然后每种表面有4到10几种不同的吸附位点,这样算下来完成Ni/Pt二元合金的HER计算,就生成了5 3 10 = 150个结构,如果考虑到整个周期表中的金属和二元合金的话,组合情况至少要多100倍,这样就是上万个结构。每个结构如果要在20核的单节点上算5小时,算下来就是上百万的核时。壕如 Gerbrand Ceder, Kristin A. Persson也不敢轻易的去尝试表面数据库的搭建,所以笔者也只是满足于自动生成了上千种POSCAR文件结构就浅尝辄止了。

如果真要做高通量,自己去手动输入mp-23 mp-5229这样的编号显然是比较笨的方法,用元素组成,然后自动去数据库里下载结构文件时更优方法:

1
struct = mpr.get_structures('Ni-Pt-', final=True)

到这里,我们已经把高通量计算最难的结构批量搭建完成了,后面就是自动生成POTCAR,KPOINTS,INCAR文件了,这些东西用VASPKIT或者pymatgen,或者自己写个小循环就能轻易搞定的了。流程如下:

VASP批量计算完成以后,还有重要的内容就是后处理了,首先就是要把后处理程序化,这一点pymatgen和VASPKIT都做得非常好,把这些数据批处理以后再导入到materials project数据库里或者自己搭建一个新得表面模型数据库。然后,怎么在浩如烟海的计算结果里日提取出规律性的东西,这就是另外一回事了。


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