计算化学公社

标题: 比SMD算溶解自由能更好的溶剂模型uESE的使用 [打印本页]

作者
Author: sobereva 时间: 2021-4-30 08:29
标题: 比SMD算溶解自由能更好的溶剂模型uESE的使用

比SMD算溶解自由能更好的溶剂模型uESE的使用

Use of uESE, a solvent model that is better than SMD for calculating solvation free energy

文/Sobereva@北京科音 2021-Apr-30

1 uESE溶剂模型简介

SMD是如今非常常用的溶剂模型，在《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）里笔者专门介绍过。在Gaussian程序中，SMD的极性部分等同于IEFPCM，而由于其非极性部分考虑较周到，因此溶解自由能计算精度较好。但是SMD算离子体系的溶解自由能误差远大于算中性体系的，通常需要用显式+隐式的杂化溶剂模型才能得到较好结果。

最近Vyboishchikov等人提出了uESE (universal Easy Solvation Energy Evaluation)溶剂模型，原文见J. Comput. Chem. (2021) DOI: 10.1002/jcc.26531。之前还有xESE溶剂模型，和uESE形式很相似，原文见Phys. Chem. Chem. Phys., 22, 14591 (2020)。uESE/xESE都像SMD一样考虑了溶质-溶剂相互作用的极性和非极性部分，但uESE/xESE和SMD有很大不同：
(1)uESE/xESE的极性部分利用的是COSMO模型。通常COSMO是基于电子密度计算的，而是uESE/xESE的COSMO则是基于CM5原子电荷计算的；更具体来说，在求解表面显著电荷时利用使用CM5电荷计算溶质孔洞表面的静电势。
(2)uESE/xESE模型并不纳入SCF迭代过程中，也不影响电子结构，不能计算带溶剂模型时的受力和Hessian，相当于是一个“后”溶剂模型。uESE/xESE仅适合基于CM5电荷算溶解自由能的目的。由于uESE/xESE中的参数是基于气相下B3LYP/def2-TZVP级别的结构和波函数计算的CM5原子电荷拟合的，所以自己用uESE/xESE的时候CM5电荷也得在这个级别下计算。
(3)uESE/xESE包含对溶剂-溶质极性作用部分的经验校正项。
(4)uESE/xESE对于不同类型溶剂使用不同的极性部分的校正项和非极性部分形式。溶剂分为四类：水，非水极性质子溶剂，极性非质子溶剂，非极性溶剂。
(5)uESE/xESE拟合的经验参数中既有依赖于元素的，也有依赖于溶剂类型的，也有依赖于溶剂本身的。目前xESE只支持水，uESE支持约100种溶剂。极性部分校正项和非极性部分只有H,C,N,O,F,P,S,Cl,Br,I元素的参数（对于其它元素只能给出COSMO极性部分的结果）。而SMD不含依赖原子半径以外元素的参数。

uESE和xESE有什么存在意义？主要意义在于算溶解自由能从统计结果上比SMD更好，特别是对于离子体系。根据uESE原文的测试，溶解自由能计算精度（平均绝对误差，MAE）有下面的关系
中性体系：xESE≈SM12>=uESE>=SMD
阳离子体系：uESE>xESE>SMD>=SM12
阴离子体系：uESE>xESE≈SMD12>>SMD
即曰，算中性体系溶解自由能，用xESE会比SMD占一点优势。算阳离子用uESE比用SMD强得多。算阴离子更是绝对不能直接用SMD，不想用杂化溶剂模型的话至少也应当用uESE。根据测试，uESE算中性体系溶解自由能的MAE小于1 kcal/mol，算阴、阳离子溶解自由能的MAE差不多，都小于3 kcal/mol。而SMD算阴离子的MAE则达到8 kcal/mol左右。注：上面的SM12和SMD是同门的，但前者和uESE一样也是基于CM5原子电荷进行计算，不过极性部分利用的是广义Born方程形式。

下图左侧是uESE原文里算大量中性和离子在非水质子溶剂中的溶解自由能（纵坐标）与实验值（横坐标）的对比，右图是SMD算的情况，误差超过4 kcal/mol的用红色符号高亮。可见uESE算中性溶质的情况比SMD稍好，而计算离子的情况能好很多。

(, 下载次数 Times of downloads: 68)

下图是计算在极性非质子溶剂中的溶解自由能的情况。可见此时SMD算离子的情况表现极差，特别是对于溶解自由能非常大的离子

(, 下载次数 Times of downloads: 54)

