计算化学公社
标题: 使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索 [打印本页]
作者Author: sobereva 时间: 2020-2-8 15:57
标题: 使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索
提醒：如果你要构象搜索的分子的可旋转的键不多，比如7个键及以内，强烈建议用gentor+molclus而不要用本文的方法！因为过程明显更省事、不用跑动力学，而且能更充分确保构象空间采样的完整。做法看《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）。仅有可旋转的键特别多的时候才适合用下文的方法。
使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索
Conformational search of Remdesivir using Molclus combined with xtb dynamics simulation
文/Sobereva@北京科音 First release: 2020-Feb-8 Last update: 2024-Nov-11
0 前言
Molclus是免费灵活易用的团簇构型与分子构象搜索程序，如果不熟悉的话，请先看官网http://www.keinsci.com/research/molclus.html和《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）中的介绍。Molclus做构型/构象搜索需要提供含有初始结构信息的多帧xyz文件（通常叫traj.xyz），如笔者写的多篇Molclus教程所示，这样的文件可以通过其它程序跑分子动力学（MD）程序生成，也可以通过Molclus自带的genmer或gentor生成。本文的目的是通过一个实例，演示用Molclus基于xtb程序通过GFN-xTB方法（一种半经验形式的DFT）跑出的分子动力学轨迹做高精度分子构象搜索。
如果你还对xtb一无所知，看《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）中的相关介绍。安装方式在里面都写明了。
本文的做法相对于以其它方式产生traj.xyz文件有以下好处：
(1)相对于用经典力学程序跑MD，用xtb跑MD极其简单，不需要经过准备拓扑文件、补参数之类麻烦的过程，初学者可以轻松上手。而相对于用量子化学程序跑从头算动力学，用xtb基于GFN-xTB理论来跑能快几个数量级，用比较烂的机子对常见大小的有机分子也能跑出足够长的轨迹。
(2)虽然Molclus的gentor也可以利用系统式扫描的方式产生含有分子各种初始构象的文件，但对于含有构象可变的环状区域的体系（比如环己烷有椅式和扭船式两种极小点），gentor没法考虑环的构象。而跑动力学的话，环状区域的各种构象是可以自然而然地被采样的。
(3)对于可旋转的键非常多的体系（比如明显超过10个键），靠gentor的话要考虑的初始构象过多。而靠动力学的话，可以通过跑的时间长度来控制最终能找到能量最低构象的概率。在参数合适的前提下，跑的时间越长，最终能找到能量最低构象的概率就越高。
注：跑动力学产生traj.xyz的缺点是如果跑的时间不够长或条件不适合，最终得到的能量最低构象会对跑动力学的初始构象有依赖性，而且有随机性（相同流程两次跑出来的结果可能不同），因为此时不像gentor系统式扫描那样可以基本保证势能面所有势阱都能被采样到。
本文下面的例子涉及的文件都可以在这里下载：http://sobereva.com/attach/532/file.rar
本文例子使用的各个程序和版本：Molclus 1.9.2，xtb 6.3pre2、Gaussian16 A.03、ORCA 4.2.1。机子是Intel 2*XEON E5-2696v3（36核），系统是CentOS 7.3 64bit。
2 被考察的体系：瑞德西韦
本文的例子是对瑞德西韦(Remdesivir)在水环境中做构象搜索，确定能量最低的和最低的一批构象。此分子是美国吉利德公司研制的一种用于抗埃博拉病毒的药物，但效果一般，在2020年1月26日被美国医生用于治疗一个患有发源于武汉的肺炎（2019-nCoV）病毒的感染者，发现效果极佳，之后吉利德公司向中国超低价慷慨援助了大量瑞德西韦。此分子的结构如下
(, 下载次数 Times of downloads: 537)
获得此体系三维结构很简单。去chemspider（https://www.chemspider.com）搜Remdesivir，然后在结构图下面点3D按钮，再点Save按钮就得到了.mol文件。里面是基于二维结构经验生成的三维结构。
之后用xtb做计算时需要基于xyz格式的文件才行。把mol文件转成xyz文件办法很多，比如Multiwfn就可以。Multiwfn可以在其官网http://sobereva.com/multiwfn免费下载。启动Multiwfn，载入本文文件包里的Remdesivir.mol，然后输入100，再选子功能2，再选择导出xyz文件的选项，然后输入导出的文件名Remdesivir.xyz，文件就马上输出在了当前目录下。
3 计算流程说明
对瑞德西韦这个体系做构象搜索流程如下。虽然过程看起来有点多，但这是保证结果高度可靠的前提下耗时最低的流程。有个几十核的像样服务器的话不超过一天时间就能搞完。
(1)用xtb在GFN0-xTB下跑动力学
解释：GFN0-xTB是GFN-xTB系列中非常便宜的一个版本，由于形式简单、不需要SCF迭代，而且xtb代码效率极高，因此比一般的半经验方法耗时还低得多得多，跑小分子的动力学极其适合。虽然GFN0-xTB计算的能量并不算太准确，但用于动力学采样的目的足够可靠了。模拟温度建议用比常温明显要高的400K，使得模拟过程构象变化较快，在有限的模拟时间内能更充分采样。模拟时间设100 ps，对于这样的小分子而言长度基本够了（但有条件的话建议跑几百ps）。
后记：xtb后来的版本支持了GFN-FF力场，耗时比GFN0-xTB更低，而计算普通有机体系的构象能量差比GFN0-xTB的精度还更好点。对于动力学过程不涉及成键/断键的情况，也很推荐用GFN-FF跑动力学。
(2)通过Molclus调用xtb，在GFN0-xTB下对动力学的每一帧进行批量优化
解释：过程(1)得到的动力学轨迹通常含有几千帧，我们需要对每一帧进行几何优化然后对比能量。但这么多帧都直接用像样的方法做几何优化是极度困难的。实际上有很多帧在优化之后会收敛到相同的结构，因此可以先用极快的GFN0-xTB对所有帧进行优化，之后用Molclus里的isostat去重后，也就剩几百个结构。
(3)通过Molclus调用xtb，对(2)产生的每一个结构用GFN2-xTB在溶剂模型下进行批量优化
解释：目前还有几百个结构，都直接用DFT优化还是太昂贵。因此先用比GFN0-xTB更可靠但也更贵的GFN2-xTB进行批量优化。虽然其耗时是GFN0-xTB的2~20倍，但仍低于DFT好几个数量级。此时应当考虑GBSA溶剂模型，这使得优化后得到的能量也是有意义的。用isostat对优化后的结果进行排序，选择能量最低的一批保留下来。
(4)通过Molclus调用Gaussian，对(3)筛选出来的结构用B3LYP-D3(BJ)/6-31G*在真空下做优化和振动分析得到较可靠的结构和自由能热校正量，再调用ORCA对每个结构在PWPB95-D3(BJ)/def2-QZVPP结合SMD模型表现的水环境下算高精度单点能，二者相加得到水环境下的高精度自由能
解释：现在的候选结构已经比较少了，可以用DFT做优化和振动分析、再用高精度级别算单点能了。做振动分析一方面是为了获得自由能的热校正量，一方面是检验优化出的结构有无虚频，有虚频的结构肯定不能要。B3LYP-D3(BJ)/6-31G*是做几何优化和振动分析的最低可接受级别了，虽然看起来略low，但对于优化和振动分析目的足够了。用ORCA里的B97-3c做优化和振动分析也可以，精度通常也够，由于用的是纯泛函，耗时还更低。这个过程不必考虑溶剂模型，因为哪怕是极性溶剂，对一般体系的几何结构影响也非常微小，而且用溶剂模型后往往几何优化收敛所需步数会增加、收敛更困难一些，耗时也高一点。对于局部呈现显著离子性的体系则优化和振动分析时也要考虑溶剂模型，否则结果可能显著和实际情况不符。之所以高精度单点部分改用ORCA程序，是因为ORCA做PWPB95-D3(BJ)/def2-QZVPP这样双杂化泛函结合很大基组的计算，结果又准确而耗时也不是特别高，而Gaussian里双杂化泛函结合这样的基组算当前体系绝对算不动。用PWPB95-D3(BJ)的理由我在此文写明了：《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）。另外，使用顶尖的普通泛函wB97M-V也完全可以，算构象能量差问题的精度整体不比PWPB95-D3(BJ)差，还能便宜一倍。
PS 1：在ORCA里用DLPNO-CCSD(T)算单点可以得到更准确的结果，但对于当前这个77个原子的不小体系来说，其结合较好带弥散函数的3-zeta基组（如jun-cc-pVTZ）或者def-QZVPP基组时的耗时和内存、硬盘消耗都很高，所以不用于本文例子。但如果你的体系明显小于本文的例子，或者对本文大小的体系你能忍较长的耗时的话也可以用DLPNO-CCSD(T)。若用NormalPNO版本的DLPNO-CCSD(T)结合cc-pVTZ的话，在笔者的Intel 36核机子上算瑞德西韦的单点的耗时是PWPB95-D3(BJ)/def2-QZVPP的两倍，为两个半小时。但这样不带弥散函数的3-zeta基组对于描述分子内弱相互作用还是差点意思，所以不太推荐用。
PS 2：虽然用Spartan、Macromodel、GMMX等一些傻瓜式的构象搜索程序基于分子力场做构象搜索仿佛看起来又简单又快，但由于主流分子力场精度太low，结果可靠性极低，根本没多大实际意义，而且也不像本文的做法那样可以用于包含几乎所有元素的任意类型的体系、能合理地考虑溶剂效应。在Int. J. Quantum Chem., 118, e25512 (2018)文中通过大规模测试发现哪怕是MMFF94这样的高精度有机分子力场用于算有机体系，算的构象相对能量也颇不理想，对能量最低一批结构的排序也没多大可信度，远不如PM7等主流半经验方法。虽然Grimme也搞了个名为crest的利用xtb跑metadynamics的构象搜索工具，但从实效上看并不比本文的做法更有任何优势，过程也不像本文做法那么透明和可控，且最终也还是得借助Molclus用更高级别来refine才行，所以我不认为crest有多大实际使用价值。
4 第一步：用xtb在GFN0-xTB下做动力学模拟
xtb程序做动力学的超级详细的介绍在北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）中的“从头算分子动力学”部分讲授，并给出了十分丰富的例子，欢迎参加，这里只是很简单说一下xtb跑动力学实现当前目的的做法。创建一个含有xtb的MD任务设置的文件md.inp，内容如下，双斜杠后面的是注释。
$md
temp= 400 //温度(K)
time= 100.0 //模拟的总时间(ps)
dump= 50.0 //每多少fs往轨迹文件里写入一次
step= 1.0 //步长(fs)
hmass=1 //氢原子的质量是实际的多少倍
shake=1 //将与氢有关的化学键距离都用SHAKE算法约束住
$end
这里之所以把hmass明确设为1，是因为默认设置下hmass是4，即氢原子质量是实际的四倍，目的是减缓氢的运动，从而能用更大步长（即跑同样的时间可以用更少的步数），但这终究会给动力学行为带来人为虚假性，一般做MD也很少这么干。步长用的是1 fs，这是很稳妥的。其实由于用了shake=1，步长设2 fs也不是不行，但考虑到当前是较高的温度下模拟，原子运动较剧烈，为了确保动力学稳定性，还是用了保守的1 fs。保存间隔设50 fs还是比较合适的，设得太小的话相邻的帧结构太相似，之后要处理的不必要的结构太多；而设得太大的话，MD过程中的很多有意义的结构又会漏过去。模拟总长是100 ps = 100000 fs，因此最终得到的轨迹包含100000/50=2000帧。
确保xtb已正常安装了，然后运行xtb Remdesivir.xyz --input md.inp --omd --gfn 0。其中--omd告诉程序先做几何优化，然后再跑MD。做几何优化是为了避免由于初始结构太烂，导致一开始某些原子受到过大的斥力使其速度太大而令动力学崩溃。--gfn 0代表用GFN0-xTB。
MD阶段一开始会显示模拟的速度：
est. speed in wall clock h for 100 ps : 0.68
即模拟每100 ps耗时0.68小时，当前我们要跑的总长度也正好是100 ps，因此预计总耗时为60*0.68=41分钟。在笔者的Intel 36核服务器上，此任务最后花了42分钟跑完。计算过程中当前跑到第多少帧了会在屏幕上直接提示，并且不断往当前目录下的xtb.trj文件中写入结构。为了监控模拟的合理性，在任务跑的中途就可以把xtb.trj拷出来并改名为xtb.xyz，然后拖到免费的可视化程序VMD（http://www.ks.uiuc.edu/Research/vmd/）的VMD main窗口里观看已跑出来的轨迹，如果发现明显不合理，就应该赶紧停掉。每经过1 ps模拟，就会在当前目录下产生scoord.为开头的文件记录了此时的结构。
将最终的xtb.trj改名为traj.xyz，用于下一阶段的计算。此文件对应于本文文件包里xtb_opt目录下的traj.xyz，如果看这个轨迹的话，就会发现分子不断大幅运动，因此这条轨迹的采样足够充分了（另外，跑完后当前目录下还会有mdrestart，记录了最后一帧的坐标和速度，如果感觉跑的还不够长，此文件可以用于续跑）。
5 第二步：用Molclus调用xtb在GFN0-xTB下做批量优化
去Molclus官网下载最新版本并解压。把traj.xyz放到Molclus目录下。把Molclus目录下settings.ini里的iprog设为4代表调用xtb，确保itask为0代表是做优化任务。把xtb_arg参数设为"--gfn 0 --chrg 0 --uhf 0"，代表这是中性单重态体系且用GFN0-xTB计算。之后输入./molclus命令启动当前目录下的Molclus程序可执行文件，程序就开始调用xtb优化traj.xyz里的每一帧了。在笔者的Intel 36核机子上优化每一帧结构的耗时也就2~5秒钟，总共花了111分钟。
（重要技巧：这个批量优化过程可以利用crest程序来降低一个数量级的耗时，见本文最后一节的说明）
任务结束后，优化后的每一帧的结构和能量就都存到当前目录下的isomers.xyz里了。使用isostat对isomers.xyz里的结构去重和做能量排序，即输入
./isostat
[回车] //处理当前目录下的isomers.xyz
0.5 //能量去重阈值(kcal/mol)
0.5 //结构去重阈值(埃)
[回车] //不计算Boltzmann分布比例
期间在屏幕上会看到下面的信息：
The lowest energy is -122.66129225 a.u.
Sorting clusters according to energy...
# 1 Count: 1 E= -122.661292 a.u. DGmin= 0.23 DE= 0.00 kcal/mol
# 2 Count: 3 E= -122.660372 a.u. DGmin= 0.28 DE= 0.58 kcal/mol
# 3 Count: 5 E= -122.660067 a.u. DGmin= 0.17 DE= 0.77 kcal/mol
# 4 Count: 1 E= -122.659769 a.u. DGmin= 0.33 DE= 0.96 kcal/mol
# 5 Count: 3 E= -122.659587 a.u. DGmin= 0.20 DE= 1.07 kcal/mol
# 6 Count: 1 E= -122.659178 a.u. DGmin= 0.22 DE= 1.33 kcal/mol
...略
# 268 Count: 2 E= -122.637199 a.u. DGmin= 0.08 DE= 15.12 kcal/mol
# 269 Count: 3 E= -122.637073 a.u. DGmin= 0.33 DE= 15.20 kcal/mol
# 270 Count: 4 E= -122.636938 a.u. DGmin= 0.44 DE= 15.28 kcal/mol
# 271 Count: 1 E= -122.636612 a.u. DGmin= 0.28 DE= 15.49 kcal/mol
# 272 Count: 1 E= -122.635534 a.u. DGmin= 0.33 DE= 16.16 kcal/mol
最后得到的共272个非重复结构产生在了当前目录下的cluster.xyz里。删掉之前的traj.xyz，然后把cluster.xyz改名为traj.xyz，使得这些去重后的结构作为接下来Molclus任务的输入结构。
6 第三步：用Molclus调用xtb在GFN2-xTB结合隐式水模型下做批量优化
把settings.ini里的xtb_arg参数设为"--gfn 2 --gbsa h2o --chrg 0 --uhf 0"，代表用GFN2-xTB级别结合GBSA模型表现的水环境进行计算。启动molclus来调用xtb对当前272个结构进行优化。优化每个构象耗时十几秒钟，整个任务耗时48分钟（这个批量优化过程可以利用crest程序来降低一个数量级的耗时，见本文最后一节的说明）。之后运行isostat进行去重和排序，做法同前，输出如下：
# 1 Count: 1 E= -128.462874 a.u. DGmin= 0.29 DE= 0.00 kcal/mol
# 2 Count: 2 E= -128.462743 a.u. DGmin= 0.25 DE= 0.08 kcal/mol
# 3 Count: 2 E= -128.460194 a.u. DGmin= 0.31 DE= 1.68 kcal/mol
# 4 Count: 1 E= -128.459966 a.u. DGmin= 0.50 DE= 1.82 kcal/mol
# 5 Count: 1 E= -128.459373 a.u. DGmin= 0.24 DE= 2.20 kcal/mol
# 6 Count: 1 E= -128.458599 a.u. DGmin= 0.24 DE= 2.68 kcal/mol
# 7 Count: 2 E= -128.457983 a.u. DGmin= 0.22 DE= 3.07 kcal/mol
# 8 Count: 2 E= -128.457949 a.u. DGmin= 0.17 DE= 3.09 kcal/mol
# 9 Count: 3 E= -128.457877 a.u. DGmin= 0.18 DE= 3.14 kcal/mol
# 10 Count: 9 E= -128.457797 a.u. DGmin= 0.36 DE= 3.19 kcal/mol
...略
目前非重复的一共有148个结构。虽然当前的能量已经有一定实际实际意义了，能量高低次序已经能一定程度体现不同构象的能量高低了，但这个体系柔性高，能量最低的一批构象的能量肯定很相近，靠GFN2-xTB的精度还是明显不足以可靠地区分能量差很小的不同构象的能量相对高低。而且xtb用的GBSA溶剂模型也远不如量子化学程序用的PCM等溶剂模型可靠。而且当前也只是考虑的电子能量，没考虑自由能。因此，之后还需要用更好的级别做几何优化、振动分析、单点计算来得到准确的溶液下的自由能，这就是接下来要做的事。
GFN2-xTB在GBSA模型下给出的相对能量虽然不算很准，但也定性正确，分辨几个kcal/mol的差异还是没问题的。我们这里就保留相对能量<=3 kcal/mol的6个结构用于下一步的精炼（这个标准是随意设的。假设对于你的情况<=3 kcal/mol的都有50个结构，还是觉得太多，那可以比如只保留能量最低的20个。保留多少要结合实际计算条件、对精度要求等方面自行权衡。为了稳妥起见，条件允许的话我建议保留所有<=4或<=5 kcal/mol的构象）。
还是将traj.xyz删掉，把cluster.xyz改名为traj.xyz。
7 第四步：用Molclus调用Gaussian和ORCA得到水环境下的高精度自由能
把settings.ini里的iprog设为1代表调用Gaussian。ngeom设为6代表只处理traj.xyz里的前6个结构（即只考虑上一步得到的能量最低的6个结构，因为isostat给出的cluster.xyz是按照能量由低到高排序的）。itask设为3代表做优化+振动分析来得到自由能。把gaussian_path和orca_path分别设为调用机子里Gaussian和ORCA的实际命令。另外建议把ibkout设为1代表把算每个体系的输出文件都进行备份，便于之后必要时进行检查。
把molclus目录下Gaussian模板文件template.gjf里的关键词设为B3LYP/6-31G* em=GD3BJ opt freq。代表在B3LYP-D3(BJ)/6-31G*下做优化和振动分析。PS：如果你是Gaussian 16的用户，建议加上int=fine，能节约很多时间。默认的int=ultrafine对于B3LYP-D3(BJ)来说太浪费了。
在molclus目录下创建一个叫template_SP.inp的ORCA输入文件的模板文件，内容如下，对应在RI-PWPB95-D3(BJ)/def2-QZVPP结合SMD模型表现的水环境下做高精度单点计算。如果对ORCA的这些关键词、计算级别不熟悉，不知道为什么这么写的话，强烈建议看《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490）。
! PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX grid4 gridx4 tightSCF noautostart miniprint nopop
%maxcore 6000
%pal nprocs 36 end
%cpcm
smd true
SMDsolvent "water"
end
* xyz 0 1
[GEOMETRY]
*
（注：从ORCA 5.0开始，以上输入文件里的grid4 gridx4应当去掉，否则会报错）
在Molclus每次调用Gaussian做完优化和振动分析任务后，当程序发现当前目录下有template_SP.inp文件时，就会自动调用ORCA做单点计算，并将所得的高精度单点能自动加到振动分析得到的自由能热校正量上作为此结构的最终能量，这也正对应于溶剂下的自由能（严格来说还需要加上1.89 kcal/mol标准态转换能量后才是真正意义的溶剂下的自由能，详见http://sobereva.com/327。但是否考虑这个不影响构象间相对能量）。
现在执行./molclus |tee out.txt，这代表运行molclus的时候也把输出信息写入当前目录下的out.txt里便于之后检查。
计算过程中，调用Gaussian计算时，正在运行的Gaussian任务的输出文件名是gau.out，故可以通过tail -f gau.out命令监控计算进程，即让gau.out里新输出的信息不断显示在屏幕上。当进入调用ORCA算单点能的时候，也可以通过tail -f orcaSP.out监控计算进程。在笔者的Intel 36核机子下，对每个结构Gaussian做opt freq任务耗时都是一个小时多一点，ORCA算每个高精度单点任务也是耗时一个小时多一点，因此每个结构共花两个多小时，6个结构一共花了15个小时。
每个结构振动分析完成后Molclus会在屏幕上提示有没有虚频，在isomers.xyz文件里每个结构的Nimag =字样后面也可以看到虚频数目。当前的计算得到的每个结构都没有虚频，所以最终结果是有意义的。如果有的有虚频，应当加上有助于消除虚频的关键词重新计算，见《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（http://sobereva.com/278）。
算完后，还是用isostat处理得到的isomers.xyz，输出信息如下。可见这次有两个结构经过优化后收敛到了相同结构，故最终剩下5个。如前所述，当前显示的能量对应的直接就是溶剂环境下的自由能。
Energy of isomer 1 : -2320.81823800 a.u., new cluster
Energy of isomer 2 : -2320.81480600 a.u., new cluster
Energy of isomer 3 : -2320.80966500 a.u., new cluster
Energy of isomer 4 : -2320.81852000 a.u., new cluster
Energy of isomer 5 : -2320.81485300 a.u., attributed to cluster 2
Energy of isomer 6 : -2320.80904200 a.u., new cluster
The lowest energy is -2320.81852000 a.u.
Sorting clusters according to energy...
# 1 Count: 1 E= -2320.818520 a.u. DGmin= 1.25 DE= 0.00 kcal/mol
# 2 Count: 1 E= -2320.818238 a.u. DGmin= 1.01 DE= 0.18 kcal/mol
# 3 Count: 2 E= -2320.814853 a.u. DGmin= 1.01 DE= 2.30 kcal/mol
# 4 Count: 1 E= -2320.809665 a.u. DGmin= 1.10 DE= 5.56 kcal/mol
# 5 Count: 1 E= -2320.809042 a.u. DGmin= 1.10 DE= 5.95 kcal/mol
可见最低的两个构象能量相当接近，都不可忽略。虽然当前用的计算自由能的做法已经很好了，但由于SMD算中性体系的溶解自由能通常都有零点几kcal/mol误差（这还只是对于小体系而言），而且即便RI-PWPB95-D3(BJ)/def2-QZVPP给出的构象相对能量差也并不能保证总能精确到零点几kcal/mol，所以当前我们应当认为我们找出来的能量最低的两个构象都算是最稳定构象，在实际水环境中出现的几率都会非常大。
接着在isostat里输入y，再输入298.15，就会得到常温下的Boltzmann分布比例，如下所示，原理详见《根据Boltzmann分布计算分子不同构象所占比例》（http://sobereva.com/165）。
# 1 Count: 1 Ratio: 56.74%
# 2 Count: 1 Ratio: 42.09%
# 3 Count: 2 Ratio: 1.17%
# 4 Count: 1 Ratio: 0.00%
# 5 Count: 1 Ratio: 0.00%
可见，常温下前两个构象共同占据绝对主导，而其它构象的出现概率就都可以忽略了。
PS：如果你要考察其它温度下的Boltzmann分布，不仅isostat里要输入相应的温度，在template.gjf里也应当加上“temperature=[温度]”关键词来让Gaussian输出的自由能热校正量对应相应的温度；或者settings.ini里用ibkchk=1，使得每个opt freq任务的chk都被保留下来，便于之后自行用freqchk工具得到其它温度下的自由能热校正量。单点能从Molclus的ibkout选项备份出来的ORCA输出文件里直接就可以读到。
下面是自由能最低的那两个构象的结构图，即对应于isostat给出的cluster.xyz的前两帧：
(, 下载次数 Times of downloads: 517)
能量最低构象的旋转动画如下，分子的各个部位都能看清楚
8 弱相互作用分析
这里顺便对上面得到的两种能量最低构象考察一下分子内相互作用。Multiwfn支持的弱相互作用分析方法极多，见《Multiwfn支持的弱相互作用的分析方法概览》（http://sobereva.com/252）。为了便于构象搜索后用Multiwfn做这些分析，个人建议Molclus调用Gaussian做opt freq任务时ibkchk总是设为1从而自动保留chk文件。
这里就用《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）、《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（http://sobereva.com/667）中介绍的已经很流行的NCI方法考察一下能量最低的两个构象的分子内弱相互作用，NCI图如下所示
(, 下载次数 Times of downloads: 483)
由上图可见，这两个构象主要差异在于碱基和苯基之间的相对位置以及分叉的烃链的位置。2的pi-pi堆积更强一些，因为苯环与碱基错开的程度更低，而1中分叉烃链可以与碱基与苯基同时发生范德华吸引作用，而且还有2所没有的O-H...N内氢键。因此二者都有各自的对自己构象有利的分子内相互作用。如果你去看其它构象的话，就会发现分子内相互作用明显不如1和2这么充分和有利，因此能量明显更高。
9 附加说明
本文介绍的构象搜索方法普适性极强，几乎什么体系都可以用，各类分子、分子团簇、原子团簇皆可。只不过要注意GFN-xTB系列方法目前只能用于<=86号元素，因此碰上锕系配合物等没法用。并且原始的6-31G*只支持前四周期元素，如果涉及更靠后的元素，最简单的办法是把基组改为def2-SVP，耗时只高出一点，对前六周期元素都有定义。
本文的构象搜索流程是否一定能找到能量最低的构象很大程度上取决于体系特征、xtb做动力学跑的时间长短和温度，以及GFN2-xTB优化后保留最低的多少个构象用于之后的DFT计算。显然体系可旋转的键越少、动力学跑的越长、保留的构象越多，则最终能找到能量最低构象的成功几率就越高。动力学用较高温度也有助于翻越较高势垒，从而更全面地采样，但温度过高可能会导致动力学不稳定、出现不该有的异构化结构、在能量较低势能面区域采样不充分等问题。对于瑞德西韦这个体系我们所找出的两个能量最低构象，从几何结构和分子内弱相互作用模式来看都极其合理，从直觉上也想象不出还有什么其它构象能有更低的能量，因此可以断定我们成功找到了真实的能量最低构象。
几何优化和振动分析调用ORCA来做也不是不行，特别是如果你没有Gaussian的话这是首选。但ORCA的优化算法不如Gaussian稳健，而且振动分析的效率低，所以虽然ORCA利用RI技术可以令每一步优化耗时较低，但总的来说除非体系特别大（如一二百个原子），否则仍然更建议用Gaussian。ORCA支持的B97-3c方法做优化又好又快（大概率好于B3LYP-D3(BJ)/6-31G*，而耗时更低），但是对于笔者撰文时的最新的4.2.1版，B97-3c算的Hessian似乎有点bug，容易出现小虚频，虽然用numfreq数值频率可以解决，但此时耗时太高，故应当慎用。
有人可能会觉得6-31G*对于弱相互作用明显影响分子构象的柔性大体系的优化可能太小了。实际上用这个就已经够了。为体现这一点，笔者对前面得到的能量最低构象在6-311G**下做了进一步优化，原先的6-31G*的结构和新得到的6-311G**的结构经过叠合后的图像如下所示，红色是6-31G*的，蓝色是6-311G**的，可见没什么明显区别。内氢键键长、pi-pi堆积的平面间距也基本没变，而用6-311G**优化此体系耗时比用6-31G*高4、5倍。不过，如果当前体系只有几十原子，你的计算资源又很富裕，为了避免跟某些较真的审稿人argue，也可以用6-311G**（而再继续提升基组就绝对没有丝毫意义了）。
(, 下载次数 Times of downloads: 512)
本例计算时在DFT做优化过程中没有考虑溶剂模型，不用担心这会造成明显问题。对前面找出来的能量最低结构，笔者之后又在IEFPCM模型表现的水环境下在相同级别做了优化。下图是气相下（红线）和水模型下（蓝线）优化后的叠合图，可见结构差异基本可以忽略。
(, 下载次数 Times of downloads: 555)
不过，如果希望更严谨一些，避免审稿人对计算流程吐槽什么的，并且能接受更高一些的耗时（高几分之一），Gaussian做DFT优化的步骤也可以用Gaussian默认的IEFPCM隐式溶剂模型。注意此时不推荐用SMD模型，因为对于优化和振动分析这种目的不是算能量的任务它不会比IEFPCM更好，而且有时会出现收敛困难以及造成小虚频。如果是用ORCA来做优化和振动分析，更不建议这个阶段用溶剂模型，因为ORCA的溶剂模型的代码不如Gaussian成熟稳健，哪怕不用SMD而用CPCM，ORCA对柔性大体系优化出来的结构也容易有小虚频。而且ORCA里对大体系开溶剂模型导致耗时的增加明显高于Gaussian。
在GFN2-xTB优化后用DFT再做优化和振动分析是绝对必要的，不要图省事直接就用GFN2-xTB优化的结构算高精度单点。对于瑞德西韦这种柔性大分子，若基于GFN2-xTB的结构再做DFT优化，通常会看到结构还有显著的变化，即第一帧和最后一帧肉眼都能看出不同。为了说明这一点，对前面我们得到的能量最低结构（蓝线）我用GFN2-xTB在真空下又做了一下优化（红线），经过叠合后的图如下所示，可见从同一个初猜结构开始，GFN2-xTB的极小点结构和DFT极小点结构相差非常明显，这样的差异是绝对不能接受的（实际上这个构象还算是相差相对较小的，本文研究过程中还发现DFT和GFN2-xTB优化出来的某些构象里烷烃链那部分位置差得老远）。
(, 下载次数 Times of downloads: 528)
笔者在《谈谈谐振频率校正因子》（http://sobereva.com/221）中介绍过频率校正因子，原理上来说计算自由能热校正量时考虑校正因子有益，但从对实际构象能量差的影响程度来看考虑不考虑都无所谓，所以本文就没考虑。要考虑的话可以在opt freq任务对应的template.gjf里加上scale=关键词，就用ZPE校正因子就可以。
ORCA算高精度单点用的PWPB95-D3(BJ)/def2-QZVPP对于比当前体系显著更大的体系可能会算得很吃力，届时可以改为def2-TZVPP。但注意QZ级别才能充分发挥双杂化泛函的潜力。对于PWPB95-D3(BJ)/def2-TZVPP都很难算得动的体系，可以改为普通泛函ωB97M-V结合def2-TZVPP，参见此文的说明《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490）。
本文的算溶剂下的自由能的做法是懒人的做法，即“真空下高精度单点能”和“溶解自由能”这两部分的计算直接通过一个任务“溶剂下高精度单点能”完成，这样做原理上没有用M05-2X/6-31G*专门结合SMD模型去算溶解自由能那么理想，详见《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）。不过这也不是什么明显问题，更何况构象能量求差的时候误差还能抵消很多。
本文的构象搜索的全套流程完全可以写成shell脚本自动完成，因为涉及的所有操作环节都能简单地通过命令行完成，Molclus也是通过命令行运行的。感兴趣的读者请自行完成这样的脚本，并不复杂。
有本文读者问及用crest做metadynamics如何，前文里已经说了，有这个没有什么实际好处。笔者也实际试了一下，非常令人失望：
crest结合GFN2-xTB真空下搜索，两个半小时结束，给出的能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性
crest结合GFN2-xTB在GBSA表现的水溶剂下搜索，两个小时结束，给出的能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性
crest结合GFN0-xTB真空下做搜索，45分钟结束，所有给出的516个结构没有任何一个结构和真实的能量最低结构有相似性
可见对于当前体系大失败。crest的原理看似sophisticated，实际表现还没高温MD效果好。
技巧：利用crest程序节约xtb批量优化大量结构阶段的耗时
在前面文章中我们看到，Molclus结合xtb可以快速对xyz文件里记录的大批量结构进行快速优化。实际上这个阶段还可以明显节约时间，也就是利用crest程序。crest可执行程序在xtb的github页面上就能下载到。Molclus调用xtb时是对traj.xyz文件里的结构一个接一个来优化的，xtb在执行时是用OMP_NUM_THREADS环境变量指定的线程数做并行计算。而crest调用xtb来优化时，会同时启动OMP_NUM_THREADS个xtb程序同时对输入的xyz里的结构进行批量优化，每个xtb都只用单线程。crest这种方式调用xtb做批量优化的总耗时明显更低。
具体来说，运行以下命令就可以在GFN0-xTB级别下用normal几何优化收敛限（这是xtb默认的收敛限）对traj.xyz里所有结构进行优化：
/sob/crest -mdopt traj.xyz -gfn0 -opt normal -niceprint
其中/sob/crest是crest可执行文件路径，-niceprint代表计算过程中显示进度条。还可以用-chrg指定体系净电荷数，-uhf指定自旋，-g [溶剂名]指定溶剂模型。
计算完毕后在当前目录下得到了crest_ensemble.xyz文件，记录了优化后每一帧的结构和能量。从Molclus 1.9.3版开始，自带的isostat工具也能处理这个文件，当发现文件名里面有crest字样的时候就会自动以crest_ensemble.xyz文件的格式来读取内容，然后进行去重、排序并得到cluster.xyz。
用此做法对本文例子里高温MD跑完后的2000个结构用此方法优化，在笔者的Intel 36核机子上只花了11分钟，比起用Molclus调用xtb优化的时间节约了一个数量级。类似地，之后用GFN2-xTB在水模型下批量优化的那个阶段也可以借用crest来显著降低耗时。
作者Author: 八月的雨季 时间: 2020-2-8 16:21
社长动作好快
作者Author: 八月的雨季 时间: 2020-2-8 17:34
如果同样想得到能量最低构象，可否直接用主流的MD软件（Amber/Gromacs）在300k，TIP3P水模型下面，跑一段长时间的轨迹，同样保证采样充分，用主成分分析的方法找到代表性的构象，再取这几个结构用DFT的办法去优化？
作者Author: wuzhiyi 时间: 2020-2-8 18:09
想请教一下，用xtb跑动力学，然后xtb优化，再用orca/Gaussian算。
和直接用xtb的crest得到能量低的结构，然后orca/Gaussian算。
有啥优劣呢？
作者Author: snljty 时间: 2020-2-8 19:05
请问老师，您的贴子通过动力学搜索构象的时候经常都是跑一个高温的动力学，这和周期性退火去各个温度最低的构象各有什么优劣呢？后者更接近极小值构象，每一步优化可能快一些，但是后者得到相同数量的构象需要的总动力学步骤多得多，还有别的差别么？谢谢！
作者Author: IIIO_OIII 时间: 2020-2-9 05:23
仔仔细细看完了，真的是很详细的Protocol！
作者Author: sobereva 时间: 2020-2-9 06:13
可以。没必要用PCA分析，可以做个簇分析，把代表性结构拿出来用于之后优化
作者Author: sobereva 时间: 2020-2-9 06:16
首先，目前版本的xtb不支持退火。
用周期性退火产生大批0 K结果作为初始构象原理上是是很恰当的，但对于本文的构象搜索的protocol来说，把高温动力学改成周期性退火，在我来看从实际效率上和结果可靠性上都没什么特别的好处。
周期性退火能有意义的情况是动力学过程很便宜（诸如分子力场），而之后优化较贵（如直接用DFT），此时靠周期性退火得到更有意义的较少量的初始构象的做法显得才有价值。而GFN0-xTB跑动力学并非像力场那么便宜，而之后经过GFN0-xTB和GFN2-xTB的两步很快速过滤，留下来的少量构象也已经很有意义了，所以周期性退火相对而言没什么好处。而且，相对于周期性退火，本文的做法有更大概率把能量最低的一批结构都获得，因为初始构象采样得更全面。
作者Author: sobereva 时间: 2020-2-9 06:25
这和8L我回复的高温MD相对于周期性退火的做法有何优劣有相似之处，都是 动力学直接得到大批结构 vs. 更复杂的做法得到较少更有意义的结构 之间的权衡。
作者Author: snljty 时间: 2020-2-9 09:12
谢谢老师！
作者Author: exity 时间: 2020-2-9 15:29
XTB果然万能，社长有机会能出一个crest的帖子吗？
作者Author: IIIO_OIII 时间: 2020-2-9 15:58
从6.2版本开始，xtb就不支持周期性退火了，被crest取代了。
作者Author: sobereva 时间: 2020-2-9 16:53
文中说了，用crest没意义，所以不可能写帖子。本文的做法就已经完美了。
作者Author: IIIO_OIII 时间: 2020-2-10 04:07
汇报一个数据点。
我用文中动力学的条件跑了一遍自己的多肽离子（四个氨基酸，44个原子），有一类相对稳定的构象没有被搜索到。把温度升到500K，就被搜索到了。个人认为是因为离子中有相对较强的氢键，妨碍了构象搜索的完整度。希望自己这么理解是对的。
作者Author: sobereva 时间: 2020-2-10 06:57
原理上，400K下跑更长时间，也会增加最终找到这个构象的概率。用500K加速这个过程也是合理的。
碰到体系内有比较难以翻越势垒、阻碍构象变化的情况，都建议进一步增加温度
作者Author: pyscf 时间: 2020-2-10 08:38
你也不想想 社长花这么多时间写这个是为了推销Grimme的crest程序吗？ 当然不是
你居然还在这边问crest的事情。。。
虽然crest真的很香 亲测。。。
作者Author: pyscf 时间: 2020-2-10 08:42
本帖最后由 pyscf 于 2020-2-10 08:46 编辑
需要注意的是 你的protocol第一步就用了Grimme的xtb 殊不知crest就是原生基于xtb的
而且所有步骤一步到位（Grimme弄成了全自动的，可能“更完美”）
另外 由于其metadynamics的CV是rmsd 因此就不存在（克服了）一定需要升温至x00度才能“有机会”access到一个构象的问题
作者Author: sobereva 时间: 2020-2-11 12:09
文中已然明确写了
名为crest的利用xtb跑metadynamics的构象搜索工具
没有什么情况是“一定升温到xxx”才能翻越的，这和模拟时间长度直接相关，想加速采样效率提升到更高温度就完了。为什么没必要用crest在帖子和回复里我都已经说了。
作者Author: sobereva 时间: 2020-2-11 17:03
更新了本文，主要是在文末加入了借用crest大幅降低xtb批量优化耗时的说明，还增加了DFT与GFN-xTB优化结构的对比等内容。
作者Author: wuzhiyi 时间: 2020-2-11 17:21
对小分子来说，确定质子化状态还是一个很蛋疼的问题，crest这方面做的还挺方便的感觉。
作者Author: sobereva 时间: 2020-2-11 20:12
这里说一下crest的构象搜索的结果（输入结构和本文中高温MD初始结构相同）：
crest结合GFN2-xTB真空下搜索，两个半小时结束，给出的能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性
crest结合GFN2-xTB在GBSA表现的水溶剂下搜索，两个小时结束，给出的能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性
crest结合GFN0-xTB真空下做搜索，45分钟结束，所有给出的516个结构没有任何一个结构和真实的能量最低结构有相似性
可见对于当前体系大失败。
作者Author: elvisng 时间: 2020-2-11 22:08
可否分析一下是什么可能的原因让 crest 的结果这么差吗？
作者Author: szzjd 时间: 2020-2-12 09:50
运行第二部优化出现如下错误：
*** Configuration 1 ***
Current date: 2020-02-11 Time: 20:40:32
Loading geometry 1 from the inputted geometry file
Generating xtb.xyz file...
Running xtb: xtb xtb.xyz --opt --gfn 0 --chrg 0 --uhf 0 > xtb.out
normal termination of xtb
Extracting energy and geometry from xtbopt.xyz
forrtl: severe (59): list-directed I/O syntax error, unit 10, file /home/gauss/molclus_1.9.2_Linux/xtbopt.xyz
Image PC Routine Line Source
molclus 0000000000419DDB Unknown Unknown Unknown
molclus 0000000000446D0B Unknown Unknown Unknown
molclus 000000000040A36B Unknown Unknown Unknown
molclus 0000000000404422 Unknown Unknown Unknown
libc-2.17.so 00007F8EDB51E505 __libc_start_main Unknown Unknown
molclus 0000000000404329 Unknown Unknown Unknown
作者Author: snljty 时间: 2020-2-12 09:54
本帖最后由 snljty 于 2020-2-12 09:56 编辑
十有八九是你的xtb版本太老，Molclus版本"太新"。
Molclus的开发日志有：
2020-Jan-19：1.9.1版发布。
• 兼容了2019年12月18日发布的xtb 6.3版预览版，对老版本xtb不再兼容。
到git-hub下载最新版xtb。
作者Author: szzjd 时间: 2020-2-12 10:11
是的，问题解决，谢谢！
作者Author: sobereva 时间: 2020-2-12 13:52
从原理以及实际结果的分析，有两个可能
(1)crest默认用的iMTD-GC方法在其流程的设计上，目的就是去找能量最低构象（也包括利用MD在能量最低构象附近进行采样），因此采样并不全面、充分。然而，GFN-xTB系列方法中哪怕精度最高的GFN2-xTB，其精度还是很有限，其势能面上能量最低构象并非是真实的能量最低构象（而且可能结构相差很大），这就导致crest搜索过程采样区域高度侧重于GFN-xTB势能面的能量最低区域，这反倒往往把真正能量最低结构区域漏过去了（这点在GFN0-xTB上尤为显著，给出来的能量最低构象凭化学直觉就知道其实是高能构象才对）。
(2)crest搜索过程默认的参数对于当前体系可能还不合理，但这个因素相对于(1)来说较为次要。
另外，crest的搜索算法有明显的历史依赖性，这导致最终结果可能对初始结构的选取敏感。
相对来说，crest的搜索方式，对于那种势垒很高（直接做高温MD都不好翻越的），而且GFN-xTB的精度又能正确区分不同构象/构型能量高低次序的时候有用。而对于瑞德西韦这种靠高温MD就足矣充分采样、遍历各种构象区域，而且GFN-xTB的精度又显得明显不足的柔性大体系，用crest不仅没意义，反倒结果可能严重误导人。
作者Author: IIIO_OIII 时间: 2020-2-12 15:12
前两天我遇到了同样的问题。谢谢你汇报这个错误。
作者Author: IIIO_OIII 时间: 2020-2-12 15:13
厉害！
作者Author: elvisng 时间: 2020-2-12 21:06
本帖最后由 elvisng 于 2020-2-12 21:47 编辑
我尝试用 crest 做批次优化，但几秒就完成，出来的 crest_ensemble.xyz 是空白的。请问有人遇到这问题吗？(换成 gfn1 就没问题, 只是 gfn0 有问题)
以下是輸出內容:
> crest -mdopt traj.xyz -gfn0 -opt normal -niceprint
-gfn0 : Use of GFN0-xTB requested.
-opt 0
=============================
# threads = 36
=============================
------------------------------
GFNn-xTB Geometry Optimization
------------------------------
Geometry successfully optimized.
=======================================
| M D O P T |
=======================================
Optimization along trajectory (or ensemble).
Input file: <traj.xyz>
Containing 2000 structures.
writing TMPCONF* Dirs from file "traj.xyz" ... done.
Starting optimization of generated structures
2000 jobs to do.
[##################################################] 100.00% finished.
done.
Now appending opt.xyz file with new structures
Optimized ensemble on file <crest_ensemble.xyz>
-----------------
Wall Time Summary
-----------------
MDOPT wall time : 0h : 0m : 3s
--------------------
Overall wall time : 0h : 0m : 3s
CREST terminated normally.
作者Author: sobereva 时间: 2020-2-13 17:36
没遇见过
先看看xtb直接优化此体系能否正常完成
另外，mdopt运行中途会在当前目录下产生临时目录，里面有自动产生的xtb输入输出信息，看看里面内容
作者Author: elvisng 时间: 2020-2-13 22:12
谢谢你的回复。
用 xtb 优化此体系，即使用 gfn0 也没有问题。但是用 crest ，输出的 crest_ensemble 就是空白的。至于产生的临时档案, 也只找到 .xtboptok (0 bytes), wbo, coord 这3个临时档案, 但都没有详细输出数据。不知道是不是系统的 bug。
附上了批处理的档案（跑分子动力得到的），是一个小分子，只有200 frames 的 trajectory。
另外有一点，我跑 md 的时候，会说找不到 .param_gfn0.xtb。我把这档案从 xtb 目录复制到当前目录下，这问题就解决了。不知道两件事有没有关系？
(, 下载次数 Times of downloads: 6)
作者Author: sobereva 时间: 2020-2-14 11:31
找不到 .param_gfn0.xtb是因为xtb没有正确安装（没有恰当设置环境变量），当前版本怎么装看
将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析
http://sobereva.com/421（http://bbs.keinsci.com/thread-10106-1-1.html）
crest结合gfn0在我CentOS 7.4上没问题，你的体系的输出
- [root@192 other]# ./crest -mdopt traj.xyz -gfn0
-
- ==============================================
- | |
- | C R E S T |
- | |
- | Conformer-Rotamer Ensemble Sampling Tool |
- | based on the GFN-xTB method |
- | P.Pracht, S.Grimme |
- | Universitaet Bonn, MCTC |
- ==============================================
- Version 2.8.1, Fri 20. Dec 13:44:46 CET 2019
- Using the xTB program.
- Compatible with XTB version 6.1 and later.
-
- Cite work conducted with this code as
- P. Pracht, F. Bohle, S. Grimme, in preparation.
- and
- S. Grimme, JCTC, 2019, 15, 2847-2862.
-
- OMP: Warning #181: OMP_STACKSIZE: ignored because KMP_STACKSIZE has been defined
- Command line input:
- > ./crest -mdopt traj.xyz -gfn0
- -gfn0 : Use of GFN0-xTB requested.
- =============================
- # threads = 1
- =============================
-
- ------------------------------
- GFNn-xTB Geometry Optimization
- ------------------------------
- Geometry successfully optimized.
-
- =======================================
- | M D O P T |
- =======================================
- Optimization along trajectory (or ensemble).
-
- Input file: <traj.xyz>
- Containing 200 structures.
- writing TMPCONF* Dirs from file "traj.xyz" ... done.
- Starting optimization of generated structures
- 200 jobs to do.
- 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
- done.
- Now appending opt.xyz file with new structures
-
- Optimized ensemble on file <crest_ensemble.xyz>
-
- -----------------
- Wall Time Summary
- -----------------
- MDOPT wall time : 0h : 2m :28s
- --------------------
- Overall wall time : 0h : 2m :28s
-
- CREST terminated normally.
作者Author: silly_cheng 时间: 2020-2-16 12:14
Sob老师您好，isostat需要进入程序之后才能设置"energy threshold"之类的参数吗？请问可以在打开程序的同时向程序传递参数吗 比如 "./isostat -e XX"
作者Author: sobereva 时间: 2020-2-16 13:01
是。不能
想通过命令行直接指定就用echo -e以管道方式运行就完了，参看Multiwfn手册5.2节
作者Author: 土拔鼠 时间: 2020-3-17 23:39
请问社长 xtb 与Gromacs 、ORCA在做MD的时候主要优势在哪里？
作者Author: sobereva 时间: 2020-3-19 22:25
相对于GROMACS，xtb不用考虑力场的问题，省事，但速度远比基于经典力场的GROMACS慢
xtb实现的GFN-xTB理论达不到ORCA里的DFT那么高精度，但快几个数量级
作者Author: flybird52888 时间: 2020-3-20 10:38
请问用xtb结合molclus算气相分子团簇可以只用GFN0-xTB吗？也就是掠过加溶剂化华那一步？
作者Author: tjuptz 时间: 2020-3-21 15:55
老师，我看GitHub上的releases只到了pre1版本啊？
作者Author: ck15620956553 时间: 2020-3-21 18:51
太牛了！
作者Author: sobereva 时间: 2020-3-22 15:15
可能是撤了。6.3 pre2是今年1月发布的
作者Author: sobereva 时间: 2020-3-22 15:18
跑MD的话用GFN0-xTB可以
作者Author: flybird52888 时间: 2020-3-23 09:53
那如果做批量优化就不能了吗?
作者Author: sobereva 时间: 2020-3-25 15:45
显然精度太低了
作者Author: 量子极限 时间: 2020-3-27 21:27
感谢老师的指导和回答，但是这样的话我也有个问题，就是用xtb在GFN0-xTB下做动力学模拟时如何选取适合的temp以及time呢？
比如先根据老师给的条件然后再结果分析后再进一步更改？请问各位大佬有没有更科学的做法
作者Author: sobereva 时间: 2020-3-28 05:14
温度一般就用文中的温度就可以，如果发现最后跑出来的轨迹由于热运动不够、没能翻越一些势垒跑到其它构象去，再增加温度
可旋转键越多的体系，原理上应该跑越长的时间，使得势能面采样充分
作者Author: mutron 时间: 2020-3-29 21:01
本帖最后由 mutron 于 2020-3-30 01:47 编辑
编辑掉。。
犯了低级理解错误
，需要“当发现文件名里面有crest字样的时候就会自动以crest_ensemble.xyz文件的格式来读取内容”，把含有crest的文件名写完整没问题了
作者Author: linqiaosong 时间: 2020-4-6 22:16
想问一下社长用xtb做半经验动力学的时候能不能像Gaussian跑admp一样设置初速度的方向？我看手册上面好像没有这个选项
作者Author: sobereva 时间: 2020-4-7 06:59
先跑一下，产生mdrestart文件，然后再篡改mdrestart文件内容，通过restart=true来读入初速度（PS：ORCA的AIMD目前也得用这种做法）
例如
(, 下载次数 Times of downloads: 104)
作者Author: Evanwill 时间: 2020-6-6 15:41
请问，第一步，用xtb在GFN0-xTB下做动力学模拟，可不可以同时或者依次对多个化合物的xyz文件进行批量动力学模拟，分别形成多个结果文件？
作者Author: sobereva 时间: 2020-6-6 18:43
如果有明确目的，可以
作者Author: scf 时间: 2020-6-13 05:23
表示构象差别的overlay diagram是用什么软件做的？谢谢
作者Author: sobereva 时间: 2020-6-13 22:07
在VMD中计算RMSD衡量两个结构间的差异以及叠合两个结构
http://sobereva.com/290（http://bbs.keinsci.com/thread-1080-1-1.html）
作者Author: scf 时间: 2020-6-14 04:10
thanks
作者Author: 欢乐多 时间: 2020-6-25 21:43
老师，您好，在用XTB对14元大环柔性体系做构象搜索时，初始由chemdraw画的结构，弄成初始结构，进行构象搜索，得到的构象优化后与单晶差别很大；400~1000K中尝试不同温度，还是没有搜到与单晶相似的构象。
因为有多个14元的大环柔性化合物，仅有其中的1个有单晶数据，因此，以有单晶数据结构的化合物为参照，检测XTB做构象搜索方法的可靠性，尝试多种条件，终未搜索到与单晶相似的构象。
请问老师，xtb有没有适合大环柔性体系构象搜索的条件？
作者Author: snljty 时间: 2020-6-25 21:48
上面整个流程计算的本身就不是凝聚相的结构，但是单晶是凝聚相。你看看有没有可能是这个原因。
作者Author: 欢乐多 时间: 2020-6-25 23:41
凝聚相是怎么回事？高斯可以体现凝聚相吗？
作者Author: 喵星大佬 时间: 2020-6-27 00:36
无论用Gaussian还是xtb等优化构型的时候都是游离的分子，即在真空/溶剂中进行的，而在晶体中，分子构象明显会受到相邻其他分子的相互作用，可能包括氢键，pi-pi堆积等等，所以晶体中的构象在游离态中未必是稳定的。搜索不到晶体的构象是很正常的。
作者Author: sobereva 时间: 2020-6-27 01:10
晶体环境对于这种高度柔性的分子构象影响是相当大的
真空下模拟不可能指望得到和晶体中完全一样的结构
作者Author: 欢乐多 时间: 2020-6-27 09:33
本帖最后由 欢乐多 于 2020-6-27 09:45 编辑
好的，老师，感谢回复！老师，如果用真空下进行团簇搜索，比如放进去2~4个分子，搜到的稳定构象与单晶团簇进行对比，是否有意义？
作者Author: 喵星大佬 时间: 2020-6-27 11:08
没有，单晶结构就是只有在无限重复的情况下成立的。只有几个分子连完整的晶胞都构建不出来。
一般来说只有像蛋白，MOF(按团簇处理)之类分子内相互作用很强的，远远强于分子间作用的物种，在溶液中的结构会和晶体基本一致。
作者Author: 欢乐多 时间: 2020-6-27 13:18
本帖最后由 欢乐多 于 2020-6-27 13:57 编辑
好的，老师，感谢回复！
那么我可以从晶胞中抠出来几个分子，计算分析分子间相互作用（Hirshfeld surface），还可以用Multiwfn去考察晶体的哪些特征能够使研究充实有意义？读文献不多，文献很少有拿晶体结构讨论的？
化合物好不容易出了单晶，怎么能充分让量化计算与实验单晶数据有机结合起来，充实一下实验内容。
作者Author: sobereva 时间: 2020-6-29 01:30
研究分子单体间的弱相互作用，参考
使用Multiwfn做基于分子力场的能量分解分析
http://sobereva.com/442（http://bbs.keinsci.com/thread-10907-1-1.html）
使用PSI4做对称匹配微扰理论(SAPT)能量分解计算
http://sobereva.com/526（http://bbs.keinsci.com/thread-15902-1-1.html）
还可以通过IGM方法可视化，与Hirshfeld surface分析有一定互补性
通过独立梯度模型(IGM)考察分子间弱相互作用
http://sobereva.com/407（http://bbs.keinsci.com/thread-9472-1-1.html）
还可以考察晶体中的孔洞
使用Multiwfn图形化展示分子动力学模拟体系中的孔洞、自由区域
http://sobereva.com/539（http://bbs.keinsci.com/thread-16491-1-1.html）
作者Author: tjuptz 时间: 2020-6-30 14:31
看了一下CREST的介绍，流程上虽然作者尽量自动化了，但感觉还是不如molclus简单粗暴。对于最后利用CREST的技巧，molclus里也可以实现一下吧，这样就可以省着再调用crest了
作者Author: 欢乐多 时间: 2020-6-30 21:43
感谢老师回复，准备充分利用Mutiwfn充实一下文章内容。
作者Author: sobereva 时间: 2020-6-30 21:46
这点我也考虑过，现在没时间弄
crest虽然精密，但实际效果上反倒经常不如来简单的方法好，而且整个流程可控性差，过于黑箱，不如用molclus放心
作者Author: ldatea 时间: 2020-7-2 14:14
sob老师，跑完动力学以后，出现以下情形是否会漏掉极小值点。
1.用gfn 0几何优化，出现“原本应该有的极小值点，在gfn 0或者gfn 2的势能面中不存在此极小值点”
2.两个极小值点的坐标接近，且能量相差不大时，把两个结构并入1个cluster的。
作者Author: sobereva 时间: 2020-7-3 02:11
1 可能会。但这不是大问题。重要的极小点几乎一定能被GFN-xTB所描述出来，要漏掉也是一些很浅、化学意义不强的极小点。
2 只要阈值不是特别宽就不用担心这个
作者Author: ldatea 时间: 2020-7-3 18:03
谢谢sob老师
作者Author: hzfish 时间: 2020-7-8 20:53
创建一个含有xtb的MD任务设置的文件md.inp，内容如下，双斜杠后面的是注释。
$md
temp= 400 //温度(K)
time= 100.0 //模拟的总时间(ps)
dump= 50.0 //每多少fs往轨迹文件里写入一次
step= 1.0 //步长(fs)
hmass=1 //氢原子的质量是实际的多少倍
shake=1 //将与氢有关的化学键距离都用SHAKE算法约束住
$end
我打算使用社长的方法搜索离子液阴阳离子对的2、4、6、8、10聚体的构象。
问题：
1、temp该如何选择？有没有判定标准？
2、time该如何选择？有没有判定标准？
谢谢！
作者Author: sobereva 时间: 2020-7-9 07:00
没有什么标准。400应该合理，先跑跑，不合适再说
模拟时间越长，构型空间采样越充分，最终得到能量最低构象的几率越高。设多大和体系，以及对结果靠谱程度的要求有关。100 ps一般够了，但对于n聚体，n越大时理应设得越大
作者Author: hzfish 时间: 2020-7-9 20:25
谢谢社长！
作者Author: yyyyss11 时间: 2020-7-12 22:31
老师，请教一个问题。我尝试在md.inp里面加入限制性的约束（例如下面的形式），对于MD过程好像并不起作用。请问有没有方法能在MD中实现这种键长键角等的限制性约束？谢谢老师。
$constrain
distance: 1,2,1.4
angle: 5,7,8,auto
dihedral: 3,4,1,7,180
$end
作者Author: sobereva 时间: 2020-7-13 10:25
没办法
作者Author: yyyyss11 时间: 2020-7-13 13:40
好的，谢谢老师
作者Author: mutron 时间: 2020-7-18 04:30
试了一下ORCA里面做MD可以固定住一个二面角，所以应该没问题可以去试一下，ORCA的MD有更丰富的选项。至少之前试过做一个小分子显式溶剂的MD可以把溶质分子完全固定住，用regions或者constraint（10^-2 pm for distances, and 10^-4 degree for angles and dihedral angles）都可以，xtb没有办法在做MD的时候exact fixing
作者Author: yyyyss11 时间: 2020-7-18 13:53
我去试试orca的MD功能，谢谢指导。
作者Author: qinzhong605 时间: 2020-7-19 13:56
社长，我拿您的文件包，xtb最新的6.3.2跑，第一步只有1042帧，不太对呢，不是应该有2000帧吗？对xtb不太熟悉，盼望指点一下，谢谢了。@sobereva
作者Author: qinzhong605 时间: 2020-7-19 16:51
知道怎么回事了，社长传的文件包里md.in应该为md.inp
，初次接触，没注意到这个问题。
作者Author: mutron 时间: 2020-7-27 00:01
1.老师，第二步用GFN0-xTB优化会导致漏掉一些能量更低的构象。如果跳过第二步直接第三步用GFN2-xTB在GBSA溶剂模型下优化会多出7个构象，能量相对于文中第三步得到的最低构象（-128.462874 a.u.）会再低0.32-1.91 kcal/mol。（我觉得第二步的必要性不是太强，能省出来的计算时间很有限，尤其是结合crest批量优化）
2.crest能够获得电子能量更低的一批构象。我看到老师21楼给出的crest结果“能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性”，是否可能只是crest得出的一批能量更低的构象没有一个和老师的结果有重复？crest得出的一批构象中最低的能量（-128.469685）比文中第三步得到的最低构象（-128.462874）低了4.26 kcal/mol
3.GFN2-xTB的能量精度还是不够，用来排序得保留相对能量比较宽的一批构象（见下），碰上如本例柔性不小的分子下一步DFT优化需要处理的结构太多算不动。所以中间得加上一步比较靠谱又便宜的能量计算，比如老师推荐的B97-3c，或者RI-PWPB95-D3(BJ)/def2-TZVP（-f），获取相对准确的能量排序再筛选
比如crest得到的能量最低构象，其GBSA下的GFN2-xTB能量：-128.469685，B97-3c结合SMD模型的能量：-2320.74009298；用B97-3c重新排序的能量最低构象分别为：-128.465825和-2320.750725。GFN2-xTB能量后者比前者高了2.42 kcal/mol，B97-3c能量则反过来低了6.67 kcal/mol（RI-PWPB95-D3(BJ)/def2-TZVP(-f)则是低了4.92 kcal/mol），误差达到7-9 kcal/mol。由crest通过后续步骤算出来能量最低的构象在GFN2-xTB排序中排第141（B97-3c排序排第2），相对默认的GFN2-xTB能量排序第1位GFN2-xTB能量高了3.78 kcal/mol，也就是说如果按默认的GFN2-xTB排序，后续DFT计算需要算200多个构象才能把最稳定构象包括进去
RI-PWPB95-D3(BJ)/def2-TZVP(-f)和RI-PWPB95-D3(BJ)/def2-QZVPP在本例相对能量相差大概零点几kcal/mol，其和B97-3c对于crest能量最低的17个构象是一样的，只是顺序有变化
4.isostat结构去重的时候，结构去重阈值(埃)选取0.5可能会有点大，导致相当一部分构象被合并了。比如crest经过B97-3c能量排序后最低的17个构象（计算过程用0.2(埃)），如果用0.5(埃)会导致其中4个构象消失，但观看结构会发现烃基扭动方式和取向明显不一样。又比如拿crest得出的所有构象和文中第三步能量最低的四个构象作对比，如果选取0.5(埃)，isostat会给出3个结构有相似，但仔细观看结构会发现差得不少。如果选取选取0.2甚至0.3(埃)则不会有相似结果
5.使用mdopt对MD或者crest结果做结构优化（normal opt或者vtight）、isostat去重后，会得到好几百甚至上千个相对能量在6kcal/mol以内的构象，哪怕用很便宜比较靠谱的方法算一次能量也耗时颇久。此外，第四步用Gaussian优化过后会出现不少相同的结构（比如文中第四步里面的第5个和第2个构象），DFT优化浪费的时间就更多了。如果在xtb优化完加上一步，使用gau_xtb计算力常数再优化一遍能去掉不少相同结构，减轻后面DFT计算的负担（我觉得通过计算力常数收敛到同一个构象，相对于提高阈值判定为相同构象会更稳妥）
6.几种方法最后得到的最低能量结构的自由能（Shermo考虑Grimme的准RRHO模型）：
-2320.82217947 0.00 kcal/mol N/A 本文结果
-2320.82816084 -3.75 2（a） 本文第三步后使用B97-3c排序
-2320.82719758 -3.15 10 MD_GFN0-xTB_500ps_500K
-2320.83092026 -5.48 10 MD_GFN2-xTB_GBSA_500ps_400K
-2320.82594269 -2.36 N/A crest_GFN2-xTB_GBSA排序
-2320.82994678 -4.87 20 crest_B97-3c排序
-2320.82986194 -4.82（MD） 20+22 crest能量最低10个结构继续用MD搜索
-2320.83013936 -4.99 20+25 crest能量最低10个结构继续用MTD搜索
（a）：GFN2-xTB+GBSA溶剂模型下得到构象、在RI-PWPB95-D3(BJ)/def2-TZVP(-f)结合SMD模型水平下获得的电子能量对比，比本文最稳定构象能量低的数量。反映最低能量区域附近的采样充分程度
从结果看来，因为老师在26楼谈到的原因，做完crest如果还需要对真正能量最低结构区域采样充分、找到真正最稳定构象，有必要再做一遍MD或者MTD，具体条件（构象附近搜索）需要摸索
7.老师可否考虑在molclus里面加上批量做MD的功能，需要做MD（MTD）的构象一多就开始变得麻烦了
8.molclus 1.9.3版本在计算itask=2、3且没有提供单点任务模板文件时会给出自由能加多一遍Gcorr的错误结果，但最新版1.9.5好像吉布斯自由能G总是=0了，不知道我哪里理解错还是程序问题
暂时找到的最稳定构象的NCI图：
(, 下载次数 Times of downloads: 124)
既有图1的内氢键，pi-pi堆积也和图2差不多
作者Author: sobereva 时间: 2020-7-27 04:50
前一阵子已经发布了molclus 1.9.6版
本文只是作为演示用，需要更严格的最小点把MD轨迹跑得更长、从而采样更充分即可，简单粗暴但很有效。虽然在细节上做各种改进可能能让当前体系搜索效率更高，但不一定有普适性。
作者Author: tjuptz 时间: 2020-8-28 15:46
请问老师
1 二进制的mdrestart用文本编辑器如notepad++的二进制插件 打开可以嘛？
2 mdrestart里不是已经记录了之前md的速度吗？
3 对于其他体系也是加到最后就可以了吗？还是跟输入坐标一致的？
作者Author: sobereva 时间: 2020-8-29 00:53
xtb的mdrestart是文本文件，直接用文本编辑器打开
是记录了，你可以篡改
里面记录的速度矢量数目得和原子数对应
作者Author: tjuptz 时间: 2020-8-29 07:48
果然是，之前想错了。谢谢老师
作者Author: naoki 时间: 2020-9-15 16:32
Sob老师，我想咨询一下xtb做动力学可以用显式溶剂吗，比如水和离子液体，如果一样跑的话，盒子里溶剂一多是不是计算耗时巨大无比
作者Author: sobereva 时间: 2020-9-16 09:07
虽然可以用，但一方面会显著增大耗时，另一方面xtb自身没法控压，也没法达到理想地模拟溶液状态的效果。非要用显式溶剂应当考虑gromacs跑
作者Author: naoki 时间: 2020-9-16 09:15
谢谢Sob老师，如果想考察溶液中的聚合反应是不是首选在lammps中用ReaxFF呢
作者Author: sobereva 时间: 2020-9-16 09:24
如果模拟尺度很大，是如此
尺度不是很大的话，可以考虑CP2K跑GFN-xTB的动力学，不过我没亲自试过
作者Author: naoki 时间: 2020-9-16 10:44
非常感谢Sob老师的指点~
作者Author: wzkchem5 时间: 2020-9-16 14:12
如果是传统自由基聚合或者某些特别快的阳离子聚合，那可以考虑
如果是可控活性自由基聚合、阴离子聚合、配位聚合，可能1毫秒都未必能聚一个单体，那样就不可能直接模拟了，可以用一条足够长的聚合物链，用metadynamics模拟聚一个单体的情况，然后推而广之到整个聚合过程
作者Author: naoki 时间: 2020-9-16 14:15
有道理，谢谢您的指点。我看到之前有文章用ReaxFF做了环氧开环聚合的加速模拟，不知道靠不靠谱
作者Author: 天道啊啊 时间: 2020-9-28 11:28
老师，关于大体系分子团簇的构象搜索，是像本文这样跑用XTB跑分子动力学合适，还是genmer+molclus合适？
作者Author: sobereva 时间: 2020-9-28 12:36
如果是柔性大分子构成的团簇，适合先用xtb跑动力学，然后再用molclus进一步算更准确的能量和筛选。genmer产生初始结构的话没法考虑大分子的柔性。
作者Author: 天道啊啊 时间: 2020-9-28 17:06
谢谢老师！
作者Author: jingetiema6112 时间: 2020-11-3 16:02
老师，如果我想得到一个化合物的合理的初始构象（CHON有机化合物，结构中有桥环也有螺环），用于输入到高斯软件并在一定基组水平上进行优化，进而进行频率分析，轨道分析以及生成热的计算，是不是可以用这个帖子的示例进行合理的初始构象的搜寻？
作者Author: wzkchem5 时间: 2020-11-3 22:14
对
只要不是电子结构特别奇怪的化合物，就可以用xtb做
作者Author: jingetiema6112 时间: 2020-11-4 12:26
好的 感谢解答
作者Author: CYQ@ 时间: 2020-11-11 10:27
社长，我想请教一下，我用xtb跑动力学然后优化的构相，在用gaussian优化和算单点能之前，能够通过vmd自己观察traj.xyz轨迹来判断筛选出想要的构相，然后在提交给Gaussian优化计算吗
作者Author: CYQ@ 时间: 2020-11-11 11:53
老师，我想请教一下，xtb跑MD可以用xtb 4e15z.xyz --input md.inp --omd --gfn 2 --gbsa h2o命令跑吗
作者Author: wzkchem5 时间: 2020-11-11 12:15
这样不严格，会引入人为偏差。最好还是用自动化的方法选择，哪怕比手动选择麻烦。
作者Author: jingetiema6112 时间: 2020-11-11 12:26
你好 请请教社长吧，我也是初学者
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |
