请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 3539|回复: 27

[综合交流] gentor:扫描方式做分子构象搜索的便捷工具

[复制链接]

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

发表于 2015-12-18 04:50:50 | 显示全部楼层 |阅读模式
gentor:扫描方式做分子构象搜索的便捷工具

文/Sobereva(3)
First release: 2015-Dec-18    Last update: 2016-Aug-16



1 前言

笔者曾开发过molclus程序用于团簇构型搜索和分子构象搜索,不熟悉的话应当先看此文《使用molclus程序做团簇构型搜索和分子构象搜索》(http://bbs.keinsci.com/forum.php?mod=viewthread&tid=577)。molclus简单来说就是读取含有多帧结构信息的traj.xyz文件,依次对每个结构调用Gaussian或MOPAC或ORCA进行优化,然后对结果用isostat子程序进行统计分析。

对很多人来说,这种通过动力学方式搜索分子构象有一个不便之处就是需要通过第三方程序分子动力学程序来产生一批初始构象,很多人不会操作。现在介绍笔者开发的gentor程序来让molclus能以扫描方式做构象搜索。gentor可以对指定的键按一定规则扭转产生一批分子初始构象,然后写入traj.xyz文件中,再使用molclus调用量化程序对之依次优化后,就能得到各个极小点结构和能量最低结构。

gentor从molclus 1.3版开始作为molclus的一个组件发布。获得gentor的方法就是到molclus主页(http://www.keinsci.com/research/molclus.html)下载最新版本molclus程序,预编译的gentor就在其中。molclus程序的免费版有原子数限制,但附带的gentor没有原子数限制。

gentor这种势能面扫描方式产生初始构象相比通过动力学来产生有优有劣。优点是省得跑动力学了,方便很多,而且对势能面采样得较为均匀全面,不容易漏掉极小点。缺点就是,产生的初猜构象的数目随着可旋转的键的增加呈几何级数增长。比如每个可旋转的键扫3个点,考虑5个键就得优化243次,有10个键就得优化高达59049次,此时已经几乎没法处理了。所以对于小分子用势能面扫描很好,但体系大了的话还是做动力学产生初猜结构为宜,因为这样得到的初猜结构大多分布在感兴趣的低能构象附近,最终找到能量最低构象的可能性较大。如果用退火方式则更好,比如gromacs中可以使用周期退火,每100ps内从0K升温到500K再返回0ps,跑50ns就能得到500个处于0K的构象,用molclus对它们优化,得到能量最低结构的几率就会很大了。


2 使用方法

gentor的使用非常容易。需要在当前目录下准备两个文件:
mol.xyz:记录分子结构的.xyz文件
gentor.ini:设定要旋转的键以及旋转方式

下面是个gentor.ini的例子,考虑了4个键的旋转,囊括了各种gentor支持的键旋转的定义方式:
5-4    //绕着5-4键令4的那一头的片段进行旋转
0 180 270    //旋转0度、180度、270度
4-6    //绕着4-6键令6的那一头的片段进行旋转
e60    //每60度旋转一次直到360度,即0,60,120,180,240,300
7-8    //绕着7-8键令8的那一头的片段进行旋转
-120 120    //两个方向各旋转120度
9-12    //绕着9-12键令12的那一头的片段进行旋转
r8    //旋转0~360度内的随机8个角度

可以定义无数多个键的旋转,旋转是相对于mol.xyz里的结构进行的。上面三个键分别旋转3、6、2、8次,因此一共会产生288个构象,并输出到当前目录下的traj.xyz里。然后就可以直接运行molclus来自动调用量化程序去优化其中的每个结构,然后用isostat进行统计了。

通过旋转键产生的构象很可能会有不合理接触,gentor会自动进行检测,如果结构中有任意两个原子间距离小于它们的共价半径和的一半就视为结构不合理而被刨去,不输出到traj.xyz当中。

默认情况下,原子间距离小于它们的共价半径和的1.15倍就被视为成键。但有的时候gentor自动判断的成键方式和我们预期的不符,导致基团没法按照预期方式转动。这种情况下可以在定义键的旋转方式后面空一行自己定义成键判断方式。例如
5-4
0 180 270
4-6
e60

critbond=1.2
4 8 1
4 9 1
5 11 0

以上内容就代表,如果原子间距离小于它们的共价半径和的1.2倍就被视为成键,并且人为加上4-8、4-9键(因为最后一个数字写的是1),而去掉5-11键(因为最后一个数字写的是0)。gentor输出信息中的Recognized bond是最终的原子间成键关系。


3 结合molclus做分子构象搜索实例:多巴胺

这里我们用一个很简单的分子多巴胺作为实例展示gentor+molclus做构象搜索的基本流程。分子初始结构如下:



1.png

对应的mol.xyz文件内容如下所示,将此文件放到molclus目录下
22
dopamine
O     -2.33158827    1.90324395    0.05204669
O     -3.45288968   -0.47124460    0.37321012
N      4.33709192    0.16551438    0.29362167
C      2.10661405   -0.10993743   -0.67966314
C      0.62079598   -0.23561182   -0.42292440
C      2.91795648    0.08557132    0.60357302
C     -0.18237422    0.90195216   -0.31692881
C      0.02151120   -1.46774074   -0.25587568
C     -1.53226977    0.81838896   -0.05200952
C     -1.34075114   -1.56264519    0.01270596
C     -2.12492987   -0.43572743    0.11734436
H      2.29355697    0.72769638   -1.34826040
H      2.47698171   -0.99578880   -1.18305280
H      2.54052562    0.95777703    1.14088437
H      2.76403205   -0.77011970    1.25210418
H      0.25802965    1.87739170   -0.45049059
H      0.60483520   -2.36786844   -0.33618849
H     -1.79443393   -2.53198365    0.13749855
H      4.88321312    0.19170936    1.13290303
H      4.54272968    1.00577810   -0.21284385
H     -1.82929975    2.69323340   -0.08852494
H     -3.74330646   -1.36791983    0.46123396

此体系我们应当扫描4-6和6-3两个键,每步120度。5-4键虽然可以旋转,但旋转它没意义,因为旋转前后是相对于环平面对称的,不会出现新的极小点结构。两个羟基的旋转这里就不考虑了,有兴趣可以也考虑进去。

按照以上讨论,应把gentor.ini设为以下内容
4-6
e120
6-3
e120

然后启动gentor,会看到如下输出
Recognized bonds:
  1 -   9   1 -  21   2 -  11   2 -  22   3 -   6   3 -  19   3 -  20   4 -   5
  4 -   6   4 -  12   4 -  13   5 -   7   5 -   8   6 -  14   6 -  15   7 -   9
  7 -  16   8 -  10   8 -  17   9 -  11  10 -  11  10 -  18

Rotation angle of torsion   4  -   6:
  0.0 120.0 240.0
Rotation angle of torsion   6  -   3:
  0.0 120.0 240.0
The total number of rotations:    9

Fragment involved in torsion   4  -   6:
      3      14      15      19      20
Fragment involved in torsion   6  -   3:
     19      20

Molecular conformations have been outputted to traj.xyz in current folder
    0 conformations did not pass distance check and thus be ignored

Recognized bonds告诉你有哪些原子被判断为了成键。这很重要,如果判断错了,旋转键的时候片段可能没法跟着一起旋转。程序默认对于两个原子间距离小于其共价半径和的1.15倍就认为是成键的。
Rotation angle of torsion后面是被旋转的键,下面是旋转的角度。The total number of rotations后面是总共扫描的点数,当前体系是3*3=9个。
Fragment involved in torsion后面告诉你在旋转各个键时,哪些原子被视为片段被旋转。
最后程序还告诉你所有结构都通过了检测,没有不合理的接触,所以9个结构最终都输出到了当前目录下的traj.xyz文件中。

我们把得到的traj.xyz拖进VMD程序,可看到以下初始结构,显然都是合理的,可以做下一步了。

init.gif


我们用molclus调用MOPAC2012在PM7级别下对各个结构进行优化。把traj.xyz放到molclus目录下,确认MOPAC的模板输入文件里的关键词是PM7 precise,确保molclus的settings.ini的文件里的iprog为2。然后启动molclus,程序就会调用MOPAC进行优化,在Intel i7-2630QM CPU下每个构象优化只用平均2秒钟,总共用了18秒就都优化完了,结果输出在了当前目录下的isomer.xyz中。运行isostat进行统计,敲三下回车,看到排序后的能量
#    1 Count:    1  E=     -0.108836 a.u.  DGmin=   0.35  DE=    0.00 kcal/mol
#    2 Count:    1  E=     -0.108481 a.u.  DGmin=   0.20  DE=    0.22 kcal/mol
#    3 Count:    1  E=     -0.108090 a.u.  DGmin=   0.21  DE=    0.47 kcal/mol
#    4 Count:    1  E=     -0.105764 a.u.  DGmin=   0.20  DE=    1.93 kcal/mol
#    5 Count:    1  E=     -0.105716 a.u.  DGmin=   0.20  DE=    1.96 kcal/mol
#    6 Count:    1  E=     -0.105189 a.u.  DGmin=   0.15  DE=    2.29 kcal/mol
#    7 Count:    1  E=     -0.105171 a.u.  DGmin=   0.15  DE=    2.30 kcal/mol
#    8 Count:    1  E=     -0.104526 a.u.  DGmin=   0.40  DE=    2.70 kcal/mol
#    9 Count:    1  E=     -0.103219 a.u.  DGmin=   0.40  DE=    3.52 kcal/mol
对能量排序后的结构都输出到了当前目录下的cluster.xyz中。

当然半经验算的能量还不够可靠,不仅定量精度低,能量高低顺序都可能不对,我们用Gaussian在更好的M06-2X/TZVP级别下对它们算单点能。删掉之前的traj.xyz,把cluster.xyz改名为traj.xyz使得刚才优化好的结构作为初始结构。然后把settings.ini里的iprog设为1,ieneonly设为1。Gaussian模板文件template.gjf里面关键词写M062X/TZVP,并设好%nproc和%mem。然后启动molclus。算完之后能量都会记录到当前目录下的isomer.xyz中。运行isostat,按三次回车对isomer.xyz进行统计,结果为
#    1 Count:    1  E=   -516.623091 a.u.  DGmin=   0.35  DE=    0.00 kcal/mol
#    2 Count:    1  E=   -516.622925 a.u.  DGmin=   0.20  DE=    0.10 kcal/mol
#    3 Count:    1  E=   -516.622135 a.u.  DGmin=   0.20  DE=    0.60 kcal/mol
#    4 Count:    1  E=   -516.622102 a.u.  DGmin=   0.20  DE=    0.62 kcal/mol
#    5 Count:    1  E=   -516.621145 a.u.  DGmin=   0.21  DE=    1.22 kcal/mol
#    6 Count:    1  E=   -516.621051 a.u.  DGmin=   0.15  DE=    1.28 kcal/mol
#    7 Count:    1  E=   -516.620880 a.u.  DGmin=   0.15  DE=    1.39 kcal/mol
#    8 Count:    1  E=   -516.620532 a.u.  DGmin=   0.40  DE=    1.61 kcal/mol
#    9 Count:    1  E=   -516.618691 a.u.  DGmin=   0.40  DE=    2.76 kcal/mol

然后若输入一个温度值比如300,可以将单点能近似为自由能得到此温度下的波尔兹曼分布:
#    1   Count:   1   Ratio:  33.58%
#    2   Count:   1   Ratio:  28.18%
#    3   Count:   1   Ratio:  12.27%
#    4   Count:   1   Ratio:  11.86%
#    5   Count:   1   Ratio:   4.33%
#    6   Count:   1   Ratio:   3.92%
#    7   Count:   1   Ratio:   3.27%
#    8   Count:   1   Ratio:   2.27%
#    9   Count:   1   Ratio:   0.33%
可见当前体系有多种构象在常温下都会有较大分布比例。

再次用VMD查看isostat产生的cluster.xyz,其第一帧,也就是M06-2X/TZVP下能量最低的结构为
2.png


仔细看isostat输出的各个构象的能量的话,会看到1、2之间非常相近。观看cluster.xyz中对应的结构会发现1、2的差异只是-CH2-NH2基团一个冲左一个冲右而已,破坏二者镜像对称的只不过是苯环上那个在间位的羟基而已,显然由于它离-CH2-NH2基团甚远,不会造成1、2之间能量出现显著差异。类似地3、4之间,5、6之间能量也很接近,能量差不到0.1kcal/mol。

如果想要更精准的结果,建议优化在M062X/def2SVP下进行,单点能用双杂化泛函结合cc-pVTZ或def2-TZVPP计算。但当体系较大、可能的构象很多时,如几百个,建议先用廉价的半经验,甚至分子力场先做一下粗略优化,以筛选出来能量较低的、有必要之后用DFT进一步优化的构象。用molclus做这种事极方便,详见《使用molclus程序做团簇构型搜索和分子构象搜索》。

评分

参与人数 9eV +40 收起 理由
mooninwhere + 2 赞!
byymem + 5 赞!
liu_tiao + 5 赞!
zsu007 + 5 赞!
卡开发发 + 5 GJ!
minagami + 4 好物!
aqhuangry + 5
dreamyeye + 4 赞!
captain + 5 好物!

查看全部评分

思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

227

帖子

0

威望

1925

eV
积分
2152

Level 5 (御坂)

发表于 2015-12-18 07:57:24 | 显示全部楼层
好物

227

帖子

2

威望

1867

eV
积分
2134

Level 5 (御坂)

发表于 2015-12-18 09:13:47 | 显示全部楼层
设置旋转参数的文件(gentor.ini),每个分子应该不大一样,如果打算批量搜索一群分子,设置这个参数文件,应该会变成一件头疼的事情。

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

 楼主| 发表于 2015-12-18 12:16:19 | 显示全部楼层
greatzdk 发表于 2015-12-18 09:13
设置旋转参数的文件(gentor.ini),每个分子应该不大一样,如果打算批量搜索一群分子,设置这个参数文件, ...


一般不会涉及那种情况的。要考察构象的分子数目不会很多,人工弄很容易。肉眼一看有几个键需要考虑旋转,立刻就写上去了。任何构象搜索程序都避免不了这个过程。如果程序自动指定,若设定不合理可能造成大量时间浪费。

评分

参与人数 1eV +3 收起 理由
greatzdk + 3 有道理

查看全部评分

思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

192

帖子

4

威望

983

eV
积分
1255

Level 4 (黑子)

发表于 2016-1-2 00:00:46 | 显示全部楼层
有没有可能自动识别可以旋转的键?

1665

帖子

19

威望

5375

eV
积分
7420

Level 6 (一方通行)

发表于 2016-1-2 04:06:23 | 显示全部楼层
本帖最后由 liyuanhe211 于 2016-1-2 04:08 编辑
smutao 发表于 2016-1-2 00:00
有没有可能自动识别可以旋转的键?

Spartan可以做这件事情(Set Torsions 之后会自动给出一系列建议猜测),但对较复杂的分子经常不靠谱,包括转不转、要转几个fold等等,需要人为调。
Snap1.png

不过有GUI确实比写文本方便一点。


评分

参与人数 1eV +5 收起 理由
smutao + 5 谢谢

查看全部评分

4

帖子

0

威望

16

eV
积分
20

Level 1 能力者

发表于 2016-5-29 23:09:58 | 显示全部楼层
老师,为什么没有生成cluster.xyz的文件,而且运行isostat,按三次回车后,软件闪了一下就关闭了,我运行的就是帖子中的例子

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

 楼主| 发表于 2016-5-30 00:49:56 | 显示全部楼层
lidanhoney 发表于 2016-5-29 23:09
老师,为什么没有生成cluster.xyz的文件,而且运行isostat,按三次回车后,软件闪了一下就关闭了,我运行的 ...


绝对是你的文件有问题,把你的isomer.xyz贴出来
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

4

帖子

0

威望

16

eV
积分
20

Level 1 能力者

发表于 2016-5-30 19:49:45 | 显示全部楼层
sobereva 发表于 2016-5-30 00:49
绝对是你的文件有问题,把你的isomer.xyz贴出来

是我的mopac的软件安装有问题,现在安装好了就可以用了,不过用gaussian优化是不是很慢,好久都没有什么反应

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

 楼主| 发表于 2016-5-31 00:38:22 | 显示全部楼层
lidanhoney 发表于 2016-5-30 19:49
是我的mopac的软件安装有问题,现在安装好了就可以用了,不过用gaussian优化是不是很慢,好久都没有什么 ...

自行监控gau.out看看算到什么情况了
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

4

帖子

0

威望

16

eV
积分
20

Level 1 能力者

发表于 2016-6-1 18:49:24 | 显示全部楼层
对应的mol.xyz文件内容如下所示,将此文件放到molclus目录下
22
dopamine
O     -2.33158827    1.90324395    0.05204669
老师,这个xyz文件中的22是什么意思 ,我发现在我提取坐标后,设置不一样的值,用VMD看到的结果完全不一样。那这个22是什么意思

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

 楼主| 发表于 2016-6-1 21:20:23 | 显示全部楼层
lidanhoney 发表于 2016-6-1 18:49
对应的mol.xyz文件内容如下所示,将此文件放到molclus目录下
22
dopamine

原子数
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

65

帖子

0

威望

147

eV
积分
212

Level 3 能力者

发表于 2016-6-21 20:13:41 | 显示全部楼层
非常感谢

9

帖子

0

威望

735

eV
积分
744

Level 4 (黑子)

发表于 2016-8-10 12:04:22 | 显示全部楼层
老师好,如果我想搜索过渡态的conformer要怎样处理呢?是把过渡态的部分固定住还是应该怎样比较好呀?谢谢老师!

9119

帖子

16

威望

1万

eV
积分
20756

管理员

公社社长

 楼主| 发表于 2016-8-10 12:34:47 | 显示全部楼层
ZCSco 发表于 2016-8-10 12:04
老师好,如果我想搜索过渡态的conformer要怎样处理呢?是把过渡态的部分固定住还是应该怎样比较好呀?谢谢 ...


先找一个构型下的过渡态,然后对这个过渡态结构用gentor产生各个二面角的各种可能情况的结构。优化的时候把过渡态的原子冻住。
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395(2000人已满),思想家公社QQ群2号:466017436(讨论计算化学为主,加入时必须注明研究方向)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

Archiver|手机版|小黑屋|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2017-5-25 22:30 , Processed in 0.114875 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表