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

计算化学公社

 找回密码
 现在注册!
查看: 21043|回复: 118

[综合交流] 使用molclus程序做团簇构型搜索和分子构象搜索

  [复制链接]

1万

帖子

25

威望

1万

eV
积分
26294

管理员

公社社长

发表于 2015-1-5 18:29:02 | 显示全部楼层 |阅读模式
重要补充说明1:molclus从1.2版开始加入了genmer工具用来产生团簇的初始构型,详见《genmer:生成团簇初始构型的超便捷工具》(http://bbs.keinsci.com/forum.php?mod=viewthread&tid=2369),比起本文例子中通过分子动力学程序产生分子团簇初猜结构方便得多得多。
重要补充说明2:molclus从1.3版开始加入了gentor工具用来通过扫描方式产生分子构象搜索用的初始构型,详见《gentor:扫描方式做分子构象搜索的便捷工具》(http://bbs.keinsci.com/forum.php?mod=viewthread&tid=2388),因此对小分子做构象搜索时就不是必须得像本文那样通过动力学来产生初始构象了,省事不少。


使用molclus程序做团簇构型搜索和分子构象搜索

文/Sobereva(3)
First release:2015-Jan-5   Last update:2017-Dec-1




前言:molclus是北京科音自然科学研究中心开发的一款用于分子团簇构型搜索和分子构象搜索的程序,需要结合第三方分子动力学和量子化学程序使用,本文将对此程序进行介绍。第1节介绍基于分子动力学做构型/构象搜索的基本原理,第2节介绍molclus程序的基本使用,3、4节通过实例进行演示。


1 构型/构象搜索原理

分子团簇的构型搜索和分子的构象搜索是计算化学经常涉及的问题。前者用来寻找分子团簇的各种构型,比如可以回答这样的问题:5个某分子组成的团簇,总共有多少种结构?能量各是多少,哪个构型最稳定?后者用来寻找特定分子的各种可及的构象,比如一个柔性分子,有很多键都可以旋转从而构成不同构象,我们想知道哪个构象能量最低,各种构象能量差是多少,由此我们还可以进一步根据Boltzmann分布估计出特定温度下不同构象出现的比率是多少。

构型/构象搜索的算法很多,比如势能面扫描、分子动力学、蒙特卡罗、基因算法、LMOD等等,适用场合不同。对于分子团簇构型搜索和分子构象搜索而言,分子动力学方法是比较好的选择。其基本过程是:在合适的方法下做一定长度的分子动力学模拟,由于热运动,结构能够反复穿越势垒到达势能面不同区域。每隔一定步数取一次结构,将这些分布在势能面的不同位置上的结构作为初猜结构分别进行优化,最终就会得到势能面上各个极小点位置的构型和能量。

对于分子团簇的构型搜索和分子的构象搜索而言,并不会涉及到键的形成和断裂,因此不需要用极耗时的从头算分子动力学,而只需要计算量极低的分子力场来跑分子动力学就够了。这样的程序很多,比如免费的gromacs、tinker、NAMD、lammps,现在amber实际上也算基本免费了(用免费的ambertools里的sander来跑就行了,不需要用收费的PMEMD部分),以及收费的Material Studio。至于结构优化和能量计算部分,可以用分子力场、半经验、DFT、高精度后HF方法等来做,这就取决于需要的精度和运算能力了。



2 molclus程序

2.1 简介

molclus有Windows、Linux和MacOS版,由Fortran编写。molclus主页为:http://www.keinsci.com/research/molclus.html。程序会不断改进、更新。molclus的免费版可处理的原子数目上限为30个。如果需要处理30个以上原子,需购买完整版,价格96元。源代码有偿提供,价格399元。付费购买2年内可免费升级,并提供技术支持,免费版不提供技术支持。购买和下载事宜请访问molclus的主页。

molclus程序可以视为是分子动力学程序与量子化学程序之间的一个接口。几乎所有动力学程序都可以用,包括自编的,只要能输出.xyz格式的轨迹或者输出的轨迹能被免费的VMD程序(http://www.ks.uiuc.edu/Research/vmd/)识别即可。量化程序目前支持Gaussian、ORCA、MOPAC(2012及之后的版本)、xtb。其中Gaussian是大家最常用的,主要适合小体系。ORCA由于RI技术,在纯泛函下比Gaussian快得多得多,因此可以用于中等体系,当然对于小体系也没问题。MOPAC专门做半经验计算,它支持的PM6-DH+通常足矣得到定性合理的结果,而速度比DFT快两个数量级,故很适合用于中、大体系的优化。xtb是Grimme开发的做GFN-xtb的程序,此方法相当于半经验的DFT,精度与PM7等半经验方法互有胜负,但普适性更好,而且速度非常快。可以说,molclus支持的Gaussian、ORCA和MOPAC把不同尺度、不同精度范围都覆盖了。如果不熟悉ORCA、RI技术和MOPAC,请务必参看《大体系弱相互作用计算的解决之道》(http://sobereva.com/214)当中的介绍。分子团簇和弱相互作用关系极为密切,必须选择合适的计算级别,如果不熟悉弱相互作用计算的话建议看看《乱谈DFT-D》(http://sobereva.com/83)。


2.2 使用流程

molclus直接解压后就能用。基本使用流程如下:

(1)用分子动力学程序模拟分子团簇或单个分子,每隔一定步数将结构写入一次轨迹。力场并不要求很准确。
(2)用免费的VMD程序将动力学程序输出的私有轨迹格式转化为标准的.xyz格式轨迹文件,命名为traj.xyz。
(3)设定好量化程序的几何优化模板文件,模板文件需放到当前目录下,以template命名。对于Gaussian程序,后缀为.gjf;对于ORCA,后缀为.inp;对于MOPAC,后缀为.mop;而xtb不需要设模板文件。molclus会将模板文件里面的[GEOMETRY]字段替换为traj.xyz里的结构并传递给量化程序执行,因此应该在模板文件设定好基组、方法和相关关键词(比如内存分配量、辅助收敛的设定等)。
(4)设定好molclus的参数文件settings.ini,详见2.3节。
(5)启动molclus。molclus会载入traj.xyz、settings.ini和模板文件,然后自动调用量子化学程序对traj.xyz的每一帧结构进行优化,并将优化好坐标和能量写入到当前目录下的isomers.xyz文件中。
(6)启动molclus中附带的isostat工具,它会载入isomers.xyz,对其中的结构根据能量进行排序并进行统计、去除能量和几何结构相似度过高的重复结构。isostat在当前目录下输出的cluster.xyz文件就是根据能量从低到高排列并且去除了重复结构后的文件,可以载入VMD来方便地观看每个结构。之后还可以输出指定温度下各个构型的波尔兹曼分布比例(原理见http://sobereva.com/165,理应用自由能进行计算,isostat用的是单点能作为其近似)。

如果想使用比结构优化更高级别的方法精确计算每个结构的能量,可以把settings.ini里的ieneonly改为1,把模板文件里的和优化有关的关键词都去掉,再把cluster.xyz改名为traj.xyz然后启动molclus,molclus就会调用量化程序对其中每个结构再做一次单点计算。

每次程序运行前,会对当前目录下.out、.arc等量化程序之前产生的临时文件进行清理,所以如果有重要的临时文件在当前目录下,执行molclus前应当先将它们挪至它处。


2.3 参数文件

程序目录下的settings.ini在molclus启动时会被载入,其中记录了运行参数,在下面进行介绍。大部分参数都不用改。

isys:1=Windows;2=Linux。在什么系统中运行molclus就设为什么
iprog:调用的量化程序。1=Gaussian;2=MOPAC;3=ORCA;4=xtb
ngeom:0=处理traj.xyz里每一帧;n=处理前n帧;n,m=处理n~m帧(比如2,5即处理2,3,4,5帧)
ieneonly:0=对traj.xyz里每一帧进行优化;1=对每一帧只计算单点能
distmax:n=如果结构中任意两个原子间的距离超过n埃则跳过这一帧。做团簇动力学的时候,有可能某个分子跑飞了,对这样的结构优化没有意义,利用这个参数就可以忽略掉此结构
ipause:1=对于Gaussian和ORCA,计算出错(如不收敛)则暂停,以便用户对输出文件进行检查,然后可以按回车继续处理下一帧结构;2=每处理完一个结构后都暂停;0=总是不暂停
ibkout:1=每处理完一帧结构,把程序的输出文件备份,以便之后检查或根据自己的需要提取其它数据。比如对于第5帧结构,Gaussian的输出文件会备份到gau00005.out,ORCA的会备份到orca00005.out,MOPAC的会备份到MOPAC00005.arc;0=不做备份,因此molclus运行完毕后将找不到量化程序中间输出的文件

----以下是Gaussian专用的参数----
gaussian_path:Gaussian的可执行文件的路径,比如"/sob/g09/g09"、"d:\study\g09w\g09.exe"。也可以只写"g09",前提是程序所在目录已经加入到了PATH环境变量里。如果路径中有空格,则必须有双引号,没空格的话无所谓,后同
[注意:如果是Windows版Gaussian,gaussian_path必须指向诸如g09.exe而不能是g09w.exe,经常有人在这里犯错。另外,运行前还必须进入“控制面板”-“高级系统设置”-“高级”-“环境变量”,添加GAUSS_EXEDIR环境变量,使之指向Gaussian目录,如d:\study\g09w,否则Gaussian没法以命令行方式运行]
igaucontinue:1=如果优化是因为超过了步数上限而引起的,则会基于另一个模板文件template2.gjf对最后一帧结构进行继续优化;0=不做这个处理。igaucontinue这个选项主要有两个主要用处。当初始结构偏离平衡结构太远时,无论Gaussian输入文件里opt=maxcyc=n设多大,实际步数都不允许超过程序内定上限(和坐标变量数成正比),因此可能因为步数上限不够而导致没收敛,因此可以基于template2.gjf继续优化。另一种情况是基于template.gjf里的关键词在优化时发生了震荡而没有收敛,若在template2.gjf里写上备用方案的关键词,诸如calcfc、GDIIS、opt(maxstep=x,notrust)等等,可能在这些关键词下继续优化的时候就能收敛了。关于辅助收敛的方法,在此贴有汇总:《量子化学计算中帮助几何优化收敛的常用方法》(http://sobereva.com/164
energyterm:molclus是从Gaussian输出文件中末尾的archive字段部分读取能量的,这个选项用于设定读取哪项。比如是HF/DFT计算,这个参数应设为HF=,即读取HF=这个字符串后面紧跟着的那个数(HF/DFT能量)。如果是MP2或双杂化泛函计算,就应当设MP2=。以此类推
ibkchk:和ibkout的用法一样,但这个控制的是如何备份.chk文件

----以下是ORCA专用的参数----
orca_path:ORCA可执行文件的路径,比如"D:\study\orca\orca.exe"
ibktrj:1=把每次优化后产生的.trj改名为orca[帧号].xyz。这是优化过程的轨迹文件,因此可以用VMD直接观看优化过程
mpioption: ORCA并行运算时传递给MPI的额外选项。默认为none,即什么也不设。此选项的用处例如,在root下运行较新版本的ORCA,不带--allow-run-as-root选项就运行不了,此时可以设定mpioption= --allow-run-as-root来解决此问题。

----以下是MOPAC专用的参数----
mopac_path:MOPAC可执行文件的路径,比如"C:\Program Files\MOPAC\MOPAC2016.exe"。



3 分子团簇构型搜索实例:(H2O)4

此例我们通过molclus结合Gromacs 4.6.7分子动力学程序以及MOPAC2012、Gaussian09 D.01量子化学程序搜索(H2O)4团簇构型,这个团簇已经被很多文章研究过了。这里用这种很小的团簇进行示例仅是为了节约时间,对无论多大的团簇计算过程基本都是一样的。

Gromacs是免费、最主流、速度快且很好用的小分子和生物大分子的分子动力学模拟程序,所以我们这里用Gromacs,编译方法参见《Gromacs 5.0与4.6.7编译方法》(http://sobereva.com/247)。大家也可以用其它自己擅长的动力学程序,只要能产生VMD能支持的轨迹就行,因为我们要靠VMD来把轨迹转为molclus能读取的.xyz轨迹格式(当然啦,如果自己有其它办法能转成.xyz轨迹格式也行)。

动力学我们用NVT系综。对于当前体系实际上不用周期边界条件也行。不过对于较大团簇的话,当模拟温度较高时,可能模拟过程中外侧的某个分子跑着跑着就飞了,越飞越远,之后的轨迹就没意义了。而如果用了周期边界条件,它飞过边界后还会跑回来。

Gromacs模拟部分所有涉及到的文件都可以在这里下载,如果有的地方不清楚直接看里面的文件就明白了。
H20_4_MD.rar.rar (55.96 KB, 下载次数: 308)

评分

参与人数 20eV +84 收起 理由
ZBC + 3 精品内容
Graphite + 5
byymem + 5 赞!
万山红遍 + 4 谢谢
greatzdk + 2 GJ!
momian + 5 赞!
十年磨练 + 4
fzmn + 4 赞!
hcxytpp@163.com + 5 赞!
chittyda + 4 GJ!
zsu007 + 3 谢谢
诸葛壹次心 + 3 谢谢
vigaryang + 3 牛!
ter20 + 10 好物!
therotyonth + 3 谢谢
Xououw + 2 好物!
captain + 5 牛!
卡开发发 + 5 价格相当公道:-)
平辉正 + 5 赞!简直不能更好
diaolanxinyu + 4 感谢分享!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog:http://sobereva.com
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明研究方向,否则一概不批。

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

117

帖子

0

威望

1088

eV
积分
1205

Level 4 (黑子)

发表于 2015-1-5 18:44:02 | 显示全部楼层
感谢sobereva老师分享!

66

帖子

0

威望

566

eV
积分
632

Level 4 (黑子)

发表于 2015-1-5 20:25:36 | 显示全部楼层
赞一个

142

帖子

0

威望

1161

eV
积分
1303

Level 4 (黑子)

发表于 2015-1-5 21:51:16 | 显示全部楼层
顶,绝对的支持~

243

帖子

2

威望

1673

eV
积分
1956

Level 5 (御坂)

发表于 2015-1-6 11:56:31 | 显示全部楼层
原理略有些简单啊。。用MD方法搜索团簇的效率还是比MC低了不少,这样很难得到团簇的global minima吧
还有一个问题就是MD所得到的结构用于量化的OPT的话很容易不收敛的,因为距离极小值的距离比较远

如果要进一步开发的话,建议放弃MD软件,用basin hopping进行数据采样,这样得到的结构是一个分子力学上的极小值,既可以节约采样的时间,也可以节约量化结构优化的时间

1万

帖子

25

威望

1万

eV
积分
26294

管理员

公社社长

 楼主| 发表于 2015-1-6 12:28:21 | 显示全部楼层
fhh2626 发表于 2015-1-6 11:56
原理略有些简单啊。。用MD方法搜索团簇的效率还是比MC低了不少,这样很难得到团簇的global minima吧
还有 ...


我开发过MC搜索团簇的程序,对于搜索原子团簇适合,但这种方法明显不适合分子团簇(本文强调的是分子团簇),或者说效率比MD明显要低,而且收敛难度反倒远比基于MD大得多,毕竟MC过程根本没有很好地考虑原子间的作用,别提几何优化收敛了,由于随机移动原子造成的构型不合理,甚至一开始就可能SCF不收敛,甚至分子团簇产生了错误的结构(比如氢原子干脆优化后发生转移了)。

分子动力学得到分子团簇全局极小点没有任何问题,是分子团簇构型搜索有效的方法,更是众多文献讨论分子团簇(包括很大团簇)的标准方法。basin hopping、metadynamics之类采样方未必效率更高,而且有局限性。

“原理简单”不是缺点反倒是优点,用简单的原理,很好地解决实际需要解决的问题,这才是好的、适合大众使用的程序。
反之原理复杂,甚至引入一堆参数,需要使用者需要大量经验积累,颇为麻烦地才能使用来解决看似不复杂的问题,甚至还得改反复改代码,那样的程序我不推崇,不是大众向的,而只是少数专家们发paper用的。
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog:http://sobereva.com
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明研究方向,否则一概不批。

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

107

帖子

0

威望

859

eV
积分
966

Level 4 (黑子)

发表于 2015-1-6 14:45:38 | 显示全部楼层
那个settings.ini中的ipause参数=3表示总不暂停?,但是文件中写的是=0表示不暂停。好像是,,

1万

帖子

25

威望

1万

eV
积分
26294

管理员

公社社长

 楼主| 发表于 2015-1-6 15:00:17 | 显示全部楼层
tjuchan 发表于 2015-1-6 14:45
那个settings.ini中的ipause参数=3表示总不暂停?,但是文件中写的是=0表示不暂停。好像是,,

帖子里写错了,已更正。
其实无所谓,只要不是1和2就行,所以0和3效果都一样
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog:http://sobereva.com
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明研究方向,否则一概不批。

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

107

帖子

0

威望

859

eV
积分
966

Level 4 (黑子)

发表于 2015-1-6 16:44:08 | 显示全部楼层
综合性有点强,很难把几类软件串在一起。

270

帖子

0

威望

2340

eV
积分
2610

Level 5 (御坂)

发表于 2015-1-6 18:08:54 | 显示全部楼层
Sob 老师强,赞一个!

49

帖子

0

威望

538

eV
积分
587

Level 4 (黑子)

发表于 2015-1-6 20:35:56 | 显示全部楼层
1)该软件的思路简单但巧妙,确能有效解决团簇构型和构象搜索问题。
2)MD我用过gromacs的比较老的版本,它的安装和上手还是有些麻烦;对于体系包含1:1两种分子(或三种)的比较复杂的体系,有时不太好建模。能否MD软件也能支持MS做的MD轨迹文件呢,它还是比较好上手和建模的;我记得作者曾经也写过个MS轨迹文件转换的程序,这应当不算太困难。
3)QM能支持gaussian、orca、mopac,这个基本满足各种精度的需要了。

1万

帖子

25

威望

1万

eV
积分
26294

管理员

公社社长

 楼主| 发表于 2015-1-6 21:24:02 | 显示全部楼层
wei 发表于 2015-1-6 20:35
1)该软件的思路简单但巧妙,确能有效解决团簇构型和构象搜索问题。
2)MD我用过gromacs的比较老的版本, ...


