石墨烯平面与水平面角度计算

0.写在前面

  本用例是在CentOS7系统上进行运行，使用的软件版本为python3.7。用以解决计算分子动力学模拟中石墨烯层的取向角，代码实现思路是通过将模型均匀分成十二层，识别每层中六元环与水平面的夹角，并求平均，作为该层的取向角。
  文末提供所有代码和测试用例的下载链接。

1.代码实现

  该代码分为三部分，分别是carbon.py找到体系中所有的六元环，segmentation.py分割体系为n层，计算角度.py逐帧计算每层角度并统计输出。

1.1 carbon.py

  将该文件和ICO.so放入文件夹中，并创建子文件A，A中放入test.lammpstrj轨迹文件，运行

python carbon.py -I A -t 1

import argparse
import ctypes
import os
from math import fabs
import numpy as np
from goto import with_goto

cpp = ctypes.CDLL(r'./ICO.so')
cpp.angleCheck.restype = ctypes.c_double
def vector(atomi, atomj):
    disx = atomi[1] - atomj[1]
    disy = atomi[2] - atomj[2]
    disz = atomi[3] - atomj[3]
    return (disx, disy, disz)
def checkBoundary(x, y, z):
    (x1, x2) = x.split(maxsplit=2)
    (y1, y2) = y.split(maxsplit=2)
    (z1, z2) = z.split(maxsplit=2)
    LX = float(x2) - float(x1)
    LY = float(y2) - float(y1)
    LZ = (float(z2) - float(z1))
    return LX, LY, LZ, x2, x1, y2, y1, z2, z1
def checkPBC(matrix, box,cutoff):
    [LX, LY, LZ, x2, x1, y2, y1, z2, z1] = box
    listcopy = []
    dO_O = cutoff
    for atomsarray in matrix:
        dx = fabs(atomsarray[0][1] - float(x2))
        dy = fabs(atomsarray[0][2] - float(y1))
        dz = fabs(atomsarray[0][3] - float(z1))
        if (dx <= dO_O or dy <= dO_O or dz <= dO_O):
            if dx <= dO_O:
                movex = -LX
                k = [atomsarray[0][0], atomsarray[0][1] + movex, atomsarray[0][2], atomsarray[0][3]]

                listcopy.append([k])
            if dy <= dO_O:
                movey = LY
                k = [atomsarray[0][0], atomsarray[0][1], atomsarray[0][2] + movey, atomsarray[0][3]]

                listcopy.append([k])
            if dz <= dO_O:
                movez = LZ
                k = [atomsarray[0][0], atomsarray[0][1], atomsarray[0][2], atomsarray[0][3] + movez]

                listcopy.append([k])
            if dx <= dO_O and dy <= dO_O:
                movex = -LX
                movey = LY
                movez = 0
                k = [atomsarray[0][0], atomsarray[0][1] + movex, atomsarray[0][2] + movey, atomsarray[0][3] + movez]

                listcopy.append([k])
            if dx <= dO_O and dz <= dO_O:
                movex = -LX
                movez = LZ
                movey = 0
                k = [atomsarray[0][0], atomsarray[0][1] + movex, atomsarray[0][2] + movey, atomsarray[0][3] + movez]

                listcopy.append([k])
            if dy <= dO_O and dz <= dO_O:
                movex = 0
                movey = LY
                movez = LZ
                k = [atomsarray[0][0], atomsarray[0][1] + movex, atomsarray[0][2] + movey, atomsarray[0][3] + movez]

                listcopy.append([k])
            if dx <= dO_O and dy <= dO_O and dz <= dO_O:
                movex = -LX
                movey = LY
                movez = LZ
                k = [atomsarray[0][0], atomsarray[0][1] + movex, atomsarray[0][2] + movey, atomsarray[0][3] + movez]

                listcopy.append([k])
        else:
            pass
    return listcopy

def data_extractor(line):
    (aid, atype, ax, ay, az) = line.split(maxsplit=4)
    return [float(aid), float(ax), float(ay), float(az)]


