计算化学公社

标题: gentor：扫描方式做分子构象搜索的便捷工具 [打印本页]

作者
Author: sobereva 时间: 2015-12-18 04:50
标题: gentor：扫描方式做分子构象搜索的便捷工具
注：每次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播放起来就会如下所示，可见上方和左方的甲基是顺时针同步旋转，而右侧的甲基是逆时针同步旋转：

(, 下载次数 Times of downloads: 1094)

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做构象搜索的基本流程。分子初始结构如下：

(, 下载次数 Times of downloads: 284)

一般情况下，分子结构用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/），可看到以下初始结构，显然都是合理的，可以做下一步了。

(, 下载次数 Times of downloads: 260)

我们用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下能量最低的结构为

(, 下载次数 Times of downloads: 244)

仔细看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里显示，如下所示：

(, 下载次数 Times of downloads: 195)

此体系可旋转的键比较多。此体系里有个类似肽键的区域，由于存在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帧，结构如下所示

(, 下载次数 Times of downloads: 189)

结构很合理，其中存在的内氢键必定是它是最优构象的主要原因之一。我们还可以进一步，按照《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）里介绍的方法，通过Multiwfn+VMD绘制出填色的RDG图，如下所示（笔者后来单独在B3LYP/TZVP级别下对这个结构算了个单点得到波函数文件用于做RDG分析目的）

(, 下载次数 Times of downloads: 194)

由箭头所示，不仅此结构里形成了显著的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）。

作者
Author: ruanyang 时间: 2015-12-18 07:57
好物

作者
Author: greatzdk 时间: 2015-12-18 09:13
设置旋转参数的文件（gentor.ini），每个分子应该不大一样，如果打算批量搜索一群分子，设置这个参数文件，应该会变成一件头疼的事情。

作者
Author: sobereva 时间: 2015-12-18 12:16

greatzdk 发表于 2015-12-18 09:13
设置旋转参数的文件（gentor.ini），每个分子应该不大一样，如果打算批量搜索一群分子，设置这个参数文件， ...

一般不会涉及那种情况的。要考察构象的分子数目不会很多，人工弄很容易。肉眼一看有几个键需要考虑旋转，立刻就写上去了。任何构象搜索程序都避免不了这个过程。如果程序自动指定，若设定不合理可能造成大量时间浪费。

作者
Author: smutao 时间: 2016-1-2 00:00
有没有可能自动识别可以旋转的键？

作者
Author: liyuanhe211 时间: 2016-1-2 04:06
本帖最后由 liyuanhe211 于 2016-1-2 04:08 编辑

smutao 发表于 2016-1-2 00:00
有没有可能自动识别可以旋转的键？

Spartan可以做这件事情（Set Torsions 之后会自动给出一系列建议猜测），但对较复杂的分子经常不靠谱，包括转不转、要转几个fold等等，需要人为调。

(, 下载次数 Times of downloads: 203)

不过有GUI确实比写文本方便一点。

作者
Author: lidanhoney 时间: 2016-5-29 23:09
老师，为什么没有生成cluster.xyz的文件，而且运行isostat，按三次回车后，软件闪了一下就关闭了，我运行的就是帖子中的例子

作者
Author: sobereva 时间: 2016-5-30 00:49

lidanhoney 发表于 2016-5-29 23:09
老师，为什么没有生成cluster.xyz的文件，而且运行isostat，按三次回车后，软件闪了一下就关闭了，我运行的 ...

绝对是你的文件有问题，把你的isomer.xyz贴出来

作者
Author: lidanhoney 时间: 2016-5-30 19:49

sobereva 发表于 2016-5-30 00:49
绝对是你的文件有问题，把你的isomer.xyz贴出来

是我的mopac的软件安装有问题，现在安装好了就可以用了，不过用gaussian优化是不是很慢，好久都没有什么反应

作者
Author: sobereva 时间: 2016-5-31 00:38

lidanhoney 发表于 2016-5-30 19:49
是我的mopac的软件安装有问题，现在安装好了就可以用了，不过用gaussian优化是不是很慢，好久都没有什么 ...

自行监控gau.out看看算到什么情况了

作者
Author: lidanhoney 时间: 2016-6-1 18:49
对应的mol.xyz文件内容如下所示，将此文件放到molclus目录下
22
dopamine
O -2.33158827 1.90324395 0.05204669
老师，这个xyz文件中的22是什么意思，我发现在我提取坐标后，设置不一样的值，用VMD看到的结果完全不一样。那这个22是什么意思

作者
Author: sobereva 时间: 2016-6-1 21:20

lidanhoney 发表于 2016-6-1 18:49
对应的mol.xyz文件内容如下所示，将此文件放到molclus目录下
22
dopamine

原子数

作者
Author: therotyonth 时间: 2016-6-21 20:13
非常感谢

