显著改进SMD溶剂模型描述Br和I元素精度的溶剂模型SMD18的介绍

sobereva

显著改进SMD溶剂模型描述Br和I元素精度的溶剂模型SMD18的介绍

Introduction of the solvent model SMD18, which significantly improves the accuracy of the SMD solvent model in describing Br and I elements

文/Sobereva@北京科音  2021-Aug-23

之前笔者在《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）中专门说过SMD溶剂模型，这是现阶段最适合算溶剂环境下体系能量的隐式溶剂模型。SMD模型需要通过原子球叠加来构建溶质的孔洞用于算溶剂-溶质相互作用的静电部分，定义原子球需要原子半径。在SMD原文里，仅对H、C、N、O、F、Si、P、S、Cl、Br专门优化了半径，对其它元素用的是Bondi范德华半径来凑合（如果Bondi范德华半径都没有，就用2.0埃凑合）。

SMD溶剂模型作者后来的一篇文章Chem. Eur. J., 24, 15983 (2018)是SMD18溶剂模型的原文。这篇文章根据涉及卤素的体系在乙腈中的反应自由能，发现对于碘用Bondi范德华半径1.98埃来凑合算的结果很不好，此半径值明显太小，于是改为了1991年SM1溶剂模型里对碘优化的半径2.74埃，这令计算结果与实验数据相符程度有了极大的提升。后来他们对Br的半径进行了专门优化，发现用2.60埃做半径的结果与实验数据最接近，这比SMD原文里专门给出的Br的半径3.06埃要小不少，这也令结果大有改进。将SMD里Br和I的半径分别改为2.60埃和2.74埃之后，就被此文称为SMD18溶剂模型。

对于溶液中的涉及卤键的形成或断裂的过程，与Br、I关系密切的化学过程（例如成键断键），或者算含有Br和I的体系的溶解自由能，都强烈建议使用SMD18代替SMD模型。

根据文中的测试，SMD18模型非常适合结合M06-2X使用，我也建议用SMD18时都用M06-2X泛函。为了最大程度地误差抵消，计算溶液中自由能（1M标准态浓度）用的流程也最好和原文一致，具体来说是：对溶质在M06-2X下做优化和振动分析，在这个过程就带着SMD18溶剂模型，直接取量子化学程序输出的自由能，然后再加上1.89 kcal/mol的1atm->1M标准态转变自由能变即可。对Cl、Br、I应当用def2-TZVPD（此时对I来说是赝势基组，需要结合相应的赝势），对其它元素用def2-TZVP。如果不清楚怎么在Gaussian中这么用基组，看《给ahlrichs的def2系列基组加弥散的方法》（http://sobereva.com/340）和《详解Gaussian中混合基组、自定义基组和赝势基组的输入》（http://sobereva.com/60）。按照SMD18原文所述，如果存在100cm^-1波数以下振动频率的情况，应当改用Grimme的准RRHO模型计算热力学量而非用Gaussian默认用的RRHO模型，这可以通过Shermo程序实现，详见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552），后文有例子。

PS：当然了，用def2-TZVP、def2-TZVPD做opt freq很昂贵，显然不是可取的做法，但无奈SMD18原文这么干的，为了和SMD18优化的半径尽可能实现误差抵消，自己计算时也只能忍着了。如果实在算着太吃力的话，也可以opt freq时候用便宜一些的基组比如def2-TZVP(-f)（这么做对于自由能热校正量和结构的影响很小），之后再用前述SMD18原文里用的基组算个单点，之后再把电子能量和自由能热校正量相加。另外，原理上考虑ZPE校正因子更好，但也为了与SMD18原文一致，就不去考虑这点了。

目前的Gaussian不支持SMD18溶剂模型，但只要在SMD溶剂模型计算时改一下Br和I的元素半径即可轻易实现。这需要写scrf(SMD,solvent=溶剂名,read)关键词，输入文件末尾空一行写modifysph，下面再空一行写以下内容即可
I 2.74
Br 2.60

下面给一个Gaussian的完整的计算例子输入文件。输入文件和Gaussian 16跑的输出文件可以在http://sobereva.com/attach/613/file.zip里得到。

