RESP拟合静电势电荷的原理以及在Multiwfn中的计算

wuzhiyi · 发表于 Post on 2019-5-29 03:14:59

本帖最后由 wuzhiyi 于 2019-5-29 05:55 编辑

谢谢sob老师的文章
有几个问题想请教一下sob老师关于氨基酸基于FEP的pKa计算一直很有问题，我觉得问题可能处在partial charge上，所以希望用更精确的partial charge来看可不可通过只更改电荷的方式获得更精确的结果。1，关于隐式溶剂的选取，文中说的是用默认的PCM而【谈谈隐式溶剂模型下溶解自由能和体系自由能的计算】则建议用SMD。想问一下，如果希望尽可能精确的话，是不是用SMD更好，还是有建议用PCM的原因？
2，关于基祖的选取，我感兴趣的系统是带负电的，在【原子电荷计算方法的对比】中图五是对比的是不带电的甲基。根据我的经验，不带电的系统对于电荷是非常tolerable的，partial charge哪怕用AM1－BCC也不会太糟，比如溶解自由能和实验结果很接近。但带电系统则有很大的问题，很多情况甚至R^2接近0，我想问一下，sob老师觉得如果用更好的泛函比如M06-2X相较于B3LYP或者更高的基祖比如aug-cc-PVTZ可以帮我们改善这个问题吗？文中说了用DFT方法来算cc-PVTZ的话效率很低，想问一下，假设我想用aug-cc-PVTZ的话，什么泛函比较好？系统<30个原子。谢谢sob老师

sobereva · 发表于 Post on 2019-5-29 06:05:12

wuzhiyi 发表于 2019-5-29 03:14
谢谢sob老师的文章
有几个问题想请教一下sob老师关于氨基酸基于FEP的pKa计算一直很有问题，我觉得问题可 ...

1 对于计算RESP电荷用哪个都可以，无所谓。因为只有溶剂模型的极性部分会影响电荷分布，而SMD的极性部分本来就是默认的IEFPCM。

2 A泛函算能量比B泛函好，绝对不代表A泛函算的电子密度分布也比B更精确。B3LYP完全没问题，用M06-2X只会白多花时间。DFT用aug-cc-pVTZ是极不划算的选择，要带弥散函数用ma-TZVP就完了，看
给ahlrichs的def2系列基组加弥散的方法
http://sobereva.com/340（http://bbs.keinsci.com/thread-3487-1-1.html）
带负电荷的体系，B3LYP/ma-TZVP妥妥的

wuzhiyi · 发表于 Post on 2019-5-29 15:30:30

本帖最后由 wuzhiyi 于 2019-5-29 15:35 编辑

sobereva 发表于 2019-5-29 06:05
1 对于计算RESP电荷用哪个都可以，无所谓。因为只有溶剂模型的极性部分会影响电荷分布，而SMD的极性部分 ...

谢谢sob老师的指导，不好意思我没有写的很清楚，之所以说PCM和SMD的问题因为
1，我读到在【谈谈隐式溶剂模型下溶解自由能和体系自由能的计算】中，1.3 节有写体系有的性质受隐式溶剂模型影响大，如激发能、HOMO-LUMO gap、偶极矩、原子电荷，所以很好奇这里原子电荷是怎么受不同隐式溶剂模型影响的。
2，因为做RESP肯定要算不同构象占比，而隐式模型SMD和PCM之间耗时区别不是很大，同在【谈谈隐式溶剂模型下溶解自由能和体系自由能的计算】中1.4节有写SMD算溶解自由能的精度是所有溶剂模型里拔尖的，所以我想更好的溶解自由能可以算更准确的占比，不知道可不可以让RESP更加反映现实。
3，同关于溶解自由能的问题，在和老板讨论后，目前方案是用B3LYP-D3(BJ)/6-311G**优化，然后算溶解自由能请问sob老师是建议B3LYP-D3(BJ)/6-311G**（IEEPCM）来保持机组泛函一致来算比较好，还是应该尊重隐式溶剂模型的定义用HF/6-31+G*（IEEPCM）或者M05（6）-2X/6-31G*（SMD）或者更好？

sobereva · 发表于 Post on 2019-5-29 16:29:57

wuzhiyi 发表于 2019-5-29 15:30
谢谢sob老师的指导，不好意思我没有写的很清楚，之所以说PCM和SMD的问题因为
1，我读到在【谈谈隐式溶剂 ...

1 不同隐式溶剂模型定义的分子孔洞形状、定义孔洞时用的原子球半径、表面显著电荷的计算方式等等很多数值细节方面不同，影响最终得到的溶质的波函数，进而影响原子电荷

