通过SMD溶剂模型描述离子液体溶剂环境的方法

sobereva

通过SMD溶剂模型描述离子液体溶剂环境的方法

The method to describe solvent environment of ionic liquid via SMD solvation model

文/Sobereva @北京科音  2018-Aug-7

离子液体已被广泛作为溶剂。常有人问离子液体溶剂环境怎么在量子化学计算中体现，本文就专门说一下。阅读本文之前务必先阅读《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）了解基础知识。


1 原理

Truhlar搞出来的SMD溶剂模型是一种比较普适的隐式溶剂模型。原理上，只要你给了溶剂的以下参数，就可以通过SMD溶剂模型描述溶剂环境：
· 静态介电常数eps(298K)
· 动态介电常数epsinf或折射率（前者是后者的平方）
· Abraham氢键酸度：来自Abraham的一些文章
· Abraham氢键碱度：来自Abraham的一些文章
· 气液界面表面张力(298K)
· 芳香度：溶剂分子中芳香碳占重原子（即非氢原子）的比例
· 卤素度：溶剂分子中卤素占重原子（即非氢原子）的比例

对于SMD溶剂模型而言，离子液体作为溶剂并没有任何特殊性，只要提供了离子液体的以上参数，就可以用SMD恰当地描述离子液体溶剂环境。这种做法的可靠性在SMD原作者的一篇文章J. Phys. Chem. B, 116, 9122 (2012)中得到了充分验证，文中测试了一批常见的离子液体体系。从测试结果看，在计算溶解自由能的计算精度上，离子液体和普通分子作为溶剂的情况差不多，尽管对于有的离子液体误差比较大，比如[EMIM][PF6]。

用SMD描述离子液体的一个困难之处是很多离子液体的SMD参数找不全。如果你感兴趣的只是溶剂-溶质之间的极性部分作用还好，只需要定义eps，有些时候也需要定义epsinf（或折射率），这就够了，而其它参数不会被利用；而如果还需要考虑溶剂的非极性部分贡献，那就很麻烦，表面张力可能还不难找，但Abraham氢键酸度和碱度没地方搞去。

为了解决离子液体SMD参数往往找不全的问题，JPCB那篇文章里提出一个模型叫SMD-GIL，其中GIL是generic ionic liquid的缩写。这个模型并不是一个新的溶剂模型的形式，它的意思只不过是把除了芳香度和卤素度以外的SMD参数都直接设为许多离子液体实验上已知的参数的平均值。SMD-GIL的参数值在JPCB文中给了：
静态介电常数=11.5
折射率=1.43（折合动态介电常数=1.43^2=2.0449）
Abraham氢键酸度=0.229
Abraham氢键碱度=0.265
表面张力=61.24 cal/mol/A^2
而离子液体的芳香度和卤素度自己根据其化学组成手动一算就知道。

根据JPCB文中的测试，使用SMD-GIL参数计算和使用实际离子液体的实验参数做SMD计算得到的结果精度半斤八两，互有胜负。因此，通过SMD溶剂模型描述你用的离子液体，有实验值的参数可以直接用实验值，而没有的就用SMD-GIL给出的离子液体的平均参数代替即可，结果是合理的；或者如果你懒得找实验参数，则所有参数都用SMD-GIL的也完全没问题。JPCB文中也给出了一些常见离子液体的实验参数，见文中表1、3。


2 实例：计算乙醇在[BMIM][NTf2]中的溶解自由能

[BMIM][NTf2]的结构如下

乙醇在[BMIM][NTf2]中的常温下的溶解自由能实验值为-3.76 kcal/mol（来自那篇JPCB文章补充材料Table S1.6，对应气相和溶剂下都是1M标准态浓度）。我们看看利用SMD-GIL的参数算出来的结果是多少。如《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》所述，我们先对结构优化，然后用M052X/6-31G*在溶剂模型下和气相下分别算个单点，求差就是常温下的溶解自由能（对应气相和溶剂下都是1M标准态浓度）。

SMD几乎是目前最流行的隐式溶剂模型，因此很多主流程序都已经支持SMD，比如Gaussian、NWChem、ORCA、GAMESS-US等（而非主流的Dmol3和ADF之流尚不支持）。这里我们用量化计算最常用的Gaussian作为例子，结构已经在B3LYP/TZVP下优化过（用什么级别优化无所谓，只要是合理级别即可，对溶解自由能计算结果影响很小）。

