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就建立好了表面模型。

练习二：切割复杂化合物表面，LiFePO₄-(001)-Fe-terminal

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

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。两种情况体相中PO₄都得以保留。

转载请注明来源，欢迎对文章中的引用来源进行考证，欢迎指出任何有错误或不够清晰的表达。

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

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

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

练习二：切割复杂化合物表面，LiFePO₄-(001)-Fe-terminal