注：本文内容对应molclus官网上最新版molclus程序包里的genmer组件
genmer：生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具
genmer: An ultra-convenient tool for generating initial configurations of clusters
First release：2015-Dec-15 Last update: 2022-Mar-3
1 前言
笔者开发过免费的molclus程序用于团簇构型搜索和分子构象搜索，不熟悉的话应必须先看此文搞明白molclus的思想：《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）。简单来说，molclus的主要执行过程就是读取含有多帧结构信息的traj.xyz文件，对每个结构调用Gaussian/MOPAC/ORCA/xtb/OpenBabel进行优化，然后对结果用isostat子程序进行筛选。
记录初始构型的traj.xyz可以通过多种方式产生，比如分子动力学。本文介绍用molclus中的genmer组件来快速方便地生成molclus所需要的初始构型，读者会发现基于genmer来使用molclus做原子和分子团簇构型搜索超级方便和高效！获得genmer的方法就是到molclus主页（http://www.keinsci.com/research/molclus.html）下载最新版本molclus程序，预编译的genmer就在其中。使用genmer/molclus发表文章时请按照molclus主页的说明进行引用。
目前已经有很多采用genmer+molclus的团簇构型搜索文章发表了，比如J. Mater. Chem. C, 7, 11005 (2019)、J. Phys. Chem. A, 120, 4560 (2016)、J. Theor. Comput. Chem., 16, 1750058 (2017)、Phys. Chem. Chem. Phys., 19, 28434 (2017)、Chem. Phys. Lett., 663, 128 (2016)等。
2 原理
genmer可以把任意种、任意数目的单体（可以是原子或者分子）放在一起通过随机方式生成指定数目的团簇初始构型。
genmer生成团簇初始构型的算法是洒家想出来的，又简单又实用，耗时又极低。就是让被加入的单体的中心从体系的原点出发向着一个随机的方向“发射”，发射过程中单体不断移动，当这个单体被移动到和体系已有的原子没有不合理接触（即被加入单体里的原子和当前团簇的原子间的距离都大于相应两个原子的半径和），这个单体就被加入当前团簇；接下来再这样加下一个单体，反复如此。最开始时的团簇就是第一个单体，随着单体不断加入，团簇就会不断长大，当指定数目的单体全都被加进去时这个团簇的生成就宣告结束。由于每个单体发射的方向都是随机的，还可以允许单体自身随机旋转，所以每次用这个算法生成的团簇构型都可以不一样。在生成团簇的时候可以通过设置一些约束条件来控制产生的团簇的大致形状，比如可以设置单体原子允许的X/Y/Z坐标范围（通过xmin/xmax、ymin/ymax、zmin/zmax指定），再比如可以设置单体中心距离体系中心的最大和最小距离（通过rmax、rmin指定），等等。如果沿着某个方向发射单体过程中，在单体移动到合适的加入位置之前就已经违背了约束条件，那么这次发射就失败了，程序会再尝试向着另一个随机朝向发射，反复如此，最大尝试的次数可以通过maxemit指定。如果超过这个尝试次数还没能成功地将单体加入体系，那么这个团簇的产生就宣告失败。
将上述团簇产生过程画图示意的话就是下面这样。青色的椭圆是已经加入团簇的单体，三个箭头示意了三种随机方向，当要被加入的单体按照箭头所示方向运动到图中黄色椭圆的位置时，由于和团簇里现有的原子已经没有不合理的接触了，此时这个单体就成功加入了，其坐标就被写入到团簇坐标里了。
genmer的这个算法既可以生成分子团簇的初始构型，也可以生成原子团簇的初始构型。对于前者，单体是分子，原子半径使用范德华半径；对于后者，单体是原子，原子半径使用共价半径。
假设设定了团簇由10个A、5个B单体构成，程序会先把10个A都加进去，再把5个B加进去，根据如上算法，得到的团簇中10个A就会倾向于在团簇内层，5个B会倾向于在团簇外层。如果想让它们的出现区域不受设定顺序的影响而是完全随机的，genmer也可以做到自动打乱顺序，这称为shuffle。
在产生团簇时每个分子都被当成是刚性的。但如果一种分子有多种可能的构象都需要考虑，可以提供包含多个构象结构的分子的.xyz文件，并且设定每种构象出现的概率。这样每次加入分子时就会按概率抽取一种构象。
虽然genmer产生的团簇初始构型可能看起来略微松散，但这完全没关系，一做优化就自动收缩变得致密了。molclus用的团簇构型搜索的关键是让初始构型分布得足够随机且全面，这样才可能经过优化后获得尽可能多的极小点构型，特别是采样到能量最低结构。只要让genmer产生的构型数目足够多，是肯定可以达到这个目的的。
3 使用方法
genmer的使用非常容易，只需要提供每种分子的几何结构.xyz文件，设定出现数目就行了。如果你不懂什么叫.xyz格式，看此文：《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。运行参数和分子设定都记录在genmer.ini文件中，直接启动genmer后程序会读取当前目录下的genmer.ini，然后把生成的指定数目的初始构型都输出到当前目录下的traj.xyz当中。之后就可以直接运行molclus来自动调用其它程序去优化其中的每个结构了。对于Windows版，双击genmer.exe图标即可启动，对于Linux版，通过输入genmer可执行文件的路径即可启动。PS：实际上设置文件也不是必须叫genmer.ini且放在当前目录下，genmer启动时如果在当前目录下找不到genmer.ini的话会让你直接输入其文件路径，因此设置文件也可以命名成容易分辨的名字。
参数文件genmer.ini的每个选项在此文件里都有注释，一看就懂，下面介绍一下。通常大家就在程序自带的genmer.ini基础上编辑即可，如果此文件里有的参数没有出现，就会用默认值。
ngen= 50 // 期望生成的团簇构型的数目。默认为1
ngenmust= 0 // 若>0，则ngen选项无效，genmer将无限循环不断产生构型，直到成功产生的构型达到ngenmust。此选项默认为0，代表不做这个要求
step= 0.1 // 设定每一步移动单体的距离。对于单体数少的情况不用改这个，但如果要产生的团簇里单体数甚巨，想节省时间可以把这个步长值设大，但出现不合理接触的可能性会增加。默认为0.1，一般不需要改
ishuffle= 1 // 0: 不打乱加入单体的顺序 1: 打乱顺序（如果当前只有一种单体的话，此选项设0还是1都不影响结果） -N：前N种单体加入顺序不变，只打乱其余种类分子加入顺序。 默认为0
imoveclust= 1 // 1: 将产生的团簇的中心挪到原点，即(0,0,0)的位置 0: 不挪动。默认为1
irandom= 1 // 0: 输入文件相同的情况下让每次运行结果都相同 1: 让每次运行结果都不同（因为此时每次运行时会使用不同的随机数种子）。默认为1
iradtype= 1 // 0: 使用CSD共价半径来检测原子间接触，适合生成原子团簇 1: 使用Bondi范德华半径检测原子间接触，适合生成分子团簇。默认为1
radscale= 1.0 // 原子半径乘的系数。显然设得越小，单体之间的原子间距离会越近；设得越大，单体之间的原子间距离会越大。默认为1.0，一般不需要改
maxemit= 300 // 生成每个团簇时尝试发射单体的最大次数。默认300，一般不用改。但当约束条件设得比较严、比较复杂而不容易成功产生团簇时，可以加大这个参数增加成功几率
---- 这一行必须有且包含"----"字符，告诉程序代表接下来的内容是每类单体的数目（可以附带选项）和.xyz结构文件，单体类型可以有无数多个。
4 //第1类单体数目
examples/H2O.xyz //第1类单体的结构文件。可以写绝对路径，也可以写这样的相对路径
3 //第2类单体数目
examples/ethanol.xyz //第2类单体的结构文件
...略
在单体数目后头可以加入以下选项来设置加入此类单体时的特殊设定：
norot：单体不允许随机旋转。默认情况下每次加入单体时都是让单体随机旋转的，以让团簇构型间尽量有差异性
noxrot、noyrot、nozrot：分别让单体不能绕X轴、绕Y轴、绕Z轴旋转。也可以同时写多个，比如同时写noxrot noyrot则体系就只能绕Z轴随机旋转。三个同时写的话效果和写norot相同
cenxyz [XYZ坐标]：定义单体的中心位置，坐标是相对于单体结构文件的坐标。默认情况下是以单体几何中心作为单体的中心位置
rmax [距离]：设定这个单体的中心与体系中心最远距离不能超过多少。默认不做限制
rmin [距离]：这个选项可以要求单体中心与体系中心距离不得小于多少，这实际上是发射单体时与体系中心的最初距离。默认为0
xmax [数值]：要求单体里的所有原子的x坐标不允许大于这个值。默认为不限制
xmin [数值]：要求单体里的所有原子的x坐标不允许小于这个值。默认为不限制
ymax [数值]：要求单体里的所有原子的y坐标不允许大于这个值。默认为不限制
ymin [数值]：要求单体里的所有原子的y坐标不允许小于这个值。默认为不限制
zmax [数值]：要求单体里的所有原子的z坐标不允许大于这个值。默认为不限制
zmin [数值]：要求单体里的所有原子的z坐标不允许小于这个值。默认为不限制
显然，这些选项对于精确控制团簇构型很有用，非常灵活。比如用cenxyz可以直接对某些单体进行定位；rmax可以令单体构成球形团簇；rmax和rmin联用可以实现构建空心团簇；zmax和zmin联用可以产生扁片状团簇；rmax、rmin、zmax、zmin一起联用可以生成环状或者桶装团簇。
以下是单体设置中加入上述选项时的格式示例
3 norot cenxyz 0 1.5 0.2
examples/CH3CH2CH2NO2.xyz
20 rmax 6
examples/H2O.xyz
在输入文件末尾还可以加入一行$distcons，然后在下面通过一行或多行，指定特定元素间必须满足的距离范围（没设的话就不做限制），格式是：[元素名1] [元素名2] [下限] [上限]，上限值可省略，具体看4.13节的例子。例如写H Mg 2.0就代表氢和镁之间距离不能小于2.0埃，Pt Pt 1.8 3.2代表铂之间距离必须在1.8~3.2埃，Na N 0 4.5代表Na和N之间距离必须小于4.5埃。若体系中任意一对原子违背了特定的元素间约束条件，则构型会被视为失败而不输出。利用这个设定，可以尽量避免已知不该成键的两种元素在优化后成键，以及避免本该成键的两种元素由于一开始离得太远而优化后难以出现成键。
在输入文件末尾还可以加入一行$modrad，然后在下面通过一行或多行，对一个或多个元素半径进行修改，比如Pt 3.5就把半径改为了3.5埃。
4 用genmer产生团簇初始构型例子
下面给出几个genmer的例子展现其生成团簇的能力，用到的单体结构文件都是genmer自带的。为了简洁，基本上所有默认的参数都舍掉了。以下图像都是把genmer产生的traj.xyz文件直接拖进VMD程序（http://www.ks.uiuc.edu/Research/vmd/）显示的。分子单体的结构都是用GaussView建的，然后再按照http://sobereva.com/477里说的方式转换为xyz格式。显然，这一节直接用genmer产生的结构不是有实际意义的团簇结构，仅仅是产生给molclus用的初猜结构而已。必须效仿第5节的例子，结合molclus进行依次优化、筛选才能得到最终有实际价值的团簇结构。
4.1 6个水团簇
用以下genmer.ini文件可以得到由6个水组成的20个团簇初始构型
ngen= 20
----
6
examples\H2O.xyz
结果如下：
4.2 3000个水团簇
用以下genmer.ini文件生成1个多达3000个水组成的团簇初始构型，并且水距离体系中心不能超过40埃。由于分子数较多，生成需要一定时间，大约十秒左右。注意rmax不能设得太小，否则会导致有的水在各个方向尝试加入均不成功，最后宣告失败（由于算法有一定随机性，可能这次没生成成功，但下次就成功了）。
ngen= 1
----
3000 rmax 40
examples\H2O.xyz
结果如下：
4.3 6个水和3个乙醇的团簇
用以下genmer.ini文件可生成6个水和3个乙醇构成的10个团簇，并且打乱乙醇和水的出现顺序
ngen= 10
ishuffle= 1
----
6
examples\H2O.xyz
3
examples\ethanol.xyz
结果如下：
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内容为
ngen= 1
----
30
examples\butadiene_cis_trans.xyz
结果如下：
4.5 多巴胺二聚体
用以下genmer.ini文件可以生成15个多巴胺二聚体初始构型，这算是稍微大一些的分子了。
ngen= 15
----
2
examples\dopamine.xyz
结果如下：
4.6 水层包围多巴胺
用以下genmer.ini文件可以生成200个水分子紧密包裹在多巴胺附近的10种构型。由于ishuffle设为了0（其实这本身也是默认的），而且多巴胺的结构文件先出现，因此多巴胺被包围在内侧。水的rmax不能太大也不能太小，为了产生接近球形且较为致密的包围多巴胺的水球，可以先尝试比较大的rmax，然后逐渐缩小，直到程序无法成功产生团簇结构为止。
ngen= 10
ishuffle= 0
----
1
examples\dopamine.xyz
200 rmax 15
examples\H2O.xyz
结果如下：
前面说过，如果ishuffle设-N，则前N类分子加入的时候还是按顺序加入，而其余的分子则随机顺序加入。这对于构建很多体系很有用，比如下面的例子令多巴胺出现在中心，乙醇和水混合出现在多巴胺附近（水和甲醇同步加入，不分先后次序）。
ngen= 10
ishuffle= -1
----
1
examples\dopamine.xyz
50 rmax 15
examples\ethanol.xyz
100 rmax 15
examples\H2O.xyz
4.7 水分子与多巴胺的氨基形成氢键二聚体
以下genmer.ini生成一个水与一个多巴胺形成的二聚体，共生成20个，要求水必须与多巴胺的氨基接触，目的是优化后可以得到水与氨基通过氢键形成的二聚体。
ngen= 20
imoveclust= 0
----
1 norot cenxyz 4.33709192 0.16551438 0.29362167
examples\dopamine.xyz
1 rmax 4
examples\H2O.xyz
注意这里对多巴胺用了norot，这样它在团簇里的朝向就和单体结构文件里相同了，而不会随机变化朝向。cenxyz后面的值正是dopamine.xyz里氨基的氮原子的坐标。由于多巴胺是第一个出现的分子，而且其中心被设定在了氮原子上，因此在产生的团簇里，多巴胺的氮原子恰好就是在团簇坐标的原点的位置。由于rmax 4让水的中心（默认是水的几何中心）不能与团簇原点距离超过4埃，所以保证了产生的团簇中水肯定在氨基附近。此例用了imoveclust= 0，使得生成的团簇的笛卡尔坐标原本是什么就在traj.xyz文件里记录什么，而自动不把整个团簇的几何中心在写入traj.xyz之前挪到(0,0,0)的位置，这使得traj.xyz中多巴胺的位置是固定的。如果你把这个例子弄懂了，就差不多算是完全理解genmer的算法了。
运行结果如下。为了看得清楚，这里把所有20个团簇结构叠加在了一起显示，可见水很好地包围了氨基，各个角度都有出现。如果你打开traj.xyz，会看到氮原子的坐标恰好是(0,0,0)。
4.8 Li6团簇
用以下genmer.ini文件可以生成6个Li组成的团簇的15个初始构型。注意这里用了共价半径
ngen= 15
iradtype= 0
----
6
examples\Li.xyz
结果如下：
4.9 环状的Li30团簇
用以下genmer.ini文件可以生成10个由30个Li组成的环状团簇的初始结构。这些Li只允许出现在径向距离r=5埃至9埃的范围内，并且要求z坐标必须在-1~1埃之间，因此得到的就是螺丝垫片形状的团簇了。由于当前设的约束条件比较多，Li只能出现在相对比较狭小的限定的空间中，为了增加产生的成功几率，这里将尝试发射单体次数maxemit设的比默认值300大得多得多（凡是像这种对单体出现的空间区域要求的比较严的情况，都建议改大maxemit。如果设得很大还是不能成功产生，那要么减少加入的单体数，要么把对空间范围的约束放宽点）
ngen= 10
iradtype= 0
maxemit= 1000
----
30 zmin -1 zmax 1 rmin 5 rmax 9
examples/Li.xyz
结果如下（为了看着更清楚，在VMD里显示的时候令间距小于特定距离的Li-Li显示出了键）
4.10 笼状的碳团簇包夹水分子
用以下genmer.ini文件可以生成100个碳组成的笼状团簇，其中包夹一个水分子。由于当前允许碳出现的空间范围很窄，即r=4~5埃，所以maxemit设得较大，但即便设为3000，还是有个别团簇产生失败，当前要求产生10个，实际成功产生了7个。如果把maxemit设得更大，成功产生的会更多。这里C.xyz是只含一个碳原子的xyz文件。
ngen= 10
iradtype= 0
maxemit= 3000
----
100 rmin 4 rmax 5
C.xyz
1
examples/H2O.xyz
结果如下：
4.11 CB6主体分子包夹客体分子乙醇
用以下genmer.ini文件可以生成环状主体分子CB6包夹乙醇分子的10个构型。CB6的中心被要求固定在0,0,0位置不动。由于CB6里的孔洞很窄，为了让生成容易成功，这里将原子半径设为了默认的0.8倍。由于乙醇如果不出现在主体分子孔洞内的话就没意义了，所以用了rmax 2让它的中心偏离CB6中心（即0,0,0）不得超过2埃，否则有可能会出现乙醇结合到CB6外部的结构。
ngen= 10
radscale= 0.8
----
1 norot cenxyz 0 0 0
CB6.xyz
1 rmax 2
examples/ethanol.xyz
结果如下
4.12 18碳环平行堆叠的团簇
18碳环是个很有意思的体系，笔者做过大量理论研究，见http://sobereva.com/carbon_ring.html。假设18碳环的xyz文件在当前目录下，且此文件中18碳环平行于XY平面，genmer.ini里通过下面的设置就可以产生一个分布接近球形的10个18碳环平行堆叠的团簇结构。
10 noxrot noyrot rmax 10
C18.xyz
得到的结构之一
4.13 4个Mg和8个H的团簇
下面这个输入文件可以产生4个Mg和8个H组成的团簇，共产生10个。这里用共价半径，但是把Mg的半径改大到了2.0埃（没特殊意图，只是作为演示）。为了避免H和H之间太近可能优化后马上变成氢气，所以把H和H最小距离设为了2埃。
ngenmust= 10
iradtype= 0
----
4 rmin 2.5 rmax 5
Mg.xyz
8 rmax 5.5
H.xyz
$modrad
Mg 2.0
$distcons
H H 2.0
得到的结构之一：
5 结合molclus做团簇构型搜索
上一节演示了用genmer得到了各种团簇的traj.xyz文件，这一节我们基于这些团簇初始构型，用molclus做构型的批量优化、获得能量，从而筛选出较能量较低，也因此较合理的结构。如果下文的例子没看懂，应当先把前述的《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）一文仔细看一遍了解molclus的基本用法。由于genmer每次产生的团簇构型可能不同（取决于irandom），因此你尝试重复下文的例子的话，结果可能和下文有所不同。
5.1 多巴胺二聚体例子：molclus结合MOPAC做半经验计算
这里我们用molclus调用MOPAC在PM6-DH+级别下优化多巴胺二聚体来得到它们的团簇构型。PM6-DH+是一种适合计算弱相互作用的半经验方法，虽然定量结果算不上精确，但是对弱相互作用体系给出定性正确的结果是足够的。
把前面4.5节得到的含有15个多巴胺二聚初始构型的traj.xyz文件放到molclus目录下，把molclus的settings.ini里的MOPAC路径配置好，将iprog设为2，并确认MOPAC的模板输入文件template.mop里的关键词是PM6-DH+ precise。启动molclus，程序自动调用MOPAC对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个结构能量很低是很显然的，因为两个氨基都各自和两个羟基同时形成了氢键。第2个结构是每个氨基和一个羟基形成氢键。第3个结构形成了明显的pi-pi堆积，而且两个氨基之间有氢键。第4个结构有两对氨基与羟基形成的氢键，而且还有一点pi氢键的特征。第5个结构只有一对氢键，但是pi氢键的特征更明显。
由于本例只考虑了15个多巴胺二聚体结构，不算很多，而二聚体极小点构型有很多，还不能确保二聚体最稳定构型就在其中，要想更稳妥应该考虑更多结构。
另外，半经验方法精度终究比较有限，获得能量接近的结构的能量差定量不准确，甚至顺序也不对。如果要得到准确的二聚体结构和能量，建议在PM6-DH+优化完之后用molclus调用Gaussian在B3LYP-D3(BJ)/6-311G*级别下对各个结构再优化一下，然后用molclus调用Gaussian在M06-2X-D3/6-311+G(2d,p)级别下（要求更精确可以用比如B2PLYP-D3(BJ)/jun-cc-pVTZ）计算能量最低的数个结构的单点能。这种“半经验预优化->DFT结合一般基组优化->高精度方法算单点能"是我很推荐的搜索分子团簇的做法，这用molclus实现极其容易，可参考《使用molclus程序做团簇构型搜索和分子构象搜索》中计算水四聚体团簇构型的例子。
5.2 乙醇+水团簇例子：molclus结合Gaussian做分子力场计算
这一节用Gaussian里的UFF力场优化4.3节得到的6个水+3个乙醇的团簇的10个初始构型。力场计算极其便宜，可以比半经验优化更大体系、更大批量构型。（由于MOPAC和xtb速度很快，其实对于当前这样不大的体系，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初始产生的结构对比，可看到经过力场优化后团簇变得明显紧凑了。之后我们可以再让molclus调用量子化学程序进一步优化、计算能量获得更准确的结果。
5.3 Li6团簇例子：molclus结合Gaussian做DFT计算
这一节我们用molclus调用Gaussian在B3LYP/6-31G*级别下寻找Li6团簇稳定构型。
把4.8节得到的包含15个Li6团簇初猜构型的traj.xyz放到molclus目录下。将Gaussian的模板输入文件template.gjf里的关键词填入B3LYP/6-31G* opt(cartesian,maxcyc=50)，备用模板文件template2.gjf里的关键词填入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中。从结构看会发现前两个结构其实是一样的，之所以能量有细微差别其实是收敛精度问题而已。所以实际上有三个唯一的团簇结构，如下所示，从左到右能量依次升高
这些团簇从图上看都有对称性，而当前优化过程由于收敛精度原因可能最后一步的构型不严格满足实际对称性，故建议自行在gview里把结构做对称化后再优化一下以得到更准确结果（把对称化后的坐标写成traj.xyz的格式，然后用molclus再跑一次就行了）。
5.4 环状Li30团簇例子：molclus结合xtb做半经验级别的DFT计算
在4.9节中提到了怎么产生环状Li30的初始构型。xtb程序（基本使用介绍见http://sobereva.com/421）速度极快，耗时比MOPAC还低，而且xtb程序实现的GFN-xTB方法比MOPAC支持的那些半经验方法普适性好得多，因此明显更适合处理电子结构相对来说比较复杂的原子团簇体系。为了对Li30构型进行快速搜索，在这一节我们用molclus调用xtb对已产生好的那10个Li30初始结构进行优化。之后，如果你想要更准确、更可靠的结果，可以再调用Gaussian或ORCA在更准确的DFT下优化和算能量。
把molclus的settings.ini里的iprog设为4，并且在xtb_arg选项里加上--gfn 1，这代表让xtb使用GFN-xTB方法的1代。这点很重要，目前的（2019-Mar-18发布的）的xtb程序默认用的是GFN2-xTB，即GFN-xTB方法的2代，虽然比1代在原理上更好，但是碰见复杂的原子团簇很容易出现SCF不收敛，因此得强制要求用1代。
启动molclus，在普通4核机子下xtb仅仅几秒钟就可以优化完一个结构，xtb速度超级快！！！然后我们照常用isostat对结果排序，其中能量最低的一个结构如下所示
这个结构非常理想，特别有序，凭直觉就能知道应该是个很稳定、很有意义的结构。此例反映出molclus+xtb的组合对于原子数较多的原子团簇的搜索也很合用，可谓黄金组合！
上面是研究Li30的环状团簇，笔者还用molclus结合xtb程序对30个Li构成的笼状团簇做了构型搜索，genmer产生的团簇初始结构和优化后的结构如下所示，细节详见http://bbs.keinsci.com/thread-12509-1-1.html，如图可见得到的结构意想不到地好！genmer一开始给出的初始结构比较乱，但是经过神器xtb一优化，就变得非常有序、漂亮了。
6 总结
从以上介绍和例子可见，用genmer+molclus的方式做团簇构型搜索虽然思路简单，没有用到什么高深的技术，但是非常实用，结果理想，操作很简便，对使用者的知识要求很低，完全不用掌握什么复杂的理论知识，强烈推荐大家在实际研究中使用，并且向同行推广！
有一些用户问让genmer产生多少初猜团簇比较合适，这取决于体系大小、特征、研究目的、可接受的计算花费。如果你想得到能量最低的团簇结构，那么显然产生的初猜结构数目越多，最终能获得能量最低结构的几率就越大。对于结构比较简单的小分子构成的二聚体团簇，想获得能量最低的构型，一般产生20个初猜结构就够了。但是如果是复杂、较大的柔性分子构成的二聚体团簇，或者简单小分子构成的比如六聚体团簇，那要产生的初始构型就得很多，起码一百个。对于原子团簇，也是原子数越多，要产生的初猜数就应该越多。受限于用户的实际计算条件，考虑的初猜结构数目得量力决定。
|
[复制链接 Copy URL]
发表于 Post on 2015-12-15 12:18:02
|
只看该作者 Only view this author
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
发表于 Post on 2015-12-15 15:22:29
楼主 Author
