计算化学公社

标题: 计算RESP电荷时Gaussian输入文件的关键字IOp求助 [打印本页]

作者
Author: chemzhang 时间: 2022-3-2 15:01
标题: 计算RESP电荷时Gaussian输入文件的关键字IOp求助
本帖最后由 chemzhang 于 2022-3-2 15:10 编辑

（这个帖子半是求助半是经验交流）
我最近在根据一篇文献的报道使用Amber模拟计算氯化胆碱，准备输入文件时需要知道胆碱一价正离子Ch+的RESP电荷以下是参考文献中给出的数据。参考文献使用的是g03与R.E.D.计算的RESP电荷
partial charge

几何构型与原子名称

于是我用g16做几何优化、静电势的计算后利用Multiwfn经one-stage fitting求出RESP电荷,如下所示

比较发现只有羟基的计算结果跟原文献相比还比较吻合，而氮原子以及与氮直接相连的碳原子的RESP电荷几乎是文献报道的300%（N 0.05 vs 0.148; C(-N) -0.13 vs -0.32)，出现了严重的偏离。
而使用antechamber计算得出的RESP与原文献符合程度比Multiwfn更好(N 0.026 C(-N) -0.15 比原文献偏离50%）

我还用RESP.sh计算了RESP电荷，得到的结果与Multiwfn是一致的。
经过仔细检查后我发现Multiwfn与antechamber之间的差异是由于g16计算esp时使用的关键字导致。我在choline_esp.gjf中设置的关键字是
#p hf/6-31g(d) pop=MK IOp(6/33=2)，
而Multiwfn中计算RESP的关键字，根据RESP.sh，我猜测是
# B3LYP/6-31g(d) em=GD3BJ scrf(solvent=water) pop=MK IOp(6/33=2,6/42=6)
于是我逐个更改关键字，发现问题在于IOp。如果我将RESP.sh中关键字IOp部分的6/42=6去掉，只留下IOp(6/33=2)，那么RESP.sh得到的结果便与antechamber吻合，也与文献报道符合较好。
我虽然可以修改RESP.sh得到想要的效果，但是不知道怎么在Multiwfn中做出相应的调整。即便我预先在g16中用IOp（6/33-2）计算得到chk与fchk文件，以此作为Multiwfn的输入文件，Multiwfn给出的RESP电荷依旧是严重偏离的那组。
也就是说，怎么样可以直接在Multiwfn中，而不需要借助RESP.sh，获得理论水平为#p hf/6-31g(d) pop=MK IOp(6/33=2)的RESP电荷呢？

作者
Author: chemzhang 时间: 2022-3-2 15:06
本帖最后由 chemzhang 于 2022-3-2 15:07 编辑

此外我还想向各位老师请教一下IOp中6/42=6究竟有什么具体作用，为什么这个参数可以对一部分原子（N，C）的RESP电荷产生那么大的影响，而对其他部分原子（如羟基）几乎没有影响。我在amber官网中查到这么一段内容，但不是很理解

作者
Author: sobereva 时间: 2022-3-2 22:43
不要认为文献里就是合理的

拟合静电势电荷计算时候需要设定分布于范德华表面附近的拟合点，原则上来说，拟合点的密度越大，拟合结果越好。6/42=N设的是MK方式设置拟合点的时候每平方埃内的点为N个。N默认为1，明显偏小，即得到的拟合静电势电荷质量差。6/42=6刻意用来提高拟合点密度以提高拟合静电势电荷的数值层面的质量。要是加这个有害无益的话，下文的脚本里自然不会加上这个
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
http://sobereva.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）

直接基于波函数文件用Multiwfn计算RESP电荷的时候也可以直接手动设拟合点分布情况，如下所示，可见默认设置下等同于IOp(6/42=6)时候Gaussian产生的拟合点分布（注：虽然拟合点密度相同，但Multiwfn产生的拟合点的坐标和Gaussian略有差异，都是合理的，这导致的结果差异可忽略不计）

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

原理上，结果会随着拟合点密度的增加趋于收敛，默认设6是在精度和耗时上的较好权衡。继续增大对结果没有明显改进，而减小则可能令结果精度打不可忽视的折扣。

另外，Multiwfn手册4.7.8节专门介绍了怎么可视化拟合点分布以及拟合误差，感兴趣可以看看，对理解原理有好处。

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