关于考虑多构象权重的RESP电荷计算过程中的疑问

不想飞的猫头鹰

我的疑惑比较多，而且涉及的帖子比较多，部分地方可能描述不够准确，感谢老师们的包容和不吝赐教！一开始我是接触到使用sobtop为gromacs生成分子itp（sobtop首页以及二茂铁的帖子）的过程中了解到分配的原子电荷最好是采用resp电荷，然后就是在跟随sob老师的用multiwfn算resp拟合静电势电荷的原理的帖子以及扩展到resp2(0.5)电荷的帖子中了解到，如果想要生成用于gromacs动力学模拟的电荷，最好是考虑气相和液相以及多构象权重，所以综合起来的做法基本就是首先用molclus调用xtb做构象搜索并且做gfn0-xtb和gfn2-xtb(gbsa h2o)粗略优化后用isostat筛选排序，并调用gaussian在b3lyp（d3bj）/6-311g**（如果算得太慢就6-31g*）下做相比于半经验方法更好精度的优化和振动分析，然后调用orca或还是gaussian在更高精度的pwpb95(d3bj)/def2qzvpp（cpcm smd water）（算得太慢就wB97MV/def2tzvpp，gaussian就m062x(d3)/def2tzvpp）计算方法下算单点最后加上热校正量获得自由能从而获得玻耳兹曼分布权重，进而获得考虑多构象权重的resp电荷结果。我这里有一堆问题想请教：
1.既然用sobtop在基于hessian矩阵计算力常数的前提是做了几何优化与振动分析（基本没有对hessian的频率校正因子，所以不加scale关键词），这个地方和构象搜索中的优化与振动分析步骤貌似是重复的，是不是可以取构象搜索时产生的chk（格式转为fchk）文件从而之后就不必对取出来的结构重新做一遍优化加振动分析了？但是isostat对最后高精度计算的结果已经进行了排序，而生成的文件名都是gau000xx.out/chk、gauSP000xx.out/chk与orcaSP000xx.out/chk，那么怎么知道排序后的cluster.xyz的每一帧的分子对应的是这些out/chk文件呢？
2.我看帖子中用gaussian做几何优化时都没加隐式水模型，而是在算单点的时候才加了隐式水模型，但是resp2(0.5)电荷的帖子中提到“几何优化使用B3LYP-D3(BJ)搭配6-311G**（或更好一点点的def-TZVP），优化时在真空下即可，但如果是带电体系或局部带显著电荷的体系则应当在溶剂模型表现实际溶剂环境下进行优化。”也就是说，如果自己的目的是想要用于在gromacs显式水溶液中跑动力学，是不是最好在最初的构象搜索的用gaussian做几何优化和振动分析阶段就加上scrf（默认pcm不用注释，单点时改成smd）？这样的话在使用sobtop里的例子时是不是最好也是自己手动加上scrf？
3.既然计算resp电荷时考虑了多构象的权重问题，为何用sobtop生成振动常数不需要考虑多构象权重问题？是不是只需要用最大权重的构象来产生力常数，而只在分配电荷的时候才考虑多构象？如果是这样的话，像是在将Confab或Frog2与Molclus联用对有机体系做构象搜索一文中最后搜索出来的最大占比的两个构象形态差异巨大但分布权重基本相同（17.88%与16.95%）的情况下，如果只考虑一个构象的振动分析结果来生成力常数而忽略了其他几个构象，结果是否偏颇？而如果同时考虑多构象权重的力常数和电荷问题，电荷问题简单，multiwfn能很方便地分配电荷，但是力常数这么多，也是按照单纯的权重来加和吗？又该如何操作呢？
4.在跟随sob老师用molclus调用xtb或openbabel的confab做构象搜索的自主尝试中（参数都一致，除了orca的grid关键词去掉了,这样做的做法来自于《谈谈orca5.0的新特性和改变》中“默认是defgrid2，据称没有特殊情况这个档次的格点质量就足够了，而对积分格点要求特高的情况（近乎Gaussian里int=superfine的情况），比如VPT2做非谐振计算的时候，才需要defgrid3。”我安装的是orca5.0.3就没法用grid4 grid4x关键词），我的瑞德西韦分子的最大两个构象和sob老师的帖子基本一致，但是最低的两个构象的E（molclus校正好的自由能）分别为-2320.825044 a.u.和-2320.824164 a.u.，相比于原帖的-2320.818520 a.u.和-2320.818238 a.u.稍低了一丁点，但是最后以0.5kcal/mol排序并计算298.15K下的玻尔兹曼分布比例时却是68.08%和26.81%，与帖子中的56.74%和42.09%有比较明显的差异；而这个情况在使用confab构象搜索的过程中更为明显：我的最低的两个能量是-1469.137491和-1469.137120，sob老师是-1469.136706和-1469.136656基本差异不大，但是在用isosat生成cluster.xyz（跟随教程的默认0.25kcal/mol）后用vmd查看时，能量最低的两帧构象差距并不大，和sob老师帖子上说的“构象差异极大”有所出入,这是因为初始构象随机性导致的最终优化和计算结果有一定差异的原因吗？如果是这样的话，之后在进行构象搜索的时候怎么确定自己搜索出来的能量最低的几帧的构象到底是不是真正的最低构象和占比最大权重的构象呢？（下面图片是我的cluster.xyz的前两帧的actos分子的截图，可以看出比较明显的变化主要是甲基的位置的挪动，与原贴差异巨大）

