关于使用Mutiwfn计算RESP的若干问题

YZH_jesse

我阅读完Sobereva老师《R原子电荷的思想以及在Multiwfn中的计算》(http://bbs.keinsci.com/thread-16190-1-1.html) 后尝试了一下甲醛分子RESP2的计算，以下为我的计算步骤及计算过程中出现的几个主要问题：

第一步：在gview中随便画一个甲醛分子的初始结构，保存为 .gjf 格式的Gaussian输入文件，修改其内容为：

1. %nprocshared=6
   
2. %mem=2GB
   
3. %chk=C:\Users\yxdhx\Desktop\Charge Distribution\Formaldehyde\opt.chk
   
4. # B3LYP/TZVP em=GD3BJ opt
   
5. 
   
6. niconiconi
   
7. 
   
8. 0 1
   
9. C                  0.52701500   -0.00001300    0.00042100
   
10. H                  1.11571800   -0.93807900   -0.00084900
    
11. H                  1.11357700    0.93925600   -0.00084800
    
12. O                 -0.67392300   -0.00013700   -0.00010400
    
13. 
    
14. --link1--
    
15. %oldchk=C:\Users\yxdhx\Desktop\Charge Distribution\Formaldehyde\opt.chk
    
16. %chk=C:\Users\yxdhx\Desktop\Charge Distribution\Formaldehyde\SP_gas.chk
    
17. # B3LYP/def2TZVP em=GD3BJ geom=allcheck
    
18. 
    
19. 
    
20. --link1--
    
21. %oldchk=C:\Users\yxdhx\Desktop\Charge Distribution\Formaldehyde\opt.chk
    
22. %chk=C:\Users\yxdhx\Desktop\Charge Distribution\Formaldehyde\SP_solv.chk
    
23. # B3LYP/def2TZVP em=GD3BJ scrf=solvent=ethanol geom=allcheck

复制代码

第二步：将Gaussian计算得到的两个 .chk 文件通过 formchk.exe 程序转为 .fch 文件，输入到Multiwfn中，依次选择 7 / 18 / 1 来计算RESP，得到的结果与Sobereva老师博文中所得到结果的对比如下：
（PS：其中 Test1 为我画的第一个初始结构的计算结果，Test2为我画出的第二个初始结构的计算结果，Test sob为我用当前Multiwfn版本并使用http://bbs.keinsci.com/thread-16190-1-1.html中给出的甲醛初始结构算得的结果，最后一个为sob老师在网站上直接给出的计算结果，意在对比初始结构对于计算结果的影响）


问题1：初始结构对于计算结果还是会产生影响的（小数点后第三位），而且由于Lammps中给出的TIP3P水模型中 O 和 H 原子的电荷均精确到小数点后第三位（-0.830/0.415），所以我觉得这个误差还是有必要去深究一下；而且很有趣的是我随意画出的两个初始结构的计算结果 和sob老师给出的初始结构的计算结果有明显差异，我不理解这个差异产生的原因是什么（三者除了初始原子坐标不同 其他条件保持严格相同，而且从 .gjf 第一部分可以看出是先对甲醛分子进行过结构优化的）

问题2：依次选择 7 / 18 / 1 计算RESP后，均提示了 Stage 2 of standard RESP fitting is skipped since no atom needs to be fitted （包括sob老师在博文中给出的 .gjf 文件内容），但是更神奇的是，在同样不加额外的约束的前提下，我自己构建的模型中H原子的电荷并不是完全一致的，但是sob老师的模型中H原子RESP的计算结果就是完全一致的（可由插图中看出），又是一个不能理解的地方。

希望各位老师们能给予我这个巨菜无比的初学者一点“神谕”

sobereva

1 从笛卡尔坐标上看，你的H2CO明显不满足C2v点群。有点群对称性的体系在计算前一定要先在GaussView里对称化成实际的点群，这样算得又快又精确又容易收敛

2 也是因为你的体系不满足C2v对称性所致

YZH_jesse

sobereva 发表于 2020-11-30 00:52
1 从笛卡尔坐标上看，你的H2CO明显不满足C2v点群。有点群对称性的体系在计算前一定要先在GaussView里对称化 ...

感谢sob老师  在GaussView中设置对称化后很轻松地解决了上述问题。
虽然计算结果的精度满足我的需求，但是和您给出的甲醛分子初始笛卡尔坐标算得的结果并不是严格相同（小数点后第四位有差异）（由相同版本的Multiwfn算得），这个差异可以忽略不计吗？

sobereva

YZH_jesse 发表于 2020-11-30 10:24
感谢sob老师  在GaussView中设置对称化后很轻松地解决了上述问题。
虽然计算结果的精度满足我的需求，但 ...

你还是没有充分对称化，两个氢的电荷还是有差异。

YZH_jesse

sobereva 发表于 2020-11-30 12:10
你还是没有充分对称化，两个氢的电荷还是有差异。

sob老师  Test 3 是我对称化之后计算的结果  H原子的电荷是一致的  （Test 1和Test 2是我之前的计算结果  我在回复中没有写清楚）
现在我在尝试计算 间苯二胺（MPD）的RESP，还是存在很多问题。

我的计算流程为：

1. 在GaussView中构建一个MPD分子，在点群对称中设置C2v对称（如果存在这个等级的话），保存为 .gjf 格式的文件，关键词为：（目的是对这一初始构象进行几何优化）

1. %mem=2GB
   
2. %nprocshared=6
   
3. %chk=C:\MPD.chk
   
4. # opt b3lyp/6-311g(d,p) em=gd3bj

复制代码

2. 将优化之后的结果分别进行四种不同的处理：

    2.1 直接将优化之后的结果放到Multiwfn中进行RESP计算 （7 / 18 / 1)

    2.2 将优化之后的结果根据关键词 #HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) 在Gaussian中重新计算，将输出结果放到Multiwfn中进行RESP计算 (7 /18 / 8 / 1)

    2.3 将优化之后的结果根据您在 （http://bbs.keinsci.com/thread-16190-1-1.html）一文中提到的RESP2的模板进行计算，将输出结果（gas/solv）分别在Multiwfn中计算取平均（7 / 18 / 1）

    2.4 将优化之后的结果放到Antechamber （Amber 20 版本）中进行计算，命令行为：antechamber -i MPD.out -fi gout -o MPD.mol2 -fo mol2 -c resp

3. 结果分析：

    B.R. 为参考文献通过AM1-BCC算得的结果 （文献：Computer simulations of water flux and salt permeability of the reverse osmosis FT-30 aromatic polyamide membrane）（之前社长说不去计较文献的结果，我想的是如果我算的结果和他相近的话肯定是我的结果更为准确，但是计算结果相差太大了）

    在大量的尝试过程中，不同的初始建模经过几何优化（在gaussian中采用相同的关键词，如1.中所述）之后得到的构象共有三种（还有一种没画图），三种构象的RESP结果均不相同，而且和文献的结果均相差甚远。

想到可能存在问题的地方就是：
1. 在几何优化之前需要得到MPD分子能量最低的几何构象
2. 需要在点群对称上继续处理一下（可我在GaussView中已经将对称性设置成C2v了，而且我尝试过在几何优化时增加symm=loose关键词，结果也没什么改观）







PS：计算过程中的全部文件我都放到了附件中  （https://pan.baidu.com/s/1V5lAycdyoHz2PLpWlwgh0A） 提取码：qotl  

sobereva

MPD没有多构象问题，直接gview对称化成C2v，拿Gaussian算得到波函数，然后Multiwfn里直接用两步式RESP计算就完了
如果Multiwfn给你的结果不严格满足对称性，加上等价性设置