脚本分享(二)——在不同表面之间迁移快速迁移反应路径

type

Post

status

Published

date

Dec 18, 2022

slug

summary

一个可以快速地将反应路径迁移到另一个表面的脚本

category

技术分享

icon

password

URL

Property

Sep 16, 2024 11:31 PM

其实看很多催化计算的文章，甚至是大部分有理论计算的文章中间，作者都会去对比几种不同材料的构效关系。特别的，对于元素有掺杂的情况，其实表面反应路径的变化不大，但是做新的路径需要的步骤却很繁琐：

导出上一个路径优化完的结构——>构建一个新的基底模型——>把这个表面分子模型移到新的表面上——>生成新的POSCAR及其他计算文件

而且这样做还会存在手动搭建的结构的固有问题：初始结构不合理，以及如果需要计算过渡态还需要重新调整原子顺序。

我之前正好在做材料的的掺杂对表面反应关键步骤的影响的计算，于是写了下面这个脚本。这个脚本的原理是分别指定两个表面活性位，然后将某一高度上的原子转移到另一个表面。这样做的好处很多：省去手动建模的繁琐；直接拿上个路径过渡态的IS和FS进行建模，可以不用调整顺序直接插点；获得一个比较合理的初始结构，节约机时等。

以下是代码：


#CONTCAR:A completed reaction pathway structure (substrate + surface molecule).
#BASECAR:New catalyst substrate model files to be constructed.
#POSCAR: A new model for catalyst surface calculations.
import os
import shutil
import pymatgen.core as mg
import pandas as pd
from pymatgen.io.vasp import Poscar

Hight_max = 0.43 #Divided height of the model surface in Cartesian coordinate
site_old_Species = 'Cu'#Location of the active site of the original structure
site_old_num = 1 #Active site number of the original structure
site_new_Species = 'Cu'#Location of the active site of the new structure
site_new_num = 1 #Active site number of the new structure
Species_df = pd.DataFrame(columns=['H','C','N','O','Cu'])#Order of possible elements
Species_element = Species_df.columns.values.tolist()

base = Poscar.from_file("BASE")
base = base.structure
Species_2 = 0
base_num = 0
base_species = 'H'

for i in range(base.num_sites):
    if base_species != str(base.species[i]):
        base_species = str(base.species[i])
        Species_df.loc['Start_index',base_species] = i
    if site_new_Species == str(base.species[i]):
        Species_2 = Species_2 + 1
        if site_new_num == Species_2:
            site_B_index = i
            break
           
for i in range(Species_df.shape[1]):
    if Species_df.iloc[0][i] != Species_df.iloc[0][i]:
        Species_df.iloc[0][i] = Species_df.iloc[0][i+1]
for i in range(Species_df.shape[1]):
    if i != Species_df.shape[1]-1:
        Species_df.iloc[0][i] = Species_df.iloc[0][i+1]
    else:
        Species_df.iloc[0][i] = base.num_sites
        
Species_df_s = Species_df.copy(deep=True)

for root, dirs, files in os.walk(os.getcwd()):
        dir_list = dirs
        break
dir_list.remove('.ipynb_checkpoints')

for dirs in dir_list:
    base = Poscar.from_file("BASE")
    base = base.structure
    Species_df = Species_df_s.copy(deep=True)    
    #Reading structure files
    old_str = Poscar.from_file(dirs+'/'+"CONTCAR")
    old_str = old_str.structure
    #Determine the displacement vector between the two active sites
    Species_1 = 0
    for i in range(old_str.num_sites):
        if site_old_Species == str(old_str.species[i]):
            Species_1 = Species_1 + 1
            if site_old_num == Species_1:
                site_A_index = i
                break
    Displacement_vector = base.frac_coords[site_B_index] - old_str.frac_coords[site_A_index]
    num = 0
    for i in range(old_str.num_sites):
        coord_1 = old_str.frac_coords[i-num]
        if coord_1[2] < Hight_max :
            old_str.remove_sites([i-num])
            num = num + 1
    #Removal of the base atoms from the original structure file
    num = 0
    for i in range(old_str.num_sites):
        coord_1 = old_str.frac_coords[i-num]
        if coord_1[2] < Hight_max :
            old_str.remove_sites([i-num])
            num = num + 1
    #Adding molecules to a new substrate file
    for i in range(old_str.num_sites):
        insert_index = Species_df.loc['Start_index',str(old_str.species[i])]
        old_spec = old_str.species[i]
        coord_new = old_str.frac_coords[i] + Displacement_vector
        base.insert(insert_index,old_spec,coord_new)
        for m in range(Species_element.index(str(old_str.species[i])),Species_df.shape[1]):
            Species_df.iloc[0,m] = Species_df.iloc[0,m] + 1
    #Exporting a new structure file
    base.to(dirs+'/'+"POSCAR")
    shutil.copyfile('INCAR',dirs+'/'+"INCAR")
    shutil.copyfile('KPOINTS',dirs+'/'+"KPOINTS")

