Python实现：已知化学分子的输入文件坐标（高斯计算输入文件为例），求其中任意三个原子确定的平面的法向量和单位法向量

计算化学在处理实际化学问题时，比如需要在某一化学平面的法向量上进行分子操作，这时最重要的是确定化学平面和求法向量，才能进行后续的操作（如下图所示），下面以高斯输入文件为例，用python代码实现该功能，包括以下两部分：

1. 确定化学平面：一般通过三个点确定，也就是三个原子的坐标确定一个化学平面；

2. 计算法向量和单位法向量

以下输入文件为例：

以下信息复制到新文件，命名为xx.gjf，具体如下：

%nprocshared=20
%mem=2000MB
%chk=xx.chk
# ub3lyp/6-31g(d) geom=connectivityTitle Card Required0 2C                 -0.07799900   -0.06595600    1.16020400C                  0.28486000    0.24116000   -0.05689700Sc                -0.57500900   -2.04502600    0.35332800Sc                 2.15465200    0.28669700    1.05119200Sc                -1.43262200    1.53516600    0.33122200C                  1.07300100   -3.58961700   -0.99271100C                  0.64679900   -2.84267000   -2.13946700C                  2.30603300   -3.23292800   -0.33821000C                  3.11036400   -2.22051000   -0.93080800C                  2.69014500   -1.51047400   -2.10623700C                  1.43348100   -1.78981800   -2.72004200C                 -0.09411100   -4.12378500   -0.30887600C                 -0.78412200   -2.83577100   -2.15917600C                  2.34170400   -3.36823900    1.07365900C                  3.87523700   -1.31082500   -0.11059700C                  3.18771700   -0.16842900   -2.03077700C                 -1.27997100   -3.64164900   -1.05489200C                  0.73221900   -0.72542800   -3.34429800C                  3.93218000   -0.03421600   -0.79817100C                 -0.07297600   -4.15203700    1.15172500C                  1.15383000   -3.76689100    1.79966100C                  3.10083400   -2.46616800    1.89316300C                 -1.47865600   -1.70099900   -2.72341900C                  3.89718200   -1.41037500    1.31997700C                 -0.70719700   -0.68212400   -3.34099200C                  1.22674700    0.62270500   -3.27453000C                  2.43358500    0.91972100   -2.57164700C                 -2.47198800   -3.17244600   -0.35393800C                  1.19710600   -3.10035100    3.08300900C                  2.38588400   -2.30649100    3.13804200C                 -1.29167200   -3.75279500    1.82824600C                  4.01126900    1.20707100   -0.07146000C                 -1.10617500    0.70225700   -3.25233900C                 -2.70364300   -1.33578400   -2.09900400C                  0.09109300    1.50307100   -3.20649800C                  4.07910900   -0.17705300    2.08089700

简单说明下该文件信息，前8行为计算信息，其中前四行为计算的说明，分别代表：

%nprocshared=20，这表示cpu使用核数为20个
        %mem=2000MB，这表示占用缓存大小为2000MB
        %chk=xxxx.chk，这会输出计算信息包括波函数等细节
        # ub3lyp/6-31g(d) geom=connectivity: 这是计算方法关键词部分，第一个字符 # 为固定样式

空一行
Title Card Required，一般为文件说明，可随意更改

空一行

0 2，这里是计算的电荷数为0，自旋多重度为2

。。。。然后就是具体的分子三维坐标

。。。。

任意三个坐标都可以确定一个平面，而读取坐标的逻辑如下：

1. 输入三个原子对应的序号，原子的序号在文件中实际上是按照从前到后依次排列，但是行号与原子序号又不能一一对应，需要进一步处理。

从该输入文件看，行号的第9行才对应第一个原子，在写代码时直接将输入的数字加8可以实现，但是为了代码的普适性，也就是任何高斯输入文件都可以使用该python代码，该输入文件计算说明部分一共是8行，如果加入其它计算的信息，行数很可能不是上文中的8行。解决方法如下：

高斯计算文件有一个固定的格式，就是关键词部分的第一个字符一定是 # ，而且第一个坐标一定是该行之后的第5行（该行前面的计算说明部分的行数可多可少），这就是高斯输入文件的固定格式，只需要确定该行的行号，再用该行号 + 4之后的那行就对应第一个坐标。

所以：输入的原子序号 + 读取 # 字符的行号 + 4 = 该原子序号对应的坐标，这部分在下面代码中有详细的体现。

2. 将原子对应的坐标转换成可操作的矩阵（向量），通过三点确定的平面，也就是两个向量进行叉乘，即可得到法向量，再将该向量转化成单位向量。这里通过导入numpy模块可轻松实现。

3. 将计算结果写入result.txt文件

最终python代码如下，复制粘贴进pycharm即可，使用前需要安装numpy、linecache库，安装方法请自行google或者baidu，简单到令人发指，time库python已自带，直接调用即可。

注意：高斯输入文件需要放入代码所在的同一文件夹下，否则不能执行：

