关于应用Gaussian计算复合物结合能筛选功能单体的方法

johngu111

请教各位老师，我想从功能单体1、2、3中筛选出对目标分子1结合能力最强的单体，一般文献都是直接计算单体与目标分子的结合能，筛选出结合能最强的单体（假设为单体1）。但这里如果还存在竞争分子2和3，虽然筛选出的单体1对目标分子1的结合能最大，但如果同时对竞争分子2或3的结合能也较大，是不是说明单体1就不是理想的选择性单体。所以如果想从1、2、3中筛选出对目标分子1具有最强选择性（假设只有2和3作为竞争物）的单体，我该怎么计算筛选才比较最合理呢，没找到这么做的文献，不知道Gaussian能否从理论上实现这一选择性的定量化。

wzkchem5

首先，除非分子极大、跑不动频率计算，否则应该用结合自由能而非结合能。
其次，把结合自由能换算成平衡常数，就可以知道单体1、2、3浓度相等，且目标分子1、竞争分子2、3浓度也相等时，各个单体和各个目标/竞争分子的结合比例。然后一般取能让单体1和目标分子1所形成复合物的平衡浓度最大的单体即可。对于结合是动力学控制的情形，还要考虑结合速率常数的影响。

johngu111

wzkchem5 发表于 2023-4-17 18:55
首先，除非分子极大、跑不动频率计算，否则应该用结合自由能而非结合能。
其次，把结合自由能换算成平衡常 ...

谢谢老师，我还有两个小问题：
1.您说的结合自由能是ΔG（ΔG=G(复合物)-G(单体)-G(目标分子)）吗？
2.ΔG=-RTln(K)，K是平衡常数，假设单体和目标/竞争分子浓度都是1：1，复合物形态也是1：1，那K=c(复合物)，分别计算每个单体对目标分子1、竞争分子2、竞争分子3的平衡常数K1、K2、K3；其中K1比例最大的单体就是对目标分子1来说的最优单体是吧，我这样理解对吗？

wzkchem5

johngu111 发表于 2023-4-17 13:00
谢谢老师，我还有两个小问题：
1.您说的结合自由能是ΔG（ΔG=G(复合物)-G(单体)-G(目标分子)）吗？
2. ...

1. 是
2. ”单体和目标/竞争分子浓度都是1：1“和”复合物形态也是1：1“不可能同时成立。应该是单体和目标/竞争分子的初始浓度是1:1

johngu111

wzkchem5 发表于 2023-4-17 20:25
1. 是
2. ”单体和目标/竞争分子浓度都是1：1“和”复合物形态也是1：1“不可能同时成立。应该是单体和 ...

这个结合比例弄得我有点晕了，既然K(复合物)=1是不可能的，单体和目标/竞争分子结合的比例就不可能是1，那我计算的时候预设的单体-目标/竞争分子为1-1的复合物构型还有意义吗？

wzkchem5

johngu111 发表于 2023-4-17 13:43
这个结合比例弄得我有点晕了，既然K(复合物)=1是不可能的，单体和目标/竞争分子结合的比例就不可能是1， ...

搞清楚平衡常数和结合比例的区别。把平衡常数的公式写出来，和结合比例的公式对比

johngu111

wzkchem5 发表于 2023-4-17 20:51
搞清楚平衡常数和结合比例的区别。把平衡常数的公式写出来，和结合比例的公式对比

哦哦，谢谢老师，是我弄错了，平衡常数的分母是单体和目标分子的平衡浓度，您说的浓度相等都是假设同一平衡浓度的情况下，然后对比K1、K2、K3的比例，取K1最大的单体就是最优单体对吧

wzkchem5

johngu111 发表于 2023-4-17 14:12
哦哦，谢谢老师，是我弄错了，平衡常数的分母是单体和目标分子的平衡浓度，您说的浓度相等都是假设同一平 ...

嗯，假设平衡浓度都一样也可以，假设初始浓度都一样也可以

johngu111

wzkchem5 发表于 2023-4-17 21:28
嗯，假设平衡浓度都一样也可以，假设初始浓度都一样也可以

我看文献里有机吸附体系中这种单体筛选的计算，一般都是用的结合能，这里为什么用结合自由能更好呢，老师这里有相关的帖子或者文献供我学习一下吗？

wzkchem5

johngu111 发表于 2023-4-17 14:57
我看文献里有机吸附体系中这种单体筛选的计算，一般都是用的结合能，这里为什么用结合自由能更好呢，老师 ...