def search_neighbor(matrix, atomcopy,cutoff):
    listj = []
    listneighboring = []
    arr, col, _ = matrix.shape
    dO_O = cutoff
    dO_O2 = cutoff ** 2

    for i in range(arr):
        idi = matrix[i][0][0]
        dist = np.square(matrix[i][0][1:] - matrix[:, 0, 1:]).sum(axis=1)
        for j in np.where(dist <= dO_O2)[0]:
            da = dist[j]
            if i == j:
                continue
            else:

                if j not in listj:
                    listj.append(j)

        whereCopy = np.where(atomcopy[:, 0, 0] == idi)
        for copyorder in whereCopy[0]:
            dist = np.square(atomcopy[copyorder][0][1:] - matrix[:, 0, 1:]).sum(axis=1)
            for order in np.where(dist <= dO_O2)[0]:
                distic = dist[order]


                if order not in listj:
                    listj.append(order)

        for copyorder in whereCopy[0]:
            dist = np.square(atomcopy[copyorder][0][1:] - atomcopy[:, 0, 1:]).sum(axis=1)
            for N in np.where(dist <= dO_O2)[0]:
                distic = dist[N]
                if N == copyorder:
                    continue
                else:

                    wherem = np.where(matrix[:, 0, 0] == atomcopy[N][0][0])
                    for order in wherem[0]:
                        if order not in listj:
                            listj.append(order)

        dist = np.square(matrix[i][0][1:] - atomcopy[:, 0, 1:]).sum(axis=1)
        for orderC2 in np.where(dist <= dO_O2)[0]:
            distci = dist[orderC2]
            if distci >= 3 * dO_O2:
                print('wrong!/n')

            wherem = np.where(atomcopy[orderC2][0][0] == matrix[:, 0, 0])
            for order in wherem[0]:
                if order not in listj:
                    listj.append(order)

        listneighboring.append(listj.copy())
        listj = []
    return listneighboring
@with_goto
def ring_identification(listmap):
    list4ring = []
    list4pair = []
    list5ring = []
    list5pair = []
    list6ring = []
    list6pair = []
    for host in range(len(listmap)):
        if len(listmap[host]) <= 2:
            pass
        else:
            for node1 in listmap[host]:
                if len(listmap[node1]) < 2:
                    continue
                else:
                    for node2 in listmap[node1]:
                        if node2 == host:
                            label.new2
                            continue
                        elif len(listmap[node2]) < 2:
                            continue
                        for node3 in listmap[node2]:
                            if (node3 in listmap[host]) and (node3 != node1):

                                goto.new2
                        for node3 in listmap[node2]:
                            if node3 == host:
                                continue
                            elif node3 == node1:
                                continue
                            elif len(listmap[node3]) < 2:
                                continue
                            for node4 in listmap[node3]:
                                if (node4 in listmap[host]) and (node4 not in [node1, node2, host]):

                                    goto.new2
                            for node4 in listmap[node3]:
                                if node4 == node1:
                                    continue
                                elif node4 == node2:
                                    continue
                                elif len(listmap[node4]) < 2:
                                    continue
                                for node5 in listmap[node4]:
                                    if (node5 in listmap[host]) and (node5 not in [node1, node2, node3, host]):
                                        sixMRing = {host, node1,
                                                    node2, node3,
                                                    node4, node5}
                                        if sixMRing not in list6ring:
                                            list6pair.append([{host, node1},
                                                              {node1, node2},
                                                              {node2, node3},
                                                              {node3, node4},
                                                              {node4, node5},
                                                              {node5, host}])
                                            list6ring.append(sixMRing)
                                    else:
                                        continue
    return list6ring, list6pair


def cage_identification(data,filename,output,cutoff):
    listRing = []
    listUnion = []
    listNeighboring = []
    listStep = []
    listSets = []
    listStat = []
    wallTime = [0, 0, 0, 0, 0, 0]
    Frame=0

    for eachData in data:
        out=os.path.join(output, str(Frame)+'.txt')
        copyDB = checkPBC(eachData[0], eachData[2],cutoff)
        listn = search_neighbor(eachData[0], np.array(copyDB),cutoff)

        print('neighbor')
        listr6, listp6 = ring_identification(listn)
        print('done')


        listr6 = np.array(listr6)
        print(len(listr6))


        Results=[eachData[0][list(r)][:, 0, 0] for r in listr6]

        np.savetxt(out, Results, delimiter=' ',fmt='%d')





        Frame+=1

    print('\n')
    for each in wallTime:
        print(each)
    return  []
def is_suffix_lmp(suffix: str):
    if suffix == '.lammpstrj':
        return True
    return True
def mainfun(args, fn,output):
    filename = fn
    PATH = os.path.abspath((os.path.abspath(filename)))
    foldername = args.i
    DIR = os.path.dirname(PATH)
    name, suffix = os.path.splitext(PATH)

    CAGEStype = []
    cagetype_total = []

    listdata = []
    listUnion = []
    listRing = []
    listNeighboring = []
    listStep = []
    listSets = []
    listStat = []
    if is_suffix_lmp(suffix):
        os.system('clear')
        print(DIR + '/' + foldername + '/' + filename)
        datafile = open(DIR + '/' + foldername + '/' + filename, 'r')
        ending = datafile.seek(0, 2)
        datafile.seek(0, 0)
        print('Loading...')
        data=datafile.readlines()
        data=np.array(data)
        frame=np.sum(data == 'ITEM: NUMBER OF ATOMS\n')

        listdata=[]
        for f in range(frame):
            listatom = []
            start=np.where(data == 'ITEM: NUMBER OF ATOMS\n')[0] + 7
            title=data[start[f]-8]

            aaa = data[start[f]-6].strip('\n')
            aaa= aaa.split()
            natoms=int(aaa[0])
            atoms = data[start[f]:start[f]+natoms-1]
            (x1, x2)= data[start[f] - 4].split(maxsplit=2)
            (y1, y2) = data[start[f] - 3].split(maxsplit=2)
            (z1, z2)= data[start[f] - 2].split(maxsplit=2)
            LX = float(x2) - float(x1)
            LY = float(y2) - float(y1)
            LZ = (float(z2) - float(z1))
            title=title.split()[0]


            for atomstr in atoms:
                (aid, atype, ax, ay, az) = atomstr.split(maxsplit=4)
                if atype == args.t:
                    O = [float(aid), float(ax), float(ay), float(az)]
                    listatom.append([O])
            atomnum = len(listatom)
            print(title,atomnum)
            listdata.append([np.array(listatom),title , [LX, LY, LZ, x2, x1, y2, y1, z2, z1],atomnum ])





        print('Done!\n\n\n')



        cage_identification(np.array(listdata),name,output,args.c)




if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description='Notice:\n' + '\n 1.The code and the folder containing the trajectories to be analyzed should be in the same directory.\n' + ' 2.trajectories must be in lammpstrj format and contain only water molecules.')
    parser.add_argument('-i', type=str, default='1101', help="Path of folder containing the trjs to be analysed")
    parser.add_argument('-t', type=str, default='1', help="Symbol of Carbon")
    parser.add_argument('-c', type=float, default=2.0, help="Cutoff")

    args = parser.parse_args()


    foldername = args.i

    PROJECT_DIR_PATH = os.path.dirname(os.path.abspath(os.path.abspath(__file__)))
    DIR_PATH = os.path.join(PROJECT_DIR_PATH, foldername)
    files = os.listdir(DIR_PATH)
    ringslist=[]
    for filename in files:
        abs=os.path.abspath(args.i)


        upper=os.path.split(abs)[0]

        output=os.path.join(upper, os.path.splitext(filename)[0])
        try:
            os.mkdir(output)
        except FileExistsError:
            pass
        else:
            pass
        mainfun(args, filename,output)

  -i 为子文件夹名字，-t为原子类型，可选变量-c为截断半径，默认为2.0   运行输出文件为lammpstrj轨迹同名的文件夹，与carbon.py在同一目录。将输出文件夹内的文件下载到本地文件夹ring中。

1.2 segmentation.py

import os

h_f= -25.0 #分层起始下限 
h_e = 25  #分层上限
n_layer = 5  #层数
trj = open("test.lammpstrj",'r')  #改为lammpstrj文件名称



h_c= (h_e-h_f)/n_layer #每层高度
if not os.path.exists(".//trajectory"):
    os.mkdir(".//trajectory")
for i in range(n_layer):
    if not os.path.exists(".//trajectory//"+str(i+1)):
        os.mkdir(".//trajectory//"+str(i+1))
name = 0
flag = 0 #判定是否开始写入
for line in trj:
    line = line.strip('\n')
    if(line == 'ITEM: TIMESTEP'):
        flag = 1
        files = {f'file{jj}': open(f".//trajectory//"+str(jj)+"//"+str(name)+".txt", "w") for jj in range(1, n_layer + 1)}

        name =name + 1
        continue
    flag = flag + 1
    if(flag > 9):
        data = line
        data = data.strip('\n')
        data = data.split(' ')
        cs = int((float(data[4])-h_f)/((h_e-h_f)/n_layer))+1
        if(cs<=n_layer and cs>=1):
            files[f'file{cs}'].writelines(line+'\n')

  修改py文件中的h_f为分层的起始高度值，h_e为分层的结束高度值，n_layer为层数，trj为lammpstrj文件名。
  程序原理是将z轴进行分层，分为均匀n_layer层。运行后会生成trajectory文件夹，该文件与ring文件夹处于同一目录下，并且和.py文件处于同一文件中。