作者
Author: ZCSco 时间: 2016-8-10 12:04
老师好，如果我想搜索过渡态的conformer要怎样处理呢？是把过渡态的部分固定住还是应该怎样比较好呀？谢谢老师！

作者
Author: sobereva 时间: 2016-8-10 12:34

ZCSco 发表于 2016-8-10 12:04
老师好，如果我想搜索过渡态的conformer要怎样处理呢？是把过渡态的部分固定住还是应该怎样比较好呀？谢谢 ...

先找一个构型下的过渡态，然后对这个过渡态结构用gentor产生各个二面角的各种可能情况的结构。优化的时候把过渡态的原子冻住。

作者
Author: ZCSco 时间: 2016-8-10 12:50

sobereva 发表于 2016-8-10 12:34
先找一个构型下的过渡态，然后对这个过渡态结构用gentor产生各个二面角的各种可能情况的结构。优化的时 ...

嗯嗯，明白了，谢谢老师！

作者
Author: ZCSco 时间: 2016-8-13 15:07
老师好，再请教一下，当生成过渡态的conformer结构时，因为过渡态原子距离比较远，程序判断为没有成键，所以并没有一起旋转。能否在定义转键的时候把要转的片段也定义进去呢？

作者
Author: sobereva 时间: 2016-8-14 04:53

ZCSco 发表于 2016-8-13 15:07
老师好，再请教一下，当生成过渡态的conformer结构时，因为过渡态原子距离比较远，程序判断为没有成键，所 ...

我更新了一下gentor程序：
(, 下载次数 Times of downloads: 27)

这个版本可以自定义连接关系、自己设判断成键的阈值，比如gentor.ini的内容：

7-13
e120
5-7
e120
1-5
critbond=1.15
3,6,1
3,8,1
5,7,0

复制代码

就代表两个原子间距离小于它们共价半径和的1.15倍就当成成键。然后再读取额外的设定，把3,6、3,8认为是成键，而把5,7取消成键。

你用用，如果能解决你的问题告诉我一声。

作者
Author: ZCSco 时间: 2016-8-15 11:13

sobereva 发表于 2016-8-14 04:53
我更新了一下gentor程序：

试用了，很好用！谢谢老师

您是用判断成键的方法，比之前想的每一个都自己定义片段的方法还要方便许多，非常感谢！

作者
Author: ZCSco 时间: 2016-11-4 09:44
老师好，想写一个类似gentor功能的code，以实现内转动势能面扫描、片段判断和标记、解离曲线扫描等，帮助我处理大批量的分子。想请教一下老师，程序的思路是怎样比较好。我有两种想法，一个是根据键长判断原子的连接关系，把所有有直接或间接联系的原子放进一个新的submolecule，然后再对它进行平移或者旋转；另一种是能不能通过openbabel把直角坐标转化成SMILE或别的分子标记格式，然后再判断是否为同一物种或片段，至于转键或拉伸键长的扫描，是不是用转换成内坐标的格式，然后再操作呢？其实用gentor就非常好了，但是我想把这个和其它后面的处理结合，所以希望能自己学着写一下。我程序方面很弱，不太清楚要怎样用代码实现对分子的相关操作，求指点，谢谢老师！

作者
Author: sobereva 时间: 2016-11-4 14:46

ZCSco 发表于 2016-11-4 09:44
老师好，想写一个类似gentor功能的code，以实现内转动势能面扫描、片段判断和标记、解离曲线扫描等，帮助我 ...

可以参考truhlar组的MStor程序包里的confgen程序的源码

作者
Author: ZCSco 时间: 2016-11-4 19:58
对哦！我咋没想到呢谢谢！！

作者
Author: 杨小狗 时间: 2017-2-23 15:03
文中gentor.inp好像应该是gentor.ini
我正在学习，先按照例子操作操作
另外我有一个问题，购买molclus程序后可否放在SecureCRT中然后再调用集群主机中的动力学程序

作者
Author: sobereva 时间: 2017-2-24 01:27

杨小狗发表于 2017-2-23 15:03
文中gentor.inp好像应该是gentor.ini
我正在学习，先按照例子操作操作
另外我有一个问题，购买molclus程 ...

molclus并不会调用动力学程序，而是把你单独用动力学程序跑好的xyz轨迹让molclus处理

作者
Author: 杨小狗 时间: 2017-2-24 08:38

sobereva 发表于 2017-2-24 01:27
molclus并不会调用动力学程序，而是把你单独用动力学程序跑好的xyz轨迹让molclus处理

好的，谢谢老师我懂啦。

作者
Author: janihyhubo 时间: 2017-3-7 11:59
Sob社长，您好，我按照您帖子里的例子正常运行了，但是一到我自己的分子就成功不了，下面，分别是对应的输入输出文件，我只是想要旋转一个C-H键，但是得到的输出文件，坐标完全一样，麻烦您帮忙看一下，谢谢

