计算化学公社

 找回密码 Forget password
 注册 Register
Views: 974|回复 Reply: 5

[Gaussian/gview] 基于python生成蛋白质-配体簇模型刚性扫描的gaussian输入文件[班门弄斧]

[复制链接 Copy URL]

47

帖子

3

威望

982

eV
积分
1089

Level 4 (黑子)

发表于 Post on 2022-8-9 19:33:18 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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库:
  • numpy
  • pandas
  • math

一些解释:
  • 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:
  • 自动判断是否使用readopt功能
  • ...

最后的最后:这只是新手小白的一个小小尝试,代码可能写的冗杂,如有错误,欢迎批评指正
  1. import numpy as np
  2. import pandas as pd
  3. import math


  4. # 初始参数设置
  5. fn1="1.xyz"
  6. fn2="2.xyz"
  7. atomidx1=163
  8. atomidx2=1
  9. stepsize=0.05
  10. step=50
  11. charge=0
  12. spin=1
  13. nproc=32
  14. mem="12GB"
  15. method="b3lyp"
  16. basis="6-31g(d,p)"
  17. #fixatoms="1,2,16,17,33,34,52,71,78,86,97,111,113,129,130,156,145"
  18. other=""
  19. #------------------------------------------------------------------------------#
  20. with open(fn1) as f1:
  21.     m1=f1.readlines()
  22. xyz1=list(filter(None,[ l.split() for l in m1[2:] if not l.startswith('\n') ]))
  23. with open(fn2) as f2:
  24.     m2=f2.readlines()
  25. xyz2=list(filter(None,[ l.split() for l in m2[2:] if not l.startswith('\n') ]))
  26. df1=pd.DataFrame(xyz1,columns=["atom","x","y","z"])
  27. df2=pd.DataFrame(xyz2,columns=["atom","x","y","z"])
  28. data1=df1[["x","y","z"]].values.astype(np.float64)
  29. data2=df2[["x","y","z"]].values.astype(np.float64)
  30. vec=data1[atomidx1-1]-data2[atomidx2-1]
  31. proc_data=data2
  32. for i in range(1,step+1):
  33.     proc_data=np.append(proc_data,data2+vec*stepsize*i/math.sqrt(np.square(vec).sum()),axis=0)  
  34. mm=pd.DataFrame(proc_data,columns=["x","y","z"])
  35. mm.insert(0,"atom",np.tile(df2["atom"].values,step+1))
  36. # -------------------------------------------------------------------------------- #
  37. # 输出轨迹
  38. fp=open("traj.xyz","w")
  39. for m in range(step+1):
  40.     print("{}".format(data1.shape[0]+data2.shape[0]),file=fp)
  41.     print("{}".format(m*stepsize),file=fp)
  42.     for i in range(data1.shape[0]):
  43.         print ("{}     {}  {} {}".format(df1.iloc[i]["atom"],df1.iloc[i]["x"],df1.iloc[i]["y"],df1.iloc[i]["z"]),file=fp)
  44.     for j in range(data2.shape[0]):
  45.         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)
  46. fp.close()
  47. #-------------------------------------------------------------------------------#
  48. # 输出gaussian文件
  49. fg=open("traj.gjf","w")
  50. for m in range(step+1):
  51.     chkname="%chk="+"traj_"+str(m)+".chk"
  52.     print("{}".format(chkname),file=fg)
  53.     print("{}".format("%nproc="+str(nproc)),file=fg)
  54.     print("{}".format("%mem="+str(mem)),file=fg)
  55.     print("{}".format("#p "+method+"/"+basis+" "+other),file=fg)
  56.     print("{}".format("            "),file=fg)
  57.     print("{}".format("traj_"+str(m)+" Structure"),file=fg)
  58.     print("{}".format("            "),file=fg)
  59.     print("{}".format(str(charge)+" "+str(spin)),file=fg)
  60.     for i in range(data1.shape[0]):
  61.         print ("{}     {}  {} {}".format(df1.iloc[i]["atom"],df1.iloc[i]["x"],df1.iloc[i]["y"],df1.iloc[i]["z"]),file=fg)
  62.     for j in range(data2.shape[0]):
  63.         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)
  64.     print("{}".format("            "),file=fg)
  65.     #print("{}".format("{}".format("notatoms="+fixatoms)),file=fg)
  66.     print("{}".format("            "),file=fg)
  67.     if m!=step:
  68.         print("{}".format("--link1--"),file=fg)
  69.     else:
  70.         continue
  71. fg.close()
复制代码
录制_2022_08_09.gif convert.rar (61.33 KB, 下载次数 Times of downloads: 22)

评分 Rate

参与人数
Participants 2
eV +15 收起 理由
Reason
yjb + 5 好物!
sobereva + 10

查看全部评分 View all ratings

347

帖子

0

威望

1387

eV
积分
1734

Level 5 (御坂)

发表于 Post on 2022-8-9 21:39:21 | 显示全部楼层 Show all
好奇你的 GIF 动图是如何上传的。

我之前上传的时候传不上去,文件也不算太大(1.6 MB)……
Sob 老师还说即使上传 GIF 了也会添加水印,就不能显示运动的效果……

47

帖子

3

威望

982

eV
积分
1089

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-8-9 21:46:56 | 显示全部楼层 Show all
乐平 发表于 2022-8-9 21:39
好奇你的 GIF 动图是如何上传的。

我之前上传的时候传不上去,文件也不算太大(1.6 MB)……

就很正常的上传,可能是因为大小的问题。gif太大,服务器响应不过来

4万

帖子

99

威望

4万

eV
积分
89976

管理员

公社社长+计算化学玩家

发表于 Post on 2022-8-10 07:03:57 | 显示全部楼层 Show all
乐平 发表于 2022-8-9 21:39
好奇你的 GIF 动图是如何上传的。

我之前上传的时候传不上去,文件也不算太大(1.6 MB)……

目前论坛对gif的限制是600KB
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

347

帖子

0

威望

1387

eV
积分
1734

Level 5 (御坂)

发表于 Post on 2022-8-10 16:05:26 | 显示全部楼层 Show all
sobereva 发表于 2022-8-10 07:03
目前论坛对gif的限制是600KB

了解,谢谢Sob老师

112

帖子

0

威望

3336

eV
积分
3448

Level 5 (御坂)

发表于 Post on 2022-9-23 03:03:57 | 显示全部楼层 Show all
Sob老师那篇博文的刚性扫描都是线性体系或者小分子体系。酶结合底物的地方从RC到TS的扫描能用刚性扫描吗?随着成键体系反应坐标的距离越来近(或者越来越远,断键的体系),底物的其他部分大概率还是会和结合口袋的残基产生重合,重合的这部分会引起势能的急剧增高,这必定会影响TS坐标的确定啊

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 04:30 , Processed in 0.646030 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list