请选择 进入手机版 | 继续访问电脑版
第9届北京科音分子动力学与GROMACS培训班将于4月17~20日于北京举办,请点击此链接查看培训详情,欢迎参加和相互转告!

计算化学公社

 找回密码
 现在注册!
查看: 100217|回复: 318

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

  [复制链接]

2万

帖子

25

威望

3万

eV
积分
65392

管理员

公社社长+计算化学玩家

发表于 2015-1-5 18:29:02 | 显示全部楼层 |阅读模式
本文更新频繁,每次出molclus新版本都会更新。文中评论区的很多讨论都是对于已过时的老版本的,不要作为最新版本的参考。

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

文/Sobereva @北京科音
First release:2015-Jan-5   Last update:2020-Dec-10

前言:molclus(主页:http://www.keinsci.com/research/molclus.html)是北京科音自然科学研究中心开发的一款免费的用于原子/分子团簇构型搜索和分子构象搜索的程序包,目前已被大量知名期刊的理论计算化学方向的研究文章所使用。本文将对此程序包的基本特征和其中的molclus组件的用法进行介绍。第1节介绍基于分子动力学做构型/构象搜索的基本原理,第2节介绍molclus程序的基本使用,第3、4节通过实例进行演示。


1 构型/构象搜索原理

团簇的构型搜索和分子的构象搜索是计算化学经常涉及的问题。团簇构型搜索的目的是找到原子团簇(如Li6)或者分子团簇(如水的五聚体)的最稳定的一批构型,包括坐标和相对能量。构象搜索是对于柔性分子而言的,柔性分子有许多可旋转的键,因此势能面上有很多极小点,构象搜索目的是找到最稳定的一批构象,得知哪个构象能量最低、构象的能量彼此间相差多少。在得到最稳定的一批团簇构型或者分子构象之后,我们还可以根据它们的自由能,进一步根据Boltzmann分布估计出特定温度下它们的分布比率,见《根据Boltzmann分布计算分子不同构象所占比例》(http://sobereva.com/165)。

构型/构象搜索的算法很多,比如势能面扫描、分子动力学、蒙特卡罗、基因算法、LMOD等等,适用场合不同。团簇构型搜索和分子构象搜索都可以基于分子动力学实现,基本过程是:在合适的方法下做一定长度的分子动力学模拟,由于热运动,体系能够反复穿越势垒到达势能面不同区域。每隔一定步数取一次结构,将这些分布在势能面的不同位置上的结构作为初猜结构分别进行几何优化,最终就会得到势能面上各个极小点位置的构型和能量,之后可以对能量进行排序取最低的一批。还可以通过分子动力学做周期性退火来实现构象搜索,效果往往更好。比如每200 ps期间让温度从0 K上升至500 K然后再平缓降到0 K,如果跑2000 ps,就会产生10个对应0 K的结构,对这些结构再进行依次优化,即可得到一批极小点结构和能量。退火过程中每次当体系温度较高时,体系就有较大动能,从而能够较容易地穿越势垒而跑到其它势阱区域,然后温度降低时,体系就会陷入势阱的极小点附近。对于分子团簇的构型搜索和分子的构象搜索而言,并不会涉及到化学键的形成和断裂,因此不需要用极耗时的从头算分子动力学,而只需要计算量极低的分子力场来跑分子动力学就够了。这样的程序很多,比如免费的GROMACS、Tinker、NAMD;Amber实际也算基本免费了(用免费的Ambertools里的Sander模块来跑就行了,不需要用收费的PMEMD部分)。用xtb基于半经验DFT层面的GFN-xTB方法跑动力学耗时也并不很高,而且没有折腾参数和弄拓扑文件的麻烦,还可以在动力学过程中描述键的形成和断裂,故也可以用于原子团簇的搜索。至于结构优化和能量计算部分,可以用分子力场、半经验、DFT、高精度后HF方法等来做,具体用什么就取决于体系特征、需要的精度和运算能力了。

另外,一种简单粗暴但很实用的做分子构象搜索的方法是“系统式搜索法”,即让分子的某些二面角按照指定规则旋转并将所产生的结构做为优化的初始结构。比如设定4个键可以旋转,每个键每120度旋转一次,那么总共就可以产生3*3*3*3=81个初始构象,对这些构象都做优化,就可以得到一批极小点结构,都进行优化后,肯定能找到能量最低构型。molclus程序包中的gentor组件就可以实现这个目的。也有一些专门的构象生成程序如Confab、Frog2,比gentor考虑的因素更多,也更傻瓜化。

对于团簇构型搜索,一种也是看起来很暴力但实际效果不错的做法是把指定类型和指定数目的分子或原子随机挨在一起,形成球形团簇,将这样产生的几十乃至几百个团簇当做初猜结构来依次进行优化,就可以得到一批团簇极小点构型,只要生成的初始团簇数目较多,其中总会有能量最低结构。molclus里的genmer组件就可以实现产生团簇初猜结构的目的。

简单来说,基于molclus做构型/构象搜索所要干的主要事情就是:用户先基于第三方分子动力学程序,或者基于molclus自带的gentor/genmer组件,或者用Confab等构象生成程序,先产生一批初始结构,然后molclus就可以调用量子化学程序或分子力学程序对这些结构依次进行优化(还可以做单点、振动分析等任务),然后用molclus自带的isostat组件对结果进行归簇、排序,从而最终得到能量最低结构和能量较低的一批结构。这个过程看起来没什么复杂的,但molclus将整个过程变得极尽方便省事,而且使用者可以对整个流程进行十分灵活的控制,是想要较低的耗时还是要较高的精度,都可以完全由使用者决定。目前市面上有不少收费的构象搜索程序,比如Spartan、GMMX、HyperChem、Sybyl、MOE、M$里的Conformer、MacroModel之类,我相信,只要你稍微花点时间用明白了molclus,之后一定就再也不想用那些程序了,又不灵活又贵又没法达到像样的精度。


2 molclus程序

2.1 简介

molclus主页为:http://www.keinsci.com/research/molclus.html,程序的可执行文件可在此免费下载。molclus会不断改进、更新。

如果molclus在你的文章中被使用,必须引用,基本格式为(可根据期刊要求修改):
Tian Lu, molclus program, Version x.x, http://www.keinsci.com/research/molclus.html  (accessed 月 日, 年)
这里的月、日、年就是你下载Molclus程序的日期。笔者在未来会专门撰写介绍Molclus的文章发表,届时请引用论文(发表后会在molclus官网和此帖中提及)。

在通过分子动力学程序,或者molclus自带的gentor/genmer组件,或者Confab等程序产生记录了大量初始结构的.xyz格式的轨迹文件之后,molclus就可以调用计算程序对其中的结构进行优化或者能量计算了。molclus目前支持调用以下程序:
• Gaussian:支持G09及之后版本。这是最最常用的做量子化学计算的程序,收费。Gaussian也支持UFF/AMBER/Dreiding分子力场计算,比半经验方法又快出两数量级乃至更多
• ORCA:网址https://orcaforum.cec.mpg.de/。免费。这是第二流行的量子化学程序。由于其RI技术利用得很好,在纯泛函、双杂化泛函下比Gaussian快得多得多,还支持又非常精确又能算得动中等大小体系的DLPNO-CCSD(T)。
• MOPAC:支持2012及之后的版本,网址http://openmopac.net。免费。专门做半经验计算,它支持的PM6-DH+或PM7对有机体系通常足矣得到定性合理的结果,而速度比DFT快两个数量级,故很适合用于中、大体系的优化
• xtb:做GFN-xTB方法计算的程序,免费。此方法相当于半经验的DFT,精度整体强于各种半经验方法,特别是普适性更好,而且xtb速度比MOPAC更快,在此文中有安装方法和基本用法介绍:《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》(http://sobereva.com/421
• Open Babel:网址http://openbabel.org/wiki/。免费。此程序主要用于转换化学领域的文件格式,但也可以实现MMFF94高精度有机分子力场以及其它一些力场的计算

可以说,molclus支持的Gaussian、ORCA、MOPAC、xtb和Open Babel程序把不同尺度、不同精度范围都覆盖了,精度&耗时的基本关系是:高级别后HF > DFT > 半经验/GFN-xTB > 分子力场,这些方法都可以由上述程序实现。如果读者不熟悉ORCA和RI,建议参看《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》(http://sobereva.com/490)当中的介绍,如果不熟悉MOPAC,可参看本文中的相关信息:《大体系弱相互作用计算的解决之道》(http://sobereva.com/214)。分子团簇和弱相互作用关系极为密切,必须选择合适的计算级别,否则团簇搜索将毫无意义。不熟悉弱相互作用计算的话建议看看《乱谈DFT-D》(http://sobereva.com/83)、《谈谈“计算时是否需要加DFT-D3色散校正?》。(http://sobereva.com/413)。


2.2 使用流程

以下是molclus使用的流程框图:

molclus.png

Molclus通常有以下使用方式:
• genmer结合molclus:用于团簇分子或原子团簇构型搜索。使用简便
• gentor结合molclus:用于分子构象搜索。使用简便,但不支持环状区域构象搜索
• Confab或Frog2构象生成程序结合molclus:用于分子构象搜索。使用最为傻瓜化,但不支持环状区域构象搜索、不支持普通有机分子以外的情况
• 经典力场分子动力学程序结合molclus:用于分子团簇搜索和分子构象搜索。使用稍繁琐,被动力学模拟的体系必须有恰当的力场
• xtb程序跑动力学结合molclus:普适性极强,用于构象搜索、原子/分子团簇构型搜索皆可,使用容易(不受力场约束,不需要准备拓扑文件)


molclus的基本使用流程如下:
(1) 产生记录所有初始结构的多帧traj.xyz文件。通常有这些做法:
● 用分子动力学程序模拟团簇或单个分子,每隔一定步数将结构写入一次轨迹;也可以做周期性退火,将对应0 K的帧写入轨迹。令轨迹名为traj.xyz。(如果程序不支持输出xyz轨迹格式,可用免费的VMD程序将其私有轨迹格式转化为.xyz格式的轨迹)。
● 用molclus程序自带的gentor产生含有一批分子初始构象的traj.xyz文件,或者用genmer工具产生含有一批团簇初始构型的traj.xyz文件。
● 用Confab直接产生含有分子的一批随机的构象的traj.xyz文件。
(3) 对于让molclus调用Gaussian、ORCA或MOPAC的情况,需提供这些量子化学程序的模板文件并放到当前目录下。对于Gaussian程序,名为template.gjf;对于ORCA,名为template.inp;对于MOPAC,名为template.mop。模板文件里要写好调用这些程序计算时要用的关键词,如理论方法、基组、相关关键词(内存分配量、辅助SCF收敛的设定、优化的关键词和选项等),但是坐标部分写成[GEOMETRY]。molclus计算时会自动将[GEOMETRY]字段替换为traj.xyz里的结构,然后传递给这些量子化学程序执行。对于调用xtb和Open Babel的情况,不需要设定模板文件,因为这俩本身就是通过命令行调用的。
(4) 设定好molclus的参数文件settings.ini指定任务类型和运行设置,详见2.3节。
(5) 启动molclus。molclus会载入traj.xyz、settings.ini和模板文件,然后自动调用量子化学或分子力学程序对traj.xyz的每一帧结构进行优化,并将优化好坐标和能量写入到当前目录下的isomers.xyz文件中(molclus很灵活,也可以对每帧结构做比如优化+振动分析+高精度单点能计算,这种情况molclus会把坐标、自由能、自由能的热校正量和虚频数目写入到isomers.xyz中,也同时输出在屏幕上)
(6) 启动molclus中附带的isostat工具,载入上一步得到的isomers.xyz,此工具会对其中的结构根据能量进行排序并进行统计、去除能量和几何结构相似度过高的重复结构。isostat在当前目录下输出的cluster.xyz文件就是根据能量从低到高排列并且去除了重复结构后的文件,可以载入VMD来方便地观看每个结构。之后还可以根据提取的各结构的能量,输出指定温度下各个构型的波尔兹曼分布比例(原理见《根据Boltzmann分布计算分子不同构象所占比例》http://sobereva.com/165)。isostat支持并行计算来大大降低耗时,在启动时可以设置并行核数。

如果想使用比结构优化更高级别的方法精确计算每个结构的能量,可以把settings.ini里的itask改为1,把模板文件改成单点任务所用的关键词,再把cluster.xyz改名为traj.xyz,然后启动molclus。此时molclus就会调用量子化学程序对其中每个结构再做一次单点计算,可以再次用isostat进行统计筛选。

程序的安装:Molclus本身不需要安装就可以直接启动。在Windows下启动molclus就是双击molclus.exe图标。在Linux下启动就是进入molclus所在目录,输入./molclus(有时需要先用chmod +x ./molclus加上可执行权限)。每次程序运行前,会对当前目录下.out、.arc等量子化学或分子力学程序之前产生的临时文件进行清理,所以如果有重要的临时文件在当前目录下,执行molclus前应当先手动将它们挪至别处。

Windows版molclus包里有isostat.exe和isostat_32bit.exe。如果你的机子是64位系统,一定要用前面的。后面那个纯粹是给仍然在用32bit系统的人用的,对内存可使用量限制较大,对巨大的isomers.xyz进行筛选可能会崩溃,而且不支持并行计算。

对于在超算等特殊情况下使用,必须需要通过提交的方式执行上述量子化学程序,有可能molclus没法自动调用那些程序。为此molclus程序包里还提供了xyz2QC工具来解决这个问题,请参看《将多帧xyz文件转化成量子化学输入文件的工具:xyz2QC》(http://bbs.keinsci.com/thread-12468-1-1.html)。如果你之前没接触过xyz文件格式,建议你看一看《谈谈记录化学体系结构的xyz文件》(http://sobereva.com/477)。另外,有molclus用户针对在集群上通过PBS作业提交方式使用molclus的方法做了介绍,见:http://bbs.keinsci.com/thread-15598-1-1.html

如果你不会装Gaussian,看http://sobereva.com/439;如果不会装ORCA,看http://sobereva.com/451;如果不会装MOPAC,看http://sobereva.com/262。如果不会装xtb,看http://sobereva.com/421。Open Babel在Windows下安装很简单,直接启动安装程序然后一直下一步即可,对于Redhat派系的Linux系统(如CentOS),用yum install epel-release添加EPEL源后用yum install openbabel即可安装。


2.3 参数文件settings.ini

molclus目录下的settings.ini在molclus启动时会被载入,其中记录了运行参数,在下面进行介绍(此文件里的注释也有描述)。实际情况需要改什么就在原文件基础上修改,注意格式别改乱了,用不到的设置也不要自己删掉。通常做一个新任务只需要改两、三个参数就够了。

iprog:molclus调用的量子化学或分子力学程序。1=Gaussian;2=MOPAC;3=ORCA;4=xtb; 5=Open Babel
ngeom:0=处理traj.xyz里每一帧;n=处理前n帧;也可以直接指定考虑的帧号,如2,5-8,12-14(如果只处理比如第6帧需要写6-6)。利用这个可以实现续算等目的
itask:对traj.xyz里每一帧所做的任务类型。0=几何优化;1=计算单点能;2=振动分析;3=优化+振动分析;-1=热力学组合方法计算(仅对Gaussian有效)
注:对于Gaussian、ORCA、MOPAC这样需要设置template模板文件的情况,itask仅决定按照什么任务的形式来从输出文件中提取数据,而所做的任务的关键词需要照常写(如itask=0的时候,模板文件里得照常写优化的关键词opt及相关选项,例如对于Gaussian可以是opt(gdiis,maxstep=5,notrust))。而对于xtb和Open Babel,并不需要设模板文件,itask直接就决定了在调用时做什么任务。
ibkout:1=每处理完一帧结构,把程序的输出文件备份,以便之后检查或根据自己的需要提取其它数据。比如对于第5帧结构,Gaussian的输出文件会备份到gau00005.out,ORCA的会备份到orca00005.out,MOPAC的会备份到MOPAC00005.arc;0=不做备份,因此molclus运行完毕后将找不到计算程序中间输出的文件;2=只备份成功的任务;3=只备份失败的任务
distmax:n=如果结构中任意两个原子间的距离超过n埃则跳过这一帧。做团簇动力学的时候,有可能某个分子跑飞了,对这样的结构优化没有意义,利用这个参数就可以忽略掉此结构
ipause:1=对于Gaussian和ORCA,计算出错(如不收敛)则暂停,以便用户对输出文件进行检查,然后可以按回车继续处理下一帧结构;2=每处理完一个结构后都暂停;0=总是不暂停
iappend:0=新产生的结构会覆盖原有的isomers.xyz文件;1=新产生的结构会续写到当前目录下之前的isomers.xyz文件(如果有的话)末尾,若没有则会新建之
freeze:设置含有优化过程的任务中原子的冻结情况。0代表不做冻结,设成诸如1,4-6,8-10代表对1,4,5,6,8,9,10原子的笛卡尔坐标进行冻结。对Gaussian/ORCA/MOPAC/xtb有效

为了避免有些人糊涂,在这里提供一个表格,明确说明molclus支持哪些任务(itask)和哪些程序的组合。绿色代表支持,棕色代表不支持。绿色里如果有文字的话,说明需要在模板文件里体现相应的关键词

combination.png

----以下是Gaussian专用的参数----
gaussian_path:Gaussian的可执行文件的路径,比如"/sob/g09/g09"、"D:\study\g16w\g16.exe"。也可以只写"g09",前提是程序所在目录已经加入到了PATH环境变量里。如果路径中有空格,则必须有双引号,没空格的话无所谓,后同
[注意:如果是Windows版Gaussian,gaussian_path必须指向诸如g09.exe而不能是g09w.exe,因为后者只是前者的一个图形界面而已。另外,运行前还必须进入“控制面板”-“高级系统设置”-“高级”-“环境变量”,添加GAUSS_EXEDIR环境变量,使之指向Gaussian目录(不是指向可执行文件名),如D:\study\g09w,否则Gaussian没法以命令行方式运行,也因此molclus将没法调用Gaussian。另外,不要用太老的G09版本,比如年代久远的G09 A.01、A.02版可能没法被molclus支持]
igaucontinue:1=如果优化失败,则会基于另一个模板文件template2.gjf对之前最后一帧结构继续优化,如果还失败了会再基于template3.gjf继续优化;0=不做这个处理。igaucontinue=1的主要用处在于给用户更多备选方案。比如可以在template.gjf里写常规的优化关键词(并且步数上限设得比较有限,如80),而在template2.gjf、template3.gjf里使用更有助于收敛的关键词,诸如calcfc、GDIIS、opt(maxstep=x,notrust)等等。这样当第一次优化不收敛的时候,就会自动利用其它模板文件里的关键词进一步优化。关于辅助收敛的方法,在此贴有汇总:《量子化学计算中帮助几何优化收敛的常用方法》(http://sobereva.com/164)。另外,有时候优化过程中途出现三个原子排成直线而报错,靠igaucontinue=1可以实现基于之前的结构继续优化,通常也能进一步优化下去。
energyterm:molclus是从Gaussian输出文件中末尾的archive字段部分读取电子能量的,这个选项用于设定读取哪项。对于HF/DFT计算,这个参数应设为HF=,即读取HF=这个字符串后面紧跟着的那个数(HF/DFT能量)。如果是MP2或双杂化泛函计算,就应当设MP2=;如果是CCSD(T),就应当设CCSD(T)=,以此类推,不懂的话仔细看http://sobereva.com/488
molclus也可以从Gaussian的热力学组合方法输出中读取指定的热力学量,当itask被设为了-1,template.gjf里的关键词比如设成了# G4,如果将energyterm定义为"G4 Free Energy="(此处的双引号必须写),则molclus就会从输出文件里面找这个字段来读取G4方法算的自由能。
对于freq (itask=2)和opt freq (itask=3)任务不需要设这个,因为总是读取的直接输出的自由能(但如果用提供了template_SP.gjf让molclus在这两种任务后自动算高精度单点来得到较准确自由能,则需要设这个来指定从单点任务中读取什么能量,见后文)。
ibkchk:和ibkout的用法一样,但这个控制的是如何备份.chk文件

----以下是ORCA专用的参数----
orca_path:ORCA可执行文件的路径,比如"D:\study\ORCA4\orca.exe"、"/sob/ORCA4/orca"
ibkgbw:和ibkout的用法一样,但这个控制的是如何备份.gbw文件
ibktrj:1=把每次优化后产生的后缀为_trj.xyz的文件改名为orca[帧号].xyz备份出来。这是优化过程的轨迹文件,因此可以用VMD直接观看优化过程
ibkhess:1=把每次振动分析或优化+振动分析后产生的.hess的文件改名为orca[帧号].hess备份出来。这是Hessian文件,以便用户之后以某些方式利用之
mpioption: ORCA并行运算时传递给MPI的额外选项。默认为none,即什么也不设。此选项有特殊的用处,例如在root下通过较新版本OpenMPI运行ORCA,不带--allow-run-as-root选项就运行不了,此时可以设定mpioption= --allow-run-as-root来解决此问题。如果参数里面有空格,应当写成比如mpioption= " '--bind-to core' "

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

----以下是xtb专用的参数----(molclus总是通过"xtb"命令调用xtb,因此不用指定其可执行文件路径)
xtb_arg:xtb程序运行时额外接的参数,要用双引号扩住,比如"--chrg 1 --uhf 1"代表体系电荷为1,alpha电子比beta电子多1个。(虽然自行使用xtb做优化时应加上--opt,但不要写在此处,因为molclus会根据任务类型自动加上)

----以下是Open Babel专用的参数----(以正常方式安装好Open Babel后molclus即可直接调用,不用指定其可执行文件obabel的路径)
obabel_FF:Open Babel中使用的力场。建议用MMFF94,是最准确的有机分子力场,但支持的元素也仅限有机体系所涉及的。目前的Open Babel里还可以用UFF、ghemical、GAFF力场
obabel_param:Open Babel运行时额外加入的参数,比如--steps 2500代表几何优化做2500步。(注意,2.4版之前的Open Babel使用默认的共轭梯度法优化时有bug,应加上--sd改为最陡下降法优化)


2.4 技巧&细节

xyz是一种格式非常简单的可用记录多帧结构的格式,可以直接用文本编辑器手动编辑,因此用户可以自由修改traj.xyz、cluster.xyz和isomers.xyz的内容,从而控制molclus处理哪些结构,或者手动提取特定结构、删除多余的结构等。molclus完全基于.xyz文件,这使得molclus非常灵活。不懂xyz格式的话可以参看《谈谈记录化学体系结构的xyz文件》(http://sobereva.com/477)。

在便宜级别下做优化和振动分析,而在高精度级别下算电子能量,二者相加来得到较准确的自由能是司空见惯的做法,看比如《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》(http://sobereva.com/327)。Multiwfn的振动分析(itask=2)和优化+振动分析(itask=3)任务直接支持这种做法。具体做法是在当前目录下提供template_SP.gjf(template_SP.inp)文件,设置Gaussian(ORCA)高精度单点计算的关键词,格式同template.gjf(template_SP.inp)。之后,当优化或优化+振动分析任务结束后,程序发现当前目录下有template_SP.gjf(template_SP.inp)时就会用其中的信息构成单点任务文件并调用Gaussian(ORCA)计算,然后自动将新算出来的单点能与前面算的自由能热校正量相加作为最终能量输出。几个细节:
(a)如果template.gjf和template_SP.inp同时存在,优先考虑template_SP.inp
(b)opt、opt freq任务用的程序与自动算高精度单点的可以不同,比如iprog=1调用Gaussian时也可以通过在当前目录下提供template_SP.inp来令程序自动调用ORCA算高精度单点
(c)利用template_SP.gjf基于Gaussian自动算高精度单点时,Gaussian用户需要在settings.ini里对energyterm参数恰当设置,用于指定如何读取这个级别的单点任务输出文件里的能量

molclus运行意外原因中断的续算方式:比如traj.xyz里总共有1000个结构要算,算到第180个结构的时候没算完就意外中断了,此时就把settings.ini里的ngeom设为180-1000,把iappend设为1,然后重新启动molclus计算,程序就会把之前没算的第180~1000个结构算完,并且把新产生的结构续写到之前产生的isomers.xyz里。这样,最后得到的isomers.xyz里就有1000帧了,和一次性全跑完得到的精确一样。

molclus可以同时运行多个以试图节约时间。大家都知道支持并行的计算程序的计算耗时不会和调用的核数线性相关,随着核数增加并行效率必然会愈发降低。比如你的机子核数很多,有56核,traj.xyz里有1000个结构需要处理,你希望降低molclus处理这个traj.xyz的总耗时,那么你可以比如把molclus文件夹复制成四份,比如分别叫molc1、molc2、molc3、molc4,分别将这四个目录下的settings.ini里的ngeom设为1,250、251,500、501,750、751,1000。然后开启四个终端,分别进入这四个目录并且运行molclus。等都运行完了,使用这条命令把四个目录下产生的isomers.xyz合并在一起:
cat molc1\isomers.xyz molc2\isomers.xyz molc3\isomers.xyz molc4\isomers.xyz > isomers.xyz
之后用isostat对这样的isomers.xyz产生的统计结果,和只跑一个molclus但是一次性考虑所有1000个结构的情况是完全一样的,而总耗时能降低不少。类似地,如果你有很多台机子可用,也可以这样把一个任务分开在不同机子上计算,然后把不同机子上的isomers.xyz合并到一起。
注:对于分割运行的数目非常多的情况,手动操作麻烦,可以考虑使用ggdh开发的辅助工具《分解合并Molclus任务的小脚本splitMC》(http://bbs.keinsci.com/thread-14361-1-1.html)。
为了灵活,molclus也可以直接指定设置文件和初始结构文件的路径。例如molclus /RAS/Layer.ini代表将/RAS/Layer.ini作为设置文件。再比如molclus lock.ini /RAS/chu2.xyz代表用当前目录下的lock.ini的设定对/RAS/chu2.xyz里的结构进行计算。因此设置文件和初始结构文件并非只能是当前目录下的settings.ini和traj.xyz。

有用户对isostat程序的运行机制感兴趣,这里具体交代一下。isostat使用能量和几何结构偏差作为判断两个结构是否相似的标准,同时小于两个阈值才算相似。几何结构在判断时是基于距离矩阵来比较的(因此对比前isostat并不会对两个结构做对齐(align),也不要求原子间序号顺序一致)。isostat会循环输入的.xyz里的每个isomer。对于输入文件里第一个isomer,或者当前isomer和之前已有的所有簇都不相似,那么这个isomer将被作为一个新的簇,此簇的能量、结构也等同于这个isomer。若当前isomer与之前某个簇相似,那么这个isomer就被认为归入了这个簇,因此这个簇的容量会+1;如果与此同时这个isomer的能量比这个簇的能量更低,那么这个isomer将被作为这个簇的代表,即使用这个isomer的能量和结构作为这个簇的能量和结构。等所有isomer都循环完之后,程序就会输出各个簇的绝对能量,相对于能量最低的簇的能量差(DE),以及这个簇与与它结构最相近的另一个簇的原子间距偏差最大值(DGmin)。显然,最终得到的簇的数目<=isomer的数目。



下面,本文将示例如何将GROMACS或AMBER分子动力学程序、molclus程序和量子化学程序相结合进行团簇构象搜索和分子构型搜索。如果你不打算利用这俩分子动力学程序产生初始结构(其实一般也没必要,因为过程较麻烦),那么可以把下面例子中的动力学部分直接略过,直接参考基于molclus做搜索的部分,之后再看以下文章分别了解怎么用gentor、genmer、xtb、Confab明显更便利地产生初始结构并结合molclus做搜索:
《gentor:扫描方式做分子构象搜索的便捷工具》(http://bbs.keinsci.com/thread-2388-1-1.html
《genmer:生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具》(http://bbs.keinsci.com/thread-2369-1-1.html
《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》(http://bbs.keinsci.com/thread-16255-1-1.html
《将Confab或Frog2与Molclus联用对有机体系做构象搜索》(http://bbs.keinsci.com/thread-20063-1-1.html
如果你虽然想基于动力学程序产生初始结构,但是又不太会用GROMACS、AMBER等基于经典力场的分子动力学程序,那么建议用上面帖子里提到的通过xtb程序跑动力学,方便得多得多。


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

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

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

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

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

评分

参与人数 51eV +204 收起 理由
PC-coming + 4 牛!
YuhangYao + 5 好物!
panernie + 4
qiqi7 + 2 谢谢
不去重蹈 + 3 精品内容
Reminder + 2 谢谢
吾宇倾城 + 3
双黄蛋 + 4
cyx98 + 4 谢谢
JontyLee + 3 谢谢分享
大威天龙 + 3 不明觉厉
dongdong + 4 谢谢
朙天儿 + 5 好物!
exity + 5 精品内容
xxu若白 + 5
kay + 5 精品内容
rice + 3 精品内容
she + 4 谢谢
16aDream + 5 精品内容
teao + 4 精品内容

查看全部评分

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436。3号:764390338,合计8000人,讨论范畴相同,可加入任意其一但不可都加入。
思想家公社的门口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!

117

帖子

0

威望

1126

eV
积分
1243

Level 4 (黑子)

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

197

帖子

0

威望

2747

eV
积分
2944

Level 5 (御坂)

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

142

帖子

0

威望

1161

eV
积分
1303

Level 4 (黑子)

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

758

帖子

6

威望

4203

eV
积分
5081

Level 6 (一方通行)

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

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

2万

帖子

25

威望

3万

eV
积分
65392

管理员

公社社长+计算化学玩家

 楼主| 发表于 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)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436。3号:764390338,合计8000人,讨论范畴相同,可加入任意其一但不可都加入。
思想家公社的门口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!

107

帖子

0

威望

860

eV
积分
967

Level 4 (黑子)

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

2万

帖子

25

威望

3万

eV
积分
65392

管理员

公社社长+计算化学玩家

 楼主| 发表于 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)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436。3号:764390338,合计8000人,讨论范畴相同,可加入任意其一但不可都加入。
思想家公社的门口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!

107

帖子

0

威望

860

eV
积分
967

Level 4 (黑子)

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

296

帖子

0

威望

3488

eV
积分
3784

Level 5 (御坂)

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

80

帖子

0

威望

1384

eV
积分
1464

Level 4 (黑子)

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

2万

帖子

25

威望

3万

eV
积分
65392

管理员

公社社长+计算化学玩家

 楼主| 发表于 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)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436。3号:764390338,合计8000人,讨论范畴相同,可加入任意其一但不可都加入。
思想家公社的门口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!

123

帖子

0

威望

2028

eV
积分
2151

Level 5 (御坂)

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

6

帖子

0

威望

26

eV
积分
32

Level 2 能力者

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

2万

帖子

25

威望

3万

eV
积分
65392

管理员

公社社长+计算化学玩家

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


同样适用。其实这和molclus本身没直接关系。对于一个体系,只要动力学程序能模拟,量化程序能算,那么就都可以用molclus做搜索。
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436。3号:764390338,合计8000人,讨论范畴相同,可加入任意其一但不可都加入。
思想家公社的门口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!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )

GMT+8, 2021-4-11 18:30 , Processed in 0.245144 second(s), 30 queries .

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