使用之前需要：1. 准备文件 2.修改脚本参数

准备文件

使用的时候需要拿出上个路径的CONTCAR，放在单独的路径文件夹里：

加上你需要的INCAR，KPOINTS：这两个可选，如果不需要记得把下面这两行注释掉


    shutil.copyfile('INCAR',dirs+'/'+"INCAR")
    shutil.copyfile('KPOINTS',dirs+'/'+"KPOINTS")

BASE:新的表面

mol_tran.ipynb：脚本文件

这是可以运行脚本的文件结构：


|____INCAR
|____KPOINTS
|____BASE
|____mol_tran.ipynb
|____09ts02_tran1
| |____CONTCAR
|____04Cu_O
| |____CONTCAR
|____03Cu_O_PO
| |____CONTCAR
|____00Cu
| |____CONTCAR
|____01Cu_O2
| |____CONTCAR
|____05Cu_O_pro
| |____CONTCAR
|____08ts01_tran2
| |____CONTCAR
|____07ts01_tran1
| |____CONTCAR
|____06Cu_PO
| |____CONTCAR
|____02Cu_O2_pro
| |____CONTCAR

修改脚本参数

这些参数需要修改：


Hight_max = 0.43 #Divided height of the model surface in Cartesian coordinate
site_old_Species = 'Cu'#Location of the active site of the original structure
site_old_num = 1 #Active site number of the original structure
site_new_Species = 'Cu'#Location of the active site of the new structure
site_new_num = 1 #Active site number of the new structure
Species_df = pd.DataFrame(columns=['H','C','N','O','Cu'])#Order of possible elements这些参数需要修改：

这些参数需要修改：

Hight_max：这个高度以上的原子会被转移到新的表面

site_old_Species：上个路径的活性位点的元素，原子会保持相对跟这个原子的相对位置，转移到另一个表面

site_old_num：上个路径活性位点序号

site_new_Species：新路径的活性位点的元素，原子会保持相对跟上个表面活性位相对位置，转移到这个原子

site_new_num：新路径活性位点序号

Species_df：请填写整个路径所有可能出现的元素，并依序排列。注意BASE的元素顺序不能打乱，例如BASE中是C，N材料，反应路径需要吸附氧气，则可以写成：O，C，N或者C，O，N或者C，N，O，而写成N，C，O等是错误的

以及修改下面这一行：


base_species = 'H'

把这里的H改成Species_df里面第一个元素

之后运行脚本即可：


|____INCAR
|____KPOINTS
|____BASE
|____mol_tran.ipynb
|____09ts02_tran1
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____04Cu_O
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____03Cu_O_PO
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____00Cu
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____01Cu_O2
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____05Cu_O_pro
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____08ts01_tran2
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____07ts01_tran1
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____06Cu_PO
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS
|____02Cu_O2_pro
| |____INCAR
| |____CONTCAR
| |____POSCAR
| |____KPOINTS

注：

由于pymatgen导出POSCAR不会重排原子顺序，如果直接导出，会出现：H，C，H，N，C这种古怪的POSCAR，所以代码写得比较长，这样导出的时候会以Species_df中的元素顺序输出POSCAR，这样过渡态插点不会乱。

由于我在用的服务器没连接外网，安不了pymatgen，所以我涉及到pymatgen和ase的代码一直是放在我的macbook上跑的。而我生成POTCAR一般会在服务器上用vaspkit，我怕混淆赝势版本，所以一直没在本地加PBE的路径。如果需要在生成的文件夹里加POTCAR，可以加上下面的代码（需要安装了vaspkit）：


#! /usr/bin/env python
# -*- coding: UTF-8 -*-
import os
import time


for root, dirs, files in os.walk(os.getcwd()):
        dir_list = dirs
        break

cdir_str = os.getcwd()
path_list = cdir_str.split("/")
cur_filename = path_list[-1]]

vaspkit_103_comm = '(echo 402)|vaspkit > vaspkit.out'

for dirs in dir_list:
    tar_dir = os.getcwd()+'/'+dirs
    os.chdir(tar_dir)
    os.popen(vaspkit_103_comm)
    time.sleep(4)

随记：

后面一周打算写类似vaspkit 501、502功能做热力学校正的脚本。之前师姐让我写下这个脚本，一直摆烂没动哈哈哈。主要是vaspkit一次只能生成一个温度的零点能和TS校正，很烦。根本没办法看能量-温度的变化，数据导出也不太行，我每次都是复制粘贴到EXCEL整理，立个flag下周来写。

xmu已经放假，看见几个老师都问我怎么还不回家哈哈，卢嘉锡基本没啥人来办公室了。

最近在搞毕设的实验，估计隔了半年没做过实验了，真的生疏了。毕设再不做毕不了业了😭，准备这两周把仪器搭完，回家安心写代码。准备溜啦~~~

2020.12.18于卢嘉锡 300

附一张下雨天我超级有氛围感的工位：