计算化学公社
标题: 将Confab或Frog2与Molclus联用对有机体系做构象搜索 [打印本页]
作者Author: sobereva 时间: 2020-10-28 19:48
标题: 将Confab或Frog2与Molclus联用对有机体系做构象搜索
将Confab或Frog2与Molclus联用对有机体系做构象搜索
Using Confab or Frog2 with Molclus to perform conformational searches for organic systems
文/Sobereva@北京科音 First release: 2020-Oct-28 Last update: 2022-Oct-8
1 前言
molclus是免费、灵活、好用的团簇构型与分子构象搜索程序，如今已经被广泛使用，请阅读官网http://www.keinsci.com/research/molclus.html中的信息了解详情。简单来说，在分子构象搜索方面，用户首先需要产生一批初始构象记录到traj.xyz文件中，然后molclus程序可以按照用户的要求调用特定计算程序对这些构象做特定任务（优化/单点/振动分析等），最后通过molclus里的isostat工具进行去重和能量排序，更多信息见《Molclus FAQ》（http://sobereva.com/730）。在产生初始构象方面，有不同可行的做法，比如可以做分子动力学产生，例子见《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）和《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）；也可以通过molclus自带的gentor以用户要求的方式系统性地产生，见《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）。本文介绍的是另一种产生方式，通过专门的构象生成程序confab或Frog2来产生。用这俩程序的好处是可以傻瓜化地对有机分子产生一批质量不错、较有意义的初始构象，而且产生的构象数目可控，与molclus相结合做有机分子构象搜索十分方便，尤为适合化学直觉差、对构象搜索理解不深的人。
下面就先简单介绍一下构象生成以及Frog2和Confab的相关知识，然后给出一个与molclus联用搜索药物分子Actos优势构象的实例。
2 简谈构象生成
一般说的构象搜索，目的是获得能量最低的一批结构。而构象生成（conformer generation）主要是生成一批在低能量~中能量的构象空间中充分采样的构象，而没有化学意义的过高能量的构象被忽略，彼此间结构相似性太高的构象也被自动去掉。如果生成方式得当的话，其中一个或多个构象应当与真实的能量最低构象比较接近，或者说经过较理想级别的几何优化后应当能有至少一个结构收敛到真实的能量最低构象。另外，这些构象中也应当有与生物活性构象（如真实的蛋白质-配体复合物中的配体构象）相接近的构象，这样的话利用这些构象做分子对接就有望得到较真实的复合物结构。
用于构象生成的程序很多。对BALLOON、Confab、Frog2和RDKIT这四种免费的具有构象生成功能的程序的对比测试和简要介绍看J. Chem. Inf. Model., 52, 1146 (2012)。文章测试发现RDKIT表现较好，速度也不错，但鉴于安装麻烦一些，使用的话还得写Python脚本，会令计算化学初学者困扰，所以不在本文里提了。BALLOON结果不佳，而且我亲测发现速度很慢，本文也不提了。Confab的结果和速度都不错，一行命令就能出结果，所以本文介绍下。Frog2虽然比Confab的结果的质量差点，但速度明显更快，而且有网页版，都不需要安装任何程序，故也有一定独特的价值，所以在本文里也说一下。
使用这些构象生成程序产生给molclus用的traj.xyz，相对于其它方式有什么优点和缺点呢？这里对比一下
• 相对于gentor的优点：无需用户自行判断哪些键是需要考虑旋转的，也不需要人为设定怎么旋转，相对于“手动挡”的gentor来说这些构象生成程序算是自动挡，因此也便于结合molclus对大批体系自动做构象搜索。而且对于可旋转的键特别多，因而构象空间极大而不可能全都考虑的情况，可以由用户大致控制总共产生多少构象（当然，产生的构象越多，最终经由molclus得到真实能量最低构象的可能性越高）。
• 相对于gentor的缺点：键怎么旋转不能由用户根据实际要求和自身的化学直觉、经验进行精细操控，没法令某些键同步旋转，而且没法用于普通有机体系以外的体系（如过渡金属配合物、连接有机基团的原子团簇等）。
• 相对于做动力学的优点：用户不必额外掌握动力学程序的使用，比做动力学明显更省时省事。而且不像动力学那样最后产生的结构可能和初始用的构象有关、需要考虑如何恰当设置模拟温度，也不必纠结于当前的动力学轨迹是否对构象空间采样得够全面。
• 相对于做动力学的缺点：没法考虑环状结构的构象（至少对于Confab和Frog2来说），没法用于普通有机体系以外的体系，没法用于分子间形成的复合物。例如环糊精这种柔性环状体系包夹小分子的构型/构象搜索就明显最适合靠动力学产生traj.xyz文件，而不能靠构象生成程序。
3 Confab程序简介
Confab算法的原文见Journal of Cheminformatics, 3, 8 (2011)。Confab已经集成在了著名的化学体系文件格式转换兼化学信息学工具OpenBabel里。Open Babel最新版可以在https://github.com/openbabel/openbabel/releases免费下载。读者应当使用OpenBabel 3.0及以后的版本，不要用更老的诸如2.x版。
Confab的运行很简单。把OpenBabel安装好之后，进入操作系统的命令行模式，输入以下命令就基于test.mol2里的结构产生一批构象并一起写入到traj.xyz里了。
obabel test.mol2 -O traj.xyz --confab --verbose --conf 10000
这里--verbose要求输出详细信息，从屏幕上可以看到有多少个键被判断为可旋转的键、每个键对应哪两个原子、每个键有多少种可能的取值、最后实际产生了多少个构象。
提醒：经常有人问怎么他的机子上的OpenBabel没法正常使用以上命令产生构象。我建议在大家使用我写这篇文章时用的OpenBabel 3.1.1，下载地址：http://sobereva.com/soft/OpenBabel-3.1.1-x64.exe。并且最好不要装到默认的C盘的目录下，免得由于操作系统可能的权限限制导致OpenBabel无法成功创建新文件。
Confab本质上用的是系统式构象生成算法，但比gentor考虑得更多。首先根据原子所处的化学环境判断键的特征，指认出各个可旋转的键，并根据内置的库判断每个键都有哪些扭转角度需要考虑。假设有四个键，分别有6、3、12、4个需要考虑的扭转角度，则构象空间中总数是6*3*12*4=864个。--conf设的是实际考虑的构象数，比如设了100，那么通过随机选取扭转角的方式会得到100个随机的结构。之后程序会将能量过高的过滤掉，也通过构象间重原子的RMSD将过于相似的构象去除，以保证产生的构象的多样性。默认情况下，构象间重原子RMSD差异要求大于0.5埃，可以通过--rcutoff选项修改。默认情况下，构象产生过程中会将高于目前能量最低构象50 kcal/mol以上的构象忽略掉，可以通过--ecutoff选项修改。显然，最终产生的构象数目是少于等于--conf所设的，而且--ecutoff设得越大，最终给出的构象数可能越多。Confab是通过MMFF94力场计算的构象能量，并且在计算构象能量前，程序会自动进行优化，先优化最中间的扭转角，然后再依次优化更靠边的扭转角，这样优化耗时比起对所有变量一起优化低得多，但精度也打了折扣。
输入给Confab的文件应当含有当前体系的三维结构，具体处于什么构象无所谓，不影响产生的构象，但键长和键角不能离谱，因为Confab产生的所有构象里键长和键角都和输入结构里是一致的。也因此，不同构象之间的能量差只由范德华作用、静电作用和扭转项决定。
Confab无法产生环的不同构象，产生的构象中环的结构和输入的结构相同。如果环的构象很关键的话，应当用前述的分子动力学的方式产生给molclus用的traj.xyz。
OpenBabel也可以基于基因算法产生构象，但没有文章详细介绍其算法，我感觉也不会比Confab方法有什么实际优势（我没有实际测试），本文就不提了，介绍见
https://open-babel.readthedocs.io/en/latest/3DStructureGen/multipleconformers.html
4 Frog2程序简介
Frog2的介绍见Nucleic Acids Research, 38, W622 (2010)，程序页面见https://bioserv.rpbs.univ-paris-diderot.fr/services/Frog2/。此程序可以下载后离线使用，也可以在在线服务器上直接提交。
此程序支持的信息形式有SMILES字符串、SDF文件和mol2文件。用户可以只提供体系的一维/二维信息，也可以直接提供三维信息（立体结构文件）。对于我们构象搜索来说，一般建议用mol2格式提供三维结构信息，这样在分子的手性方面不会有含糊性。
生成构象的过程中，Frog2根据化学环境判断可旋转的键及其有代表性的角度，然后对每个产生的构象通过蒙特卡罗方式确定各个扭转角取的值。如果某个构象能量高于能量最低构象E window（默认50 kcal/mol）以上，则此构象会被忽略。如果构象的绝对能量高于E max（默认100 kcal/mol），则这个构象也会被忽略。这里说的能量是Frog2根据MMFF94力场计算的，但只考虑范德华作用部分。Frog2也根据构象间RMSD来去除相似性过高的构象。Frog2里由用户自定义#conf设置起初生成多少个构象，最终保留下来的数目是那些通过上述方式筛选的。
可以要求Frog2对产生的构象通过AMMOS工具做能量极小化，但这不会考虑静电作用。实测发现做不做这个极小化对耗时影响不大，原理上做了有好处。和Confab一样，对于提供三维结构作为输入的时候，Frog2也无法考虑环的不同构象。
Frog2的在线提交页面是https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::Frog2。如果想得到用于molclus做构象搜索的文件，在页面中，Input type选3D，Input drug description通常选mol2，点击upload标签页，点击“浏览”，然后选择分子的mol2文件。Minimize建议设为Yes，#conf根据期望得到的构象多少恰当设置，其它设置就用默认即可（另：界面里面有个Produce选项，默认是Multi，代表构象生成任务，如果切换为Single，就只给出一个能量最低构象。Disambiguate选项是对于输入一维/二维结构的时候是否区分手性异构体而言的，对于输入三维结构的情况没有用处）。点击Run提交后，过一会儿就会出现结果页面，里面若显示诸如455 conformation(s) calculated就是说最终产生了455个构象。在Frog.mol2旁边的小磁盘按钮上选“另存为”就可以把包含所有最终构象的多帧mol2文件保存下来了，之后可以放到GaussView或VMD里观看所有结构进行检查。
注意如果你的分子的mol2文件是用GaussView创建的，在提交前应当用文本编辑器打开，把里面的“Molecule Name”改成“Molecule_Name”或者随便什么名字，反正里面不能有空格才行，否则Frog2无法运行。
如果想把如上方式产生的包含所有构象的Frog.mol2转化成给molclus用的traj.xyz文件，需要运行molclus文件包里的mol2_xyz程序（从molclus 1.9.9版开始才有），输入.mol2文件路径，即会在当前目录下产生traj.xyz。
5 结合molclus做构象搜索实例：Actos
5.0 前言
这里通过Actos这个治疗糖尿病的药物在真空下的构象搜索展示一下Confab与molclus的联用的实际应用。没看过《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）的读者务必先看一下了解molclus的基本使用常识。本例相关的文件都可以在本文的文件包里下载：http://sobereva.com/attach/575/file.zip。
文件包里的Actos.pdb是我用GaussView直接画的结构保存出的文件，如下图所示
(, 下载次数 Times of downloads: 241)
可见此体系可旋转的键非常多，因此能充分体现出用Confab构象生成工具产生初始构象的价值。这里没有用Frog2产生初始构象，这是因为笔者发现至少对于这个体系，Frog2产生的初始构象不如Confab。下图左侧和右侧分别是Confab和Frog2产生的构象的叠加图（都产生了400多个的情况），对各个结构相对于吡啶部分消除了平动和转动。可见Frog2产生的初始构象遍及的构象形态明显不如Confab的那么充分，因此有漏掉一些重要构象的风险。
(, 下载次数 Times of downloads: 228)
此例的计算流程分为三步：
(1)用Confab产生几百个初始构象
(2)用crest调用xtb在GFN2-xTB方法下做批量优化，然后用isostat筛选出能量最低的一批留到下一步
(3)用molclus调用Gaussian在B3LYP-D3(BJ)/6-31G*下做优化和振动分析，并且在M06-2X/def2-TZVPP下算单点，从而得到较高精度的气相自由能
此例涉及的程序包括OpenBabel 3.0.0，molclus 1.9.9，xtb 6.3.3，crest 2.10.2，Gaussian 16 A.03。可视化使用VMD 1.9.3（http://www.ks.uiuc.edu/Research/vmd/）。
以上步骤显然还有改进余地使结果更准确。比如单点可以用ORCA在明显更高精度的DLPNO-CCSD(T)/cc-pVQZ或便宜一些的RI-PWPB95-D3(BJ)/def2-QZVPP下进行，自由能热校正量可以通过笔者的Shermo程序结合恰当的频率校正因子和准RRHO方法得到更理想的结果。相关知识见《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490）和《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）。如果想在溶剂环境下做构象搜索，xtb做优化和Gaussian/ORCA做单点计算的时候应带着溶剂模型，例如水溶剂下的构象搜索在此文里给了具体例子：《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）。
大家也可以根据实际量子化学研究经验，实际体系特点、对精度的要求、计算资源等因素对搜索流程进行改造，molclus的设计令构象搜索流程极其灵活，可控性极强。下面几节依次进行上述三个步骤。
5.1 用Confab产生初始构象
安装OpenBabel后，进入操作系统的命令行模式，把Actos.pdb放到当前目录下，然后运行
obabel Actos.pdb -O traj.xyz --confab --verbose --conf 10000
瞬间就运行完了，产生的traj.xyz在本文文件包的obabel目录下。在屏幕上会看到以下信息
..Input format = ent
..Output format = xyz
..RMSD cutoff = 0.5
..Energy cutoff = 50
..Conformer cutoff = 10000
..Write input conformation? False
..Verbose? True
**Molecule 1
..title = Actos.pdb
..number of rotatable bonds = 7
....rotor 1 from 11 to 12 has 3 values
....rotor 2 from 10 to 11 has 12 values
....rotor 3 from 9 to 10 has 3 values
....rotor 4 from 9 to 6 has 12 values
....rotor 5 from 15 to 18 has 6 values
....rotor 6 from 19 to 18 has 3 values
....rotor 7 from 3 to 7 has 12 values
..tot conformations = 279936
....Using a cutoff of 10000 we will only explore 3.6% of these
..tot confs tested = 10000
..below energy threshold = 2131
....tree size = 998 confs = 812
....new tree size = 2020 confs = 591
..generated 591 conformers
由上可见，当前体系有7个键被判断为了可旋转的键，每个键对应哪两个原子成键、有几种可能的取值都注明了。构象空间总数是所有键可取值的乘积，即3*12*3*12*6*3*12=279936，显然全都计算的话对一般条件来说肯定算不动，因此当前我们用了--conf 10000，即在这构象空间中随机生成10000个构象，占总构象空间10000/279936*100%=3.6%。显然--conf设得越大越有可能最终找到能量最低构象，但当前经过能量和相似度筛选处理后，最终产生591个构象也不算少了。大家感兴趣的话可以把当前产生的traj.xyz拖到VMD里看看，可以看到各种样子的构象都有所出现。因此，对于当前这个体系来说，基于这591个初始构象，应当已经有足够几率最终找到能量最低和最低的一批构象了。
5.2 使用xtb做初步优化和能量计算
xtb是一个能够实现GFN-xTB理论的免费程序，这相当于半经验版本的DFT方法，这一节我们用它来对上述591个构象批量做快速的预优化并得到粗略的能量。xtb的安装方法和基本用法在《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）都已经给出了。如果你的机子核数较多的话，为了用最短的时间优化完这591个构象，建议按照《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）文末介绍的技巧，通过免费、下载即用的crest程序调用xtb同时对大量体系并行优化（虽然用molclus调用xtb一个一个优化这些结构当然也是可以的，但耗时明显更高）。
假设xtb已经照常装好了，crest可执行文件是/sob/crest并且已赋予了可执行权限，Confab产生的traj.xyz已放在了当前目录下，我们现在执行下面的命令
/sob/crest -mdopt traj.xyz -niceprint
这会令crest调用xtb在默认的GFN2-xTB方法下对traj.xyz里所有的结构进行优化，优化后的所有结构连同能量会写入当前目录下的crest_ensemble.xyz文件中。在笔者的Intel双路服务器下，利用36核并行，花了7分多钟就完成了全部优化。（顺带一提，改用名义上更便宜的GFN1-xTB并不会在当前步骤节约耗时。如果使用比如-opt tight选项把默认的收敛限从verytight放宽到tight，虽然可以节约一些时间，但可能导致最终得到一些相似度非常高但又有一定差异的结构。鉴于本身这一步耗时很低，所以没必要放宽收敛限）。
现在启动molclus目录下的isostat程序，输入crest_ensemble.xyz的路径，然后能量和几何阈值都输入0.5，程序经过归簇后得到110个结构，根据能量由低到高的顺序输出到了当前目录下的cluster.xyz文件中。此文件在本文文件包里的crest-xtb目录下提供了。
isostat显示的信息如下
# 1 Count: 4 E= -72.233099 a.u. DGmin= 0.25 DE= 0.00 kcal/mol
# 2 Count: 7 E= -72.232944 a.u. DGmin= 0.30 DE= 0.10 kcal/mol
# 3 Count: 3 E= -72.228248 a.u. DGmin= 0.25 DE= 3.04 kcal/mol
...略
# 14 Count: 4 E= -72.221337 a.u. DGmin= 0.35 DE= 7.38 kcal/mol
# 15 Count: 8 E= -72.221030 a.u. DGmin= 0.13 DE= 7.57 kcal/mol
# 16 Count: 4 E= -72.220773 a.u. DGmin= 0.16 DE= 7.73 kcal/mol
# 17 Count: 2 E= -72.220672 a.u. DGmin= 0.20 DE= 7.80 kcal/mol
# 18 Count: 3 E= -72.220460 a.u. DGmin= 0.28 DE= 7.93 kcal/mol
...略
毕竟GFN2-xTB的精度非常有限，上面的能量次序和相对能量并不够可靠，必须进一步refine，这里我们就把前15个构象留到下一步，这足够多了。
5.3 使用DFT优化和算单点得到可靠的自由能
这一节我们通过molclus调用Gaussian对上一步筛选出的能量最低的15个构象在B3LYP-D3(BJ)/6-31G*下进一步做优化和振动分析，一方面令结构明显更精确，另一方面获得自由能热校正量，并且在M06-2X/def2-TZVPP下做精度较好的单点计算，将之与自由能热校正量相加就是较准确的气相自由能。
我们把上一节得到的cluster.xyz改名为traj.xyz，将molclus目录下的settings.ini里的iprog设为1（调用Gaussian），itask设为3（做优化+振动分析），ngeom设为15（处理traj.xyz的前15个），其它保持默认不变。然后编辑molclus目录下的Gaussian模板文件template.gjf，将关键词改为# B3LYP/6-31G* em=GD3BJ opt freq int=fine。然后在当前目录下创建一个名为template_SP.gjf的Gaussian做单点任务的模板文件，内容如下（记得末尾要有两个空行）。
# M062X/def2TZVPP
Template file
0 1
[GEOMETRY]
之后启动molclus，它会调用Gaussian对traj.xyz的前15个结构依次做优化和振动分析，并且由于当前目录下存在template_SP.gjf，molclus还会基于优化完的结构自动算单点能。在Intel 36核机子上6个小时算完，并在当前目录下得到了isomers.xyz文件。其中的结构是优化后的结构，用文本编辑器打开还可以看到每个结构对应的自由能（G）和自由能的热校正量（Gcorr）。这里的自由能就是B3LYP-D3(BJ)/6-31G*计算的标况下的自由能热校正量和M06-2X/def2-TZVPP单点能加和得到的。
运行isostat，按回车载入当前目录下的isomers.xyz，再按两次回车用默认的能量和结构的归簇阈值，然后就会看到下面的排序后的构象能量。在当前情况下，由于isomers.xyz里记录的是自由能，所以isostat显示的E对应的也是自由能。
# 1 Count: 1 E= -1469.136706 a.u. DGmin= 0.42 DE= 0.00 kcal/mol
# 2 Count: 1 E= -1469.136656 a.u. DGmin= 0.37 DE= 0.03 kcal/mol
# 3 Count: 1 E= -1469.136367 a.u. DGmin= 0.16 DE= 0.21 kcal/mol
# 4 Count: 1 E= -1469.136220 a.u. DGmin= 1.07 DE= 0.30 kcal/mol
# 5 Count: 1 E= -1469.136218 a.u. DGmin= 0.28 DE= 0.31 kcal/mol
# 6 Count: 1 E= -1469.135917 a.u. DGmin= 1.07 DE= 0.50 kcal/mol
# 7 Count: 1 E= -1469.135857 a.u. DGmin= 0.37 DE= 0.53 kcal/mol
# 8 Count: 1 E= -1469.135641 a.u. DGmin= 0.14 DE= 0.67 kcal/mol
# 9 Count: 1 E= -1469.135435 a.u. DGmin= 0.51 DE= 0.80 kcal/mol
# 10 Count: 1 E= -1469.134939 a.u. DGmin= 0.42 DE= 1.11 kcal/mol
# 11 Count: 1 E= -1469.134080 a.u. DGmin= 0.28 DE= 1.65 kcal/mol
# 12 Count: 1 E= -1469.133714 a.u. DGmin= 0.16 DE= 1.88 kcal/mol
# 13 Count: 1 E= -1469.133464 a.u. DGmin= 0.43 DE= 2.03 kcal/mol
# 14 Count: 1 E= -1469.133416 a.u. DGmin= 0.38 DE= 2.06 kcal/mol
# 15 Count: 1 E= -1469.132133 a.u. DGmin= 0.14 DE= 2.87 kcal/mol
在上面这部分信息之前还可以看到Nimag information was found, no isomer has imaginary frequency的提示，说明这15个结构的振动分析都没有虚频，所以上面的数据都是有意义的。
之后可以输入比如298.15让isostat计算这个温度下的Boltzmann分布，结果如下
# 1 Count: 1 Ratio: 17.88%
# 2 Count: 1 Ratio: 16.95%
# 3 Count: 1 Ratio: 12.48%
# 4 Count: 1 Ratio: 10.68%
# 5 Count: 1 Ratio: 10.66%
# 6 Count: 1 Ratio: 7.75%
# 7 Count: 1 Ratio: 7.27%
# 8 Count: 1 Ratio: 5.79%
# 9 Count: 1 Ratio: 4.65%
# 10 Count: 1 Ratio: 2.75%
# 11 Count: 1 Ratio: 1.11%
# 12 Count: 1 Ratio: 0.75%
# 13 Count: 1 Ratio: 0.58%
# 14 Count: 1 Ratio: 0.55%
# 15 Count: 1 Ratio: 0.14%
可见，当前这个体系并没有唯一的优势构象，有很多构象都有不可忽视的出现比率，这是Actos这个体系自身结构特点所决定的。自由能最低的两个构象，也即cluster.xyz的前两帧在下图给出了，虽然二者看起来构象差异极大，但从上面isostat给出的数据可以看到在当前计算级别下二者的自由能是几乎一样的。
(, 下载次数 Times of downloads: 206)
值得一提的是，这上面这些构象中，第13、14号构象都存在明显的N-H...N内氢键，虽然由于这点看似它俩应该是最稳定的构象，但是它们的自由能都比最低构象高出约2 kcal/mol，这是因为在这样的结构下存在明显的结构扭曲。例如第13号构象为了形成内氢键而导致的两处显著偏离平面的扭曲的地方在下图用箭头标注了，这带来的能量升高因素大幅抵消了内氢键对构象的稳定化作用。
(, 下载次数 Times of downloads: 212)
如果大家查看本文文件包里crest-xtb目录下的cluster.xyz，会发现第一帧就是存在内氢键的结构，即GFN2-xTB理论将这种实际上不是特别稳定的结构误当成了最稳定结构，因此仅仅靠xtb做构象搜索在精度方面是明显不够的，尤其是对于Actos这种优势构象彼此间能量差异很小的情况（涉及到溶剂的情况此问题更严重，因为xtb里的GBSA模型是极其粗糙的）。
另：对本例，建议运行molclus前将settings.ini里的ibkout设为1，这样molclus运行过程中会把每个结构的优化+振动分析以及单点任务的输出文件都备份出来，一方面便于你之后检查（比如发现有虚频的话可以看虚频模式如何），另一方面便于你之后通过Shermo程序在更好的热力学模型或其它温度/压力下批量计算这15个构象的自由能热校正量、得到Boltzmann分布比率以及权重平均的热力学数据，这仅仅需要给Shermo提供一个文件列表即可。详见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）。
6 总结
本文介绍了Confab和Frog2两种有机体系构象生成程序的原理和用法，并且结合一个高度柔性分子Actos的构象搜索演示了如何将之与molclus联用来便利地搜索优势构象。这种联用是构象空间很大、可旋转的键很多的有机分子构象搜索的一个非常理想的途径，鼓励大家在未来广泛使用。
作者Author: zsu007 时间: 2020-10-29 18:38
让大家轻松使用构象搜索，社长威武！
作者Author: winnerwill 时间: 2020-10-29 19:11
谢谢社长分享！
请教Open babel或者Confab 有并行机制吗？需要怎么设置环境变量吗？
发现只有单核运行，虽然小点的分子还挺快的。
作者Author: sobereva 时间: 2020-10-31 09:29
对此任务并行不了
作者Author: linq5203 时间: 2020-11-9 01:11
Sob老师，对于构象生成是否应该强调一下溶剂分子的作用？例如在跑动力学时加上几个水？不然到最后费时费力用DFT＋隐式溶剂筛出来的构象不仅没用还可能有严重误导性。个人被坑过，对此有些敏感。
作者Author: liyuanhe211 时间: 2020-11-9 04:19
个人经验确实如此，对一些比较柔而尺寸稍大的有机分子，不加显式溶剂的后果常常是分子过分“蜷曲”
作者Author: sobereva 时间: 2020-11-10 03:21
Confab之类构象生成程序没法考虑溶剂分子
对于分子内存在极性基团相互作用，通过跑动力学来产生初始构象，是可以如此
作者Author: xxzj 时间: 2020-12-11 21:11
老师，我看见贴子中有这样一句话“ 相对于做动力学的缺点：没法考虑环状结构的构象（至少对于Confab和Frog2来说），没法用于普通有机体系以外的体系，没法用于分子间形成的复合物。例如环糊精这种柔性环状体系包夹小分子的构型/构象搜索就明显最适合靠动力学产生traj.xyz文件，而不能靠构象生成程序。”对于两个分子构成的复合物寻找最稳定的构型是不是就不可以用了，例如卟啉二聚体，还是类似环糊精这种柔性环状体系不可以使用，谢谢老师
作者Author: sobereva 时间: 2020-12-13 04:34
confab这种程序是对于单分子而言的，显然不能用于产生二聚体
如果你要考虑环糊精的环的构象，也不能用confab
作者Author: hunterpyj 时间: 2020-12-14 14:55
请问： centos下安装完open babel 后 运行obabel Actos.pdb -O traj.xyz --confab --verbose --conf 10000
只显示 1 molecule converted，并没有调用confab 运算，问题出在哪里呢？
作者Author: sobereva 时间: 2020-12-15 02:28
Openbabel版本太老。自行装个新的。不方便装就用Windows版
作者Author: 新时代农民工 时间: 2021-2-23 22:08
请教社长，可旋转的键太多，conf设置了10000，然后命令行显示 using cutoff of 10000 we will only explore 0.0%of these ,怎么办？
还有就是如果分子本身是helix的构象可以用confab 么？
之前用薛定谔搜索过，不知道薛定谔的macromodel 社长觉得怎样？
作者Author: sobereva 时间: 2021-2-24 12:10
能算得动多少就算多少
我不清楚具体是什么情况，最好给个结构图
macromodel对于构象搜索不会比用此文介绍的方式更好
作者Author: Notoschord 时间: 2021-2-24 23:16
请教老师，在Open Babel GUI下应该如何使用confab呢？
我点选"Confab, the diverse conformer generator"之后，
再点选"CONVERT"，总是显示"0 molecules converted"。
而且，GUI下也不知道在哪里调整confab的参数……
作者Author: sobereva 时间: 2021-2-26 14:56
我从不通过GUI运行此功能，没必要，纯粹自找麻烦
作者Author: 苏玖染 时间: 2021-5-13 14:14
sob老师，Confab或Frog2，这两种方法对于阴离子体系，是不是无法处理？
作者Author: sobereva 时间: 2021-5-13 19:18
没试过，试试便知，理应没这个限制
作者Author: 苏玖染 时间: 2021-5-14 09:33
试了一下，确实不能算阴离子，默认会补氢变成中性分子计算
作者Author: 花花向日葵 时间: 2021-7-2 10:18
你好，请问运行molclus前将settings.ini里的ibkout设为1，发现molclus运行过程中只输出了每个结构的优化+振动分析文件，但没有单点任务的输出文件
作者Author: sobereva 时间: 2021-7-3 04:15
本来就是这样
作者Author: 花花向日葵 时间: 2021-7-3 11:06
本帖最后由 花花向日葵 于 2021-7-3 11:11 编辑
好的，谢谢
作者Author: qingqing 时间: 2021-7-21 09:58
两个归簇阈值设置不同会不会被审稿人diss。我现在的情况是要对大量分子构象搜索，xtb计算的能量归簇阈值设小0.5 kcal/mol，0.25 angstrom的话，后续g16和orca计算大量结构耗时非常高，所以目前设置的第一次归簇阈值为0.5 kcal/mol，0.4 angstrom, 第二次归簇为0.5 kcal/mol，0.25 angstrom. 但是若审稿人问起，要怎样回复比较好呢？谢谢
作者Author: sobereva 时间: 2021-7-22 00:46
两个阈值物理意义、单位都完全不同，没有任何理由要相同
作者Author: qcn1211 时间: 2021-8-11 11:58
Frog2 构象搜索后，在Frog.mol2旁边的小磁盘按钮上选“另存为”显示如下图所示界面，请问如何得到多个独立的帧文件？
作者Author: sobereva 时间: 2021-8-12 01:55
不要说什么小磁盘，都不知道你用的什么文本编辑器
如果给出的是多帧的mol2文件，自己写个脚本拆分
作者Author: zhouwenhao 时间: 2021-8-16 16:28
我和你遇到了同样的情况，现在也不晓得是什么原因，请问你解决了嘛
作者Author: 丁越 时间: 2021-12-16 19:33
本帖最后由 丁越 于 2021-12-16 19:34 编辑
这是因为你编译open babel的时候没有编译eigen库。这一点官网上就没有提及，我也是查了一下才知道的。
https://www.cnblogs.com/wq242424/p/8231600.html
作者Author: 无敌菠萝派 时间: 2022-5-28 19:30
请教一下老师，为什么说这些构象程序不适合用于过渡金属配合物呢，对于一些未知配位键长键角甚至金属离子与配体上具体哪个原子配位都不确定的情况下，好像gentor也不太适用，请问这种时候选择其他哪种构象程序比较好呢
作者Author: sobereva 时间: 2022-5-29 01:25
你这种问题明显有更好的办法
拿到一个配体结构，通常通过基本化学常识，或者其它已知的含有此配体的化合物，就知道哪里可能配位，因此可以手动摆出来大致有意义的初猜结构。而且实在不好确定时，还可以用Multiwfn对配体表面计算静电势极值点、ALIE极值点辅助判断，看：
使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用
http://sobereva.com/159
使用Multiwfn和VMD绘制平均局部离子化能(ALIE)着色的分子表面图（含视频演示）
http://sobereva.com/514（http://bbs.keinsci.com/thread-14632-1-1.html）
使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频演示）
http://sobereva.com/443（http://bbs.keinsci.com/thread-11080-1-1.html）
通常根本没必要靠随机产生初始结构方式来搜索过渡金属配合物结构，效率太低，而且这么搞出来的初始结构往往SCF收敛很困难。
作者Author: 无敌菠萝派 时间: 2022-5-29 10:31
谢谢sob老师，我的体系就是因为金属离子的配位能力强，配体中又有多个可能与其配位的原子，所以我没法直接判断，我会仔细研究一下您贴出的链接内容。
作者Author: 无敌菠萝派 时间: 2022-6-2 17:54
本帖最后由 无敌菠萝派 于 2022-6-2 18:01 编辑
sob老师，我还想问一个问题，我用文章中ALIE绘制出了分子表面图和极值点，但是我还想知道极值点和里面某个原子的距离该怎么操作呢
作者Author: sobereva 时间: 2022-6-3 03:36
Mouse - Label - Bonds，点击极值点和一个原子的中心进行测量
作者Author: 无敌菠萝派 时间: 2022-6-3 10:22
多谢sob老师
作者Author: laplacesuanzi 时间: 2022-10-7 22:46
本人新手，这个程序怎么打不开呢，求各位赐教，谢谢
作者Author: sobereva 时间: 2022-10-8 17:12
卸了当前的Openbabel，尝试用我目前用的这个版本：http://sobereva.com/soft/OpenBabel-3.1.1-x64.exe
并且把程序装到C盘以外的目录下避免由于权限原因无法创建新文件
作者Author: laplacesuanzi 时间: 2022-10-18 07:48
问题已解决，谢谢老师。
作者Author: anlancx 时间: 2022-10-18 17:03
老师，请问在5.3 使用DFT优化和算单点得到可靠的自由能里，是否有办法能够让15个结构同时优化？
作者Author: sobereva 时间: 2022-10-18 17:26
把molclus目录复制多份，settings.ini里ngeom指定处理traj.xyz里的帧号范围，同时运行多个molclus，最后把各个isomers.xyz合并即可
作者Author: qcn1211 时间: 2022-12-31 14:40
sob老师，请问我用confab进行构象搜索出现如图片所示，是什么问题呢？
作者Author: sobereva 时间: 2022-12-31 17:27
zt.mol2不在当前目录下
输入dir看当前目录都有什么
作者Author: 蟹黄宝 时间: 2023-3-15 10:38
本帖最后由 蟹黄宝 于 2023-3-15 10:40 编辑
我使用 obabel 0.pdb -O traj.xyz --confab --verbose --conf 10000 处理一个86个原子（12聚合度）的全氟聚醚分子，出现这种报错：*** Open Babel Error in operator()
memory limit exceeded...
我找网上和Open Babel 教程上没找到内存设置的地方，
请教老师，这该怎么解决？
作者Author: 护卫天使 时间: 2023-5-27 13:21
老师您好，我用Confab搜索了图中单体的构象，其中该单体是考虑N的去质子化，因为文献报道该聚合物-NH的pKa为6.6左右，我用xtb程序的-gfn2 -gbsa h2o -chrg -1 方法进行预优化，结果经过Molclus筛选得到的低能构象中羟基的氢会迁移到N原子上，请问老师我后面要如何处理？
作者Author: sobereva 时间: 2023-5-28 02:04
先取一个结构，用wB97XD/6-311G**结合IEFPCM描述的水环境看看会优化成什么样，是否质子会转移。如果不会的话，考虑xtb优化时加上约束避免质子转移，或者尝试调用MOPAC在PM6-D3H4结合隐式溶剂模型下做优化看看
作者Author: 美泽虎 时间: 2023-6-24 00:54
请教老师， 我运行Confab产生了1249个conformers，但最下面显示 0 molecules converted （见图片附件），这种是正常的吗？用的是最新的windows版。 3.1.1 x64那个
(, 下载次数 Times of downloads: 85)
作者Author: sobereva 时间: 2023-6-24 01:30
不正常。先尝试瑞德西韦，看是不是程序本身的毛病
作者Author: 美泽虎 时间: 2023-6-24 02:01
本帖最后由 美泽虎 于 2023-6-24 06:50 编辑
Sob老师，我尝试了一下用文中的Actos.pdb这个例子。除了最后一行多了一个0 molecules converted，其他的输出信息都和文中一致（见下图）。C:\Users\14047\OneDrive\Desktop\Weixin Image_20230623135650.png，之后产生的xyz文件也可以正常在VMD里面load出591个conformer（见下图）C:\Users\14047\OneDrive\Desktop\Weixin Image_20230623140028.png
作者Author: sobereva 时间: 2023-6-24 04:55
贴图方式不对，仔细看指定的社员必读贴
如果用相同的步骤，actos产生的xyz文件里有正常的构象，你的体系没有产生哪怕一个，那就是bug，建议找开发者解决，或者用别的方式产生初始构象，如gentor、frog2、跑动力学等
作者Author: xiaokexxk 时间: 2023-12-5 09:28
本帖最后由 xiaokexxk 于 2023-12-5 16:27 编辑
老师你好，请问运行molclus前将settings.ini里的ibkchk设为1，发现molclus运行过程中只输出了每个结构的优化+振动的chk文件，没有单点能任务的chk文件，是默认如此吗？
作者Author: sobereva 时间: 2023-12-6 07:28
是
作者Author: ldsserver 时间: 2024-1-24 10:49
请教老师，confab的MMFF94没有B，有办法改成UFF吗，obabel -L confab里边没有ForceField参数
作者Author: sobereva 时间: 2024-1-25 03:55
因为硼的成键形式复杂、情况太多，导致支持硼的力场极少。虽然UFF支持硼，但也不意味着用UFF就能优化出合理结果（虽然你可以用Gaussian做UFF优化一个构象看看情况）。建议至少用半经验做构象筛选
作者Author: ldsserver 时间: 2024-1-31 10:36
感谢老师！
作者Author: 护卫天使 时间: 2024-3-30 21:05
老师好，安装了Confab，使用时报错：Missing or unknown output file or format spec or possibly a misplaced option.但是不知道如何解决，麻烦老师指教。
作者Author: sobereva 时间: 2024-3-30 23:32
先尝试用帖子里相同的输入文件看是否正常
如果正常，把输入文件改成mol2格式再试
作者Author: 护卫天使 时间: 2024-3-31 10:45
老师好，尝试了博文中的相同的输入文件、重新安装、换电脑、.mol2为输入格式，都是一样的报错。改变输出文件的格式为.pdb或者.mol2，能运行，但是显然得不到结果，请问一下老师，这怎么解决呢？
作者Author: sobereva 时间: 2024-3-31 15:01
把pdb轨迹载入VMD，然后保存成xyz轨迹
作者Author: 2877321934 时间: 2024-8-23 14:16
本帖最后由 2877321934 于 2024-8-23 15:52 编辑
请问老师，我将xtb优化的结构用isostat筛选后，都用mopac优化了一遍,用isostat再筛选一遍后，再取分布排名前15的结构进行dft优化，和直接取xtb生成并用isostat筛选后分布前15的结构进行dft优化比较。发现两种方法得到的优势构象完全不一样，请问用mopac多筛选一步会使结果更可靠吗，哪一种方法得到的构象更可靠？
作者Author: sobereva 时间: 2024-8-24 06:24
此时用MOPAC完全多余
本来GFN2-xTB的平均精度就比MOPAC支持的任何半经验方法要高
作者Author: 2877321934 时间: 2024-8-24 11:48
好的，谢谢老师
作者Author: zhangjiacheng 时间: 2024-8-27 09:37
我也遇到了同样的问题，请问你最后解决了吗？
作者Author: 2877321934 时间: 2024-8-31 00:12
追问老师一个问题，我将xtb优化的结构用isostat筛选后，用mopac优化了一遍，再用dft优化，得到了两个构象，分别占50%和40%，但是，直接取xtb生成并用isostat筛选后分布前15的结构进行dft优化，只有一个优势构象，占100%，且我发现我使用后者方法得到的占100%的优势构象就是前者方法中占40%的构象，这是否说明我使用后者方法漏掉了一种优势构象，前者方法得到的占50%的构象也应当在研究中考虑
作者Author: sobereva 时间: 2024-8-31 06:02
xtb跑过一遍后就没有任何意义用MOPAC
要么用MOPAC，要么用xtb，没有任何结合使用的意义。通常首选molclus结合xtb做GFN2-xTB的计算，这是半经验级别里整体来说相对最靠谱的。
作者Author: 2877321934 时间: 2024-8-31 12:04
好的，谢谢老师
作者Author: 2877321934 时间: 2024-9-3 00:25
本帖最后由 2877321934 于 2024-9-3 00:55 编辑
请问老师，我使用confab生成了100多个结构，再按照您说的用xtb优化，但是我的研究体系是在二氯甲烷溶液中的，因此我按照您的方法在xtb中优化后，取了前15个构象进行dft优化，此时我使用iefpcm溶剂模型，但是后来想到，这样做不一定合理，因为前面在xtb中优化时没考虑溶剂模型，而溶剂模型对构象分布的影响很大，已知我的体系对溶剂效应十分敏感，我将xtb在气相下优化好的结构分别在气相和二氯甲烷中做dft，在气相中占100%的结构到了二氯甲烷中只占4%，因此我的体系可能在二氯甲烷中的实际优势构象在xtb气相优化后不被排到前15名，因此我想在xtb优化时就使用溶剂模型，请问老师是否有必要在xtb优化时考虑溶剂
作者Author: sobereva 时间: 2024-9-3 11:40
优化过程本身不用，但如果你直接对xtb算的能量在isostat里排序来做筛选的话，那不考虑溶剂效应可能明显不合理。此时要么优化时就带着溶剂模型，要么对优化后每一帧再让molclus调用xtb带着溶剂模型算单点。
作者Author: 2877321934 时间: 2024-9-3 15:54
好的，谢谢老师，倘若在优化过程加溶剂应该如何操作呢
作者Author: sobereva 时间: 2024-9-4 11:14
按xtb手册里说的加上溶剂模型设定就完了
作者Author: 2877321934 时间: 2024-9-4 15:55
ok 感谢
作者Author: 2877321934 时间: 2024-9-5 03:25
我查阅了站内信息得到gfn-2-xtb支持的两种溶剂模型gbsa和alpb都很粗糙，但聊胜于无，请问您觉得哪种在研究构象分布中更为可靠呢
作者Author: sobereva 时间: 2024-9-5 10:56
原理上来说ALPB更先进，对于低极性溶剂原理上更好。
作者Author: qingran 时间: 2025-6-17 19:30
本帖最后由 qingran 于 2025-6-17 21:20 编辑
谢谢老师！
| 欢迎光临 计算化学公社 (http://ccc.keinsci.com/) |
Powered by Discuz! X3.3 |