需要注意的是以上都是统计数据，大家也不要指望算每个中性体系xESE都比SMD好、算每个离子体系uESE都一定比SMD好。笔者写此文的目的主要是提醒大家算离子体系的溶解自由能的时候，如果懒得用杂化溶剂模型，则至少要用uESE，切勿拿SMD凑合（除非有uESE不支持的元素）。

2 uESE和xESE程序的使用

uESE和xESE方法有同名的计算程序，可以在http://iqcc.udg.edu/~vybo/ESE/直接下载，Windows和Linux版可执行文件都提供了，不开源。这两个程序用的输入文件相同，uESE可以指定溶剂，而xESE只支持水溶剂。
2023-Sep-27注：我发现此时以上链接已无法访问。笔者写此文时用的uESE和xESE程序可以在http://sobereva.com/attach/593/uESE,xESE.zip下载。

uESE和xESE都需要基于命令行使用。uESE在Windows下的使用格式是：uESE.exe [输入文件路径] -solvent [溶剂名]。可以用的溶剂见https://github.com/vyboishchikov/ESE/blob/main/solvent-list.md。

uESE/xESE的输入文件里需要提供各原子的元素在周期表里的序号、原子坐标，以及原子的CM5电荷。笔者开发的Multiwfn波函数分析程序（主页&免费下载：http://sobereva.com/multiwfn）直接就能计算CM5电荷。为了令用户能尽可能方便地产生uESE/xESE的输入文件，从2021-Apr-30更新的Multiwfn开始只要将settings.ini里的uESEinp参数设为1，则Multiwfn计算完CM5电荷后就会问你是否产生uESE/xESE的输入文件。使用uESE/xESE算溶解自由能用的输入文件若是Multiwfn产生的，请注意在你发表的文章中提及，并按照Multiwfn启动时的提示恰当引用Multiwfn。

Multiwfn计算CM5电荷需要提供含有波函数的信息作为输入文件，诸如.wfn/.wfx/.fch/.molden/.mwfn等等，产生方式见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。不了解Multiwfn的话参看《Multiwfn FAQ》（http://sobereva.com/452）。

3 使用uESE计算PhO-在乙腈中的溶解自由能实例

下面就通过一个例子演示uESE程序的使用，是PhO-（苯酚的羟基的质子解离掉后的阴离子）在乙腈中的溶解自由能计算。此例所有相关文件可以在http://sobereva.com/attach/593/file.rar下载。本例Gaussian使用G16 A.03版，Multiwfn是2021-Apr-30更新的3.8(dev)版，操作系统是Win 10 64bit。

首先用Gaussian对PhO-进行优化。根据uESE原文，应当用B3LYP/def2-TZVP在真空下进行。输入文件如下

%chk=C:\PhO-_optfreq.chk
#p B3LYP/def2TZVP opt freq
[空行]
Title Card Required
[空行]
-1 1
C                0.00000000 0.00000000 -1.80531390
C                0.00000000 -1.20821014 -1.10754617
C                0.00000000 -1.20820376 0.28716582
C                0.00000000 0.00000000 0.98475086
C                -0.00000000 1.20820376 0.28716582
C                -0.00000000 1.20821014 -1.10754617
H                0.00000000 0.00000000 -2.90492418
H                0.00000000 -2.16044049 -1.65754504
H                0.00000000 -2.16042940 0.83722286
H                -0.00000000 2.16042940 0.83722286
H                -0.00000000 2.16044049 -1.65754504
O                0.00000000 0.00000000 2.41475076

计算完成后，检查输出文件确认没有虚频。然后用formchk将PhO-_optfreq.chk转换为PhO-_optfreq.fch。formchk的使用在《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）里说了。

将Multiwfn目录下的settings.ini里的uESEinp设为1，然后启动Multiwfn，载入PhO-_optfreq.fch，之后依次输入
7  //布居分析与原子电荷计算
16  //计算CM5电荷
1  //使用内置的球对称化的自由原子密度
[按回车]  //产生的uESE输入文件为当前目录下的PhO-_optfreq.inp
然后就可以把Multiwfn窗口关了。

下载uESE的Windows版可执行文件uESE.exe到某处，把PhO-_optfreq.inp放到与此文件相同的目录下，然后按Win+R键，输入cmd进入命令行界面，通过cd命令进入此目录下（这种计算机基本常识不会的话自行Google），运行uESE PhO-_optfreq.inp -solvent acetonitrile。仅需一秒钟就能运行完，在输出信息中会看到
COSMO electrostatic energy =    -66.946 kcal/mol
Correction term =       12.393 kcal/mol
Total solvation free energy =    -54.554 kcal/mol