2 从溶解自由能的计算准确度、进而影响构象平均的角度来说，用SMD原理上更合适

3 HF/6-31+G*（IEEPCM）完全没必要考虑。应当在SMD原文专门测试过的M05-2X/6-31G*级别下算

wuzhiyi · 发表于 Post on 2019-5-30 18:11:58

本帖最后由 wuzhiyi 于 2019-5-30 18:39 编辑

sobereva 发表于 2019-5-29 16:29
1 不同隐式溶剂模型定义的分子孔洞形状、定义孔洞时用的原子球半径、表面显著电荷的计算方式等等很多数值 ...

谢谢sob老师的耐心回复
现在准备
第一步用M05-2X/6-31G* scrf(SMD) opt freq 优化并获得液相自由能
第二步把有第一步优化好的结构再在气象下优化 (B3LYP-D3(BJ)/6-311G** opt freq) 得热校正量再把优化好的结构M06-2X/def2-TZVP算单点得气象自由能
第三步然后把第一步优化好的液相结构在B3LYP-D3(BJ)/6-311G** scrf(SMD)算电荷最后组合算RESP

还有另一个问题想请教一下sob老师
对于scrf(SMD,solvent=water)来说eps和epsinf都是真正的水
但对分子动力学来说，一般的tip3p和tip4p水都有不一样的eps和epsinf，请问一下sob老师觉得需不需要自定义eps和epsinf以求取得更加契合tip3p和tip4p水的RESP？
谢谢

sobereva · 发表于 Post on 2019-5-31 07:20:20

wuzhiyi 发表于 2019-5-30 18:11
谢谢sob老师的耐心回复
现在准备
第一步用M05-2X/6-31G* scrf(SMD) opt freq 优化并获得液相自由能

你的步骤完全莫名其妙，意义不明

M05-2X/6-31G* scrf(SMD) opt freq这种极差的档次怎么可能算准液相自由能。该怎么算这里写明了
谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
http://sobereva.com/327（http://bbs.keinsci.com/thread-3345-1-1.html）

不需要，纯粹想多了。本来TIP3P/4P水模型的原子电荷也不是RESP算的

wuzhiyi · 发表于 Post on 2019-5-31 16:51:38

sobereva 发表于 2019-5-31 07:20
你的步骤完全莫名其妙，意义不明

M05-2X/6-31G* scrf(SMD) opt freq这种极差的档次怎么可能算准液相自 ...

不好意思，我是第一次做考虑多构象的RESP，还不是很清楚，我再仔细看了一下，不知道我的理解对不对。

因为我的系统带电，而且一端带正电另一边带负电，在液相下则有很多稳定结构，而气相下很多液相下稳定的结构会变形出现头尾相连。而因为要在tip3p的环境下模拟，所以在考虑多构象的时候应该用溶剂下的Boltzmann权重。

所以要算每个构象的自由能=气象自由能+溶解自由能

在【谈谈隐式溶剂模型下溶解自由能和体系自由能的计算】的2 溶解自由能的计算方法中说
溶解自由能=溶剂单点能－气象单点能
因为用SMD溶剂，所以溶剂单点能必须要在M05-2X/6-31G* scrf(SMD)下算，所以B3LYP-D3(BJ)/6-311G** opt scrf(SMD)优化，然后M05-2X/6-31G* scrf(SMD)算单点能。
然后算气象单点能，把第一部优化好的结构按进M06-2X/def2-TZVP算单点。
溶解自由能 = G（M05-2X/6-31G* scrf(SMD)）－G（M06-2X/def2-TZVP）

然后算气象自由能
气象自由能=热校正量+单点能
按照【RESP拟合静电势电荷的原理以及在Multiwfn中的计算】3.2的说法，我应该在气象下优化 (B3LYP-D3(BJ)/6-311G** opt freq) 得热校正量再把优化好的结构M06-2X/def2-TZVP算单点能

但在液相下的一些稳定结构在气象下不稳定所以我不能在气象下优化结构，但如果不是同泛函／基组的下的优化好的结构的话无法用freq算热矫正能。
想问一下sob老师的建议是什麽？

另不知道是不是我理解错了，我觉得
因为液相自由能=气象自由能+溶解自由能
而气象自由能=热校正量+气象单点能
溶解自由能=溶剂单点能－气象单点能
所以液相自由能=（热校正量+气象单点能）+（液相单点能－气象单点能）
=热校正量+溶剂单点能
所以其实没有必要算气象单点能能？可以直接B3LYP-D3(BJ)/6-311G** opt freq scrf(SMD)获得优化结构和热校正量，然后M05-2X/6-31G* scrf(SMD)算液相单点能，这两个相加就是液相自由能？