1.3 计算角度

import numpy as np
import os

framenum = 2 #帧数
n_layer = 5 #切片层数

def calculate_angle_between_triangle_and_horizontal_plane(x, y, z):
    # 计算三角形的法向量
    normal_vector = np.cross(y - x, z - x)

    # 水平面的法向量是(0, 0, 1)
    horizontal_vector = np.array([0, 0, 1])

    # 计算两个向量的夹角
    dot_product = np.dot(normal_vector, horizontal_vector)
    magnitudes = np.linalg.norm(normal_vector) * np.linalg.norm(horizontal_vector)

    # 使用 arccos 计算夹角（弧度）
    angle_radians = np.arccos(dot_product / magnitudes)

    # 将弧度转换为度
    angle_degrees = np.degrees(angle_radians)

    return angle_degrees



if not os.path.exists(".//angle"):
    os.mkdir(".//angle")

files = {f'file{jj}': open(f".//angle//"+str(jj)+".txt", "w") for jj in range(1, n_layer + 1)}
all = open('.//angle//all.txt','w')
for i in range(framenum):

    rings = []
    source_ring = open('ring//'+str(i)+'.txt')
    for line in source_ring:
        list_ring = []
        line = line.strip('\n')
        data = line.split(' ')
        list_ring.append(int(data[0]))
        list_ring.append(int(data[1]))
        list_ring.append(int(data[2]))
        list_ring.append(int(data[3]))
        list_ring.append(int(data[4]))
        list_ring.append(int(data[5]))
        rings.append(list_ring)
    source_ring.close()
    all.writelines(str(i)+' ')
    for j in range(n_layer):
        angles = []  # 保存角度
        layerC = open('trajectory//'+str(j+1)+'//'+str(i)+'.txt')
        xh=[]
        Cx=[]
        Cy=[]
        Cz=[]
        setC = set()
        for line in layerC:
            line = line.strip('\n')
            data = line.split(' ')
            xh.append(int(data[0]))
            setC.add(int(data[0]))
            Cx.append(float(data[2]))
            Cy.append(float(data[3]))
            Cz.append(float(data[4]))
        layerC.close()
        for k in range(len(rings)):
            if((rings[k][0] in setC) and (rings[k][1] in setC) and (rings[k][2] in setC) and (rings[k][3] in setC) and (rings[k][4] in setC) and (rings[k][5] in setC)):
                a0 = -1
                a1 = -1
                a2 = -1
                for kk in range(len(xh)):
                    if(rings[k][0] == xh[kk]):
                        a0 = kk
                    elif(rings[k][2] == xh[kk]):
                        a1 = kk
                    elif (rings[k][4] == xh[kk]):
                        a2 = kk
                a0_list = []
                a1_list = []
                a2_list = []
                a0_list.append(Cx[a0])
                a0_list.append(Cy[a0])
                a0_list.append(Cz[a0])
                a1_list.append(Cx[a1])
                a1_list.append(Cy[a1])
                a1_list.append(Cz[a1])
                a2_list.append(Cx[a2])
                a2_list.append(Cy[a2])
                a2_list.append(Cz[a2])
                a00 = np.array(a0_list)
                a11 = np.array(a1_list)
                a22 = np.array(a2_list)
                jd = calculate_angle_between_triangle_and_horizontal_plane(a00, a11, a22)
                angles.append(jd)
        for ll in range(len(angles)):
            if(angles[ll]>90):
                angles[ll] = abs(angles[ll]-180)
        sum = 0
        for ll in range(len(angles)):
            sum = sum + angles[ll]
        ave = 0
        if(len(angles) != 0):
            ave = sum / len(angles)

        files[f'file{j+1}'].writelines(str(ave) + '\n')

        all.writelines(str(ave) + ' ')
    all.writelines('\n')
    print(str(i)+'已完成')

  修改framenum为帧数，n_layer为分层层数（需要和1.2中分层数一致）。运行后输出文件夹angle，angle文件夹下有all.txt文件，存储角度信息。竖列表示帧，横行表示层。

2.代码下载

测试用例
carbon.py
ICO.so
segmentation.py
计算角度.py