即曰溶解自由能是-54.554 kcal/mol。uESE/xESE给出的溶解前后都是1M浓度的溶解自由能。从SMD溶剂模型原文的补充材料里可以查到实验值是-55.1 kcal/mol（也是溶解前后都是1M的情况），可见此例uESE表现得极为理想，误差仅有0.55 kcal/mol，不过这很大程度上是运气好。

下面也用SMD模型算一下溶解自由能，做法在《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）里专门讲过。基于上面B3LYP/def2-TZVP气相下优化的结构，在M05-2X/6-31G*级别下算一次气相单点能和SMD模型表现的乙腈环境下的单点能，结果为
# M052X/6-31G* scrf(SMD,solvent=acetonitrile)：-306.941786445 Hartree
# M052X/6-31G*：-306.847792107 Hartree
因此溶解自由能为627.51*(-306.941786445+306.847792107)=-58.98 kcal/mol。相对于实验的误差达到-3.9 kcal/mol，明显大于uESE的，这体现出了uESE的显著优势。

注意uESE模型在测试、拟合参数的时候对离子体系都是只考虑了带+1、-1电荷的离子。高价离子与溶剂相互作用显著强于一价离子，笔者估计用uESE肯定不会得到较好结果，必须用杂化溶剂模型。

作者
Author: abdoman 时间: 2021-4-30 12:31
PhO-感觉大概率是在训练集里面。
不知道其他化合物的效果如何，值得期待。

作者
Author: wzkchem5 时间: 2021-4-30 16:20
如果这个溶剂模型对SCF迭代没有影响的话，大概率对高价阴离子不靠谱吧，有的高价阴离子在真空中会自发电离，但是在水之类的强极性溶剂里还比较稳定，这种体系用uESE肯定定性错误

作者
Author: sobereva 时间: 2021-4-30 17:02

wzkchem5 发表于 2021-4-30 16:20
如果这个溶剂模型对SCF迭代没有影响的话，大概率对高价阴离子不靠谱吧，有的高价阴离子在真空中会自发电离 ...

此模型只考虑+1和-1电荷溶质的情况，明尼苏达溶解自由能数据库也只包含这种情况
高价离子与溶剂分子作用很强，不用显式溶剂模型肯定没戏

作者
Author: winnerwill 时间: 2021-4-30 17:05

wzkchem5 发表于 2021-4-30 16:20
如果这个溶剂模型对SCF迭代没有影响的话，大概率对高价阴离子不靠谱吧，有的高价阴离子在真空中会自发电离 ...

这种情况是不是应该在SMD溶剂模型几何优化后，再在气相下B3LYP/def2-TZVP级别算个单点，再用按本文方法用uESE溶剂模型算溶解自由能？

作者
Author: wuzhiyi 时间: 2021-4-30 17:06
想问一下
1. 是用B3LYP/def2TZVP做几何优化，而不是B3LYP-D3(BJ)/def2TZVP?
2. 这个溶剂模型和COSMO-RS系列比，哪个强啊？

作者
Author: sobereva 时间: 2021-4-30 18:58

wuzhiyi 发表于 2021-4-30 17:06
想问一下
1. 是用B3LYP/def2TZVP做几何优化，而不是B3LYP-D3(BJ)/def2TZVP?
2. 这个溶剂模型和COSMO-RS ...

1 作者没加D3。而且考虑的都是小分子，加不加D3影响可忽略不计
2 没测试过。算中性的应该半斤八两，算离子应该这个强

作者
Author: wuzhiyi 时间: 2021-5-3 05:38
想问一下，如果对关键区域做杂化模型，SMD和uESE/xESE哪个比较好啊？

感觉uESE/xESE的话，应该就是在真空中优化水和感兴趣的分子，然后再算原子电荷，然后再用uESE/xESE计算溶解自由能。

作者
Author: sobereva 时间: 2021-5-3 07:02

wuzhiyi 发表于 2021-5-3 05:38
想问一下，如果对关键区域做杂化模型，SMD和uESE/xESE哪个比较好啊？

感觉uESE/xESE的话，应该就是在真 ...

如果是考虑显式溶剂的话，那就没必要用uESE/xESE了，而且大概率不会比SMD好，毕竟经验性更强。

作者
Author: wuzhiyi 时间: 2021-5-3 22:25
感谢sob老师的讲解，我下载 Version 3.8(dev), release date: 2021-Apr-30版的Multiwfn
读取chk文件
然后
7 // Population analysis and atomic charges
16 // CM5 atomic charge
1 // Use build-in sphericalized atomic densities in free-states (more convenient)