作者
Author: sobereva 时间: 2017-3-7 14:09

janihyhubo 发表于 2017-3-7 11:59
Sob社长，您好，我按照您帖子里的例子正常运行了，但是一到我自己的分子就成功不了，下面，分别是对应的输 ...

H又没连着其它基团，你转动C-H，显然结构看不出变化

作者
Author: mooninwhere 时间: 2017-4-25 20:17
感谢大神的介绍！

作者
Author: 维维最爱小黄人 时间: 2017-5-28 09:58
老师，我在最后调用gaussian进行优化的时候，molclus只一下就闪退，生成的gau.gif里只有一个结构的坐标，是gaussian的原因么？

作者
Author: sobereva 时间: 2017-5-28 12:38

维维最爱小黄人发表于 2017-5-28 09:58
老师，我在最后调用gaussian进行优化的时候，molclus只一下就闪退，生成的gau.gif里只有一个结构的坐标，是 ...

估计你没正确写模板文件，导致也因此没法正确产生高斯输入文件

作者
Author: 小苹果 时间: 2017-7-9 08:34
运行gentor for linux 版，出现如下错误，而运行win版本正常得到结果，请问为什么？谢谢！

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

At line 70 of file gentor.f90 (unit = 10, file = 'gentor.ini')

Fortran runtime error: End of file

作者
Author: sobereva 时间: 2017-7-9 09:06

小苹果发表于 2017-7-9 08:34
运行gentor for linux 版，出现如下错误，而运行win版本正常得到结果，请问为什么？谢谢！

Criterion f ...

gentor.ini末尾加几个空行再试

作者
Author: 小苹果 时间: 2017-7-9 09:31
一个空行，解决问题，谢谢！

作者
Author: 茶新味 时间: 2017-8-15 08:35
Sob老师，读不懂您说的什么意思，文中“我们用molclus调用MOPAC2012在PM7级别下对各个结构进行优化。把traj.xyz放到molclus目录下，确认MOPAC的模板输入文件里的关键词是PM7 precise，确保molclus的settings.ini的文件里的iprog为2。然后启动molclus，程序就会调用MOPAC进行优化“那我是放在molclus-MOPAC2012下，还是molclus-1.3.5-Win-free下面呢？软件中的genmer又有什么作用呢？有没有直接调用Gauss的方法。

作者
Author: sobereva 时间: 2017-8-15 12:07

茶新味发表于 2017-8-15 08:35
Sob老师，读不懂您说的什么意思，文中“我们用molclus调用MOPAC2012在PM7级别下对各个结构进行优化。把traj ...

molclus帖子一开头就给出了genmer的帖子
genmer：生成团簇初始构型的超便捷工具
http://bbs.keinsci.com/forum.php?mod=viewthread&tid=2369

怎么调用Gaussian在贴子里说得很详细了

文中已经说了，traj.xyz必须放在molclus目录下，从来没说要放到MOPAC目录下

作者
Author: 茶新味 时间: 2017-9-4 10:54
Sob老师，您好。我想请教一下，您帖子里说先molclus调用MOPAC2012在PM7级别下对各个结构进行优化，删掉之前的traj.xyz，把cluster.xyz改名为traj.xyz使得刚才优化好的结构作为初始结构，再调用Gaussian，那可不可以直接把gentor输出的traj.xyz直接用molclus调用Gaussian来优化呢？

作者
Author: sobereva 时间: 2017-9-4 11:19

茶新味发表于 2017-9-4 10:54
Sob老师，您好。我想请教一下，您帖子里说先molclus调用MOPAC2012在PM7级别下对各个结构进行优化，删掉之前 ...

可

作者
Author: 茶新味 时间: 2017-9-4 14:49
Sob老师，用Gaussian程序的人比较多，是不是可以分开写分别用MOPAC2012和Gaussian的，而不是混合写在一起，这样可能看起来会不会更简明一些呢。我按照帖子算例，把traj.xyz放到molclus目录下，把settings里面改成1，然后启用molclus程序，闪一下就关闭了，没有显示调动起来Gaussian程序，也没有结果输出在isostat里，请问怎么回事呢？是哪里有问题吗

作者
Author: 茶新味 时间: 2017-9-4 21:00
Sob老师，按照帖子算例，直接调用高斯程序，不用MOPAC2012其中setting中的ieneonly设为1还是0呢，直接调用高斯总会出现图一的错误，isostat里面没有输出文件，图二到图四分别是我在setting，template1和template2设置的关键字，请问Sob老师，是什么原因算不过去呢？

作者
Author: sobereva 时间: 2017-9-4 22:56