单点任务的输入文件如下：
# M052X/6-31G*

B3LYP/TZVP opted

0 1
C                  1.17643900   -0.39950500    0.00000000
H                  1.15236000   -1.03828400    0.88463100
H                  2.11738900    0.15557800    0.00000000
H                  1.15236000   -1.03828400   -0.88463100
C                  0.00000000    0.55417200    0.00000000
H                  0.04041900    1.20079700   -0.88644000
H                  0.04041900    1.20079700    0.88644000
O                 -1.20073500   -0.22291400    0.00000000
H                 -1.95570500    0.37470800    0.00000000

溶剂模型下的单点任务的输入文件如下。这里我们自己完全自定义了一个新溶剂，除了芳香度和卤素度以外的溶剂参数就是上面提到的SMD-GIL的平均离子液体溶剂参数。
# M052X/6-31G* SCRF(SMD,read,solvent=generic)

B3LYP/TZVP opted

0 1
C                  1.17643900   -0.39950500    0.00000000
H                  1.15236000   -1.03828400    0.88463100
H                  2.11738900    0.15557800    0.00000000
H                  1.15236000   -1.03828400   -0.88463100
C                  0.00000000    0.55417200    0.00000000
H                  0.04041900    1.20079700   -0.88644000
H                  0.04041900    1.20079700    0.88644000
O                 -1.20073500   -0.22291400    0.00000000
H                 -1.95570500    0.37470800    0.00000000

eps=11.5
epsinf=2.0449
HBondAcidity=0.229
HBondBasicity=0.265
SurfaceTensionAtInterface=61.24
CarbonAromaticity=0.12
ElectronegativeHalogenicity=0.24

上面文件中设的CarbonAromaticity就是芳香度，ElectronegativeHalogenicity是卤素度。[BMIM][NTf2]的化学组成为C10H15F6N3O4S2，重原子数（非氢原子数）为25，此体系中咪唑环上的碳属于芳香碳，数目为3，因此芳香度为3/25=0.12。体系中卤素有6个，因此6/25=0.24就是卤素度。

对两个文件用Gaussian16 A.03版进行计算，气相单点能为-154.9998251，溶剂下单点能为-155.0061949，因此溶解自由能为627.51*(-155.0061949+154.9998251)=-3.997 kcal/mol，和实验值-3.76 kcal/mol相符不错。

苏玖染

sob老师，我想请问下，如果计算的是有机物钠盐，比如十二烷基苯磺酸钠，计算溶剂自由能的时候也把钠离子考虑进去一块算溶剂下的和真空中的能量么？

sobereva

苏玖染 发表于 2018-8-7 09:05
sob老师，我想请问下，如果计算的是有机物钠盐，比如十二烷基苯磺酸钠，计算溶剂自由能的时候也把钠离子考 ...

是

Dyyyzz

sob老师，请问如果在定义了溶剂模型下还要给def2基组添加弥散函数，使用《给def2以ma-方式加弥散函数的Gaussian格式的基组定义文件（含所有def2支持的元素）》中提到的ma-方式，此时“@C:\ma-SVP.txt”应写在哪里呢？

sobereva

Dyyyzz 发表于 2022-4-19 10:13
sob老师，请问如果在定义了溶剂模型下还要给def2基组添加弥散函数，使用《给def2以ma-方式加弥散函数的Gaus ...

看http://sobereva.com/g09/m_input.htm里关于从末尾读取不同字段的顺序关系，按顺序写，每个字段之间有一个空行

Dyyyzz

sobereva 发表于 2022-4-19 18:17
看http://sobereva.com/g09/m_input.htm里关于从末尾读取不同字段的顺序关系，按顺序写，每个字段之间有 ...

我按照http://sobereva.com/g09/m_input.htm里面说的gen在前，scrf=read在后的顺序写了输入文件并相应添加了空行，但还是不能成功计算，输出文件中提示没有定义eps，请老师指教。
输入文件：
%chk=FeCl3.chk
%nprocshared=20
%mem=5GB
#p m062x/gen opt freq em=gd3 scrf(pcm,read,solvent=generic)