这之后得到的是

Final atomic charges, after normalization to actual number of electrons
Atom 1(O ): -0.61679616
Atom 2(C ):    0.04381511
Atom 3(C ): -0.27180443
Atom 4(C ): -0.27182572
Atom 5(C ): -0.27182720
Atom 6(H ):    0.04837694
Atom 7(H ):    0.03270410
Atom 8(H ):    0.04839308
Atom 9(H ):    0.04839674
Atom 10(H ):    0.04838862
Atom 11(H ):    0.03270885
Atom 12(H ):    0.04839246
Atom 13(H ):    0.03269543
Atom 14(H ):    0.04838218

Calculation took up    6 seconds wall clock time

If outputting atom coordinates with charges to Deprotonated.chg in current folder? (y/n)

y
Result have been saved to Deprotonated.chg in current folder
Columns 1 to 5 are name,X,Y,Z,charge respectively, unit is Angstrom

而不是

[按回车] //产生的uESE输入文件为当前目录下的PhO-_optfreq.inp

作者
Author: sobereva 时间: 2021-5-3 22:33

wuzhiyi 发表于 2021-5-3 22:25
感谢sob老师的讲解，我下载 Version 3.8(dev), release date: 2021-Apr-30版的Multiwfn
读取chk文件
然后 ...

uESEinp设成1了么？启动时是否提示settings.ini没找到？

作者
Author: wuzhiyi 时间: 2021-5-3 23:01

sobereva 发表于 2021-5-3 22:33
uESEinp设成1了么？启动时是否提示settings.ini没找到？

非常抱歉，忘记看的改setting这一句话了。改完之后就没有问题了。

作者
Author: yangch17 时间: 2021-5-13 11:28
请问一下，用这个uESE模型怎么进行温度和压力校正啊

作者
Author: sobereva 时间: 2021-5-13 19:16

yangch17 发表于 2021-5-13 11:28
请问一下，用这个uESE模型怎么进行温度和压力校正啊

只能算标况的。其它常见溶剂模型也都如此

作者
Author: yuclark 时间: 2021-5-23 11:41
是否适用离子液体溶剂环境

作者
Author: sobereva 时间: 2021-5-23 16:47

yuclark 发表于 2021-5-23 11:41
是否适用离子液体溶剂环境

不。
就算此模型或许适合，类似下文，但至少作者目前没有专门考虑过这种情况
通过SMD溶剂模型描述离子液体溶剂环境的方法
http://sobereva.com/431（http://bbs.keinsci.com/thread-10607-1-1.html）

作者
Author: yuclark 时间: 2021-6-19 16:08

sobereva 发表于 2021-5-23 16:47
不。
就算此模型或许适合，类似下文，但至少作者目前没有专门考虑过这种情况
通过SMD溶剂模型描述离子 ...

多谢！

作者
Author: fatpig 时间: 2021-6-23 16:57
用uESE测试了几个含过渡金属(Fe/Co)的中性分子及一价阳离子，对照体系是和实验值比较接近的SMD(M05-2X/6-31g*)数据。

对中性分子，得出的结果基本和SMD一致，可以用。

但对一价阳离子，uESE在某些solvent出现0.5eV级别的误差，但在某些solvent又比较准，所以这部分计算可能是有问题的，最好不要使用在含金属离子上。

作者
Author: sobereva 时间: 2021-6-23 16:59

