文章背景

在科研工作中，我的工作需要接触大量的石墨烯团簇结构。对结构掺入一个氮时，氮的分布位置可以分为三类：在团簇内部，团簇边缘，团簇空位附近。对于后两种，一般会形成吡啶氮，但掺氮后我们得到的结构其吡啶氮往往还包括了氢，即有个氢与吡啶氮相连。这是不合理的，因此需要将与吡啶氮相连的氢删除掉。
常规操作是逐次打开结构，然后手动删除相应的氢。但如果有几千个这样的结构，手动删除需要花费很长的时间。我一开始是手动删除，执行了一个早上跟晚上也只删除了六百多个，还有两千多个。这样的重复工作让我不禁思考，怎么样可以让程序自动识别出吡啶氮的位置，然后找到相连的氢，自动把它删除呢？
一切自动执行的任务，其实现的核心是程序员对任务过程每个步骤的规则都烂熟于心，然后采用数据与逻辑的方法将每个步骤程序化。结合循环便可以实现批量自动执行。
因此，我针对这个任务的特点，重新剖析了每个实现步骤，尝试使用简单的数学方法结合for循环、if条件结构将步骤程序化。

任务内容

打开石墨烯团簇结构，删除与吡啶氮相连的氢原子
在这里插入图片描述

程序实现思路

自动打开石墨烯团簇结构，识别出与吡啶氮相连的氢原子，自动删除
在这里插入图片描述

实现代码

建立标准结构中边缘碳与氢的位置差值标准数据集

def get_str_NH_varyposition(str_path=r'D:\software output files\initial_str_addH',str_file='POSCAR0'):
    '''
    该函数实现自动识别吡啶氮，找出与吡啶氮相连的氢的index,并进行删除
    完整设计思路为：
    1.首先读入一个完整的结构，提取出所有吡啶氮与氢的位置信息
    2.计算吡啶氮与对应氢之间的位置差，存入列表，作为标准数据
    3.读入待删除含有吡啶氮相连氢的结构，提取所有氮、氢的位置
    4.使用迭代，计算氮与氢的相对位置，与标准数据进行比对，如果完全一致，则识别出与吡啶氮相连氢的序号
    5.删除对应的氢，将新结构存入新文件路径
    本函数实现步骤1，2
    str_path：为标准结构所在文件夹的路径
    str_file：为标准结构文件的名字
    '''
    #读入结构
    str_atom=read(os.path.join(str_path,str_file),format='vasp')#读入结构信息，转为atoms object
    #获取结构位置信息
    str_position=str_atom.get_positions()
    #提取结构中所有的N
    N2=str_atom[[atom.index for atom in str_atom if atom.symbol=='N']]
    #提取结构中所有的H
    H27=str_atom[[atom.index for atom in str_atom if atom.symbol=='H']]
    #获取N,H的位置
    H27_position=H27.get_positions()
    N2_position=N2.get_positions()
    #边缘氮、空位氮与相连氢的坐标差
    Npev_index=[73,118,57,45,29,81,92,83,89,93,117,64,65,59,72,31,32,25,70,61,27,53,54,99,102,103,50]
    Hev_index=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,119,120,121]
    N_H_varyx_L=[]
    N_H_varyy_L=[]
    for i in range(0,len(Npev_index),1):
        Npev_position=str_position[Npev_index[i]]
        Hev_position=str_position[Hev_index[i]]
        N_H_varyx=Npev_position[0]-Hev_position[0]
        N_H_varyy=Npev_position[1]-Hev_position[1]
        N_H_varyx_L.append(N_H_varyx)
        N_H_varyy_L.append(N_H_varyy)
    #print("N_H_varyx_L",N_H_varyx_L)
    #print("N_H_varyy_L",N_H_varyy_L)
    return N_H_varyx_L,N_H_varyy_L

读入待修改结构，识别氮与氢位置差值是否存在标准数据集

def find_NpH_index(str_folder,save_path):
    '''
    该函数实现自动识别吡啶氮，找出与吡啶氮相连的氢的index,并进行删除
    完整设计思路为：
    1.首先读入一个完整的结构，提取出所有吡啶氮与氢的位置信息
    2.计算吡啶氮与对应氢之间的位置差，存入列表，作为标准数据
    3.读入待删除含有吡啶氮相连氢的结构，提取所有氮、氢的位置
    4.使用迭代，计算氮与氢的相对位置，与标准数据进行比对，如果完全一致，则识别出与吡啶氮相连氢的序号
    5.删除对应的氢，将新结构存入新文件路径
    本函数实现步骤3，4，5
    str_folder：为待修改结构文件所在文件夹路径
    save_path：为存储修改后的结构文件的存储文件夹路径
    '''
    str_files=os.listdir(str_folder)
    for file in str_files:
        str_PATH=os.path.join(str_folder,file)
        str_atom=read(str_PATH,format='vasp')#读入结构信息，转为atoms object
        #获取结构位置信息
        str_position=str_atom.get_positions()
        #提取结构中所有的N
        N2=str_atom[[atom.index for atom in str_atom if atom.symbol=='N']]
        #提取结构中所有的H
        H27=str_atom[[atom.index for atom in str_atom if atom.symbol=='H']]
        #获取N,H的位置
        H27_position=H27.get_positions()
        N2_position=N2.get_positions()
        #获取吡啶氮与氢位置相差的标准数据
        N_H_varyx_L,N_H_varyy_L=get_str_NH_varyposition()
        #识别读入的氢与氮的相对位置差
        H_delete_index_L=[]#建立新列表存储识别出的氢列表
        for i in range(0,len(N2_position),1):
            for j in range(0,len(H27_position),1):
                N2_H27_varyx=N2_position[i][0]-H27_position[j][0]#计算x坐标差值
                #print("N2_H27_varyx",N2_H27_varyx)
                N2_H27_varyy=N2_position[i][1]-H27_position[j][1]#计算y坐标差值
                #print("N2_H27_varyy",N2_H27_varyy)
                z_L=[z for z in N_H_varyx_L if N2_H27_varyx==z]#识别坐标差值是否符合标准数据
                k_L=[k for k in N_H_varyy_L if N2_H27_varyy==k]
                if not z_L == []:
                    if not k_L == []:
                        str_position_L=str_position.tolist()#np.array转为list
                        H27_position_j=H27_position[j].tolist()
                        H_delete_index=str_position_L.index(H27_position_j)
                        H_delete_index_L.append(H_delete_index)
        #print("H_delete_index_L",H_delete_index_L)
        del str_atom[[i for i in H_delete_index_L]]#删除识别出来的与吡啶氮相连的氢原子
        write(os.path.join(save_path,file),str_atom,format='vasp')#将修改后的结构保存

函数调用：

str_folder=r'D:\software output files\initial_str_addH'
save_path=r'D:\software output files\initial_str_auto_deleteH'
find_NpH_index(str_folder,save_path)