【ASE+python学习】批量实现将含有空位的结构进行氢饱和(实际就是在某个位置添加氢原子)

在搬砖过程中，我经常要对大批量的含有空位的结构进行氢饱和处理。如果一个个打开结构，然后进行加氢，我估计要累死，对颈椎及其不友好。而且这种简单且没有技术含量的工作，我想可以用更科学的方式进行实现。

任务大纲

本文以含有一个空位的石墨烯团簇为例，通过ASE与python程序对该结构进行加氢。单个结构加氢的任务大纲见下，如果是多个任务，结合for循环执行即可。

任务思路

实现给结构加氢的任务思路为：

1. 首先需要基于已有的已经进行加氢的结构，计算在空位处氢原子与三个碳原子的相对位置信息，包括氢原子偏离结构平面的距离、与碳相连的键长、以及氢原子与结构平面的夹角、氢原子与碳相连的键在结构平面的投影在x方向与y方向的夹角(大概计算就行，后期会进行结构优化，优化氢原子的位置)
2. 接着使用ASE读取结构信息，获取空位处三个碳原子的坐标信息
3. 基于三个碳原子的坐标信息，计算三个氢原子的坐标信息
4. 然后通过ASE将三个氢原子逐次添加到结构上
5. 将新结构写入新的文件中

氢原子相对位置计算

这里计算的参数大概就行，只要能让氢与相应的碳原子连接在一起，后续进行结构优化会优化到实际位置

代码实现-单个结构加氢

导入所需库/包

import os
import math
import shutil
from ase.io import read,write
from ase.visualize import view
from ase.build import molecule
from ase.build import add_adsorbate

计算氢原子坐标

def addHsite(C95Vac1_path=r'D:\software output files\quntum_experiment\POSCAR\non_periodic_CN\C96_molecule\4-C95Vac1_str_1e-2',C95Vac1_file='POSCAR9'):
    '''
    该函数用于获取空位附近三个碳原子的位置，并计算出三个氢的位置
    思路设计：
    读入含有空位的结构
    获取空位处三个碳原子的位置
    在MS打开含有空位且已进行氢饱和的结构，量碳氢键长，计算该键与xy平面的夹角，以及在xy平面的投影
    接着根据空位的排布，计算碳碳氢二面角，从而计算碳氢键长在xy平面的投影在x方向与y方向的分量
    基于空位三个碳原子的坐标，计算三个氢原子的坐标
    '''
    C95Vac1=read(os.path.join(C95Vac1_path,C95Vac1_file),format='vasp')#读入结构信息，转为atoms object
    view(C95Vac1)
    #获取空位碳原子的位置,仅提取x,y方向
    position_62=C95Vac1[62].position[:2]
    position_61=C95Vac1[61].position[:2]
    position_114=C95Vac1[114].position[:2]
    #计算x,z方向偏移量
    z=1.06*math.sin(49*math.pi/180)#1.06为氢键的长度，这里计算的是氢原子偏离结构平面的高度
    x=1.06*math.cos(49*math.pi/180)#这里计算的是氢原子与碳原子所成的键，在xy平面的投影长度
    #分别计算与C62,C61,C114相连氢原子的位置
    #C62-H位置计算
    position_62_newx=[i+x for i in position_62 if i ==position_62[0]]#计算在x方向氢原子与C62位置的偏移量
    position_62_new=[]#建立空列表用于存储氢原子在平面的坐标
    position_62_new.append(position_62_newx[0])
    position_62_new.append(position_62[1])
    #计算C61与C114的氢原子在x,y方向的偏移量
    xx=x*math.cos(71*math.pi/180)#71为二面角
    xy=x*math.sin(71*math.pi/180)
    #C61-H位置计算
    position_61_newx=[i-xx for i in position_61 if i ==position_61[0]]
    position_61_newy=[i+xy for i in position_61 if i ==position_61[1]]
    position_61_new=[]
    position_61_new.append(position_61_newx[0])
    position_61_new.append(position_61_newy[0])
    #C114-H位置计算
    position_114_newx=[i-xx for i in position_114 if i ==position_114[0]]
    position_114_newy=[i-xy for i in position_114 if i ==position_114[1]]
    position_114_new=[]
    position_114_new.append(position_114_newx[0])
    position_114_new.append(position_114_newy[0])
    return x,z,position_62_new,position_61_new,position_114_new

给结构加氢

def str_addH(C95Vac1_path=r'D:\software output files\quntum_experiment\POSCAR\non_periodic_CN\C96_molecule\4-C95Vac1_str_1e-2',C95Vac1_file='POSCAR9'):
    '''
    该函数用于给结构的空位位置附件添加氢原子
    思路：基于函数addHsite()所得到的三个氢原子的坐标，依次添加到结构表面，接着将修改后的结构写入新的POSCAR文件中
    '''
   x,z,position_62_new,position_61_new,position_114_new=addHsite()
    H=molecule('H')
    C95Vac1=read(os.path.join(C95Vac1_path,C95Vac1_file),format='vasp')#读入结构信息，转为atoms object
    add_adsorbate(C95Vac1,H,position=position_62_new,height = z)#在结构上依次添加氢原子
    add_adsorbate(C95Vac1,H,position=position_61_new,height = z)
    add_adsorbate(C95Vac1,H,position=position_114_new,height = z)
    write(os.path.join(C95Vac1_path,C95Vac1_file+'_addH'),C95Vac1,format='vasp')#将添加了氢原子的结构写入文件
    view(C95Vac1)

运行结果

调用函数运行

str_addH()

运行结果

代码实现-批量结构加氢

批量给结构加氢，实际上只需要在读入结构那部分代码改为使用for循环遍历文件夹中的结构即可。

导入相关库/包

import os
import math
import shutil
from ase.io import read,write
from ase.visualize import view
from ase.build import molecule
from ase.build import add_adsorbate

遍历文件夹内所有结构，并逐次给结构加氢

def batchstr_addH(file_path,save_str,save_folder):
    '''
    该函数实现批量读入含有空位的结构，并将三个氢原子添加到空位附近，实现氢饱和，并将新结构写入文件，存入新文件夹
    file_path:结构文件所处的文件夹路径
    save_str:添加氢原子后的结构存储文件夹所在路径
    save_folder：添加氢原子后的结构存储文件夹的名字
    '''
    save_path=os.path.join(save_str,save_folder)
    str_files=os.listdir(file_path)
    x,z,position_62_new,position_61_new,position_114_new=addHsite()
    H=molecule('H')
    if os.path.exists(save_path):
        shutil.rmtree(save_path)
    if not os.path.exists(save_path):
        os.makedirs(save_path)
    for str_file in str_files:#遍历结构文件
        structure=read(os.path.join(file_path,str_file),format='vasp')
        add_adsorbate(structure,H,position=position_62_new,height = z)#在结构上依次添加氢原子
        add_adsorbate(structure,H,position=position_61_new,height = z)
        add_adsorbate(structure,H,position=position_114_new,height = z)
        write(os.path.join(save_path,str_file),structure,format='vasp')

file_path=r'D:\software output files\quntum_experiment\POSCAR\non_periodic_CN\C96_molecule\5-C94Vac1N1_str_1e-2(based_on_poscar9)_ASE\Npe'
save_str=r'D:\software output files\quntum_experiment\POSCAR\non_periodic_CN\C96_molecule\5-C94Vac1N1_str_1e-2(based_on_poscar9)_ASE'
save_folder='Npe_addH'
batchstr_addH(file_path,save_str,save_folder)

批量实现的逻辑与结果

如此，便实现了自动批量加氢的操作。