"""
求法向量
"""
import time
import numpy
import linecache# 将向量转化为单位向量，如果向量的0向量，返回原向量
def UnitVect(v):norm = numpy.linalg.norm(v)if norm == 0:return vreturn v / norm
# 传入三个坐标，求三点确定的面对应的法向量
def NormVect(array1, array2, array3):a = numpy.array(coord_1)b = numpy.array(coord_2)c = numpy.array(coord_3)global hh = numpy.cross(a-b,a-c)print(f"法向量是：{h}")x = UnitVect(h)return x# 输入文件名和坐标序号
filename = input("请输入文件名：")
i = int(input("请输入原子1的序号："))
j = int(input("请输入原子2的序号："))
k = int(input("请输入原子3的序号："))
# 记录开始运行程序的时间
start = time.time()
# 将输入的数字全部加 n，n指的是当第m行识别到有 # 字符之后的第5行为第一个计数的坐标，即：n = m + 4，当遇到有 # 字符时中断for循环
# range中的10可以改变，一般 # 出现在前10行以内，为了减小计算量限制在10以内（如果文件很大，又没有 # 出现会造成计算资源浪费）
for m in range(1,10):mark_line = linecache.getline(f'{filename}', m)mark_x = mark_line[0]if mark_x == '#':n = m + 4break
# 这里获取第一个坐标
cood1_line = linecache.getline(f'{filename}', i+n)
coord1_line = cood1_line.split( )
coord_name1, *coord_temp = coord1_line
coord_1 = [float(x) for x in coord_temp]
print(f"第一个原子是：{coord_name1}, 坐标是：", coord_1)
# 这里获取第二个坐标
cood2_line = linecache.getline(f'{filename}', j+n)
coord2_line = cood2_line.split( )
coord_name2, *coord_temp = coord2_line
coord_2 = [float(x) for x in coord_temp]
print(f"第二个原子是：{coord_name2}, 坐标是：", coord_2)
# 这里获取第三个坐标
cood3_line = linecache.getline(f'{filename}', k+n)
coord3_line = cood3_line.split( )
coord_name3, *coord_temp = coord3_line
coord_3 = [float(x) for x in coord_temp]
print(f"第三个原子是：{coord_name3}, 坐标是：", coord_3)
# 接下来对坐标进行操作
NormVect_val = NormVect(coord_1,coord_2,coord_3)
print(f"单位法向量是：{NormVect_val}")# 写入一个文件，用于放最后结果
with open("result.txt", mode = "w", encoding = "utf-8") as f:f.write(f"第一个原子是：{coord_name1}, 坐标是：{coord_1}\n")f.write(f"第一个原子是：{coord_name2}, 坐标是：{coord_2}\n")f.write(f"第一个原子是：{coord_name3}, 坐标是：{coord_3}\n")f.write(f"法向量是：{h}\n")f.write(f"单位法向量是：{NormVect_val}")# 记录结束运行程序的时间
end = time.time()
print(end - start)

Python实现：已知化学分子的输入文件坐标（高斯计算输入文件为例），求其中任意三个原子确定的平面的法向量和单位法向量相关推荐

Python实现“已知三角形两个直角边，求斜边”
用Python实现"已知三角形两个直角边,求斜边" 要求:用户输入两个直角边(数值为浮点类型),若非浮点类型,则提示用户,继续输入. 思路:伪代码描述下步骤 1.-input a ...
python已知两条直角边求斜边,Python实现“已知三角形两个直角边，求斜边”
用Python实现"已知三角形两个直角边,求斜边" 要求:用户输入两个直角边(数值为浮点类型),若非浮点类型,则提示用户,继续输入. 思路:伪代码描述下步骤 1.-input a ...
已知两点的经度和纬度，计算两点间的距离(php，javascript)
php代码:转载 http://www.cnblogs.com/caichenghui/p/5977431.html 1 /** 2 * 求两个已知经纬度之间的距离,单位为米 3 * 4 * @par ...
已知三角形三边求面积的c语言程序,已知三角形三边分别为4，5，6，求三角形的面积。用c语言编写程序...
已知三角形三边分别为4,5,6,求三角形的面积.用c语言编写程序关注:114 答案:6 mip版解决时间 2021-01-18 16:33 提问者谁把流年搁浅 2021-01-17 23:52 ...
旋转矩阵的应用：已知旋转前后的点坐标计算旋转中心坐标
点A(x1,y1)绕圆心C(x,y)旋转角为点B(x2,y2),已知点A.B坐标以及角,求解圆心C(x,y)坐标. 根据二维旋转公式: 根据此问题将A相对C的坐标和B相对C的坐标代入以上公式: 因为, ...
已知某个班有 30 个学生，学习 5 门课程，已知所有学生的各科成绩。计算每个学生的平均成绩，并输出
已知某个班有 30 个学生,学习 5 门课程,已知所有学生的各科成绩.计算每个学生的平均成绩,并输出注意: 定义一个二维数组 A,用于存放 30 个学生的 5 门成绩. 定义一个一维数组 B,用于 ...
已知abc+cba=1333，其中abc均为一位数，求出符合条件的abc的值
已知abc+cba=1333,其中abc均为一位数,求出符合条件的abc的值 #include<stdio.h> int main(){int a,b,c;for(a=1;a<=9; ...
python怎么根据点来拟合曲线_2019_nCoV_利用python根据已知点求拟合曲线及简单预测(无实际意义)...
前言本文仅做根据已知点求拟合曲线的几种方法的python实现,无任何实际意义数据来源(另一篇博文) 利用Python爬取新冠肺炎疫情实时数据,Pyecharts画2019-nCoV疫情地图参考 ...
【Python】已知一张图片中的框图坐标，切割出目标框图（单个）
Target:目标检测已知框的坐标,将框中的图像从原图片中分割出来做了一下午都要做自闭了,到晚上终于切割出了想要的那部分图片(我是真的菜,下午直接反省了一遍自己的大学生活QAQ) 话不多说,大家肯定 ...

Python实现：已知化学分子的输入文件坐标（高斯计算输入文件为例），求其中任意三个原子确定的平面的法向量和单位法向量

Python实现：已知化学分子的输入文件坐标（高斯计算输入文件为例），求其中任意三个原子确定的平面的法向量和单位法向量相关推荐

最新文章

热门文章