sobereva · 发表于 Post on 2019-6-1 09:25:51

wuzhiyi 发表于 2019-5-31 16:51
不好意思，我是第一次做考虑多构象的RESP，还不是很清楚，我再仔细看了一下，不知道我的理解对不对。

...

我文中已经说了（我这里引用的是目前最新的版本）

(3)对绝大多数中性且局部不显离子性的分子，气相和溶剂下结构差异甚微，此时结构优化及之后的任务用气相下优化的结构就够了。但如果溶剂效应对当前溶质结构影响可能较大，或者想追求稳妥，结构优化、振动分析获得热力学校正量都应当在隐式溶剂模型下进行，并且之后算溶解自由能和高精度单点能的时候也用这个结构。

当然得算气相单点能，不同步骤用的级别不同，显然不能简单抵消

wuzhiyi · 发表于 Post on 2019-6-1 17:29:56

sobereva 发表于 2019-6-1 09:25
我文中已经说了（我这里引用的是目前最新的版本）
当然得算气相单点能，不同步骤用的级别不同，显然不能 ...

不好意思，我还是不太懂，所以是应该
每个构象的自由能=液相热校正量+气象单点能+溶解自由能
第一步B3LYP-D3(BJ)/6-311G** opt freq scrf(SMD)获得优化结构和液相热校正量。
第二步算溶解自由能=液相单点能－气象单点能，就是把第一步优化的结构算M05-2X/6-31G* scrf(SMD)的液相单点能减去M05-2X/6-31G*的气象单点能。
第三步用M06-2X/def2-TZVP算气象单点能
最后把第一步液相热校正量，加上M05-2X/6-31G*下液相与气象单点能差，最后加上M06-2X/def2-TZVP下气象单点能就是最后的每个构象的自由能是么

sobereva · 发表于 Post on 2019-6-1 20:07:16

wuzhiyi 发表于 2019-6-1 17:29
不好意思，我还是不太懂，所以是应该
每个构象的自由能=液相热校正量+气象单点能+溶解自由能
第一步B3L ...

是

wuzhiyi · 发表于 Post on 2019-6-2 17:26:43

sobereva 发表于 2019-6-1 20:07
是

谢谢 sob老师

sobereva · 发表于 Post on 2019-10-19 13:03:25

今日更新了官网上的Multiwfn，一个重要改进是计算RESP电荷的功能支持了根据结构对称性自动设置等价性约束，只需要把某个区域的原子序号输入，Multiwfn就能自动判断点群并将等价的原子写入等价性约束设置文件，使得局部或整体有对称性的体系的RESP电荷计算明显更为便利。与此同时，此版本还令两步式标准RESP电荷的第一步支持了自定义等价性和电荷约束（以前这两种约束只对一步式拟合有效），这使得RESP电荷拟合更为灵活。

已更新本文，特别是添加了3.5节，通过例子充分体现了本次的改进。

wuzhiyi · 发表于 Post on 2019-10-21 22:18:23

想問一下sob老師，我有一個磷脂，我像計算他的電荷，磷脂的頭部應該是在水里，所以eps=80，而磷脂的尾部則是在eps=2的環境里，想請教一下sob老師這種情況應該如果做。我感覺應該不能對一個分子的不同部位用不同的隱式模型。

sobereva · 发表于 Post on 2019-10-21 22:23:39

wuzhiyi 发表于 2019-10-21 22:18
想問一下sob老師，我有一個磷脂，我像計算他的電荷，磷脂的頭部應該是在水里，所以eps=80，而磷脂的尾部則 ...

其实都用隐式水模型就行，因为烃链部分极性很低，本身受到隐式水环境的极化程度也很低
更周到的话，也可以考虑用电荷约束，烃链部分约束总电荷为0，其余部分也约束成相应的整数。在水模型下拟合一次，得到烃链以外的部分的电荷，再在隐式溶剂模型表现烷烃环境下拟合一次，得到烃链的电荷。实际上很多磷脂力场都直接给出了烃链上CH2、CH3的原子电荷，尾部其实也不用自己拟合，就用现成的就行了，这样还更合理点。

wuzhiyi · 发表于 Post on 2019-10-22 20:08:42

sobereva 发表于 2019-10-21 22:23
其实都用隐式水模型就行，因为烃链部分极性很低，本身受到隐式水环境的极化程度也很低
更周到的话，也可 ...

謝謝sob老師

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[Multiwfn资源与经验] RESP拟合静电势电荷的原理以及在Multiwfn中的计算