gromacs还好。而且也有第三方编译的windows版,上来直接就能用。至于上手,其实效仿文章的例子,很快就能弄会,参数文件也就改改模拟时间、保存频率、温度这些,没多少需要变动的。
建模的话强烈推荐packmol(http://www.ime.unicamp.br/~martinez/packmol/),想建立什么样的都能建立,十分方便灵活。
MS模拟的轨迹可以通过perl脚本导出成xyz(http://sobereva.com/143),直接就能给molclus使用,不需要对molclus做任何修改。但我觉得与其学MS,还不如就用gmx或amber,MS价格不菲。
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog:http://sobereva.com
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明研究方向,否则一概不批。

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

118

帖子

0

威望

1512

eV
积分
1630

Level 5 (御坂)

发表于 2015-1-7 10:01:50 | 显示全部楼层
强就一个字

6

帖子

0

威望

26

eV
积分
32

Level 2 能力者

发表于 2015-1-7 18:08:23 | 显示全部楼层
SOB威武!这个软件对阴离子、阳离子构成的化合物是否适用?

1万

帖子

25

威望

1万

eV
积分
26294

管理员

公社社长

 楼主| 发表于 2015-1-7 18:48:52 | 显示全部楼层
花非花 发表于 2015-1-7 18:08
SOB威武!这个软件对阴离子、阳离子构成的化合物是否适用?


同样适用。其实这和molclus本身没直接关系。对于一个体系,只要动力学程序能模拟,量化程序能算,那么就都可以用molclus做搜索。
北京科音自然科学研究中心:http://www.keinsci.com  不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog:http://sobereva.com
Multiwfn量子化学波函数分析程序主页:http://sobereva.com/multiwfn
计算化学公社论坛:http://bbs.keinsci.com
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明研究方向,否则一概不批。

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

本版积分规则

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

GMT+8, 2017-12-15 12:34 , Processed in 0.141738 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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