计算化学公社

标题: 关于计算RESP2电荷的一个疑问 [打印本页]

作者
Author: a9471163 时间: 2021-12-24 18:09
标题: 关于计算RESP2电荷的一个疑问
本帖最后由 a9471163 于 2021-12-24 18:18 编辑

各位老师，我想做氧化石墨烯（GO）在OPLS-AA力场下的分子动力学模拟，由于羧基会去质子化形成带电体系，我选择RESP2（0.5）电荷作为原子电荷。
首先，我用xtb在extreme和GFN2-xTB这种优化级别下，施加隐式水溶剂模型完成了几何优化，随后使用B3LYP-D3(BJ)/def2-SVP级别，使用ORCA计算两个单点能（即SMD溶剂模型下和真空下的单点能。原子数约300个，模型有2个羧基，模型如图。体系较大，所以优化结构使用了xtb）
随后，使用Multiwfn计算RESP电荷。我先计算羧基去质子化的GO片层在水溶剂模型下的RESP电荷，使用电荷约束使得羧基的三个原子（C、O、O）的总带电量为-1，并且使用等价性约束使得羧基的两个O原子带电量相同，这步计算没有遇到问题。这步计算，我的模型所带RESP电荷的最终值是-2（如图）。
随后，我计算羧基没有去质子化的GO片层在真空中的RESP电荷（这是因为在真空中羧基的H不会解离，所以我必须带着羧基的H去计算），没有施加等价性约束和电荷约束，计算结果如图，可见真空中的所有原子RESP电荷之和（除了羧基的2个H原子）是-0.77277。如果加上2个羧基中的H原子，则真空中的RESP电荷之和是0。由于分子动力学模拟中，两个羧基的氢原子是要去掉的，所以我的模型所带RESP电荷的最终值就是-0.77277。
随后，取上述两步计算所得RESP电荷的平均值，作为RESP2(0.5），如图，可见我的模型中所有原子的RESP2(0.5）之和为-1.386
但-1.386这个值在分子动力学模拟时带来了问题，因为不是整数，加抗衡的Na离子后体系还是带电，想请教一下这个问题该怎么办呢？

我自己想了2套解决问题的方案，还请大佬们把把关，这样做合理吗：
1、我计算没有去质子化的GO片层在真空中的RESP电荷时，施加一个电荷约束，让2个羧基的H原子所带电荷之和为+2，则此时其余原子所带电荷刚好是-2，模型中所有原子的RESP2(0.5）之和将是-2。
2、分子动力学模拟时，GO片层带-1.386电荷就已经可以认为是合理的，只需编辑力场文件，改变Na离子的带电量，使其小于1，恰好中和掉-1.386的带电量即可，不知这么做是否合理？
(, 下载次数 Times of downloads: 16) (, 下载次数 Times of downloads: 19)

作者
Author: 喵星大佬 时间: 2021-12-24 18:24
一定哪里操作出问题了，不管用不用溶剂模型，RESP电荷总数肯定是你的净电荷

作者
Author: a9471163 时间: 2021-12-24 18:26
本帖最后由 a9471163 于 2021-12-24 18:30 编辑

喵星大佬发表于 2021-12-24 18:24
一定哪里操作出问题了，不管用不用溶剂模型，RESP电荷总数肯定是你的净电荷

您好。我没加溶剂模型时，是带着羧基上的H原子去计算了RESP电荷，这时候羧基的H原子化学价是+1，但实际的RESP电荷肯定是比+1小的，算出来就是+0.36左右，这里不是整数，然后整个体系的RESP电荷是0，是整数，减去羧基H原子的RESP电荷以后，整个体系的RESP电荷就不是整数了，问题就是在这一步出现的。

我知道只需加个电荷约束，让羧基的H原子带的电量是+1就能解决问题，但合理吗，羧基的H原子化学价是+1，但实际的RESP电荷肯定是比+1小的。我规定成+1，感觉确实有点勉强

作者
Author: eagletyr 时间: 2021-12-24 19:03
你在MD的时候是去质子化的，拟合RESP电荷的时候肯定就要用去质子化的模型来计算

作者
Author: a9471163 时间: 2021-12-24 19:18

eagletyr 发表于 2021-12-24 19:03
你在MD的时候是去质子化的，拟合RESP电荷的时候肯定就要用去质子化的模型来计算

感谢大佬解惑，计算真空里RESP电荷时模型选错了，应该改为与水溶液里使用的模型完全一致，并且施加-2的总电荷

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