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

计算化学公社

 找回密码
 现在注册!
查看: 4663|回复: 41

[综合交流] genmer:生成团簇初始构型的超便捷工具

[复制链接]

9451

帖子

16

威望

1万

eV
积分
21571

管理员

公社社长

发表于 2015-12-15 12:18:02 | 显示全部楼层 |阅读模式
genmer:生成团簇初始构型的超便捷工具

文/Sobereva(3)  2015-Dec-15



1 前言

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

对很多人来说,molclus搜索团簇构型有一个不便之处就是需要通过第三方程序(一般是分子动力学程序)来产生一批初始构型。此文介绍笔者开发的genmer程序来快速方便地生成molclus所需要的初始构型。genmer作为molclus 1.2版及之后版本当中的一个组件来发布,genmer的加入令molclus搜索团簇构型比之前版本方便得多得多,可以说是significant程度的改进。

获得genmer的方法就是到molclus主页(http://www.keinsci.com/research/molclus.html)下载最新版本molclus程序,预编译的genmer就在其中。molclus程序的免费版有原子数限制,但附带的genmer没有原子数限制。


2 原理

genmer可以把任意种、任意数目的分子放在一起随机生成指定数目的团簇初始构型。

genmer生成初始构型的算法是洒家想出来的,很简单又高效实用。是让被加入的分子先随意旋转一下,然后从团簇几何中心出发,向着任意一个方向不断移动,直到被加入分子和团簇没有不合理接触为止(即被加入分子和团簇的原子间的距离都大于原子的范德华半径和)。最开始时团簇就是第一个分子,让分子不断加入,团簇就会不断长大,直到指定数目的分子都被加进去为止。由于分子自身的旋转,还有加入的朝向都是随机的,所以每次生成的团簇构型都会不一样。

假设设定了团簇由10个A、5个B分子构成,程序会先把10个A都加进去,再把5个B加进去,根据如上算法,得到的团簇中10个A就会倾向于在团簇内层,5个B会倾向于在团簇外层。如果想让它们的出现区域不受设定顺序的影响而是完全随机的,genmer也可以做到自动打乱顺序。

在产生团簇时每个分子都被当成是刚性的。但如果一种分子有多种可能的构象都需要考虑,可以提供包含多个构象结构的分子的.xyz文件,并且设定每种构象出现的概率。这样每次加入分子时就会按概率抽取一种构象。

genmer也可以生成原子团簇的初始构型。算法同上,只不过是把分子改成单个原子,并改用共价半径来判断接触。

虽然genmer产生的团簇初始构型可能看起来略微松散,但这完全没关系,一做优化就自动收缩变得致密了。molclus用的团簇构型搜索的关键是让初始构型分布得足够随机且全面,这样才可能经过优化后获得尽可能多的极小点构型,特别是采样到能量最低结构。只要让genmer产生的构型数目足够多,是可以达到这个目的的。


3 使用方法

genmer的使用非常容易,只需要提供每种分子的几何结构.xyz文件,设定出现数目就行了。运行参数和分子设定都记录在genmer.ini文件中,直接启动genmer后程序会读取当前目录下的genmer.ini,然后把生成的指定数目的初始构型都输出到当前目录下的traj.xyz当中。然后就可以直接运行molclus来自动调用量化程序去优化其中的每个结构了。

下面通过一个例子介绍genmer的参数文件genmer.ini的格式:

50  // 要生成的团簇构型的数目
0.1  // 设定每一步移动分子的距离。对于原子数少的情况不用改这个,但如果要产生的团簇原子数甚巨,想节省时间可以把这个步长值设大,但出现不合理接触的可能性会增加。
1    // 0: 不打乱分子加入顺序  1: 打乱顺序
1    // 0: 使用CSD共价半径来检测原子间接触  1: 使用bondi范德华半径检测原子间接触
---- 以下是每类分子的数目和结构文件,可以写无数多个 ----
4    //第1类分子数目
examples\H2O.xyz   //第1类分子的结构文件
3    //第2类分子数目
examples\ethanol.xyz   //第2类分子的结构文件


4 例子

下面给出几个genmer的例子展现其生成团簇的能力。用到的分子结构文件都是genmer自带的。

4.1 6个水团簇

用以下genmer.ini文件可以得到由6个水组成的20个团簇初始构型
20
0.1
1
1
----
6
examples\H2O.xyz

结果如下:
6H2O.gif


4.2 3000个水团簇

用以下genmer.ini文件生成1个多达3000个水组成的团簇初始构型,在笔者i7-2630QM CPU上只需13秒
1
0.1
1
1
----
3000
examples\H2O.xyz

结果如下:
3000H2O.png

4.3 6个水和3个乙醇的团簇

用以下genmer.ini文件可生成6个水和3个乙醇构成的10个团簇,并且打乱乙醇和水的出现顺序
10
0.1
1
1
----
6
examples\H2O.xyz
3
examples\ethanol.xyz

结果如下:
6H2O_3ethanol.gif

4.4 顺式和反式1,3丁二烯混合团簇

看过上面的例子,肯定知道怎么获得顺式和反式1,3丁二烯的混合团簇了,只要提供一个顺式的.xyz文件,一个反式的.xyz文件,设定好各自的数目即可。但如前所述,也可以把多构象的分子的各个构象合一起作为单个.xyz文件提供。这里产生一个由30个1,3丁二烯构成的团簇,并让顺式、反式出现概率为30%和70%。

1,3丁二烯.xyz文件结构如下
10
0.3 cis   //顺式的出现比例为0.3(30%)
[顺式的坐标]
10
0.7 trans   //反式的出现比例为0.7(70%)
[反式的坐标]

genmer.ini内容为
1
0.1
1
1
----
30
examples\butadiene_cis_trans.xyz

结果如下:
30butadiene_cis_trans.png

4.5 多巴胺二聚体

用以下genmer.ini文件可以生成15个多巴胺二聚体初始构型,这算是稍微大一些的分子了。
15
0.1
1
1
----
2
examples\dopamine.xyz

结果如下:
dopamine_dimer.gif

4.6 水层包围多巴胺

用以下genmer.ini文件可以生成200个水分子包裹在多巴胺附近的10种构型
10
0.1
0   //不打乱顺序
1
----
1
examples\dopamine.xyz
200
examples\H2O.xyz

结果如下:
dopamine_H2O.gif

4.7 Li6团簇

用以下genmer.ini文件可以生成6个Li组成的团簇的15个初始构型。
15
0.1
1
0   //注意这里用了共价半径
----
6
examples\Li.xyz

结果如下:
Li6.gif


5 结合molclus做团簇构型搜索

上一节我们得到了团簇的traj.xyz文件,这里我们基于这些团簇初始构型,用molclus做构型的批量优化、获得能量。我们用多巴胺二聚体、乙醇+水团簇和Li6团簇这三个体系作为例子。如果有没看懂,应当先把前述的《使用molclus程序做团簇构型搜索和分子构象搜索》一文仔细看一遍了解molclus的用法。

5.1 多巴胺二聚体

这里我们用molclus调用MOPAC2012在PM6-DH+级别下优化多巴胺二聚体来得到它们的团簇构型。

把前面4.5节得到的含有15个多巴胺二聚初始构型的traj.xyz文件放到molclus目录下,把molclus的settings.ini里的MOPAC路径配置好,确认MOPAC的模板输入文件template.mop里的关键词是PM6-DH+ precise,并且iprog设为了2。启动molclus,程序自动对15个结构依次优化,很快就算完了,结果输出到当前目录下的isomer.xyz里。然后启动isostat,连敲回车三次,即看到按照由低到高排序的能量值,这里只列出前5个
#    1 Count:    1  E=     -0.243275 a.u.  DGmin=   0.28  DE=    0.00 kcal/mol
#    2 Count:    1  E=     -0.242410 a.u.  DGmin=   0.26  DE=    0.54 kcal/mol
#    3 Count:    1  E=     -0.241462 a.u.  DGmin=   0.40  DE=    1.14 kcal/mol
#    4 Count:    1  E=     -0.240707 a.u.  DGmin=   0.27  DE=    1.61 kcal/mol
#    5 Count:    1  E=     -0.239548 a.u.  DGmin=   0.36  DE=    2.34 kcal/mol
然后可以按回车退出,也可以输入温度,比如300,得到各种构象的波尔兹曼分布比例(如果不知道原理,看此文http://sobereva.com/165。isostat是把单点能近似作为自由能代入公式计算的)
#    1   Count:   1   Ratio:  60.30%
#    2   Count:   1   Ratio:  24.28%
#    3   Count:   1   Ratio:   8.95%
#    4   Count:   1   Ratio:   4.04%
#    5   Count:   1   Ratio:   1.19%

按照能量排序后的结构都输出到了当前目录下的cluster.xyz当中。拖入VMD观看,前5个结构图如下,与上面列出的5个依次对应
1.png

为什么第1个结构能量很低是很显然的,因为两个氨基都各自和两个羟基同时形成了氢键。第2个结构是每个氨基和一个羟基形成氢键。第3个结构形成了明显的pi-pi堆积,而且两个氨基之间有氢键。第4个结构有两对氨基与羟基形成的氢键,而且还有一点pi氢键的特征。第5个结构只有一对氢键,但是pi氢键的特征更明显。

由于本例只考虑了15个多巴胺二聚体结构,不算很多,而二聚体极小点构型有很多,还不能确保二聚体最稳定构型就在其中,要想更稳妥应该考虑更多结构。另外,半经验方法精度终究比较有限,获得能量接近的结构的能量差定量不准却,甚至顺序也不对。如果要得到准确的二聚体结构和能量,建议在PM6-DH+优化完之后用molclus调用Gaussian在M06-2X-D3/6-311G*级别下对各个结构再优化一下,然后用M06-2X-D3/6-311+G**级别(要求更精确可以用B2PLYP-D3/jul-cc-pVTZ)计算能量最低的数个结构的单点能。这种“半经验预优化->DFT结合一般基组优化->高精度方法算单点能"是我很推荐的研究分子团簇的做法,这用molclus实现极其容易,可参考《使用molclus程序做团簇构型搜索和分子构象搜索》中计算水四聚体团簇构型的例子。


5.2 乙醇+水团簇

这一节用Gaussian里的UFF力场优化4.3节得到的6个水+3个乙醇的团簇的10个初始构型。力场计算极其便宜,可以比半经验优化更大体系、更大批量构型。(由于MOPAC2012速度很快,对于当前这样不大的体系,Gaussian做UFF计算比之快不了多少)。

把4.3节的traj.xyz拷到molclus目录下,把template.gjf里的关键词设为UFF=qeq opt,然后确认settings.ini文件里iprog设为了1。启动molclus,约一分钟就把所有构型都优化完了。用isostat统计结果如下
#    1 Count:    1  E=     -0.063488 a.u.  DGmin=   0.50  DE=    0.00 kcal/mol
#    2 Count:    1  E=     -0.063090 a.u.  DGmin=   0.38  DE=    0.25 kcal/mol
#    3 Count:    1  E=     -0.062488 a.u.  DGmin=   0.50  DE=    0.63 kcal/mol
#    4 Count:    1  E=     -0.062053 a.u.  DGmin=   0.38  DE=    0.90 kcal/mol
#    5 Count:    1  E=     -0.058900 a.u.  DGmin=   0.34  DE=    2.88 kcal/mol
#    6 Count:    1  E=     -0.058547 a.u.  DGmin=   0.34  DE=    3.10 kcal/mol
#    7 Count:    1  E=     -0.057159 a.u.  DGmin=   0.53  DE=    3.97 kcal/mol
#    8 Count:    1  E=     -0.056104 a.u.  DGmin=   1.38  DE=    4.63 kcal/mol
#    9 Count:    1  E=     -0.055407 a.u.  DGmin=   2.19  DE=    5.07 kcal/mol
#   10 Count:    1  E=     -0.049182 a.u.  DGmin=   0.39  DE=    8.98 kcal/mol

300K下的分布比例
#    1   Count:   1   Ratio:  44.57%
#    2   Count:   1   Ratio:  29.34%
#    3   Count:   1   Ratio:  15.56%
#    4   Count:   1   Ratio:   9.85%
#    5   Count:   1   Ratio:   0.36%
#    6   Count:   1   Ratio:   0.25%
#    7   Count:   1   Ratio:   0.06%
#    8   Count:   1   Ratio:   0.02%
#    9   Count:   1   Ratio:   0.01%
#   10   Count:   1   Ratio:   0.00%

按能量排序后的结构如下。和4.3节中genmer初始产生的结构对比,可看到经过力场优化后团簇变得明显紧凑了
6H2O_3ethanol_UFF.gif


5.3 Li6团簇

这一节我们用molclus调用Gaussian在B3LYP/6-31G*级别下寻找Li6团簇稳定构型。

把4.7节得到的包含15个Li6团簇初猜构型的traj.xyz放到molclus目录下。将Gaussian的模板输入文件template.inp里的关键词填入B3LYP/6-31G* opt(cartesian,maxcyc=50),备用模板文件template2.inp里的关键词填入B3LYP/6-31G* opt(cartesian,maxstep=5,notrust,maxcyc=50,gdiis)。将iprog设为1,igaucontinue设为1。含义就是说,对每个结构先用第一个模板文件的设定优化,如果不收敛或报错,则取最后结构用第二个模板的设定继续优化(若不明白为什么这么用关键词看此文http://sobereva.com/164)。模板文件步数改为50是因为对当前体系Gaussian默认的步数上限太少,所以适当增大。优化过程中可能会由于构成键角、二面角的冗余内坐标原子间排成直线而报错,所以这里用笛卡尔坐标优化来避免。

启动molclus,程序对15个结构依次优化,结果输出到isomer.xyz里。然后启动isostat,连敲回车三次,看到有4种唯一的结构
#    1 Count:    2  E=    -45.110127 a.u.  DGmin=   0.02  DE=    0.00 kcal/mol
#    2 Count:    3  E=    -45.109176 a.u.  DGmin=   0.02  DE=    0.60 kcal/mol
#    3 Count:    7  E=    -45.105556 a.u.  DGmin=   1.55  DE=    2.87 kcal/mol
#    4 Count:    2  E=    -45.103271 a.u.  DGmin=   2.03  DE=    4.30 kcal/mol
相应的结构输出到了cluster.xyz中。从结构看会发现前两个结构其实是一样的,之所以能量有细微差别其实是收敛精度问题而已。所以实际上有三个唯一的团簇结构,如下所示,从左到右能量依次升高
2.png

这些团簇从图上看都有对称性,而当前优化过程由于收敛精度原因可能最后一步的构型不严格满足实际对称性,故建议自行在gview里把结构做对称化后再优化一下以得到更准确结果(把对称化后的坐标写成traj.xyz的格式,然后用molclus再跑一次就行了)。

评分

参与人数 14eV +75 收起 理由
一颗赛艇 + 4 赞!
archer + 4 とてもいい!
tjuchan + 5 好物!
zidu113 + 5 牛!
zhou + 5 赞!
captain + 5 好物!
幻七熏 + 5 赞!
zsu007 + 5
aqhuangry + 5
chrischen1128 + 12 好物!
978142355 + 5 GJ!
hcxytpp@163.com + 5 谢谢
zhanfei + 5 赞!
greatzdk + 5 牛!

查看全部评分

思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。共计近4000人,讨论计算化学为主,加入时必须注明研究方向,否则一律不批。

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

2066

帖子

8

威望

2440

eV
积分
4666

Level 6 (一方通行)

首席卖萌官

发表于 2015-12-15 15:22:29 | 显示全部楼层
沙发是娃娃的!!

164

帖子

0

威望

1433

eV
积分
1597

Level 5 (御坂)

Miss X

发表于 2015-12-15 19:36:57 | 显示全部楼层
谢谢老师介绍这么方便使用的软件!!!
真诚相待,淡淡相处,方天久地长^^

236

帖子

0

威望

1988

eV
积分
2224

Level 5 (御坂)

发表于 2015-12-15 22:52:54 | 显示全部楼层
让我想到了packmol !  Sob老师

64

帖子

0

威望

1250

eV
积分
1314

Level 4 (黑子)

发表于 2015-12-16 09:56:45 | 显示全部楼层
不用跑分子动力学了,减少了不少麻烦啊。

20

帖子

0

威望

471

eV
积分
491

科音成员

发表于 2015-12-16 14:21:44 | 显示全部楼层
大赞Sob!

32

帖子

0

威望

544

eV
积分
576

Level 4 (黑子)

发表于 2015-12-22 13:10:38 | 显示全部楼层
这软件真好,赞!

30

帖子

0

威望

217

eV
积分
247

Level 3 能力者

发表于 2015-12-23 07:19:35 | 显示全部楼层
做团簇的福音啊!大赞!

64

帖子

0

威望

1250

eV
积分
1314

Level 4 (黑子)

发表于 2015-12-28 17:05:52 | 显示全部楼层
对于电中性体系非常不错,但是带电荷的体系就不太好,有什么办法来计算带电体系吗?

13

帖子

0

威望

97

eV
积分
110

Level 2 能力者

发表于 2015-12-30 17:04:40 | 显示全部楼层
      之前都是用Packmol来构建团簇的,那个软件的逻辑是你自己设定一个三维空间的尺寸,然后再设定该空间内的原子数目。所以用Packmol经常还要自己预先估算一下多大空间适合放多少分子,有时候数量填多了会报错,填少了又太疏了,总之是用起来不太方便。这个算法只要输入数量就可以自动生成团簇了,确实简单高效啊。
      不知道用该软件可否构建适用于分子动力学模拟的初始构型?如果可以的话,那么用该软件所构筑的团簇的Box尺寸有一起输出吗,此外除了正方体的Box,可否构建适用于其他更复杂形状box的初始构型,比如长方体或截角八面体?我觉得再多加一些限定条件应该就行了。

9451

帖子

16

威望

1万

eV
积分
21571

管理员

公社社长

 楼主| 发表于 2015-12-31 00:40:31 | 显示全部楼层
幻七熏 发表于 2015-12-28 17:05
对于电中性体系非常不错,但是带电荷的体系就不太好,有什么办法来计算带电体系吗?

如果是正、负电荷离子混杂,需要额外修改程序,增加一些规则,尽量让两类交替出现。
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。共计近4000人,讨论计算化学为主,加入时必须注明研究方向,否则一律不批。

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

9451

帖子

16

威望

1万

eV
积分
21571

管理员

公社社长

 楼主| 发表于 2015-12-31 00:41:50 | 显示全部楼层
tzweir 发表于 2015-12-30 17:04
之前都是用Packmol来构建团簇的,那个软件的逻辑是你自己设定一个三维空间的尺寸,然后再设定该空间 ...


此程序目的不在于建立动力学模拟的初始构型,只是用于产生团簇构型搜索的初始构型,还没打算在建模的功能设置上做扩展。
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。共计近4000人,讨论计算化学为主,加入时必须注明研究方向,否则一律不批。

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

10

帖子

0

威望

19

eV
积分
29

Level 2 能力者

发表于 2016-1-5 14:05:01 | 显示全部楼层
多谢分享,对于计算化学初学者相当有用。

24

帖子

0

威望

776

eV
积分
800

Level 4 (黑子)

发表于 2016-1-15 11:44:48 | 显示全部楼层
虽然genmer产生的团簇初始构型可能看起来略微松散,但这完全没关系,一做优化就自动收缩变得致密了。

这样是说结构优化可以起到分子动力学一样的效果吗?我用genmer产生了几个分子团簇,用高斯结构优化,gview看优化的过程觉得结构优化的过程与分子动力学不一样,这样的未经处理的随机初始构型用来结构优化真的可以做到构象搜索吗?
2015毕业季

9451

帖子

16

威望

1万

eV
积分
21571

管理员

公社社长

 楼主| 发表于 2016-1-15 21:03:40 | 显示全部楼层
zhou 发表于 2016-1-15 11:44
这样是说结构优化可以起到分子动力学一样的效果吗?我用genmer产生了几个分子团簇,用高斯结构优化,gvie ...

无论是genmer这样随意堆积,还是分子动力学,目的都是产生随机分布的构型(团簇要说构型而非构象),在势能面上随机采样,这样优化后就能找到最近的极小点。这种做法搜索分子团簇没有任何问题,又方便效果又好。
思想家公社的门口Blog:http://sobereva.com
北京科音自然科学研究中心:http://www.keinsci.com
Multiwfn主页:http://sobereva.com/multiwfn
计算化学公社:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。共计近4000人,讨论计算化学为主,加入时必须注明研究方向,否则一律不批。

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

本版积分规则

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

GMT+8, 2017-6-24 21:43 , Processed in 0.126901 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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