转 溶剂化
2017-01-04 18:35阅读：
谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
文/Sobereva(3)
First release: 2016-May-24 Last
update: 2016-Oct-26
量子化学中，怎么在隐式溶剂模型下计算体系的溶解自由能和体系在溶剂下的自由能是频繁被问及的问题，其实计算方式非常简单，但很多人还是没有搞清楚或认识上有误区，这里就说说如何正确地计算。首先在第1节先简要说几个和隐式溶剂模型有关的常识，这几个常识不知道的话是不可能正确进行计算的。
1 几个与隐式溶剂模型有关的常识
1.1 什么是隐式溶剂模型
简单来说，隐式溶剂模型就是不具体描述溶质附近的溶剂分子的具体结构和分布，而是把溶剂环境简单地当成可极化的连续介质来考虑的。这种考虑溶剂效应的好处是可以表现溶剂的平均效应而不需要像显式溶剂模型那样需要考虑各种可能的溶剂层分子的排布方式，而且不至于令计算耗时增加很高，因此被广泛用于量子化学和分子模拟领域。在分子模拟方面的应用笔者以前专门写过一个博文有兴趣可以看看：《谈谈分子模拟中的隐式溶剂模型与GB模型》（http://sobereva.com/42）。在量化领域，隐式溶剂模型包括Onsager、PCM、CPCM、IPCM、SCIPCM、COSMO、SMD、SMx系列（如SM12）等等。隐式溶剂模型的缺点是无法表现溶剂与溶质之间的强相互作用，如氢键。而且对于溶质是离子的情况，隐式溶剂模型计算溶解自由能的精度明显低于溶质是中性分子的情况。
1.2 溶剂的极性和非极性部分
对于隐
式溶剂模型，溶剂效应可分为极性和非极性部分。极性部分体现出溶剂-溶质间的静电相互作用，也包括溶剂对溶质电子分布的极化，这是隐式溶剂模型的主体，不合理考虑这部分结果会定性错误。非极性部分相对次要，体现出溶质-溶剂间的各种非静电作用，成份较复杂，包括溶质挪进溶剂构成孔洞时排挤开溶剂所需要做的功、溶质对溶剂分子熵效应的影响、溶质-溶剂分子间的交换互斥和色散吸引作用。只有把非极性部分也考虑时溶解自由能才能定量计算准确。
Onsager、PCM、CPCM、IPCM、SCIPCM、COSMO等模型只明确定义了溶剂的极性部分怎么算，但是非极性部分的计算方式并没有在溶剂模型中明确给出，算的方式视具体程序而定。而SMD、SMx系列则明确定义了怎么计算非极性部分，并且在不同计算级别下专门拟合了参数，只要恰当使用，整体精度肯定是优于PCM、CPCM、COSMO之类的。
1.3 隐式溶剂模型对体系性质的影响
隐式溶剂模型会改变体系的势能面，因此，和势能面直接相关的，比如单点能、极小点和过渡态结构、IRC、振动频率、不同构象分布比例、激发能等也都会受到影响。隐式溶剂模型还会影响体系电子结构，所以gap、偶极矩、键级等等性质也一律受到影响。隐式溶剂模型还会间接地影响NMR等体系其它各种性质的计算结果。
体系有的性质受隐式溶剂模型影响大，如激发能、gap；有的受影响小，如几何结构、NMR、振动频率。但这里说的只是多数情况，具体还是看实际体系。对于绝大多数中性分子，是否考虑隐式溶剂模型对几何结构以及振动频率的影响很小。另外，溶质和溶剂的极性越大，相互作用就越强，溶剂效应越明显。
1.4 Gaussian中用什么隐式溶剂模型合适？
Gaussian里用隐式溶剂模型就写scrf关键词即可，详见手册。此时输出的单点能可以认为就是考虑了溶剂效应后的单点能（这么说并不确切，但姑且这么理解问题倒也不大）。
Gaussian09默认用的是曾经最主流的PCM模型（具体来说用的是IEFPCM形式），但直接使用的话，默认并不计算非极性部分。要计算非极性部分的话，得在scrf里面写read，并且坐标末尾空一行写
Dis
Rep
Cav
这代表将构成溶剂的非极性部分的值都计算出来并纳入到输出的单点能里。
Gaussian09支持的SMD是目前几乎最好的隐式溶剂模型，笔者建议若无特殊情况一律用SMD模型，即写scrf=SMD。此时给出的单点能已经同时包含了极性和非极性部分的贡献了。
1.5 如何在Gaussian中选择溶剂？
scrf里写solvent=xxx即可，xxx是溶剂名字，Gaussian手册scrf关键词部分有罗列。
如果想用的溶剂在Gaussian自带的列表里面没有，若溶剂模型不是用的SMD的话，可以自定义溶剂，在scrf里写solvent=generic,read，然后输入文件末尾空一行对溶剂信息进行定义，比如
eps=23.0
epsinf=3.3
这就代表将溶剂静态和动态介电常数分别定义为了23.0和3.3，这两个参数是隐式溶剂模型最最最重要的两个参数，溶剂效应的极性部分的表现只取决于这两个。（可能大家以前也看过一些网上的Gaussian的例子里面定义了溶剂半径之类，但其实对于Gaussian09那是完全多余的，因为G09默认用的定义溶质孔洞的方式是将原子球叠加得到，原子半径用的是UFF力场的vdW半径乘以1.1。所以溶质孔洞根本不依赖于溶剂半径的定义，除非你自己专门让Gaussian使用溶剂可及表面方式定义溶质的孔洞）
如果只是做基态计算，只需要定义eps就够了。如果做TDDFT计算、含频极化率计算之类，还得同时定义epsinf。振动分析、NMR等计算也得定义epsinf，但数值可以随意设，不影响结果。
对于混合溶剂，应该将几种溶剂分子的eps和epsinf按照体积比例进行混合来定义自定义溶剂。
SMD模型不适合自定义溶剂。因为此模型计算非极性部分的公式依赖于溶剂的多种物性参数，光定义eps、epsinf是不够的，而SMD依赖的溶剂的一些其它属性（如酸度）又没法往Gaussian里定义。
1.6 关于外迭代
溶剂的极性部分对溶质的影响以反应场来描述，反应场本身也决定于溶质的电荷分布。在隐式溶剂模型下做HF计算时，溶剂效应通过向Fock算符里添加与反应场相应的外势算符来表现。在SCF迭代收敛时，不光溶质分子的能量和波函数的变化相对上一轮迭代已经非常小，溶质电子密度的分布与反应场之间也达成了自洽，此时溶质的波函数就是被溶剂充分极化了的波函数，而且溶剂环境也被溶质电荷分布所充分极化，这正是自洽反应场(SCRF)的含义。
大家知道，后HF计算时，是先做SCF计算得到收敛的HF波函数后，再做电子相关计算得到相关能，二者加和就是后HF能量。后HF方法结合隐式溶剂模型下进行计算，实际上溶剂只是对溶质的HF密度充分进行了响应，而并没有充分对当前的后HF密度进行响应。如果想把溶剂效应描述得更准确些，就要让溶剂对后HF密度充分响应，需要基于后HF密度构建反应场并纳入体系哈密顿，通过迭代让溶剂的反应场与后HF密度间达到自洽，这叫外迭代(External
iteration)。外迭代过程往往得做几轮或十几轮，每一轮都相当于做一次HF自洽场计算+计算后HF密度的耗时，所以耗时很高。
对于DFT来说，做TDDFT某种意义上也相当于后HF计算，因此也涉及到外迭代，这点专门在此文专门说了：《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）。
用外迭代的话在scrf里写externaliteration即可。显然，外迭代仅对于后HF、TDDFT等计算有意义，而对于本身就是SCF计算的过程，包括半经验、HF、DFT、MCSCF，此时写上externaliteration没有丝毫意义，反应场都已经和溶质密度自洽了还做外迭代想迭代什么？很多人，都没搞懂什么叫外迭代，就胡乱使用externaliteration，这点是要严厉批评的！
2 溶解自由能的计算方法
溶解自由能(solvation free
energy)也叫溶剂化自由能，计算方法非常简单：在溶剂模型下计算的单点能减去气相下计算的单点能
例如，计算某分子在乙醇中的溶解自由能，就用比如# M052X/6-31G*
scrf(SMD,solvent=ethanol)计算得到的单点能，减去#
M052X/6-31G*计算得到的单点能。这种计算溶解自由能的方式是绝对合理的，那些专门搞溶剂模型的人计算溶解自由能也是用这种方式计算，并将得到的结果与实验结果相对比来拟合参数。如果对这点还心存疑虑的话，一定仔细去看J.
Phys. Chem. A, 114, 13442
(2010)，包括其中补充材料部分，文章专门说了这个问题。
有几个细节问题需要专门说一下：
(1)标准态问题
化学上，气相标准态的浓度是1atm（常温下相当于1/24.46
mol/L），液相标准态的浓度是1M(即1mol/L)。教科书以及实验上对溶解自由能的定义是这个过程的自由能变：把1atm气相状态的溶质分子挪入到溶剂中成为1M溶剂化状态。然而，按照上面的方式计算的溶解自由能则对应的是从1M气相状态挪入到溶剂中成为1M溶剂化状态过程的自由能变。也就是说，隐式溶剂模型算的溶解自由能对应的初态浓度和一般化学上定义的溶解自由能的初态的浓度不同。因此，使用隐式溶剂模型得到的溶解自由能还要加上对应于气态下1atm->1M浓度改变造成的自由能变，然后才能和一般化学上说的溶解自由能数值相对比。
在298.15K下，气体1atm->1M浓度改变对应的自由能变为1.89kcal/mol。
(2)计算级别的选择
隐式溶剂模型在优化啊、振动分析啊、光谱计算啊等等目的时候用什么级别都行，但是算溶解自由能的时候，一定一定一定要用溶剂模型参数化的级别。比如PCM溶剂模型可以结合不同原子半径来定义溶质的孔洞，UAHF半径是专门在HF/6-31G*级别（对阴离子来说是6-31+G*）参数化的，也就是说，开发者调节UAHF半径的目的是让HF/6-31G*级别通过PCM计算的溶解自由能和实验测定值最相符。因此，哪怕用更昂贵的理论方法或基组，比如在CCSD/aug-cc-pVTZ级别下结合PCM/UAHF计算的溶解自由能，也一定会比在HF/6-31G*下计算的更差！！！对于SMD模型，原文对几种级别分别进行了参数化，最终发现M05-2X/6-31G*级别下溶解自由能算得最好。所以，在SMD模型下计算溶解自由能，一定要用M05-2X/6-31G*。而至于在计算溶解自由能之前的优化过程，则爱用什么级别就用什么。（顺带一提，SMD提出来的时候M06-2X刚提出来，所以SMD文中里没考虑如今流行的M06-2X的情况。鉴于M06-2X和M05-2X有很大相似性而且HF成份几乎一样，所以笔者认为用M06-2X/6-31G*结合SMD算溶解自由能也没问题，而且在JPCB,115,14556(2011)中Truhlar他们也确实用M06-2X来算的。
(3)应当在什么结构下计算溶解自由能？
前面提到，绝大部分中性分子的结构在气相下和溶剂环境下是没有多大区别的，所以计算溶解自由能之前，一般在气相下优化就够了。也因为这个原因，诸如SMD等模型在拟合溶剂模型参数的时候用的也都是分子在气相下的结构，所以在一定程度上，用气相优化的结构比用溶剂下优化的结构更适合计算溶解自由能。不过，对于某些情况气相和溶剂下结构差异甚大，甚至不用溶剂模型都优化不出某些结构的情况，还是应当在溶剂模型下优化结构再计算溶解自由能。
顺带一提，曾经在一篇名为Comment on the Correct Use of Continuum Solvent
Models（JPCA,114,13442）的文章中，包括COSMO模型的提出者Klamt在内的三个人鼓吹计算溶解自由能的时候必须使用气相下的结构，然后次年被Cramer和Truhlar在JPCB,115,14556(2011)中打脸，文中的观点就是我上面写的。
(4)怎么得到不同温度下的溶解自由能？
没戏。隐式溶剂模型的参数是冲着标况下的实验溶解自由能来拟合的。
(5)怎么得到溶解焓？
没戏。原因同上。要怪就怪隐式溶剂模型的提出者没有冲着溶解焓来拟合参数。溶解熵什么的同样也得不到。
3 溶剂环境下溶质自由能的计算方法
溶剂环境下溶质的自由能就是溶质在气相下的自由能加上溶解自由能。为了避免歧义，这里说得更明确、严谨一些：
溶质在溶剂环境下标况（298.15K、1M）时的自由能 =
溶质在1atm气相下的自由能 + 直接通过隐式溶剂模型计算的溶解自由能（始末都是1M） +
1.89kcal/mol（标况下1atm->1M自由能变）
有几点要注意：
(1)默认设定下，通过freq关键词或者热力学组合方法得到的气相下的自由能就是标况（298.15K、1atm）状态下的。
(2)用什么方式、什么级别算气相下的自由能，跟用什么级别算溶解自由能没有丝毫关系！这两部分都分别尽量算准了，溶剂下溶质的自由能才能算准。
(3)鉴于对绝大多数中性分子气相和溶剂下结构差异甚微，所以算溶解自由能和做振动分析用的结构就都用气相下优化的结构即可。
怎么计算标况下气相下自由能，初学者总是没完没了犯错，故有必要再提一下。计算方式分为以下三种情况，结果准确度依次提升（这里B3LYP/6-31G*仅用来泛指中低计算级别，最适合用什么中低计算级别视具体体系而定）：
(1)巨懒无敌懒人：# B3LYP/6-31G* opt freq，读Sum of electronic and
thermal Free Energies=。
(2)庸人：# B3LYP/6-31G* opt freq，读取吉布斯热校正量（Thermal correction to
Gibbs Free Energy=）。然后在此结构下计算高精度单点（如#
B2PLYP/def2TZVP）。二者相加。
(3)讲究的人：先用B3LYP/6-31G*优化体系，然后在考虑B3LYP/6-31G*级别下ZPE、焓升温、熵的频率校正因子的情况下计算出吉布斯热校正量，加到高精度单点上。如果不知道具体怎么搞，仔仔仔细细细阅读此文：《谈谈谐振频率校正因子》（http://sobereva.com/221）。
注：对于(2)、(3)情况，当极个别时候气相和溶剂下结构差异大，计算气相自由能时应当先在溶剂模型下优化结构，然后带着溶剂模型做振动分析得到吉布斯热校正量，并加到气相高精度单点下作为气相自由能。当然，之后计算溶解自由能的时候也是用溶剂下优化的结构。
如果算的体系本身很小（热力学组合方法能算得动），且溶剂效应对结构影响很小时，吐血推荐直接用热力学组合方法算气相下自由能，不仅超省事而且精度高。比如关键词就写个G4，其它什么都不写，输出文件最后给出的G4
Free
Energy=就是要取的气相下自由能的值了。比较值得使用的热力学组合方法的精度&耗时关系是：W1>CBS-APNO>=G4>G4(MP2)>CBS-QB3>G3(MP2B3)>>CBS-4M。
另外，一种粗略但比较省事的计算溶剂下焓或自由能的方法如下：
(1)用opt
freq优化结构同时获得焓或自由能的热校正量（获得焓的热校正量时建议用ZPE校正因子）。溶剂模型是否需要用看情况
(2)在(1)的结构下结合溶剂模型用高级别计算单点
然后将1和2的结果相加即可。此种做法步骤少，但相当于是在高级别下计算了溶解自由能，会比在溶剂模型专门参数化的级别下算的溶解自由能精度低。
http://ccc.keinsci.com/forum.php?mod=viewthread&tid=3345&highlight≠̸��ʽ�ܼ�ģ��
