注：每次molclus程序包里的gentor更新时此文也都会同步更新。本文内容始终对应于molclus官网上的最新版本的情况。
gentor：扫描方式做分子构象搜索的便捷工具 gentor: A convenient tool for molecular conformation search using scanning method
文/Sobereva @北京科音 First release: 2015-Dec-18 Last update: 2020-Nov-11
1 前言
笔者开发过molclus程序用于团簇构型搜索和分子构象搜索，不熟悉的话请务必先看程序官网http://www.keinsci.com/research/molclus.html以及《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html），否则无法理解本文的内容。molclus简单来说就是读取含有多帧结构信息的traj.xyz文件，依次对每个结构调用量子化学或分子力学程序进行优化，然后对结果用isostat子程序进行统计、处理。
给molclus用的traj.xyz可以通过第三方程序如xtb、GROMACS、Amber等做分子动力学产生，然而很多人不会操作。本文介绍笔者开发的gentor程序来让molclus能以扫描方式做构象搜索。gentor可以对指定的化学键按一定规则扭转产生一批分子初始构象，然后写入traj.xyz文件中，再使用molclus调用计算程序如Gaussian对之依次优化后，就能得到各个极小点结构和能量最低结构。完整读完本文后读者会发现以这种方式做构象搜索使用起来相当方便、可靠，而且精度和耗时之间的关系高度可控。gentor是molclus的一个组件，去molclus主页（http://www.keinsci.com/research/molclus.html）下载最新版本molclus程序，预编译的gentor就在其中。
gentor这种系统式扫描方式产生初始构象相比通过动力学来产生有优有劣。优点是省得跑动力学了，方便很多，而且对势能面采样得较为均匀全面，不容易漏掉极小点。缺点就是，产生的初猜构象的数目随着可旋转的键的增加呈几何级数增长。比如每个可旋转的键扫3个点，考虑5个键就得优化243次，有10个键就得优化高达59049次，此时已经几乎没法处理了。所以对于小分子用系统式扫描很好，但体系大了的话还是做动力学产生初猜结构为宜，因为这样得到的初猜结构大多分布在感兴趣的低能构象附近，最终找到能量最低构象的可能性较大。如果用退火方式则更好，比如GROMACS中可以使用周期退火，每100 ps内从0 K升温到500K再返回0 ps，跑50 ns就能得到500个处于0K的构象，用molclus对它们优化，得到能量最低结构的几率就会很大了。
后记：值得一提的是笔者后来又另写了一篇构象搜索文章《将Confab或Frog2与Molclus联用对有机体系做构象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）。如果你缺乏化学直觉，不知道哪些键应当以何种方式考虑旋转；或者gentor产生的初始构象太多，全都考虑的话太昂贵，而你又不会用或不想用动力学程序产生初始构象，那么用Confab或Frog2构象生成程序来产生初始构象是不错的选择，比较傻瓜化，而且可以一定程度上控制产生多少初始构象。但用这种构象生成程序只适合有机体系，而且没法做到像gentor那样精确控制构象生成细节。
2 使用方法
gentor的使用非常容易。需要在当前目录下准备两个文件：
mol.xyz：记录分子结构的.xyz文件
gentor.ini：设定要旋转的键以及旋转方式
然后启动gentor就会在当前目录下产生traj.xyz了。
下面是个gentor.ini的例子，考虑了4个键的旋转，囊括了各种gentor支持的键旋转的定义方式（旋转的度数是相对于输入文件里结构里的二面角而言的，0度相当于不旋转）
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。也可以用比如e-120来依次旋转0,-120,-240度。e是every的缩写
15-16 //绕着15-16键令16的那一头的片段进行旋转
e5 70 90 //从70度开始每5度转一次，即依次旋转70,75,80,85,90。也可以用比如e-10 120 90来依次旋转120,110,100,90度
7-8 //绕着7-8键令8的那一头的片段进行旋转
-120 120 //两个方向各旋转120度
9-12 //绕着9-12键令12的那一头的片段进行旋转
r8 //旋转0~360度内的随机8个角度。r是random的缩写
可以定义无数多个键的旋转。初始产生的构象数目是所有指定的键的旋转次数的乘积。比如有四个键分别旋转3、6、2、8次，一共就会产生3*6*2*8=288个构象。
在gentor.ini里还可以定义同步旋转。比如
2-24
s 2
6-10
e30
4-17
s -2
代表绕着6-10键每30度旋转一次，每次旋转时2-24键也与之旋转相同的度数，4-17则与之旋转相反的度数，因此总共产生360/30=12个结构。s是synchronous的缩写，后面接的是要同步旋转的键在整个.ini文件里出现的序号。由于6-10是整个文件里第二个键，因此2-24的下面s后面是2。由于4-17下面的s后面接的是负数，所以是与6-10相反方向同步旋转的。对于三乙基苯，如果这三个键对应的是乙基与苯环相连的C-C键，那么得到的traj.xyz播放起来就会如下所示，可见上方和左方的甲基是顺时针同步旋转，而右侧的甲基是逆时针同步旋转：
gentor.ini里可以设定无数个独立和同步旋转。比如可以让B依赖与A同步旋转，且让D与C同步旋转，那么总共产生的构象数目是A的旋转次数乘以C的旋转次数。
产生的初始构象会被输出到当前目录下的traj.xyz里。然后就可以直接运行molclus来自动调用量化程序去优化其中的每个结构，然后用isostat进行统计了。
默认情况下，原子间距离小于它们的共价半径和的1.15倍就被视为成键。但有的时候gentor自动判断的成键方式和我们预期的不符，导致基团没法按照预期方式转动。这种情况下可以在定义键的旋转方式后面空一行写上critbond来设置阈值，之后还可以再进一步自定义特定原子间成键与否。例如
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是最终的原子间成键关系。
通过旋转键产生的构象很可能会有不合理接触，因此gentor会对产生的结构进行检测，如果有任意两个原子间的距离小于二者的共价半径和的0.75倍，则这个构象将不写入最终的traj.xyz。但0.75这个阈值也未必总是最合适，因此gentor也允许用户自己改这个阈值，也就是在设置键的旋转方式最后空一行写上比如critskip=0.7（如果也设了critbond则必须在critbond之前出现），则默认的0.75倍就被改为了0.7倍了。
再强调一下，要旋转的键一定要结合化学直觉精挑细选，旋转的方式也要恰当设定，毕竟总计算量会随着被指定为可旋转的键的数目增加而呈几何级数增加，因此瞎设的话，产生的初始构型数目会超级多（动辄几百万、几千万），此时根本没法算。只有那些真正可以旋转的键才应当写进gentor.ini里，而诸如环当中的键或者双键，根本转不动，一转要么分子被撕裂，要么发生构型变化而不是原先的分子，因此都不能被作为可旋转的键。而且，有很多键转不转对体系构象影响很小，特别是比如甲基，是否旋转都看不出差异，因此显然不应当被设为可旋转的键。通常来说，可旋转的键的数目最多不超过8个，如果确实有更多可旋转的键，那么建议读者还是用分子动力学程序通过做模拟退火来产生交给molclus处理的traj.xyz文件。
3 分子构象搜索实例1：多巴胺（利用MOPAC+Gaussian）
这里我们用一个很简单的分子多巴胺作为实例展示gentor+molclus做构象搜索的基本流程。分子初始结构如下：
一般情况下，分子结构用GaussView自己画就可以，然后按照此文说的方式转成xyz格式：《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。此体系对应的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，会看到如下输出
Criterion for skipping structure: 0.8000
Criterion for determining bonding: 1.1500
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
9 conformations were originally generated
0 conformations did not pass distance check and thus be ignored
Final number of conformations: 9
Molecular conformations have been outputted to traj.xyz in current folder
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程序（http://www.ks.uiuc.edu/Research/vmd/），可看到以下初始结构，显然都是合理的，可以做下一步了。
我们用molclus调用MOPAC在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，itask设为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下能量最低的结构为
仔细看isostat输出的各个构象的能量的话，会看到1、2之间非常相近。观看cluster.xyz中对应的结构会发现1、2的差异只是-CH2-NH2基团一个冲左一个冲右而已，破坏二者镜像对称的只不过是苯环上那个在间位的羟基而已，显然由于它离-CH2-NH2基团甚远，不会造成1、2之间能量出现显著差异。类似地3、4之间，5、6之间能量也很接近，能量差不到0.1kcal/mol。
如果想要更精准的结果，建议优化在M06-2X/def2-SVP下进行，单点能用双杂化泛函结合cc-pVTZ或def2-TZVPP基组计算。但当体系较大、可能的构象很多时，如几百个甚至上千个，建议先用廉价的半经验，甚至分子力场先做一下粗略优化（通过OpenBabel用MMFF94力场或Gaussian用UFF力场即可），以筛选出来能量较低的、有必要之后用DFT进一步优化的构象。用molclus做这种事极方便，详见《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）。
4 分子构象搜索实例2：阿斯巴甜（利用xtb+ORCA）
注：这个例子用的是gentor较老的版本，对于目前最新版gentor，中间数据有所不同，但最终得到的关键性数据是一样的。
此例的目的是演示怎么用gentor+molclus对一个构象非常多的偏大的柔性体系做构象搜索。本例将要用到xtb和ORCA程序，xtb是2019-Mar-18发布的版本，ORCA是4.1.1。本文涉及的所有输入、输出文件、molclus配置文件（对应1.8.1版）都可以在此下载：http://bbs.keinsci.com/kein/Aspartame.rar
其实用MOPAC和Gaussian来实现此例的目的也完全可以，只不过获得同样可靠程度的结果要花的时间比用xtb和ORCA的组合要高。因为xtb支持的GFN-xTB方法比MOPAC支持的那些半经验方法更可靠，而且xtb运行速度还更快。xtb的基本使用和安装看此文http://sobereva.com/421。而ORCA则支持不少Gaussian不支持的又好又快的理论方法，而且在ORCA里利用RI加速技术，可以使很多理论方法的耗时显著低于Gaussian，ORCA的安装方法看此文《量子化学程序ORCA的安装方法》（http://sobereva.com/451）。
下图是阿斯巴甜（Aspartame）的初始结构。这个结构是直接通过OpenBabel基于经验方法快速产生的，也就是在wiki的阿斯巴甜页面上查到此分子的SMILES字符串O=C(O)C[C@H](N)C(=O)N[C@H](C(=O)OC)Cc1ccccc1，将之写入一个文本文件里并且命名为比如a.smi，然后运行obabel a.smi -O Aspartame.xyz --gen3D，此时产生的Aspartame.xyz就是经验方式给出的很粗糙的阿斯巴甜结构。在GaussView里显示，如下所示：
此体系可旋转的键比较多。此体系里有个类似肽键的区域，由于存在pi共轭，所以不需要考虑C8-N10的旋转。我们把Aspartame.xyz拷到gentor目录下并改名为mol.xyz，把gentor.ini内容写为
2-4
e120
4-5
e120
5-7
e120
5-8
e120
10-11
e120
11-13
e120
11-17
e120
17-18
e180
13-15
e180
运行gentor，程序提示
8748 conformations were originally generated
972 conformations did not pass distance check and thus be ignored
Final number of conformations: 7776
因此最终得到的traj.xyz文件里有7776个结构。
由于要优化的结构极多，所以我们自然而然会想到先用分子力学方法去初步筛选，选出能量最低的几百个结构再用半经验之类方法去进一步优化。不过，用过xtb程序的人都知道，虽然xtb程序实现的GFN-xTB是属于量子化学的方法（某种程度上类似于基于DFT的半经验方法），但是xtb运行速度如飞。实测一下就会发现，其实对当前这个原子数不算特别多的体系，用molclus调用xtb优化每个结构的耗时比起调用OpenBabel做MMFF94力场的优化或者调用Gaussian做UFF力场的优化也高不到哪去，而可靠程度则要高得多。因此，此例我们索性直接用xtb来对这多达7776个结构批量做优化，你会发现耗时绝没有想象得那么恐怖。
将molclus的settings.ini里的iprog设为4设得被调用的程序对应于xtb，把traj.xyz放入molclus文件夹，并确认itask目前为0、机子里已经装好了xtb，然后启动molclus。在笔者的Intel 36核服务器上，平均优化一个结构只需要两三秒钟，几个小时就把这7776个结构全都优化完了。如果你有多台机子可用，想降低耗时的话也可以地把当前任务拆分放在多台机子上分别跑，然后再把结果合并到一起，实现起来极其容易，详见《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）的2.4节。
运行isostat，按回车载入当前目录下刚产生的isomers.xyz，然后输入0.5作为归簇的能量阈值，再输入0.2作为归簇的几何阈值（二者设得比默认的略宽，是因为xtb在默认设定下优化出的结构有很多很相近的，由于收敛精度关系又不是重合得很理想，用稍微宽一点的阈值可以把这些重复的结构归到一个簇里）。之后程序提示当前归簇后簇的数目特别多，问是否只保留能量最低的一批。这里输入50，这样只有能量最低的50个簇的信息会被输出出来，如下所示，并且产生的cluster.xyz里也只有这50个簇对应的结构。
# 1 Count: 26 E= -65.365051 a.u. DGmin= 0.35 DE= 0.00 kcal/mol
# 2 Count: 33 E= -65.363811 a.u. DGmin= 0.21 DE= 0.78 kcal/mol
# 3 Count: 4 E= -65.363779 a.u. DGmin= 0.20 DE= 0.80 kcal/mol
# 4 Count: 23 E= -65.363776 a.u. DGmin= 0.19 DE= 0.80 kcal/mol
# 5 Count: 4 E= -65.363714 a.u. DGmin= 0.21 DE= 0.84 kcal/mol
# 6 Count: 12 E= -65.363444 a.u. DGmin= 0.24 DE= 1.01 kcal/mol
# 7 Count: 2 E= -65.362493 a.u. DGmin= 0.25 DE= 1.61 kcal/mol
# 8 Count: 32 E= -65.362246 a.u. DGmin= 0.26 DE= 1.76 kcal/mol
# 9 Count: 4 E= -65.362225 a.u. DGmin= 0.27 DE= 1.77 kcal/mol
# 10 Count: 2 E= -65.362157 a.u. DGmin= 0.26 DE= 1.82 kcal/mol
...略
# 50 Count: 19 E= -65.359549 a.u. DGmin= 0.22 DE= 3.45 kcal/mol
xtb优化的结构、算出的能量终究比较粗糙，若想严格地找出来能量最低结构，肯定还得在更好的级别下优化几何结构，一般用B3LYP-D3(BJ)/6-311G*这样的档次就够了。但ORCA里有个更好选择，是B97-3c方法，此方法在ORCA里耗时很低（与RI-BLYP/def2-SVP相近），而精度达到B3LYP-D3(BJ)/6-311G*这样的档次不成问题。B97-3c起码对于优化分子结构已经绝对够用了。
我们将cluster.xyz改名为traj.xyz并覆盖掉原先的。上面归簇出来的50个结构没必要都用B97-3c优化，xtb算出来能量特别高的构象改用B97-3c优化后肯定还是对应能量很高的构象。所以优化靠前的（能量最低的）的一批，比如前7个就够了，因此我们将settings.ini里的ngeom设为7。如果你想更稳妥些，优化比如前15个也可以。把settings.ini里的iprog改为3以对应调用ORCA，确认orca_path设的确实是当前机子里ORCA可执行文件的实际路径。然后把ORCA任务的模板文件template.inp内容改成下面的样子（实际调用的核数nprocs、每个ORCA进程调用的内存量上限maxcore请根据机子实际情况设置。笔者用的是36核72线程256GB内存的机子）
! B97-3c opt miniprint nopop noautostart
%maxcore 4000
%pal nprocs 36 end
* xyz 0 1
[GEOMETRY]
*
然后启动molclus，在笔者的机子上几分钟就优化完一个。之后再次启动isostat，还是分别输入0.5和0.2来自定义归簇的能量阈值和归簇的几何阈值，显示的信息如下
# 1 Count: 1 E= -1029.559646 a.u. DGmin= 0.29 DE= 0.00 kcal/mol
# 2 Count: 1 E= -1029.558436 a.u. DGmin= 0.33 DE= 0.76 kcal/mol
# 3 Count: 2 E= -1029.557705 a.u. DGmin= 0.33 DE= 1.22 kcal/mol
# 4 Count: 1 E= -1029.557276 a.u. DGmin= 0.48 DE= 1.49 kcal/mol
# 5 Count: 2 E= -1029.555755 a.u. DGmin= 0.29 DE= 2.44 kcal/mol
虽然B97-3c优化出的几何结构已经绝对足够理想了，但是其能量计算精度还只能说是不错，要说几乎完美还不够。因此，我们对这5个结构再用更好的级别RI-PWPB95-D3(BJ)/def2-QZVPP去算单点能。把刚产生的cluster.xyz改名为traj.xyz并覆盖掉原先的，将settings.ini里的ngeom复原成默认值0，把itask设为1代表只计算单点能，然后把template.inp内容改成下面这样（注：如果是ORCA 5.0及之后的版本，去掉gridx4 grid4，否则会报错）
! RI-PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX nopop gridx4 grid4 tightscf noautostart
noautostart
%maxcore 4000
%pal nprocs 36 end
* xyz 0 1
[GEOMETRY]
*
启动molclus运行，在笔者的机子上大约8分钟算完一个结构。启动isostat，每个提示都直接敲回车，这5个结构按照高精度的PWPB95-D3(BJ)/def2-QZVPP能量排序后的结果就输出出来了：
# 1 Count: 1 E= -1029.844292 a.u. DGmin= 0.29 DE= 0.00 kcal/mol
# 2 Count: 1 E= -1029.842590 a.u. DGmin= 0.33 DE= 1.07 kcal/mol
# 3 Count: 1 E= -1029.841414 a.u. DGmin= 0.48 DE= 1.81 kcal/mol
# 4 Count: 1 E= -1029.840902 a.u. DGmin= 0.33 DE= 2.13 kcal/mol
# 5 Count: 1 E= -1029.839924 a.u. DGmin= 0.29 DE= 2.74 kcal/mol
如果接着输入298.15，就基于单点能近似地根据Boltzmann公式估计出了常温下的各个构象的分布比率，可见第一个构象是占了绝对的主导地位的。
# 1 Count: 1 Ratio: 80.01%
# 2 Count: 1 Ratio: 13.20%
# 3 Count: 1 Ratio: 3.80%
# 4 Count: 1 Ratio: 2.21%
# 5 Count: 1 Ratio: 0.78%
看看其中能量最低的那个结构是什么样。用VMD载入cluster.xyz，切换到第0帧，结构如下所示
结构很合理，其中存在的内氢键必定是它是最优构象的主要原因之一。我们还可以进一步，按照《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）里介绍的方法，通过Multiwfn+VMD绘制出填色的RDG图，如下所示（笔者后来单独在B3LYP/TZVP级别下对这个结构算了个单点得到波函数文件用于做RDG分析目的）
由箭头所示，不仅此结构里形成了显著的O-H...O氢键（因为对应的RDG等值面颜色非常蓝），氨基的地方还形成了一定程度的氢键（对应的等值面颜色淡蓝），而且分子内还形成了大量色散作用区域，这些因素都对降低这个构象能量产生了重要贡献。
其实更严格判断分子的最稳定构象应该去计算自由能。可以对当前的traj.xyz里每个结构在B97-3c下用freq关键词通过ORCA计算自由能的热校正量，然后手动加到isomers.xyz里的各个能量值上，然后再用isostat去处理此文件，结果就相当于按照自由能排序了，然后最后给出的Boltzmann分布比例也是按照自由能来算的了。
当前例子最耗时的步骤是调用xtb去优化那7000多个构象。想节约时间的话，其实有些可旋转的键可以把每120度转一次改成180度转一次，可以使产生出的构象数目减少N倍。虽然这会增加优化不到能量最低结构的风险，但是考虑到molclus调用xtb产生的isomers.xyz里经过isostat统计后，能量最低结构的簇里面的构象数目多达26个，能量次高的簇里含有33个，因此即便初始产生的构象少一些也不会有什么问题，照样至少会有一个初始结构最后能优化到能量最低构象上去。
值得一提的是molclus的settings.ini里的itask如果设为-2的话，则调用xtb的时候运行的是优化+振动分析任务，而且读取的不是电子能量而是自由能。以自由能来对构象/构型排序比用本例的电子能量排序明显更有意义，但相应地需要额外付出做振动分析的耗时，因此用不用itask= -2由你。
5 基于gentor+molclus的完美构象搜索流程总结
对于较大体系的精确的构象搜索，如果你是比较讲究，有一定水平的人，笔者强烈建议使用以下流程，可谓最完美的做法（上面的例子里的流程都相当于以下过程的简化版）：
• 1 用gentor或分子动力学程序产生初始构象
• 2 对上一步的结构，用molclus结合OpenBabel通过MMFF94力场或结合Gaussian通过UFF分子力场（后者关键词为UFF=qeq opt）快速优化，取几百个能量最低构象。
• 3 对上一步筛出的结构，用molclus结合xtb程序通过GFN-xTB方法或结合MOPAC程序通过PM7半经验方法进一步优化，筛出十几或几十个能量最低构象。这一步用Gaussian做PM7半经验也可以，但明显比MOPAC和xtb慢。
注：如果体系不是特别大，比如还没到50原子，倒不如直接就让molclus调用xtb，不比第1步慢太多，这样可以略过基于分子力场的初筛的过程。
• 4 对上一步剩下的结构，用molclus结合Gaussian或ORCA用普通泛函比如B3LYP-D3(BJ)对这些构象进行优化，基组用6-311G*就够了，取能量最低的数个或十几个结构。之后对这些结构做振动分析获得自由能热校正量（也可以在这一步用molclus的时候直接给模板文件里加上做振动分析的关键词，并且settings.ini里设ibkout= 1让molclus把每一步的输出文件都备份出来，就省得之后再手动做振动分析了）。
注：如果想既便宜结果又比较准确，建议用ORCA里的B97-3c方法做这一步。
• 5 对上一步选出的结构，用双杂化泛函（或者用ORCA里的更精准但更昂贵的DLPNO-CCSD(T)方法）计算它们的单点能，基组起码用may-cc-pVTZ或ma-TZVPP（是分别给cc-pVTZ和def2-TZVPP加上最低限度的弥散函数的基组，详见http://sobereva.com/119和http://sobereva.com/340），若能用到cc-pVQZ或def2-QZVPP基组则更佳（此时一般得利用ORCA里的RI来加速，靠Gaussian很难算得动）。如果这样的级别算着吃力，可改用便宜得多但精度也低一些的M06-2X-D3/ma-TZVP。
• 6 最后，将第5步得到的高精度单点能与第4步产生的自由能热校正量相加，得到这些构象的精确自由能。显然，此时自由能最低的就是热力学角度上最稳定的构象了，同时也可以按照http://sobereva.com/165介绍的方法得到这些构象的玻尔兹曼分布比例。
总的来说，如果你只会用Gaussian，完全可以顺利地结合molclus构象搜索。但如果你还会用xtb、ORCA程序，那么你可以用更短的时间得到更准确的结果。
PS：如果对ORCA、xtb、GFN-xTB、B97-3c这些东西不熟悉，可以看这些了解常识：《盘点Grimme迄今对理论化学的贡献》（http://sobereva.com/388）、《大体系弱相互作用计算的解决之道》（http://sobereva.com/214）、《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）。
有些初学者可能觉得，用Spartan之类的纯图形界面的程序做构象搜索仿佛更容易，连傻子都会用。但实际上，你只要稍微花一丁点时间了解一下molclus的用法，做构象搜索的耗时就会远比使用那些带图形界面的无脑的程序低得多，而精度则比用无脑的程序高得多。molclus的设计已经几乎达到了灵活和易用的极限。有不少人做构象搜索的目的是为了计算ECD（参考比如《使用Multiwfn绘制构象权重平均的光谱》http://sobereva.com/383），这对构象分布比例的准确度要求是非常高的，靠Spartan等封闭的、支持的计算级别很low的构象搜索程序基本不可能满足要求。
本文的例子是真空下的，对于溶剂下的构象搜索，量子化学计算时千万别忘了考虑隐式溶剂模型，因为这会严重影响构象分布，比如笔者的这篇文章Struct. Chem., 25, 1521 (2014)的表2就清楚地体现了这点。关于隐式溶剂模型的考虑见《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）。
顺带一提，使用gentor还可以很容易地结合Gaussian做二面角的刚性扫描获得势能曲线，对实际研究非常有用，详见《详谈使用Gaussian做势能面扫描》（http://sobereva.com/474）。
| 
[复制链接 Copy URL]
发表于 Post on 2015-12-18 04:50:50
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
发表于 Post on 2015-12-18 07:57:24
楼主 Author