#P M062X/genecp SCRF(SMD,solvent=acetonitrile,read) opt freq
[空行]
Title
[空行]
0 1
N                  2.53946200   -1.07717200   -0.01319200
C                  1.74543500    0.00000700   -0.00027500
N                  2.53948000    1.07717200    0.01885100
C                  3.85512000   -0.67530300   -0.00660900
H                  4.66037300   -1.38842900   -0.01940300
C                  3.85513000    0.67528200    0.00619900
H                  4.66038900    1.38840200    0.01901300
I                 -0.36464300    0.00000000    0.00004200
C                  2.13385200    2.47671900   -0.00046400
H                  1.05119700    2.53713200    0.26086500
H                  2.74614800    3.02347800    0.75566000
H                  2.31269200    2.89578200   -1.02444600
C                  2.13383000   -2.47671400   -0.00044100
H                  1.05117500   -2.53711700   -0.25836000
H                  2.74216300   -3.02400000   -0.75519100
H                  2.31661900   -2.89526100    1.02491500
Br                -3.43164200    0.00000200    0.00045300
[空行]
@/sob/SMD18.gbs/N
[空行]
@/sob/def2-ECP.txt/N
[空行]
modifysph
[空行]
I 2.74
Br 2.60
[空行]
[空行]

被引入的SMD18.gbs是BSE基组数据库的def2-TZVP基组文件，但是把Cl、Br、I替换成了def2-TZVPD之后的情况。def2-ECP.txt是def2系列基组标配的Stuttgart赝势基组的定义。即便你的体系里没有I，用上面的输入文件直接替换坐标部分后也可以算，只不过读取赝势的设定以及修改I的半径不生效而已。这两个文件可以在http://sobereva.com/attach/613/SMD18.gbs和http://sobereva.com/attach/613/def2-ECP.txt下载，上例放到了/sob目录下。

对于Gaussian 09，特别要注意对元素半径的修改没法在opt freq任务中传递给freq部分，因此必须opt和freq分别做，这是个bug，在Gaussian 16中得到了修正。

上面的例子的最低频率只有44.24 cm^-1，因此不建议直接读取输出文件里的基于RRHO模型算出来的Sum of electronic and thermal Free Energies，否则低频对熵的贡献很不准确，而应当用Shermo程序基于准RRHO模型算。把Shermo的settings.ini里的ilowfreq=后面的值设为2，启动Shermo并载入上例的输出文件，可看到结果Sum of electronic energy and thermal correction to G:       -3176.6090365 a.u.。将之加上标准态转换自由能变1.89/627.51=0.00301 a.u.后即得到1M浓度下的乙腈中的自由能-3176.60602 a.u.。

SMD18模型体现出隐式溶剂模型计算能量的准确度对原子半径是多么的敏感。顺带一提，J. Phys. Chem. A, 123, 9498 (2019)提出的SMD-B模型将SMD定义的原子半径都改为了Bondi范德华半径，发现对离子体系（特别是含硫的）的溶解自由能得到一定改进。既可以如上通过modifysph来对各个元素一一修改半径，也可以在scrf里写read的同时在输入文件末尾后面空一行直接加上Radii=Bondi来实现，当然后者省事得多。

lanthanum

请问，如果想综合smd18和smd-b的成果，应该怎么写输入文件？是写成下面这样么？
Radii=Bondi
modifysph
[空行]
I 2.74
Br 2.60

sobereva

lanthanum 发表于 2021-8-24 02:18
请问，如果想综合smd18和smd-b的成果，应该怎么写输入文件？是写成下面这样么？
Radii=Bondi
modifysph

大抵是这样
计算开始后输出文件里会给出各个原子实际用的半径，你看看是不是你期望的

lanthanum

sobereva 发表于 2021-8-24 06:20
大抵是这样
计算开始后输出文件里会给出各个原子实际用的半径，你看看是不是你期望的

谢谢社长。

lanthanum

