计算化学公社

标题: 在使用隐式溶剂模型优化结构并计算单点能时遇到的问题 [打印本页]

作者
Author: TUT 时间: 2025-5-30 16:34
标题: 在使用隐式溶剂模型优化结构并计算单点能时遇到的问题
本帖最后由 TUT 于 2025-5-30 18:25 编辑

近期在使用Gaussian隐式溶剂模型时遇到了一个难题，自己钻研很久仍未找到合适的解法，希望得到各位老师的指点！具体如下：

我想计算同一物质的两种构象a和b在溶剂中（以乙腈为例）的相对能量，实验已知b是更稳定的构象，所以模拟预期得到的结果应该是b构象相对于a构象的能量更低。我先对a和b真空状态下的构象使用B3LYP-D3BJ/6-31G+(d,p)进行了结构优化，再用M062x-D3/6-31G+(d,p)计算单点能，此时与预期一致得到了b的能量更低的结果。

但当我进一步计算它们在乙腈隐式溶剂模型中的相对能量时，尽管进行了多种尝试，得到的结果却都不符合预期：

第一次尝试：首先我对真空优化后的a和b在乙腈隐式溶剂模型IEFPCM中使用B3LYP-D3BJ/6-31G+(d,p)进行结构优化，并在SMD中使用M062x-D3/6-31G+(d,p)计算能量，可得到的结果却是a的能量更低，与预期不符。

第二次尝试：经过观察后发现，在IEFPCM中优化a的结构，使a的结构相对于真空时发生了肉眼能看出的微小变化，猜测这可能是结果不合理的原因，于是我决定在结构优化这一步也使用SMD模型。但SMD模型优化a的结构时，因为遇到了Inv3 failed in PCMMkU (L502)的溶剂孔洞报错，所以只能在不加弥散并使用SAS表面的情况下才用B3LYP-D3BJ/6-31G (d,p)优化成功（b结构不加SAS也不会报错，但为保持计算水平严格一致也加了SAS。虽然乙腈中加弥散可以计算成功，但在有需要的另一种溶剂乙酸乙酯中加弥散就会报错），并再用SES计算单点能，依然得到了a的能量更低的结果。

第三次尝试：考虑到隐式溶剂模型对结构优化影响不大，不再用隐式溶剂模型做结构优化，直接拿真空优化出的结构在乙腈的SMD中用M062x-D3/6-31G+(d,p)计算能量，结果还是a的能量更低。

从实验结果来在各溶剂中看应该是b比a更稳定，可我的计算结果却一直无法匹配（总是a的能量更低），想请问各位老师我的计算过程是不是哪里出了问题，以及我这个体系需不需要加弥散呢？

附上使用SMD做结构优化时的gif，几次尝试的初始gif结构都是一样的
(, 下载次数 Times of downloads: 8)
(, 下载次数 Times of downloads: 3)
(, 下载次数 Times of downloads: 60)

作者
Author: chands 时间: 2025-5-30 20:44
本帖最后由 chands 于 2025-5-30 20:47 编辑

我用你附件里的结构、级别和溶剂模型直接计算单点能
A=-1196.620869 Hartree
B=-1196.621324 Hartree
B低了11 kJ/mol 左右吧。不过严格来说要比较自由能。

作者
Author: TUT 时间: 2025-5-31 11:53

chands 发表于 2025-5-30 20:44
我用你附件里的结构、级别和溶剂模型直接计算单点能
A=-1196.620869 Hartree
B=-1196.621324 Hartree

原来还要计算自由能，谢谢老师指导！附件里的gif是真空优化出的结构，请问老师是在m062x-D3/6-31G+(d,p)下计算的能量吗？不经过溶剂模型做结构优化直接用溶剂模型计算单点能的结果也是可用的吗？

作者
Author: chands 时间: 2025-5-31 12:15

TUT 发表于 2025-5-31 11:53
原来还要计算自由能，谢谢老师指导！附件里的gif是真空优化出的结构，请问老师是在m062x-D3/6-31G+(d,p) ...

要计算自由能差值，因为自由能差值（等温等压下是Gibbs自由能）与平衡系数直接相关。
我按照你gjf的级别算的，没动过。
除非溶质分子与溶剂分子有相当强的作用（比如形成氢键），真空优化的结构和加溶剂模型优化的结构应该差不多。你用真空优化也没啥问题。考虑到你的分子有对称性，优化和振动分析时加上对称性会节省时间。

作者
Author: TUT 时间: 2025-5-31 18:06

chands 发表于 2025-5-31 12:15
要计算自由能差值，因为自由能差值（等温等压下是Gibbs自由能）与平衡系数直接相关。
我按照你gjf的级别 ...

好的，那我再算一下试试，感谢老师的帮助

作者
Author: sai77 时间: 2025-5-31 20:51
本帖最后由 sai77 于 2025-5-31 20:55 编辑

为什么单点能的基组跟优化保持一样的？至少应该提升一到两个档次吧
另外结构a应该是C2v对称的，你给的结果只优化到了C2，还可以进一步调整对称性继续优化

作者
Author: chands 时间: 2025-6-1 00:50

sai77 发表于 2025-5-31 20:51
为什么单点能的基组跟优化保持一样的？至少应该提升一到两个档次吧
另外结构a应该是C2v对称的，你给的结果 ...

单点能的级别是要高点，我同意。
对称性方面，Gaussian默认的对称性的阈值比较严，差一丁点都不行，我的经验是Gaussian经常判低了对称性，除非更改Tolerance。

作者
Author: TUT 时间: 2025-6-3 20:38
本帖最后由 TUT 于 2025-6-3 20:40 编辑

sai77 发表于 2025-5-31 20:51
为什么单点能的基组跟优化保持一样的？至少应该提升一到两个档次吧
另外结构a应该是C2v对称的，你给的结果 ...

其实我尝试过用mo62x/def2TZVP算单点能，只是结果还是背离预期的a的能量比b更低，所以就没有继续用了。请问老师调整a的对称性应该怎么操作呢？

作者
Author: sai77 时间: 2025-6-4 07:39

TUT 发表于 2025-6-3 20:38
其实我尝试过用mo62x/def2TZVP算单点能，只是结果还是背离预期的a的能量比b更低，所以就没有继续用了。请 ...

Tools-Point Group里调整

作者
Author: TUT 时间: 2025-6-4 10:54

sai77 发表于 2025-6-4 07:39
Tools-Point Group里调整

好的，谢谢老师的帮助

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