茶新味发表于 2017-9-4 14:49
Sob老师，用Gaussian程序的人比较多，是不是可以分开写分别用MOPAC2012和Gaussian的，而不是混合写在一起， ...

完全没有必要分开，不会造成任何使用上的麻烦，只会增加便利

作者
Author: sobereva 时间: 2017-9-4 22:57

茶新味发表于 2017-9-4 21:00
Sob老师，按照帖子算例，直接调用高斯程序，不用MOPAC2012其中setting中的ieneonly设为1还是0呢，直接调用 ...

molclus帖子里写得很明确：
gaussian_path：Gaussian的可执行文件的路径，比如"/sob/g09/g09"、"d:\study\g09w\g09.exe"。也可以只写"g09"，前提是程序所在目录已经加入到了PATH环境变量里。如果路径中有空格，则必须有双引号，没空格的话无所谓，后同

不要指向g09w.exe

作者
Author: 茶新味 时间: 2017-9-5 08:06
Sob老师，您好。

路径问题按照您说的改好了改好了，可是算的过程中出现了图一的错误，算出的gau.gif都打不开，isostat里面没有输出文件。Sob老师的帖子是在MOPAC2012计算的基础上得到可靠的数据再用高斯，我直接用Gaussian顺序是不是应该先优化ieneonly设为0，在计单点ieneonly设为1呢？

作者
Author: sobereva 时间: 2017-9-5 10:59

茶新味发表于 2017-9-5 08:06
Sob老师，您好。路径问题按照您说的改好了改好了，可是算的过程中出现了图一的错误，算出的gau.gif都打 ...

单点还是优化随意，没有顺序，molclus的执行流程是超级灵活的
当前没设GAUSS_EXEDIR所致

作者
Author: 茶新味 时间: 2017-9-5 14:03
谢谢Sob老师，非常受用。

作者
Author: ggdh 时间: 2018-6-25 13:50
本帖最后由 ggdh 于 2018-6-26 09:46 编辑

今天使用gentor,目的是扫描一个二面角，产生大量xyz文件后，转成gaussian09输入文件，然后计算某些性质，在sob帮助下解决了，分享在此：
1.使用脚本产生大量角度
假设想转60度，每0.3度转一次，那么一共会产生60/0.3=200个结构，这样手动输入比较麻烦
在linux下可以输入
for i in `seq 200`; do echo -n `echo "$i*0.3" | bc`" " ; done >> gentor.ini
在window下可以
用excel生成个数列，然后把换行符替换成空格填进gentor.ini

2.把traj.xyz分割成多个文件，需要安装obabel
obabel -i xyz traj.xyz -o xyz -m -O confor.xyz
这时会产生confor1.xyz confor2.xyz ....... confor200.xyz 每个xyz中一个结构。
然后用脚本把这些xyz转换为gjf文件，或者orca的输入文件，就可以为所欲为了

作者
Author: sobereva 时间: 2018-6-25 19:32

ggdh 发表于 2018-6-25 13:50
今天使用gentor,目的是扫描一个二面角，产生大量xyz文件后，转成gaussian09输入文件，然后计算某些性质，在 ...

刚更新了molclus 1.4.3版，现在里面的gentor可以用比如e5 120 140来让片段依次旋转120,125,130,135,140度，也可以用e-5 140 120来反过来扫。

作者
Author: Daniel_Arndt 时间: 2018-11-3 05:51
我想问一下，gentor或者MSTor中的ConfGen能否应用于一些大环化合物的构象搜索？

作者
Author: sobereva 时间: 2018-11-3 07:49

Daniel_Arndt 发表于 2018-11-3 05:51
我想问一下，gentor或者MSTor中的ConfGen能否应用于一些大环化合物的构象搜索？

不能。因为旋转环当中的一个键会导致环的破坏

作者
Author: Daniel_Arndt 时间: 2018-11-6 03:35

sobereva 发表于 2018-11-3 07:49
不能。因为旋转环当中的一个键会导致环的破坏

多谢。看来我只能尝试分子动力学的方法了。

作者
Author: TMULYU 时间: 2019-1-24 10:56
我是一个研一小菜鸟，请问您，我目前在处理一个柔性比较大的链状分子，先做构象搜索，然后优化，打算找全局的最优构象，但是需要考虑的键10-12个，整个分子有86个原子，其他软件也试过，系统搜索什么的，但是构象太多，最后以失败告终，请问可以用您的这个软件吗，您可以分享给我一些思路吗

作者
Author: 黎明前夜 时间: 2019-1-27 16:47
你好老师，我想问一下出现图一上面的找不到*.out是什么意思，我哪里弄错了么？始终无法攻克这一步，望指点

作者
Author: sobereva 时间: 2019-1-28 01:12

黎明前夜发表于 2019-1-27 16:47
你好老师，我想问一下出现图一上面的找不到*.out是什么意思，我哪里弄错了么？始终无法攻克这一步， ...