因为实验看结合强度用的是K，dG = -RTlnK，所以dG和K直接相关。但你没办法把结合能dE表示成K的函数，所以dG显然比dE更适合解释实验的结合强度。
有些文章用结合能是因为：（1）结合能好算，不需要算频率也能拿到一个数据（但这个数据是不严格的，见下）；（2）当用RRHO计算自由能时，自由能受数值误差影响较大，容易算不准（N. B. Andrea, W. Steven, ChemRxiv 2019）。但是如果不算频率的话，第一没法做零点能校正，第二（也是更重要的一点）没法确认这个结构没有虚频，所以严格来说算出来的结合能本身都不是严格的。而RRHO计算自由能不准的问题也可以通过改用Shermo（http://sobereva.com/552）读取高斯的输出文件重新算自由能，或者改用ORCA做计算的方法来解决。所以某些高斯用户用结合能来做筛选的做法，要么是图计算快而牺牲了严谨性，要么是不会用Shermo或ORCA，脱离时代了。

johngu111

wzkchem5 发表于 2023-4-17 22:56
因为实验看结合强度用的是K，dG = -RTlnK，所以dG和K直接相关。但你没办法把结合能dE表示成K的函数，所以 ...

谢谢老师提供的资料，但我记得之前用高斯算结合能的时候也是用opt+freq计算的，优化出的构型不也是无虚频的最低能量构型吗？跟这个自由能计算模式有什么本质的区别？

wzkchem5

johngu111 发表于 2023-4-17 16:19
谢谢老师提供的资料，但我记得之前用高斯算结合能的时候也是用opt+freq计算的，优化出的构型不也是无虚频 ...

读取的量不一样，能量和自由能不是从同一个地方读的。
另外还有一个问题：一般建议在更高理论级别下算单点能，加上opt+freq级别下的自由能校正得到自由能。而不是在opt+freq级别下直接读取自由能

johngu111

wzkchem5 发表于 2023-4-17 23:21
读取的量不一样，能量和自由能不是从同一个地方读的。
另外还有一个问题：一般建议在更高理论级别下算单 ...

根据您的建议，我想尝试一下计算吸附剂的季胺单元与硝酸根离子的结合自由能，列出了如下方案：
一.使用B3LYP/6-31G*对季胺单元（A)、硝酸离子（B）和它们的复合物（A-B）进行结构优化和频率计算，然后使用PWPB95-D3(BJ)/def2-QZVPP结合SMD溶剂模型对优化的结构进行单点能计算。根据公式ΔE=E(A-B)-E(A)-E(B)计算结合能。
二.计算结合自由能，将PWPB95-D3(BJ)/def2-QZVPP级别下得到的单点能（E）加上B3LYP/6-31G*级别下的Gcorr得到G。根据公式ΔG=G(A-B)-G(A)-G(B)得到ΔG，用ΔG换算平衡常数。
1.这个计算步骤合理吗？
2.构建吸附剂单元的时候，因为计算资源有限，所以简化成一个苯环+功能单体的结构（图1），不知道合不合理，还是直接用功能单体的结构更好（图2）？
3.PWPB95-D3(BJ)/def2-QZVPP泛函Gaussian09中没有找到，因为计算资源不富裕，不知道这个泛函基组搭配性价比怎么样，还能选择精度更低一点的泛函和基组吗？

wzkchem5

johngu111 发表于 2023-4-18 04:55
根据您的建议，我想尝试一下计算吸附剂的季胺单元与硝酸根离子的结合自由能，列出了如下方案：
一.使用B ...

有几个问题：（1）必须做构型搜索，初始构型不能随便用手摆；（2）必须加色散校正；（3）如果溶剂是不太弱的氢键给体（例如至少有水、醇那么强），必须对NO3-做显式溶剂化；（4）opt和freq的基组可以考虑大一点，例如6-311G**，毕竟体系小，算得动；氧最好加弥散；（5）NO3-是平面的，初始结构应该画成平面。此外注意Gcorr最好用Shermo算（或者用orca做opt和freq计算），而非读高斯输出的值，因为高斯的RRHO方法对于弱相互作用体系不准。
吸附剂整体有多大？一般来说，与其截断，不如用QM/MM、ONIOM等方法近似考虑吸附剂其他基团的影响。不然你截断引入的误差可能都远大于计算方法误差了，虽然用PWPB95但是精度远达不到PWPB95的实际精度。
用ORCA算，不算太慢。如果选择精度低的泛函和基组，你得做好说服审稿人的准备，万一审稿人问你为什么不用PWPB95等双杂化泛函，你得能回答才行，又不能以”高斯里没有PWPB95“、”高斯算双杂化泛函的功能做得不好“为理由，毕竟orca是免费的，你没道理不用。如果说你实际试了一个双杂化泛函计算，发现确实跑不动，这个倒可以算一个理由，但这就意味着你必须做至少一个双杂化泛函计算。

johngu111

wzkchem5 发表于 2023-4-18 15:29
有几个问题：（1）必须做构型搜索，初始构型不能随便用手摆；（2）必须加色散校正；（3）如果溶剂是不太 ...

谢谢老师的耐心解答，我大概知道怎么做了，我自己继续研究研究。