5.按照sob老师 《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》开头的方法安装了预编译的xtb6.6.0,在编辑环境变量的过程中发现有export MKL_NUM_THREADS=N这一句，搜了搜发现MKL库是英特尔的数学函数库，但是我的机器是AMD的，而且在《Grimme的xtb程序的编译方法》手动源码编译安装的一文中有说需要需要ifort和icc，那么AMD的机器比如锐龙zen2,zen3,zen4产品以及霄龙系列能够使用英特尔的编译器吗？使用了的话相比于同级别产品会有负优化的问题吗（网上众说纷纭，而且测试的年代都太久远）？我搜到sob老师的B站上装vasp的视频（VASP 5.4.4极简安装方法（CentOS 7.6+ifort 19））中是在图形页面下装好的，那么有没有像sob老师的教程一样简洁通用的传统的纯文字远程终端安装教程可以分享呀？另外我并没有没装过什么数学函数库如blas、lapack和mkl库，为什么xtb能够照样运行，而不管它的话会对什么有影响？在此处有没有适用于AMD处理器环境下的环境变量可声明？
6.学习sobtop的初衷是解决很多时候gromacs没法用pdb2gmx的问题，在sobtop教程中，有些例子是不做振动分析而直接查表式地用gaff加uff生成拓扑参数，这样的拓扑会对分子动力学的结果产生多大影响？我之前所做的工作都只是在同一个力场下的模拟，这里我第一次接触到同时使用amber加gaff的参数，有些疑惑，以蛋白质配体模拟为例，假如蛋白质能用pdb2gmx直接生成拓扑，配体用以上的构想搜索加振动分析加resp2（0.5）电荷计算的流程生成拓扑，得到结果后在模拟时只需要把生成的小分子的itp中的atomtype提取出来跟蛋白质的放在一起并且最后include一下配体itp就能跑了，也就是说，sobtop生成的拓扑已经包含了全部所需要的gaff力场参数，与蛋白质的amber力场参数互不影响（就如同Justin.A的膜蛋白教程一样，lipid.itp中包含了与gromos53a6兼容的原子名、原子类型与成键参数的描述方式），而amber的原子名是大写，gaff是小写，力场是怎么认出来蛋白质与配体之间的相互作用的？（顺便想问一下，看多数人推荐是用amber14sb_parmbsc1，这时有过度生成螺旋的倾向，而精度不高的tip3p水模型正好相弥补，结果还不错，而amber19sb还没发布适用于gromacs的版本，听说用opc3水最好，而且看了水模型的横测，在其他力场里是否可以用opc3水来代替spc/e水？）


提前感谢热心网友们的答疑解惑，这些问题我郁闷了很久，也希望能给同样在疑惑的网友提供参考，感谢大家的指教呀

sobereva

1 将cluster.xyz里的能量和备份的输出文件里的能量去匹配

2 只要不是局部带显著电荷的体系，opt freq加不加隐式溶剂模型无所谓。懒得判断、不在乎多花点耗时就始终都加

3 需要基于Hessian计算力常数的场合都是刚性或至少局部是刚性的体系，这种情况本来就没多构象的事（至少对刚性局部而言）

4 这不算显著的结构差异。对于那么大的体系局部甲基朝向相差一点对自由能的影响微乎其微，目前再准确的计算级别也基本不可能可靠地算出来这点差异带来的影响。骨架结构有大幅差异那才是值得说事的。

5 AMD机子上用Intel的东西的时候别开专供Intel CPU的优化选项就行了。静态库是程序编译的时候就连接并且纳入到可执行文件里的，不需要自己的机子上装相应的库。

6 这是GROMACS常识知识。基于AMBER力场的[defaults]字段的设置，不同原子类型之间LJ参数会自动基于原子类型自身的LJ参数通过混合规则产生，包括AMBER-AMBER、AMBER-GAFF、GAFF-GAFF原子类型之间的。

不要以为一个水模型对于纯水的特征模拟有长处就认为代替蛋白质力场原配水模型会更好，本身蛋白质力场在拟合参数就和原配的水模型会有误差抵消，换了别的水模型就没这种抵消效应，到时候审稿人还可能问你怎么证明你用的水模型和那个蛋白质力场的兼容性。

不想飞的猫头鹰

sobereva 发表于 2023-3-23 07:24
1 将cluster.xyz里的能量和备份的输出文件里的能量去匹配

2 只要不是局部带显著电荷的体系，opt freq加 ...

十分感谢卢老师，受启发很大，感谢您的耐心回答想要再问下第四点，也就是说，构象搜索出来的结果有比较大的差异才值得专门详细谈，而由于初始结构的随机性，最后高精度优化的结果有所不同是可以理解的吧，在这种情况下如果想要尽可能有说服力就需要选取更多帧的结构进行最后的计算或者做多次模拟可以吗？

sobereva

不想飞的猫头鹰 发表于 2023-3-23 14:19
十分感谢卢老师，受启发很大，感谢您的耐心回答想要再问下第四点，也就是说，构象搜索出来的结果有比 ...

如果想尽可能增加最终得到能量最低结构的概率，就增加初始traj.xyz的帧数（轨迹跑得更长来实现）、初筛后也保留较多的帧。并不需要多次模拟

不想飞的猫头鹰

sobereva 发表于 2023-3-24 00:26
如果想尽可能增加最终得到能量最低结构的概率，就增加初始traj.xyz的帧数（轨迹跑得更长来实现）、初筛后 ...

嗯嗯，十分感谢！