不是报错，不用管。本来当前就没有那些文件

作者
Author: sobereva 时间: 2019-1-28 01:13

TMULYU 发表于 2019-1-24 10:56
我是一个研一小菜鸟，请问您，我目前在处理一个柔性比较大的链状分子，先做构象搜索，然后优化，打算找全局 ...

用gromacs做周期性模拟退火，把每次退火循环中对应0K的帧写入traj.xyz（用VMD可以导出此文件），比如这样搞出100帧或更多帧，再用molclus处理

作者
Author: 黎明前夜 时间: 2019-1-28 22:43
本帖最后由黎明前夜于 2019-1-28 22:46 编辑

SOB老师，调用高斯的时候出现图1这样的情况，该如何解决呀，试了好多次没有攻克。先行感谢！

作者
Author: sobereva 时间: 2019-1-28 23:06

黎明前夜发表于 2019-1-28 22:43
SOB老师，调用高斯的时候出现图1这样的情况，该如何解决呀，试了好多次没有攻克。先行感谢！

没有叫TZVP*的基组

作者
Author: 黎明前夜 时间: 2019-1-28 23:38
可改成如图2那样子也会出现像图1一样的报错。。。

作者
Author: sobereva 时间: 2019-1-29 00:08

黎明前夜发表于 2019-1-28 23:38
可改成如图2那样子也会出现像图1一样的报错。。。

如果你是Windows 32bit版Gaussian，显然%mem和%nproc不能设这么大，这里提了
Gaussian的安装方法及运行时的相关问题
http://sobereva.com/439（http://bbs.keinsci.com/thread-10814-1-1.html）

作者
Author: sobereva 时间: 2019-2-13 12:34
molclus 1.6已发布，见http://bbs.keinsci.com/thread-12134-1-1.html，其中对gentor做了更新，加入了扭转项间的同步旋转设定。本文也做了相应更新、添加了新的例子

作者
Author: 胡说 时间: 2019-3-5 20:52
老师你好请教一个我在使用gentor和molclus时碰到的问题：
我用gentor旋转某些键产生了上千个结构，xtb优化后，使用isostat默认的能量和结构阈值可以得到三百多个结构，发现能量间相差很小，分布很密集，最大的能量差也小于2kcal/mol。用VMD观察，发现很多结构肉眼无法分辨差别。我不想只取前几或几十个结构DFT优化，但此步产生结构数量实在是多。
这种情况我想设大 isostat 阈值，能量设为0.5 kcal/mol 结构设为0.25 angstorm。这样最后产生的结构就只有几十个结构。用VMD观察产生的结构，也无觉得不妥。
请问这种情况下的这种做法是否安全可取？
谢谢。

作者
Author: sobereva 时间: 2019-3-6 07:14

胡说发表于 2019-3-5 20:52
老师你好请教一个我在使用gentor和molclus时碰到的问题：
我用gentor旋转某些键产生了上千个结构，xtb ...

可取

作者
Author: 胡说 时间: 2019-3-6 09:40

sobereva 发表于 2019-3-6 07:14
可取

多谢老师

作者
Author: cep 时间: 2019-3-18 16:20
老师，像包结物AB，假设A分子是固定的，要想知道B分子位于A分子哪个位置比较好可以采用构象搜索吗？包结物AB之间没有键连接

作者
Author: sobereva 时间: 2019-3-18 17:03

cep 发表于 2019-3-18 16:20
老师，像包结物AB，假设A分子是固定的，要想知道B分子位于A分子哪个位置比较好可以采用构象搜索吗？包结物A ...

这属于团簇构型搜索问题，应当用genmer，看
genmer：生成团簇初始构型的超便捷工具
http://bbs.keinsci.com/thread-2369-1-1.html

作者
Author: cep 时间: 2019-3-18 21:13

sobereva 发表于 2019-3-18 17:03
这属于团簇构型搜索问题，应当用genmer，看
genmer：生成团簇初始构型的超便捷工具
http://bbs.keinsci ...

好的，谢谢老师

作者
Author: sobereva 时间: 2019-3-30 07:19
对本文进行了重要更新，增加了第4节，是xtb结合ORCA搜索高度柔性的尺寸不算小的分子构象的例子

作者
Author: happyknighthawk 时间: 2019-3-30 16:28
Sob老师，有两个小问题想请教您：
1、对于存在柔性环系（如六元环、七元环等）的化合物，在构象搜索中需要涉及到类似船式、椅式这样的翻转，这个时候在gentor.ini中应如何设置，才能使环上的原子flip从而产生翻转的效果？
2、如果想比较客观的去比较和评价几款构象搜索软件在构象筛选过程中的表现（主要针对分子量大概在300—1000之间的、主要含CHON等原子的化合物），对于评价方法您有什么建议？
非常感谢您