Title Card Required

0 2
Fe                 0.00000000    0.00000000   -0.03726500
Cl                 0.00000000    0.00000000   -1.97726500
Cl                 0.00000000    1.77120634    0.78922668
Cl                 0.00000000   -1.77120634    0.78922668

@/home/admin/txls/ma-TZVP.txt/N

eps=11.5
epsinf=2.0449


输出文件末：
Using the following non-standard input for PCM:
--- end of non-standard input.
Must specify Eps for generic solvent.
Error termination via Lnk1e in /home/admin/software/gaussian/g09/l301.exe at Wed Apr 20 10:34:37 2022.
Job cpu time:       0 days  0 hours  0 minutes 55.3 seconds.
File lengths (MBytes):  RWF=      5 Int=      0 D2E=      0 Chk=      7 Scr=      1

sobereva

Dyyyzz 发表于 2022-4-20 10:41
我按照http://sobereva.com/g09/m_input.htm里面说的gen在前，scrf=read在后的顺序写了输入文件并相应添 ...

1. #p pbe1pbe/gen opt freq em=gd3 scrf(pcm,read,solvent=generic)
   
2. 
   
3. Title Card Required
   
4. 
   
5. 0 2
   
6. Fe                 0.00000000    0.00000000   -0.03726500
   
7. Cl                 0.00000000    0.00000000   -1.97726500
   
8. Cl                 0.00000000    1.77120634    0.78922668
   
9. Cl                 0.00000000   -1.77120634    0.78922668
   
10. 
    
11. Fe
    
12. def2SVP
    
13. ****
    
14. Cl
    
15. def2SVP
    
16. ****
    
17. 
    
18. eps=11.5
    
19. epsinf=2.0449
    
20. 
    

复制代码

没任何问题
注意引用的基组文件末尾是否有空行，多空哪怕一行就没法读入溶剂设置

另外，绝对不要用M06-2X优化过渡金属配合物，此文明确说了
简谈量子化学计算中DFT泛函的选择
http://sobereva.com/272（http://bbs.keinsci.com/thread-536-1-1.html）

MiaoZuo

sobereva 发表于 2022-4-21 06:24
没任何问题
注意引用的基组文件末尾是否有空行，多空哪怕一行就没法读入溶剂设置

三价铁是不是应该用高自选多重度6来计算啊社长？

chands

MiaoZuo 发表于 2022-11-27 10:29
三价铁是不是应该用高自选多重度6来计算啊社长？

你的问题严重歪楼，不应该发到这里。
是自旋多重度不是自选多重度

看情况，如果强场配体可能是低自旋。

MiaoZuo

chands 发表于 2022-11-27 10:54
你的问题严重歪楼，不应该发到这里。
是自旋多重度不是自选多重度

不好意思打错字了，就是在这边看到了才问一下。那是否需要分别计算2和6的单点能才可以判断？

chands

MiaoZuo 发表于 2022-11-27 11:17
不好意思打错字了，就是在这边看到了才问一下。那是否需要分别计算2和6的单点能才可以判断？

如果不确定就算一下。其实很多配合物的电子排布是清楚的，或者有实验磁矩数据，就不用算两个。某些物质的高自旋态和低自旋态的能量相差很小，就两个都要考虑。

MiaoZuo

chands 发表于 2022-11-27 11:24
如果不确定就算一下。其实很多配合物的电子排布是清楚的，或者有实验磁矩数据，就不用算两个。某些物质的 ...

谢谢大佬~~

Lemon97z

老师，我想请问一下这里离子液体的参数能适用于IEFPCM溶剂模型吗？

sobereva

Lemon97z 发表于 2022-11-30 20:14
老师，我想请问一下这里离子液体的参数能适用于IEFPCM溶剂模型吗？

这是两码事。IEFPCM溶剂模型没直接定义非极性部分，自然也不可能用文中的参数

Lemon97z

sobereva 发表于 2022-12-1 05:27
这是两码事。IEFPCM溶剂模型没直接定义非极性部分，自然也不可能用文中的参数

老师那如果只是针对极性部分呢，比如 在IEFPCM中我可以假设离子液体的介电常数为11.5吗？