Python实现：超分子化学的建模------如何操控客体分子穿过主体分子和计算该过程能量变化（高斯(Gauss)输入文件为例，一键批量处理）
2024-05-28 12:19:36
注：该Python代码可以实现所有高斯（Gauss）计算输入文件的处理（特别是针对超分子体系的路径建模非常有用，也可以作为处理其它体系的参考）。
正文如下：
超分子建模过程中，如果需要模拟客体分子穿过主体分子孔洞的过程（如小分子穿过大环分子）的能量变化，或者相互作用力的变化，具体实例如图一所示，球状富勒烯C60穿过一个分子环的过程。
图一
假设大环是由10个苯环通过σ键连接而成，形成[10]CPP，中文名称10环对苯撑，富勒烯C60分子作为客体分子，与大环会存在相互作用力，显而易见，在尺寸上这两者非常匹配，类似于在富勒烯球外面套上一根腰带，而这类分子是典型的超分子结构，理论计算对于研究这类体系尤为重要，特别是研究主体分子和客体分子间的相互作用力，扯远了，主客体相互作用不是本帖子讨论的重点。。。。。
回归正题，实现计算客体分子穿过主体分子这一过程的能量变化，设计这个模型的逻辑如下：
1. 总体思路是，富勒烯从一端到另一端的路径上设置多个点，可密可稀，对路径上的每一个点对应的超分子结构进行平面内的优化，求单点能做成二维图即可（也可以是刚性，默认中心的每一个点就是在该平面上的能量最低点，这样不用再优化结构，大大降低计算耗时）。
2. 这个路径是在[10]CPP大环的环平面对应垂直方向上，也就是说需要确定一个[10]CPP的环平面，可以默认为苯环之间的碳碳σ键均处在大平面上，也就是说，只要在这个位置找三个点就可以确定大环的化学平面（如C13-C4-C43三个碳原子对应的坐标点），再求三个坐标形成平面的法向量，再将得到的法向量转化为单位向量，以方便后续对坐标进行操作。
3. 在给定坐标基础上，只需要将富勒烯分子所有的坐标值在这个法向量上进行数学加减，即可代表在这个方向上的移动，移动的位置通过实际距离长度乘以单位法向量，即可得实际移动距离的具体坐标。
4. 通过距离批量命名输入文件名，得到对应高斯计算的输入文件。
5. 通过Gauss软件(主流是g09或者g16版本)计算得到每一个点能量，提取能量作图（非本文讨论内容）。
具体实现如下：
球形富勒烯C60穿过大环的过程的模拟，第一步是先得到优化后的超分子结构，如图二所示：
图二
用Gaussview打开该分子，存储为高斯输入文件（后缀为.gjf文件），这里需要注意，坐标顺序必须保持主体分子在前面，客体分子在后面，如图二中的数字编号，[10]CPP在前，C60在后。具体坐标如下所示：
%nprocshared=12
%mem=2500mb
%rwf=1,2000mb,2,2000mb,3,2000mb,4,2000mb,5,2000mb,6,2000mb,7,2000mb,8,2000mb,9,-1
%chk=xxxx.chk
# b3lyp/3-21g* opt c60-cycloparaphenylene0 1C 1.73732300 -6.65945200 0.02809900 C 1.02755200 -7.09867900 -1.09886500 C -0.35923100 -7.15166600 -1.09301500 C -1.09214200 -6.76784500 0.03963400 C -0.37754700 -6.49784200 1.21181700 C 1.01002500 -6.44545800 1.20553100 C -2.51785300 -6.37048700 -0.04431700 C -2.97016100 -5.74036500 -1.21000000 C -4.13845500 -4.99392400 -1.21062100 C -4.90312800 -4.84814500 -0.04645400 C -4.53591800 -5.60946900 1.07118000 C -3.36411000 -6.35602500 1.07264900 C -5.85139600 -3.71116200 0.03786900 C -5.90440600 -2.95677700 1.21537900 C -6.40004400 -1.65956700 1.21155600 C -6.85703100 -1.06756500 0.02933800 C -6.97558800 -1.88698300 -1.10270300 C -6.48354100 -3.18455800 -1.09801100 C -6.90248000 0.41140400 -0.06167900 C -6.40817300 1.02325000 -1.21997500 C -6.03985300 2.36013100 -1.21989300 C -6.14496300 3.13798100 -0.06066500 C -6.78453300 2.56780500 1.04839600 C -7.15754000 1.22949300 1.04716500 C -4.64005800 4.65248500 1.21041300 C -3.55537200 5.52027200 1.21350300 C -3.12994100 6.14885600 0.03711600 C -3.94425100 6.01574000 -1.09686700 C -5.02692900 5.14797700 -1.10110900 C -1.74018800 6.65969700 -0.04980100 C -1.00956300 6.44183900 -1.22441800 C 0.37809400 6.49263500 -1.22623300 C 1.08927100 6.76636900 -0.05284800 C 0.35328500 7.15304900 1.07679200 C -1.03349400 7.10088000 1.07840200 C 2.51547300 6.37190200 0.03499900 C 2.96709500 5.74379200 1.20198200 C 4.13620000 4.99840600 1.20469200 C 4.90210900 4.85217600 0.04139700 C 4.53622800 5.61330700 -1.07684800 C 3.36365400 6.35842800 -1.08049400 C -5.34104400 4.38075600 0.03042300 C 5.84993700 3.71487500 -0.04318300 C 5.90485500 2.96328900 -1.22243000 C 6.40060600 1.66622900 -1.22114400 C 6.85599400 1.07116200 -0.03981500 C 6.97232200 1.88772500 1.09456900 C 6.47981400 3.18517000 1.09247600 C 6.90179400 -0.40810500 0.04702900 C 6.41044100 -1.02363400 1.20472100 C 6.04105800 -2.36014200 1.20121100 C 6.14065700 -3.13376900 0.03855600 C 6.77748100 -2.56040100 -1.07033100 C 7.15289300 -1.22272400 -1.06519000 C 5.33527100 -4.37545100 -0.05446300 C 5.02398600 -5.14574800 1.07577100 C 3.94233000 -6.01463200 1.07152800 C 3.12606900 -6.14628300 -0.06125000 C 3.54840200 -5.51412100 -1.23685100 C 4.63187500 -4.64464500 -1.23359300 H 1.56541000 -7.34672800 -2.00967900 H -0.88545100 -7.43854900 -1.99928800 H -0.91060600 -6.18229600 2.10468700 H 1.52018700 -6.08573800 2.09444000 H -2.33707200 -5.72521100 -2.09203900 H -4.38431200 -4.40772900 -2.09088400 H -5.14544900 -5.57352400 1.96995500 H -3.07485500 -6.89299300 1.97196600 H -5.41551900 -3.32690100 2.11219600 H -6.28727900 -1.05339100 2.10620900 H -7.39935800 -1.48118400 -2.01720700 H -6.53116200 -3.77542700 -2.00849700 H -6.17975000 0.41790900 -2.09207500 H -5.53838000 2.76614300 -2.09331200 H -6.94304100 3.16412200 1.94303300 H -7.60200200 0.79948200 1.94069200 H -4.84870400 4.07052900 2.10418900 H -2.94650100 5.58820300 2.11037400 H -3.69015100 6.55329400 -2.00614400 H -5.60153600 5.01837500 -2.01422100 H -1.51756600 6.08167300 -2.11440800 H 0.91392000 6.17460700 -2.11655600 H 0.87701800 7.44229900 1.98377200 H -1.57408600 7.35169700 1.98681500 H 2.33260300 5.72835100 2.08305300 H 4.38139700 4.41318100 2.08579100 H 5.14722600 5.57782700 -1.97462900 H 3.07511300 6.89474600 -1.98042300 H 5.41762500 3.33561100 -2.11920300 H 6.28935500 1.06262500 -2.11768600 H 7.39476000 1.47994100 2.00878700 H 6.52580000 3.77355800 2.00465100 H 6.18445400 -0.42114700 2.07940700 H 5.54299500 -2.76899400 2.07524500 H 6.93213400 -3.15342900 -1.96782700 H 7.59563600 -0.79052700 -1.95849300 H 5.60050600 -5.01828500 1.98796400 H 3.69078200 -6.55483000 1.97993100 H 2.93823900 -5.58114600 -2.13297100 H 4.83796800 -4.06014100 -2.12631500 C 3.08629400 1.24678300 -1.18663000 C 2.65597700 2.32370100 -0.31319800 C 1.56703000 3.11047900 -0.66430900 C 0.85900600 2.85163200 -1.90693600 C 1.27520600 1.82434600 -2.74314300 C 2.41654400 1.00552200 -2.37623000 C -0.55709300 3.06239900 -1.67116800 C -1.49567400 2.24295600 -2.28288400 C -2.64937100 1.77240700 -1.53568200 C -2.81079900 2.14125400 -0.20797900 C -1.82534800 2.99500300 0.43065900 C -0.72662300 3.44561000 -0.28147900 C -3.26136200 1.16590200 0.76842500 C -2.55086700 1.41863900 2.01139700 C -2.13865100 0.35711200 2.80514200 C -2.41373400 -1.00859700 2.39275400 C -3.08356400 -1.24993700 1.20322000 C -3.52064900 -0.13964000 0.37401000 C -1.27238900 -1.82738300 2.75972000 C -0.85626900 -2.85461400 1.92351300 C 0.55986500 -3.06524800 1.68783300 C 1.49858600 -2.24587200 2.29955800 C 1.06222000 -1.17289700 3.17262700 C -0.29223400 -0.96870100 3.39765200 C 2.92698600 -0.41057000 1.96370400 C 3.35478400 0.52527600 1.03081100 C 3.52319500 0.13647400 -0.35745300 C 3.26430100 -1.16913700 -0.75185200 C 2.81349900 -2.14436900 0.22468200 C 2.65226100 -1.77547900 1.55242000 C 2.14154100 -0.36018800 -2.78858800 C 2.55380600 -1.42178700 -1.99480700 C 1.82818700 -2.99805800 -0.41386800 C 0.72937800 -3.44822400 0.29824800 C -1.56436700 -3.11370700 0.68080900 C -2.65326800 -2.32684700 0.32968100 C -3.35198300 -0.52834000 -1.01420900 C -2.92402900 0.40755000 -1.94701700 C -1.05929900 1.16995800 -3.15594000 C 0.29511600 0.96568500 -3.38110800 C 1.94365300 -0.03842200 2.96525300 C -0.82734600 0.38054400 3.42607800 C -1.66690800 2.55081600 1.80312300 C 0.58416400 3.47687300 0.34014400 C 2.81871800 1.87494100 1.05829200 C 0.83027500 -0.38358800 -3.40955200 C -1.94069100 0.03545200 -2.94851500 C -2.81592100 -1.87796100 -1.04170800 C -0.58141700 -3.47979100 -0.32355700 C 1.66973200 -2.55380600 -1.78644100 C -1.87700800 -2.23334000 -2.00088800 C -1.42953700 -1.25416300 -2.97509800 C -0.73573500 -3.05340800 -1.63463700 C 0.41602800 -2.57836600 -2.38169300 C -0.01295600 -1.46817700 -3.21040300 C -0.41306600 2.57533700 2.39827700 C 0.73858100 3.05044500 1.65128300 C 1.87994000 2.23039500 2.01752600 C 1.43248600 1.25119500 2.99177500 C 0.01591200 1.46514200 3.22701500
最终python代码如下，复制粘贴进pycharm即可，使用前需要安装numpy、linecache库，安装方法请自行google或者baidu，简单到令人发指，time库和decimal库python已自带，直接调用即可。
注意：高斯输入文件需要放入代码所在的同一文件夹下，否则不能执行：
"""
通过三个化学点求法向量和单位法向量,并实现客体分子从主体分子平面的垂直方向穿过，并并批量产生该路径上的高斯计算的输入文件
"""
import time
import numpy
import linecache
from decimal import Decimal# 将向量转化为单位向量，如果向量是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# 小数的区间变化函数定义
def frange(x, y, jump):while Decimal(str(x)) < Decimal(str(y)):yield xx += float(Decimal(str(jump)))filename = input("请输入文件名：")
atoms_num_all = int(input("请输入原子数："))
atoms_num_plane = int(input("请输入主体分子原子数："))
atom_1, atom_2, atom_3 = eval(input("三个原子确定的化学平面(请输入三个原子序号 e.g. 1,2,3)："))
initial_pos = float(input("请输入起始位置（e.g. -5，表示距离平面-5.0 Å）："))
final_pos = float(input("请输入终点位置（e.g. 5，表示距离平面5.0 Å）："))
CalGrid = float(input("请输入计算格点密度（e.g. 0.5表示每0.5 Å取一个点）："))# 记录开始运行程序的时间
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
# 这里获取第一个坐标
coord1_line = linecache.getline(f'{filename}', atom_1+n)
coord1_line = coord1_line.split( )
coord_name1, *coord_temp = coord1_line
coord_1 = [float(x) for x in coord_temp]
print(f"第一个原子是：{coord_name1}, 坐标是：", coord_1)
# 这里获取第二个坐标
coord2_line = linecache.getline(f'{filename}', atom_2+n)
coord2_line = coord2_line.split( )
coord_name2, *coord_temp = coord2_line
coord_2 = [float(x) for x in coord_temp]
print(f"第二个原子是：{coord_name2}, 坐标是：", coord_2)
# 这里获取第三个坐标
coord3_line = linecache.getline(f'{filename}', atom_3+n)
coord3_line = coord3_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}")for i_distance in frange(initial_pos, final_pos, CalGrid):i_distance_decimal = Decimal(i_distance).quantize(Decimal('0.00'))# 读取文件每一行数据with open(f"{i_distance_decimal}.gjf", mode = "w", encoding = "utf-8") as f:for m in range(1, 10):mark_line = linecache.getline(f'{filename}', m)mark_x = mark_line[0]if mark_x == '#':n = m + 4breakfor x in range(1,n + 1):data_line = linecache.getline(f'{filename}', x)f.write(f"{data_line}")for y in range(n + 1,n + atoms_num_plane + 1):data_line = linecache.getline(f'{filename}', y)f.write(f"{data_line}")for z in range(n + atoms_num_plane+1,n + atoms_num_all+1):coord_line = linecache.getline(f'{filename}', z)coord_line = coord_line.split( )coord_name, *coord_temp = coord_linecoord = [float(p) for p in coord_temp]a = numpy.array(coord)b = numpy.array(NormVect_val)x = float(str(Decimal(i_distance).quantize(Decimal('0.00'))))c = a+b*xd = [str(i) for i in c]e = ' '.join(d)f.write(f"{coord_name} {e}\n")f.write("\n")# 记录结束运行程序的时间
end = time.time()
print(end - start)上面代码运行后，按照下面提示，输入相应内容即可：
至此输入文件制作完成，将得到的高斯输入文件传入高斯软件进行计算，最后统计得到的能量，origin做二维图即可展现该路径上的能量变化。
Python实现：超分子化学的建模------如何操控客体分子穿过主体分子和计算该过程能量变化（高斯(Gauss)输入文件为例，一键批量处理）相关推荐
- python参数化建模加工图_基于Python的ABAQUS层压板参数化建模
唐维 康泽毓 杨婷 曾凤 蒋莉 摘要:为了提高层压板在ABAQUS仿真中建模的效率与准确性,提出利用Python语言对ABAQUS二次开发进行层压板参数化建模的方法.基于ABAQUS有限元软件,采用P ...
- Python实现：已知化学分子的输入文件坐标（高斯计算输入文件为例），求其中任意三个原子确定的平面的法向量和单位法向量
计算化学在处理实际化学问题时,比如需要在某一化学平面的法向量上进行分子操作,这时最重要的是确定化学平面和求法向量,才能进行后续的操作(如下图所示),下面以高斯输入文件为例,用python代码实现该功能 ...
- python 白噪声检验-利用python实现平稳时间序列的建模方式
假如某个观察值序列通过序列预处理可以判定为平稳非白噪声序列,就可以利用ARMA模型对该序列进行建模.建模的基本步骤如下: (1)求出该观察值序列的样本自相关系数(ACF)和样本偏自相关系数(PACF) ...
- 《Python金融大数据风控建模实战》 第6章 变量分箱方法
<Python金融大数据风控建模实战> 第6章 变量分箱方法 本章引言 Python代码实现及注释 本章引言 变量分箱是一种特征工程方法,意在增强变量的可解释性与预测能力.变量分箱方法主要 ...
- 时间序列分析 | Python实现时间序列统计分析与建模
时间序列分析 | Python实现时间序列统计分析与建模 目录 时间序列分析 | Python实现时间序列统计分析与建模 基本介绍 数据准备 程序设计 学习总结 参考资料 基本介绍 本文介绍时间序列P ...
- python实现gauss-seidel迭代公式_python实现高斯(Gauss)迭代法的例子
python实现高斯(Gauss)迭代法的例子 我就废话不多说了,直接上代码大家一起看吧! #Gauss迭代法 输入系数矩阵mx.值矩阵mr.迭代次数n(以list模拟矩阵 行优先) def Gaus ...
- python使用statsmodels包中的robust.mad函数以及pandas的apply函数计算dataframe中所有数据列的中位数绝对偏差(MAD)
python使用statsmodels包中的robust.mad函数以及pandas的apply函数计算dataframe中所有数据列的中位数绝对偏差(MAD.Median Absolute Devi ...
- python使用matplotlib可视化跨年数值指标中位数变化率、使用pct_change函数计算变化率、年环比变化率（pct_change function）
python使用matplotlib可视化跨年数值指标中位数变化率.使用pct_change函数计算变化率.年环比变化率(visualize year change rate with pct_cha ...
- Python之pyecharts：利用pyecharts绘制2020年11月16日微博话题热度排行榜实时变化
Python之pyecharts:利用pyecharts绘制2020年11月16日微博话题热度排行榜实时变化 目录 利用pyecharts绘制2020年11月16日微博话题热度排行榜实时变化 Bar( ...
- CV：计算机视觉技术之图像基础知识(一)—以python的cv2库来了解计算机视觉图像基础(傅里叶变换-频域-时域/各种滤波器-线性-非线性-均值-中值-高斯-双边)
CV:计算机视觉技术之图像基础知识(一)-以python的cv2库来了解计算机视觉图像基础(傅里叶变换-频域-时域/各种滤波器-线性-非线性-均值-中值-高斯-双边) 目录 一.图像中的傅里叶变换 1 ...
最新文章
- jvm 内存结构默写
- Ubuntu 设置NAT共享网络（命令行方法）
- Apache Flink 简介和编程模型
- 算法题：水洼有多少（C++）
- 最短路计数(spfa)
- 合并k个有序链表 python_[LeetCode] 23. Merge k Sorted Lists 合并k个有序链表
- 【机器视觉】 repeat算子
- 游戏设计与计算机,RPG游戏设计与实现-数学与计算机系.doc
- linux 串口 lsr 0xc9,串口发送0x0D后，从串口接收到数据被转换成了0x0A
- maven的Windows环境下安装配置
- get post put delete在vue中传参方式
- python实现关键词提取
- django 1366, “Incorrect string value: for column ‘‘ at row
- 使用Android Studio 开发APP入门经验
- available: expected at least 1 bean which qualifies as autowire candidate
- T470P笔记本安装固态以及固态中安装系统_完整步骤
- 线程有哪些状态？创建、就绪、运行、阻塞和死亡
- 还在相信男女之间真的有纯友谊？太傻太天真！
- 已知某校有以下老师及教授课程，1） 使用一个Map，以老师的名字作为键，以老师教授的课程名作为值，表示上述 课程安排。
- iso转cue mac_mac如何播放cue文件？