作者
Author: sobereva 时间: 2019-3-31 05:00

happyknighthawk 发表于 2019-3-30 16:28
Sob老师，有两个小问题想请教您：
1、对于存在柔性环系（如六元环、七元环等）的化合物，在构象搜索中需要 ...

1 这种情况应当做动力学，不适合gentor。如此文里的MD+molclus的做法
使用molclus程序做团簇构型搜索和分子构象搜索
http://bbs.keinsci.com/thread-577-1-1.html

2 分子量是次要的，关键是可旋转的键的数目
评价标准：
1 费用
2 能支持的理论方法种类、精度
3 确保找到能量最低构象花费的总耗时（这个最关键，但很不好严格比）
4 使用方便程度、支持的操作系统种类
5 使用流程的灵活性
6 对体系种类的普适性、体系里能存在的元素种类
等等

作者
Author: happyknighthawk 时间: 2019-3-31 10:06

sobereva 发表于 2019-3-31 05:00
1 这种情况应当做动力学，不适合gentor。如此文里的MD+molclus的做法
使用molclus程序做团簇构型搜索和 ...

非常感谢Sob老师这么细致的回复，谢谢。

作者
Author: captain 时间: 2019-4-1 11:48
请问大神
在这个ORCA计算中，
! RI-PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX nopop gridx4 grid4 tightscf noautostart
辅助基组 def2-QZVPP/C 是用在相关部分计算中的吗？
不加 def2-QZVPP/C，也可以吧？

作者
Author: sobereva 时间: 2019-4-1 11:53

captain 发表于 2019-4-1 11:48
请问大神
在这个ORCA计算中，
! RI-PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX nopop gridx4 gri ...

是
得加，要不然MP2部分RI没法用

作者
Author: captain 时间: 2019-4-1 12:09

sobereva 发表于 2019-4-1 11:53
是
得加，要不然MP2部分RI没法用

明白了谢大神

作者
Author: asl1994 时间: 2019-5-11 16:54
老师，我参考实例1做了构想搜索，有个审稿人问在计算PM7结构时设置的energy cut-off是多少，但是好像没有特别设置，请问要如何回答呢？

作者
Author: sobereva 时间: 2019-5-11 17:14

asl1994 发表于 2019-5-11 16:54
老师，我参考实例1做了构想搜索，有个审稿人问在计算PM7结构时设置的energy cut-off是多少，但是好像没有特 ...

我不知道你说的cutoff具体指什么
如果是指把相对于最低构象多少能量以下的取出来用于下一步精确计算，这个标准是你自己定的

作者
Author: sobereva 时间: 2019-6-12 20:31
今日更新了molclus 1.8.4，其中gentor增加了critskip参数，可以设定初始结构中两个原子距离小于二者共价半径和的多少倍就被刨去。原先程序默认值是0.5，现在用户可以自己改了。

作者
Author: pcxlm1 时间: 2019-6-26 20:58
sob老师，请问我运行gentor以后输出的traj.xyz文件是空的，linux和windows版本都一样，这可能是什么原因？

作者
Author: sobereva 时间: 2019-7-1 10:56

pcxlm1 发表于 2019-6-26 20:58
sob老师，请问我运行gentor以后输出的traj.xyz文件是空的，linux和windows版本都一样，这可能是什么原因？

仔细看gentor运行过程的提示，每个词都仔细看。如果正常运行完了，而且至少有一个是成功生成的，不可能是空的。

如果说不明白，把你的gentor.ini和相关的分子xyz文件都打包传上来

作者
Author: pcxlm1 时间: 2019-7-5 10:49

sobereva 发表于 2019-7-1 10:56
仔细看gentor运行过程的提示，每个词都仔细看。如果正常运行完了，而且至少有一个是成功生成的，不可能是 ...

压缩包内是我的mol.xyz和gentor.ini，结果显示应该289个构型生成，但是traj.xyz是空文件，请sob老师多指教。

作者
Author: sobereva 时间: 2019-7-5 16:07

pcxlm1 发表于 2019-7-5 10:49
压缩包内是我的mol.xyz和gentor.ini，结果显示应该289个构型生成，但是traj.xyz是空文件，请sob老师多指 ...

诸如e5 180 90应当写为e5 90 180，你的扫描是每步增大，不能初值比终值还大

作者
Author: pcxlm1 时间: 2019-7-6 16:56

sobereva 发表于 2019-7-5 16:07
诸如e5 180 90应当写为e5 90 180，你的扫描是每步增大，不能初值比终值还大

谢谢老师，已经解决了

作者
Author: 欢乐多 时间: 2019-8-4 09:49
本帖最后由欢乐多于 2019-8-4 12:57 编辑

sob老师，将PM6-DH改为PM7后，将不能输出isomers文件，或者输出的isomers文件是空白的？这是MOPAC没有装好吗？怎么办？什么时候能把操作步骤做成视频或图像啊，这语言叙述的对于初学者很抽象。有些时候有一点不一样的，或者我们无意识做的不一样的，您的结果就出不来，而有时候都要折腾很长时间，心力交瘁，万念俱灰啊

作者
Author: sobereva 时间: 2019-8-5 01:51

欢乐多发表于 2019-8-4 09:49
sob老师，将PM6-DH改为PM7后，将不能输出isomers文件，或者输出的isomers文件是空白的？这是MOPAC没有装好 ...

仔细查看输出文件，8成是关键词没写对

文中描述得非常清楚，在描述上已经极力避免歧义。注意理解文章每一步的命令的含义，这样有点自己理解不清楚的稍微一动脑就能弄明白

作者
Author: 欢乐多 时间: 2019-8-5 14:09

sobereva 发表于 2019-8-5 01:51
仔细查看输出文件，8成是关键词没写对

文中描述得非常清楚，在描述上已经极力避免歧义。注意理解文章 ...

老师，安照博文，这一步过不去了。怎么回事，是高斯输入文件有问题吗？请老师指教

作者
Author: 欢乐多 时间: 2019-8-5 15:04

sobereva 发表于 2019-8-5 01:51
仔细查看输出文件，8成是关键词没写对

文中描述得非常清楚，在描述上已经极力避免歧义。注意理解文章 ...

这一步MOPAC优化中将PM7 precise 改为M06-2X/def2-SVP,就不能优化了，烦请老师指教。

作者
Author: sobereva 时间: 2019-8-6 07:46

欢乐多发表于 2019-8-5 14:09
老师，安照博文，这一步过不去了。怎么回事，是高斯输入文件有问题吗？请老师指教

文中明确写了关键词是M062X，你写成了M06-2X，不要随意自行发挥

作者
Author: sobereva 时间: 2019-8-6 07:47

欢乐多发表于 2019-8-5 15:04
这一步MOPAC优化中将PM7 precise 改为M06-2X/def2-SVP,就不能优化了，烦请老师指教。

文中写明了，这步调用的是MOPAC，关键词也说了是PM7 precise，你怎么能把Gaussian里的关键词写到MOPAC模板文件里？

作者
Author: 欢乐多 时间: 2019-8-6 17:08

sobereva 发表于 2019-8-6 07:46
文中明确写了关键词是M062X，你写成了M06-2X，不要随意自行发挥

老师，经过多次尝试，深入研究发现setting.int中itask改为0优化，计算才得以顺利进行。开心。

作者
Author: 欢乐多 时间: 2019-8-6 17:26
本帖最后由欢乐多于 2019-8-6 17:34 编辑

sobereva 发表于 2019-8-6 07:47
文中写明了，这步调用的是MOPAC，关键词也说了是PM7 precise，你怎么能把Gaussian里的关键词写到MOPAC模 ...

好的，老师，多次琢磨尝试，setting.int中energyterm= MP2=改为energyterm= HF=,计算得以实现，此时，开心地发现余已初步掌握molclus的基本使用方法，内心开始清明，对诸多量化软件的使用逐步领悟，多日折腾与摸索，不懈追求与尝试开始出现了一点进步，但是对量化的认识还不是很清楚，还不知道符号背后的故事，希望参加即将到来的量化班培训，能够使自己对量化认识的更加清楚，同时，更希望能够学以致用，让博大精深的量化理论配合切实精妙的实践来提高我博士研究课题的质量，攻克实际研究过程中的壁垒，完美化解实验过程遇到的问题。

作者
Author: 欢乐多 时间: 2019-10-1 15:53
本帖最后由欢乐多于 2019-10-1 16:14 编辑

sobereva 发表于 2019-3-31 05:00
1 这种情况应当做动力学，不适合gentor。如此文里的MD+molclus的做法
使用molclus程序做团簇构型搜索和 ...

老师，对于附图中这样的分子，进行构象搜索时，用gentor可以吗?还是用您说的MD程序？
我认为这样：这种类型的分子虽然存在船椅式的不同，不管用gentor还是MD，搜索到的构象进行优化，最终优化的结果区别应该不是很大。
当然，如果这样，gentor进行构象搜索太方便了。

作者
Author: sobereva 时间: 2019-10-3 03:26

欢乐多发表于 2019-10-1 15:53
老师，对于附图中这样的分子，进行构象搜索时，用gentor可以吗?还是用您说的MD程序？
我认为这样：这种 ...

若想考虑船式、椅式的差别，只能用MD产生初始的traj.xyz。xtb做基于GFN-xTB理论的MD还是挺快挺方便的，都不用管力场的事。实际上我之前就打算写个基于xtb做MD结合molclus做环状体系构象搜索的帖子，暂时没时间，这个月找时间写出来