还有三个疑问。谢谢老师解答。
（1）SMD18相对于SMD，其实只是改了Br和I的半径，其他参数并没有改变，因此，算溶解自由能的时候，Br和I应当用def2-TZVPD，而其它元素还应该用6-31g(d)。这样的话，在Br和I数量为零的分子里，SMD成为SMD18的一个特例，是自洽的。我这样理解正确否？
（2）计算阴离子时，选smd模型，那按照【高精度单点（带弥散的大基组，气相）+溶解自由能（6-31g(d)，溶剂-气相）】这个流程算的准，还是按【高精度单点（带弥散的大基组，溶剂）】这样算的准？
（3）计算阴离子时，如果选smd-b模型，按哪个流程算的准呢？

sobereva

lanthanum 发表于 2021-8-25 09:18
还有三个疑问。谢谢老师解答。
（1）SMD18相对于SMD，其实只是改了Br和I的半径，其他参数并没有改变，因此 ...

1 你这样做也不是不可以，有一定道理，但怎么做才最理想需要测试。
2 看具体什么泛函。M05-2X的话，从统计上来说应当前者更准，但后者不一定与之有多大差距，而且也有可能碰到个别情况反倒后者更准。反正都不如uESE模型准。
3 看SMD-B原文。我倒觉得SMD-B不是特别重要，算离子都建议用uESE

lanthanum

sobereva 发表于 2021-8-25 11:29
1 你这样做也不是不可以，有一定道理，但怎么做才最理想需要测试。
2 看具体什么泛函。M05-2X的话，从统 ...

谢谢社长。

喵星大佬

这么看来要是做一种按照某种规则动态调整原子半径的隐式溶剂模型似乎有明显提高精度的潜力

wzkchem5

喵星大佬 发表于 2021-8-25 08:50
这么看来要是做一种按照某种规则动态调整原子半径的隐式溶剂模型似乎有明显提高精度的潜力

有道理，像DFT-D4那样，按原子电荷来决定原子半径

喵星大佬

wzkchem5 发表于 2021-8-25 16:12
有道理，像DFT-D4那样，按原子电荷来决定原子半径

那估计得来回迭代调整原子半径啥的，而且导数可能还不好办

sobereva

喵星大佬 发表于 2021-8-25 15:50
这么看来要是做一种按照某种规则动态调整原子半径的隐式溶剂模型似乎有明显提高精度的潜力

SCIPCM利用电子密度等值面构造孔洞，在某种程度上也有类似这种adaptive的效果，不过从结果上来看没什么实际优点

喵星大佬

sobereva 发表于 2021-8-26 06:42
SCIPCM利用电子密度等值面构造孔洞，在某种程度上也有类似这种adaptive的效果，不过从结果上来看没什么实 ...

也可能孔洞构造方式主要对非极性部分的影响比较大，可以看看在SCIPCM下面这些体系里面对应原子的孔洞半径和IEFPCM的区别，要是这样的话要结合非极性部分的计算方式才能体现出差别吧

sobereva

喵星大佬 发表于 2021-8-26 09:47
也可能孔洞构造方式主要对非极性部分的影响比较大，可以看看在SCIPCM下面这些体系里面对应原子的孔洞半径 ...

SMD18原文里说的是Coulomb radius，强调的是计算静电部分用的孔洞

苏玖染

Sob老师，我需要计算的体系为铁配合物和IO4-配位时，看哪个自旋态较为稳定，对比不同自旋多重度下在水中的自由能
请问这种情况下，（1）结构优化和计算频率时需要用SMD18隐式溶剂模型么？
（2）我用的时TPSSh/def2-TZVP+D3BJ优化的结构，可以么？
（3）在计算自由能时，高精度真空中单点用的是PWPB95 D3 def2-QZVPPD，然后用优化后的结构用M06-2X结合SMD18溶剂模型计算溶剂化自由能，这样合理么？

sobereva

苏玖染 发表于 2023-7-18 23:26
Sob老师，我需要计算的体系为铁配合物和IO4-配位时，看哪个自旋态较为稳定，对比不同自旋多重度下在水中的 ...

1 不需要特意用。Gaussian默认的IEFPCM就可以，还不容易出虚频
2 可以
3 可以。建议用D4