|
|
本帖最后由 casea 于 2022-8-11 09:19 编辑
**2022/8/11更新：附件中提供了gui和命令行版本的py文件**
最近在学习复现一篇使用簇模型研究转氨酶反应机理的文献https://pubs.rsc.org/en/content/articlelanding/2015/ob/c5ob00690b。
该文章将活性区域抠出来，同时对一些骨架碳原子进行fix，从而保持蛋白的结构。该文章优化时使用的是`B3LYP/6-31G(d,p)`。在优化初始结构后，就需要根据反应的机理寻找其中间体和过渡态。
因为文章已经分析的就很详细了，在此就不再赘述。寻找过渡态的时候就需要进行扫描。在拜读了社长的文章http://sobereva.com/474之后，决定使用刚性扫描。我的这个体系因为是蛋白-配体体系，所以高斯自带的scan是无法实现刚性扫描功能。但万能的社长给我们提供了解决方法，具体方法见上方引用链接。
本来，问题到这就都可以解决了。但是突然心血来潮之下，想用python写一个简单的程序实现。所用的文件见附件。
采用的方法特别简单，就是先计算两个原子之间的方向向量，然后根据步长和步数生成一些列的初始结构，同时生成gaussian输入文件
使用的python库：
一些解释：
- fn1，fn2分别代表着蛋白和配体的xyz轨迹
- atomidx1，2分别代表着蛋白和配体上选定的原子index(注意：此处的原子需要都是从1开始）
- stepsize：步长也就是两个原子之间直线距离
- step：运行多少步
- charge：总体系电荷
- spin: 总体系自旋多重度
- nproc：gaussian计算所用核数
- mem：分配内存
- method：使用方法
- basis:基组
- fixatom:在计算中需要固定的原子序号(用于维持蛋白骨架）
- other:可以添加一些额外的计算选项
文件最后会输出两个文件：
- traj.xyz 输出的轨迹文件
- traj.gjf 输出的Gaussian输入文件
To do：最后的最后：这只是新手小白的一个小小尝试，代码可能写的冗杂，如有错误，欢迎批评指正- import numpy as np
- import pandas as pd
- import math
- # 初始参数设置
- fn1="1.xyz"
- fn2="2.xyz"
- atomidx1=163
- atomidx2=1
- stepsize=0.05
- step=50
- charge=0
- spin=1
- nproc=32
- mem="12GB"
- method="b3lyp"
- basis="6-31g(d,p)"
- #fixatoms="1,2,16,17,33,34,52,71,78,86,97,111,113,129,130,156,145"
- other=""
- #------------------------------------------------------------------------------#
- with open(fn1) as f1:
- m1=f1.readlines()
- xyz1=list(filter(None,[ l.split() for l in m1[2:] if not l.startswith('\n') ]))
- with open(fn2) as f2:
- m2=f2.readlines()
- xyz2=list(filter(None,[ l.split() for l in m2[2:] if not l.startswith('\n') ]))
- df1=pd.DataFrame(xyz1,columns=["atom","x","y","z"])
- df2=pd.DataFrame(xyz2,columns=["atom","x","y","z"])
- data1=df1[["x","y","z"]].values.astype(np.float64)
- data2=df2[["x","y","z"]].values.astype(np.float64)
- vec=data1[atomidx1-1]-data2[atomidx2-1]
- proc_data=data2
- for i in range(1,step+1):
- proc_data=np.append(proc_data,data2+vec*stepsize*i/math.sqrt(np.square(vec).sum()),axis=0)
- mm=pd.DataFrame(proc_data,columns=["x","y","z"])
- mm.insert(0,"atom",np.tile(df2["atom"].values,step+1))
- # -------------------------------------------------------------------------------- #
- # 输出轨迹
- fp=open("traj.xyz","w")
- for m in range(step+1):
- print("{}".format(data1.shape[0]+data2.shape[0]),file=fp)
- print("{}".format(m*stepsize),file=fp)
- for i in range(data1.shape[0]):
- print ("{} {} {} {}".format(df1.iloc[i]["atom"],df1.iloc[i]["x"],df1.iloc[i]["y"],df1.iloc[i]["z"]),file=fp)
- for j in range(data2.shape[0]):
- print ("{} {:.8f} {:.8f} {:.8f}".format(mm.iloc[(m*(data2.shape[0])+j)]["atom"],mm.iloc[(m*(data2.shape[0])+j)]["x"],mm.iloc[(m*(data2.shape[0])+j)]["y"],mm.iloc[(m*(data2.shape[0])+j)]["z"]),file=fp)
- fp.close()
- #-------------------------------------------------------------------------------#
- # 输出gaussian文件
- fg=open("traj.gjf","w")
- for m in range(step+1):
- chkname="%chk="+"traj_"+str(m)+".chk"
- print("{}".format(chkname),file=fg)
- print("{}".format("%nproc="+str(nproc)),file=fg)
- print("{}".format("%mem="+str(mem)),file=fg)
- print("{}".format("#p "+method+"/"+basis+" "+other),file=fg)
- print("{}".format(" "),file=fg)
- print("{}".format("traj_"+str(m)+" Structure"),file=fg)
- print("{}".format(" "),file=fg)
- print("{}".format(str(charge)+" "+str(spin)),file=fg)
- for i in range(data1.shape[0]):
- print ("{} {} {} {}".format(df1.iloc[i]["atom"],df1.iloc[i]["x"],df1.iloc[i]["y"],df1.iloc[i]["z"]),file=fg)
- for j in range(data2.shape[0]):
- print ("{} {:.8f} {:.8f} {:.8f}".format(mm.iloc[(m*(data2.shape[0])+j)]["atom"],mm.iloc[(m*(data2.shape[0])+j)]["x"],mm.iloc[(m*(data2.shape[0])+j)]["y"],mm.iloc[(m*(data2.shape[0])+j)]["z"]),file=fg)
- print("{}".format(" "),file=fg)
- #print("{}".format("{}".format("notatoms="+fixatoms)),file=fg)
- print("{}".format(" "),file=fg)
- if m!=step:
- print("{}".format("--link1--"),file=fg)
- else:
- continue
- fg.close()
convert.rar
(61.33 KB, 下载次数 Times of downloads: 22)
|
评分 Rate
-
查看全部评分 View all ratings
|
发表于 Post on 2022-8-9 19:33:18
楼主 Author
发表于 Post on 2022-8-10 07:03:57