作者
Author: 欢乐多 时间: 2019-10-4 18:34

sobereva 发表于 2019-10-3 03:26
若想考虑船式、椅式的差别，只能用MD产生初始的traj.xyz。xtb做基于GFN-xTB理论的MD还是挺快挺方便的，都 ...

好的，老师，时刻关注xtb做分子模拟搜索构象的教程。

作者
Author: 欢乐多 时间: 2019-10-25 23:45
本帖最后由欢乐多于 2019-10-26 00:22 编辑

sobereva 发表于 2019-10-3 03:26
若想考虑船式、椅式的差别，只能用MD产生初始的traj.xyz。xtb做基于GFN-xTB理论的MD还是挺快挺方便的，都 ...

老师，gentor构象搜索后，用xtb进行优化，优化结果能量排序，较低能量结构显示有些碳原子连在一起了，多出一个碳键（图1），还有一些结构显示很怪异（图2），这个怎办？
注：用的是molclus1.8.9；附件中有分子的初始结构，及构象搜索后用xtb优化的结果

作者
Author: sobereva 时间: 2019-10-26 05:00

欢乐多发表于 2019-10-25 23:45
老师，gentor构象搜索后，用xtb进行优化，优化结果能量排序，较低能量结构显示有些碳原子连在一起了，多 ...

图2的情况通常是整个isomers.xyz里有些结构和第一帧的成键关系不同所致。对于一个稳定的有机分子，应当批量优化出来的结构里成键关系都是相同的。
你的问题，有可能是对键的旋转设置有问题，导致初始结构里就有一些不甚合理，应当仔细检查traj.xyz；也有可能是当前计算级别对这个体系描述得不太理想。

作者
Author: IIIO_OIII 时间: 2020-2-18 01:57
本帖最后由 IIIO_OIII 于 2020-2-18 02:46 编辑

社长，我用gentor产生构象的时候，氨基酸侧链末端的NH3不跟着相连的C一起走。用了windows版本molclus1.9.3。

找到原因了。初始构象里C-NH3键太长了, 所以没有被识别。打扰了。

作者
Author: Lx315 时间: 2020-5-25 17:07
老师，我在做分子构象搜索实例1：多巴胺（利用MOPAC+Gaussian）的时候出现一个问题，使用VMD查看traj.xyz初始结构是没有问题的，有九个构象，之后用MOPAC优化、运行isostat存在重复构象，排序出来的cluster.xyz就只有四个构象了，我不知道是哪出了问题，是用MOPAC优化出了问题吗？

作者
Author: snljty 时间: 2020-5-25 18:26

Lx315 发表于 2020-5-25 17:07
老师，我在做分子构象搜索实例1：多巴胺（利用MOPAC+Gaussian）的时候出现一个问题，使用VMD查看traj.xyz初 ...

这不是再正常不过了么...

作者
Author: sobereva 时间: 2020-5-26 16:14

Lx315 发表于 2020-5-25 17:07
老师，我在做分子构象搜索实例1：多巴胺（利用MOPAC+Gaussian）的时候出现一个问题，使用VMD查看traj.xyz初 ...

去除重复的可不就减少了么

作者
Author: Lx315 时间: 2020-5-26 17:43

sobereva 发表于 2020-5-26 16:14
去除重复的可不就减少了么

可是为什么我的用MOPAC优化后会有重复构象，我问了同学们都没有这个问题，我用gaussview建立了好多次初始结构都还是有这个问题，用同学给我的数据就没有这个问题了。

作者
Author: sobereva 时间: 2020-5-27 05:07

Lx315 发表于 2020-5-26 17:43
可是为什么我的用MOPAC优化后会有重复构象，我问了同学们都没有这个问题，我用gaussview建立了好多次初始 ...

显然取决于初始结构，两个初始结构都在同一个势能面极小点附近自然会收敛到同一个结构取。仔细看molclus的教程
倘若不可能有重复构象的话isostat去重的功能还有什么意义

作者
Author: 赤铜彩皮兽 时间: 2020-8-14 09:51
老师好，我用gentor进行化学键旋转，但是所有的构型都没有通过distance check，我的结构实际上是一个自由基，请问出现这样的结果应该如何解决好呢？
216 conformations were originally generated
216 conformations did not pass distance check and thus be ignored
Final number of conformations: 0

作者
Author: sobereva 时间: 2020-8-15 09:57

赤铜彩皮兽发表于 2020-8-14 09:51
老师好，我用gentor进行化学键旋转，但是所有的构型都没有通过distance check，我的结构实际上是一个自由基 ...

我得有你的具体输入文件才能判断

欢迎光临计算化学公社 (http://bbs.keinsci.com/)