fatpig 发表于 2021-6-23 16:57
用uESE测试了几个含过渡金属(Fe/Co)的中性分子及一价阳离子，对照体系是和实验值比较接近的SMD(M05-2X/6-31 ...

作者的训练集里面没有考虑过渡金属的

作者
Author: fatpig 时间: 2021-6-23 19:39

sobereva 发表于 2021-6-23 16:59
作者的训练集里面没有考虑过渡金属的

对，所以对比的对象是SMD，SMD的训练好像也没有考虑过渡金属吧。

根据测试，最好的还是直接用M05-2X/6-31g*，使用M05-2X/6-31g*/SDD或M05-2X/def2-TZVP的效果都不如M05-2X/6-31g*。

虽然样本不多，但看起来是SMD对中性和正一价的含金属化合物都勉强够用，uESE则部分正一价数据偏差较大，原因不明。

作者
Author: sobereva 时间: 2021-6-24 01:08

fatpig 发表于 2021-6-23 19:39
对，所以对比的对象是SMD，SMD的训练好像也没有考虑过渡金属吧。

根据测试，最好的还是直接用M05-2X/6 ...

uESE经验性比SMD高得多，对于与训练集明显不同的体系，SMD更好是情理之中

作者
Author: Xian 时间: 2021-7-19 21:03
Sob老师，请问一下，如果是想优化-2价的配体分子的结构，不需要算溶解自由能，那么可以用这个uESE溶剂模型吗？谢谢。

作者
Author: wzkchem5 时间: 2021-7-19 21:10

Xian 发表于 2021-7-19 14:03
Sob老师，请问一下，如果是想优化-2价的配体分子的结构，不需要算溶解自由能，那么可以用这个uESE溶剂模型 ...

不行，因为uESE没有做梯度。
对于结构优化PCM一般就够了

作者
Author: Xian 时间: 2021-7-19 22:00

wzkchem5 发表于 2021-7-19 21:10
不行，因为uESE没有做梯度。
对于结构优化PCM一般就够了

好的呢，那我去比较一下真空和PCM条件下优化的结构，谢谢。

作者
Author: sobereva 时间: 2021-7-19 22:11

Xian 发表于 2021-7-19 21:03
Sob老师，请问一下，如果是想优化-2价的配体分子的结构，不需要算溶解自由能，那么可以用这个uESE溶剂模型 ...

(2)uESE/xESE模型并不纳入SCF迭代过程中，也不影响电子结构，不能计算带溶剂模型时的受力和Hessian，相当于是一个“后”溶剂模型。uESE/xESE仅适合基于CM5电荷算溶解自由能的目的。由于uESE/xESE中的参数是基于气相下B3LYP/def2-TZVP级别的结构和波函数计算的CM5原子电荷拟合的，所以自己用uESE/xESE的时候CM5电荷也得在这个级别下计算。

作者
Author: Xian 时间: 2021-7-19 22:22

sobereva 发表于 2021-7-19 22:11

谢谢老师，我再去看看帖子，加深印象

作者
Author: slope 时间: 2021-7-20 18:07
老师，请问我可以在M06-2X/6-31g(d,p)下优化结构，再在B3LYP/def2-TZVP下进行单点能计算，再用这个结构算uESE算溶解自由能吗？

作者
Author: sobereva 时间: 2021-7-20 18:50

slope 发表于 2021-7-20 18:07
老师，请问我可以在M06-2X/6-31g(d,p)下优化结构，再在B3LYP/def2-TZVP下进行单点能计算，再用这个结构算uE ...

可

作者
Author: lanthanum 时间: 2022-12-11 18:51
请问，这个uESE或者SM12溶剂模型，有没有办法外调，用于高斯带溶剂的opt ? 谢谢。

作者
Author: chands 时间: 2022-12-11 19:19

lanthanum 发表于 2022-12-11 18:51
请问，这个uESE或者SM12溶剂模型，有没有办法外调，用于高斯带溶剂的opt ? 谢谢。

没办法，见23楼。

作者
Author: linqiaosong 时间: 2022-12-12 01:55
对于cp2k，是不是也可以用B3LYP/def2-TZVP优化以后产生molden文件，再用Multiwfn算CM5后产生uESE的输入文件

作者
Author: sobereva 时间: 2022-12-14 18:14

linqiaosong 发表于 2022-12-12 01:55
对于cp2k，是不是也可以用B3LYP/def2-TZVP优化以后产生molden文件，再用Multiwfn算CM5后产生uESE的输入文件

可以是可以
不过uESE只针对分子体系，用CP2K计算效率很低

作者
Author: BUCTtyl 时间: 2022-12-30 17:22
请问一下直接用uESE计算Cl-的溶解自由能，但是输出文件没报错也没计算，是什么原因？

作者
Author: sobereva 时间: 2022-12-31 00:48

BUCTtyl 发表于 2022-12-30 17:22
请问一下直接用uESE计算Cl-的溶解自由能，但是输出文件没报错也没计算，是什么原因？

具体问开发者
我只能说算这个没意义，单原子离子的溶解自由能不依靠显式溶剂模型肯定算不准

作者
Author: Accelerator 时间: 2023-9-26 18:07
现在文中提到的下载链接、以及原作者在其发表文章中提到的.es结尾的链接都打不开了，是否有人可以分享一份？与此同时我也在向原作者询问，如有结果也将分享给大家。

作者
Author: sobereva 时间: 2023-9-27 04:05

Accelerator 发表于 2023-9-26 18:07
现在文中提到的下载链接、以及原作者在其发表文章中提到的.es结尾的链接都打不开了，是否有人可以分享一份 ...

http://sobereva.com/attach/593/uESE,xESE.zip
这是我写帖子时候的程序版本

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