计算二面角随时间变化情况

最新推荐文章于 2024-09-17 23:15:58 发布

写代码的信哥

最新推荐文章于 2024-09-17 23:15:58 发布

阅读量322

点赞数

分类专栏：科研实例文章标签：几何学 python 算法

本文链接：https://blog.csdn.net/m0_51706126/article/details/125343170

版权

科研实例专栏收录该内容

2 篇文章 0 订阅

订阅专栏

该博客介绍了一种计算咪唑环分子中咪唑环与x0y面二面角随时间变化的方法。通过读取包含6000帧512个分子数据的输入文件，选取每个分子的N1，N2，C4三个原子坐标，计算面上向量、法向量，并利用向量夹角公式得出二面角。计算结果以文本文件形式保存，便于后续分析。整个过程实现了高效的计算，对理解分子动力学有重要意义。

摘要由CSDN通过智能技术生成

目标

计算包含咪唑环分子的咪唑环与x0y面形成的二面角的随时间变化情况。咪唑环的结构相对稳定，取咪唑环上的三个原子，N1，N2,C4，作为形成面的三个点。

输入文件

输入文件中包括6000帧时间的数据，每一帧是512个分子的数据。
在本代码运行中，我只取了一个分子中的三个原子。也就是每一帧有1536行的数据。
每一行的格式为

N 13.74656 12.65898 11.56485

除此之外，为了区分每一帧，每一帧1536行数据上还有两行提示信息。+

方法

计算面上向量

首先根据三个原子的位置信息，也就是三个点的坐标，计算出面上两个向量的位移值。

def calc_vec(poix, poiy, poiz): #通过三点的位置，计算两个向量
    #poix,poiy,poiy都是有三个元素的数组，三个元素分别是点的x,y,z的值
    vec1 = [0] * 3
    vec2 = [0] * 3
    vec1[0] = poix[0] - poiy[0]
    vec1[1] = poix[1] - poiy[1]
    vec1[2] = poix[2] - poiy[2]
    vec2[0] = poiz[0] - poiy[0]
    vec2[1] = poiz[1] - poiy[1]
    vec2[2] = poiz[2] - poiy[2]
    return [vec1,vec2]

计算法向量

根据两个向量的外积，计算法向量的值。

def calc_normal_vector(vec1, vec2):  #通过面上的两个向量计算法向量
    norm_vec = [0] * 3
    norm_vec[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1]
    norm_vec[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2]
    norm_vec[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0]
    return norm_vec

根据公式计算法向量与x0y面的法向量的夹角

学习写latex公式
根据下面的求向量的夹角的公式，计算得到这两个法向量的夹角。
$\cos \theta = \frac{\vec {n1} \times \vec{n2}}{|\vec {n1}||\vec {n1}|}$

def clac_dihedral(norm1, norm2): # 通过两个法向量计算二面角的夹角
    cos_theta = 0
    abs_norm1 = math.sqrt(norm1[0] ** 2 + norm1[1] ** 2 + norm1[2] ** 2)
    abs_norm2 = math.sqrt(norm2[0] ** 2 + norm2[1] ** 2 + norm2[2] ** 2)
    cos_theta = (norm1[0] * norm2[0] + norm1[1] * norm2[1] + norm1[2] * norm2[2])/(abs_norm1 * abs_norm2)
    theta_pi = math.acos(cos_theta)
    theta = theta_pi * (180/math.pi)
    return theta

计算

计算结构保存，因为数据量比较大，使用数组的方式保存是行不通的，所以我将结果输出到文件中。然后就可以直接读取分析数据了。

def calc(filepath,out):
    cfile = open(out, "w")
    for j in range(6000): # 文件中每一帧数据前面有两行说明文字，要跳过。总共6000帧
        for i in range(512):  # 0-511
            index_N1 = 1 + 2 + 3 * i + 1538 * j # N1原子的位置信息在文件中对应的行数 3,6,9，...，512*3
            index_N2 = 2 + 2 + 3 * i + 1538 * j
            index_C = 3 + 2 + 3 * i + 1538 * j
            lineN1 = linecache.getline(filepath, index_N1)
            lineN2 = linecache.getline(filepath, index_N2)
            lineC = linecache.getline(filepath, index_C)
            listN1 = lineN1.split()
            listN2 = lineN2.split()
            listC = lineC.split()
            pox = [float(listN1[-3]), float(listN1[-2]), float(listN1[-1])]
            poy = [float(listN2[-3]), float(listN2[-2]), float(listN2[-1])]
            poz = [float(listC[-3]), float(listC[-2]), float(listC[-1])]
            [vec1, vec2] = calc_vec(pox, poy, poz)
            norm1 = calc_normal_vector(vec1, vec2)
            norm2 = [0, 0, 1]
            theta = clac_dihedral(norm1, norm2)
            cfile.write(str(theta) + "\t") # 将一帧的512个咪唑环与x0y形成的二面的数据写成一行
        cfile.write("\n") # 每一帧占一行，得到一个6000*512 的数据文件
    print("success!")

结果

#inputfilepath是保存文件的路径。
if __name__ == '__main__':
    filename = "filename"
    outname = "outname"
    filepath = os.path.join(inputfilepath,filename)
    out = os.path.join(inputfilepath,outname)
    calc(filepath,out)

最终成功计算，计算的速度还可以，不是很慢。总体还是很满意了。