谈谈隐式溶剂模型下溶解自由能和体系自由能的计算

量子化学中，怎么在隐式溶剂模型下计算体系的溶解自由能和体系在溶剂下的自由能是频繁被问及的问题，其实计算方式非常简单，但很多人还是没有搞清楚或认识上有误区，这里就说说如何正确地计算。首先在第1节简要说一些和隐式溶剂模型有关的常识，这几个常识不知道的话是不可能正确进行计算的。关于溶剂模型更深入、系统的知识以及各种实际应用如pKa和氧化还原势的计算等，笔者会在北京科音中级量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）里非常全面具体地讲解，欢迎参加。

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等体系其它各种性质的计算结果。

体系有的性质受隐式溶剂模型影响大，如激发能、HOMO-LUMO gap、偶极矩、原子电荷；有的受影响小，如几何结构、NMR、振动频率。但这里说的只是多数情况，具体还是看实际体系。另外，溶质和溶剂的极性越大，静电相互作用就越强，溶剂效应越明显。

关于溶剂模型对结构的影响值得特别说一下。对于绝大多数中性分子，而且局部也不显离子性的情况，是否考虑隐式溶剂模型对几何结构以及振动频率的影响很小，用不用溶剂模型优化出来的结构的差异肉眼上甚至都看不出来。但是对于局部带很显著电荷的情况，不考虑溶剂模型优化出的结构可能与实际溶剂环境下其结构相差甚远、结果定性错误。如果你不清楚结构优化（以及随后的振动分析）时该不该加隐式溶剂模型，建议你一律加上，反正耗时也就增加一点（对于普通泛函做DFT而言），还可以避免对某些体系得到定性错误结果的风险，也避免一些的审稿人的吐槽。对几何结构的影响基本上都是溶剂的极性部分所致，是否考虑非极性部分对结构优化和振动分析总是很小的。

还值得一提的是，有些反应，在溶剂下和气相下反应路径上的势能面差异极大。比如Menschutkin反应，在溶剂环境下才有过渡态，在气相下连过渡态都没有，这是因为这个反应的两个产物都显离子性，溶剂环境对二者稳定化程度十分显著。所以，研究溶液下化学反应过程，如果你心里没数，所有计算都应当加上隐式溶剂模型，必要时候再同时考虑显式溶剂模型。
 

1.4 Gaussian中用什么隐式溶剂模型合适？

Gaussian里用隐式溶剂模型就写scrf关键词即可，详见手册。此时输出的单点能可以认为就是考虑了溶剂效应后的单点能（这么说并不确切，但姑且这么理解问题倒也不大）。

Gaussian09/16默认用的是曾经最主流的PCM模型（具体来说用的是IEFPCM第二版形式），但直接使用的话，默认并不计算非极性部分。要计算非极性部分的话，得在scrf里面写read，并且坐标末尾空一行写
Dis
Rep
Cav
这代表将构成溶剂的非极性部分的值都计算出来并纳入到输出的单点能里。

Gaussian09/16支持的SMD (Solvation Model Based on Density)是目前几乎最好的隐式溶剂模型。之所以好，关键在于其非极性部分的函数形式好，而且参数化过程做得较好，因此算溶解自由能的精度是所有溶剂模型里拔尖的。如果没有特殊情况，建议一律使用SMD，写scrf=SMD即可使用，此时给出的体系能量已经同时包含了溶剂效应的极性和非极性部分的贡献了，而且非极性部分的贡献还会同时单独输出出来。但是，在优化和振动分析过程中使用SMD时，个别情况下可能造成本来是高对称性的结构出现虚频、优化收敛变得困难等问题，而默认的IEFPCM则没这个问题。因此使用SMD优化和振动分析时碰到这种情况时可以尝试改用默认的IEFPCM再试，对于几何优化和振动分析来说SMD和IEFPCM在结果上几乎不会有什么差异，因为对这类任务产生影响的主要是极性部分。鉴于此，我个人经验和建议：做优化、振动分析时总是用IEFPCM，而算单点能的时候则改用SMD。
 

1.5 如何在Gaussian中设置溶剂？

1.5.1 使用内置溶剂

scrf里写solvent=xxx即可，xxx是溶剂名字。比如如果想用SMD溶剂模型表现乙醇溶剂，就写scrf(SMD,solvent=ethanol)。支持的溶剂在Gaussian手册scrf关键词部分的末尾有罗列，不要随意自行发挥去写溶剂名。有的溶剂支持多种写法，可以挑简单、易记的形式写，比如DiMethylSulfoxide可以简写为DMSO。一些溶剂名的等价写法笔者列在了此处：http://bbs.keinsci.com/thread-10624-1-1.html。

1.5.2 自定义溶剂及混合溶剂（只考虑极性部分贡献时）

如果想用的溶剂在Gaussian自带的列表里面没有，而且不需要考虑非极性部分的话，可以在scrf里写read，然后输入文件末尾空一行对溶剂的介电常数进行定义，比如
eps=23.0
epsinf=3.3
这就代表将溶剂静态和动态介电常数分别定义为了23.0和3.3，这两个参数是隐式溶剂模型最最最重要的两个参数，溶剂效应的极性部分的表现只取决于这两个。（可能大家以前也看过一些网上的Gaussian的例子里面定义了溶剂半径之类，但其实对于Gaussian09/16那是完全多余的，因为G09/16默认用的定义溶质孔洞的方式是将原子球叠加得到，原子半径用的是UFF力场的vdW半径乘以1.1。所以溶质孔洞根本不依赖于溶剂半径的定义，除非你自己专门让Gaussian使用溶剂可及表面方式定义溶质的孔洞）

如果只是做一般基态计算，只需要定义eps就够了。如果做TDDFT计算、含频(超)极化率计算之类，还得同时定义epsinf，这会影响结果。而振动分析、NMR等计算也得定义epsinf，但数值可以随意设，不影响结果。epsinf一般查不到，可以用折射率的平方估算，折射率比较好查。如果连折射率也没有，那么对于非极性溶剂，可以假定epsinf的数值等于eps；对于极性溶剂，epsinf可以姑且设1.9（各种溶剂的epsinf都在这附近）。

对于混合溶剂，应该将几种溶剂的eps和epsinf按照体积比例进行混合来定义自定义溶剂。

1.5.3 自定义溶剂及混合溶剂（同时考虑极性和非极性部分贡献时）

在SMD溶剂模型下计算时如果也要考虑溶剂的非极性部分，那么必须把溶剂的完整的SMD参数全都进行定义，即写scrf(read,SMD,solvent=generic)，然后在输入文件末尾空一行写比如
eps=11.5
epsinf=2.0449
HBondAcidity=0.229
HBondBasicity=0.265
SurfaceTensionAtInterface=61.24
CarbonAromaticity=0.12
ElectronegativeHalogenicity=0.24

其中SurfaceTensionAtInterface是表面张力，单位是cal/mol/A^2。CarbonAromaticity是芳香度，是芳环上的碳原子数占所有非氢原子数的比例，ElectronegativeHalogenicity是卤素度，是卤原子占所有非氢原子数的比例。HBondAcidity和HBondBasicity分别是Abraham酸度和碱度，程序内置溶剂的这俩参数来自于Abraham的文章，对于一种新的溶剂这是没法查到的。所以，简单来说，对于一种全新溶剂，是没有严格办法考虑非极性部分参数的。在笔者来看，碰上这种情况时的相对最佳的做法是写比如scrf(read,SMD,solvent=AAA)，然后照常定义eps和epsinf，这里AAA是与你当前要用的溶剂特征最接近的自带的溶剂，这样非极性部分的参数就会借用AAA溶剂的凑合着。虽然这样不严格，但总比完全不考虑非极性部分强一点。

通过自定义SMD溶剂模型参数做SMD计算的实例见《通过SMD溶剂模型描述离子液体溶剂环境的方法》（http://sobereva.com/431）。

如果是几种内置溶剂混合的情况，就把SMD溶剂参数相应地按照体积比混合来定义新溶剂即可。内置溶剂的SMD参数在明尼苏达溶剂描述符数据库http://comp.chem.umn.edu/solvation/mnsddb.pdf里可以查到。

注意水在SMD中是特殊的，没有直接对应的非极性参数，有人问如果溶剂是水和其它溶剂分子混合构成的怎么去定义混合溶剂。在我来看没有绝对理想的办法去定义这种混合溶剂。对于SMD算这种混合溶剂下的溶解自由能，我认为一个可以接受的做法是在水中和在另一个纯溶剂中分别计算溶解自由能，然后按照水和其它溶剂分子的体积比来权重得到混合溶剂下的溶解自由能。

如果你不是用的SMD，而是其它溶剂模型结合dis cav rep选项来考虑溶剂非极性部分的贡献，那更没法严格地自定义新溶剂或者定义混合溶剂。

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，这点是要严厉批评的！

1.7 绝对不要在Gaussian里用CPCM或COSMO！

COSMO是PCM的近似形式，它先把溶剂假设是介电常数无穷大的导体以简化计算过程，然后再近似修正回成实际介电常数溶剂的情况。COSMO写进Gaussian后就被叫做CPCM。至于Dmol3之流只支持COSMO的原因，千万不要以为是因为COSMO好，而是因为这些程序开发者水平不足或者懒，才只实现了比PCM形式更简单的COSMO。

不知道一些人是从哪里学来的坏毛病，在Gaussian里老是想用CPCM、贴出来的关键词里带着CPCM，笔者每当看见这种情况就很不爽。CPCM本身是PCM的近似，而IEFPCM又是PCM的最佳实现形式。显然放着Gaussian里是默认的而且是更好的IEFPCM不用，而刻意去用CPCM，完全就是自取其辱！！！另外，CPCM在Gaussian里对于中性体系用的参数是错的，这点在J. Chem. Theory Comput., 11, 4220 (2015)里专门提到了，结果一定程度上不可靠。所以绝对不要在Gaussian里用CPCM！！！

  

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*。即便是阴离子体系，也仍然要用6-31G*而非带弥散函数的6-31+G*。至于在计算溶解自由能之前的优化过程，则爱用什么级别就用什么。（顺带一提，SMD提出来的时候M06-2X刚提出来，所以SMD文中里没考虑如今流行的M06-2X的情况。鉴于M06-2X和M05-2X有很大相似性而且HF成份几乎一样，所以笔者认为用M06-2X/6-31G*结合SMD算溶解自由能也没问题，而且在JPCB,115,14556(2011)中Truhlar他们也确实用M06-2X来算的。

常有人问过渡金属配合物应该在什么级别下算溶解自由能。SMD原文在对不同级别时，只考虑了有机体系，但我建议依然用M05-2X/6-31G*算配合物的溶解自由能，原因如下。虽然M05-2X算含过渡金属类的体系一般不怎么样，毕竟此泛函是给主族元素参数化的，但由于在一般的过渡金属配合物中有机配体原子数远远多于过渡金属，而且分子暴露在溶剂的区域也主要是配体的原子，过渡金属原子的考虑相对次要，所以结合SMD算溶解自由能时仍应当用M05-2X/6-31G*。如果6-31G*对当前算的过渡金属没有定义的话，用其它恰当的赝势基组也可以，诸如SDD、Lanl2TZ(f)等，见《谈谈赝势基组的选用》（http://sobereva.com/373）。

鉴于老有初学者（不知为什么）老是犯糊涂，在这里再郑重强调一下：M05-2X/6-31G*这个很low的级别，除了结合SMD算溶解自由能很适合外，没有任何其它实际用处！

(3)应当在什么结构下计算溶解自由能？

前面提到，绝大部分中性分子的结构在气相下和溶剂环境下是没有多大区别的，所以计算溶解自由能之前，一般在气相下优化就够了。也因为这个原因，诸如SMD等模型在拟合溶剂模型参数的时候用的也都是分子在气相下的结构，所以在一定程度上，用气相优化的结构比用溶剂下优化的结构更适合计算溶解自由能。不过，如前所述，对于某些情况气相和溶剂下结构差异甚大，甚至不用溶剂模型都优化不出某些结构的情况，还是应当在溶剂模型下优化结构再计算溶解自由能。如果你搞不清楚，就索性优化时候总是带着溶剂模型。

顺带一提，曾经在一篇名为Comment on the Correct Use of Continuum Solvent Models（JPCA,114,13442）的文章中，包括COSMO模型的提出者Klamt在内的三个人鼓吹计算溶解自由能的时候必须使用气相下的结构，然后次年被Cramer和Truhlar在JPCB,115,14556(2011)中打脸，文中的观点就是我上面写的。

(4)怎么得到不同温度下的溶解自由能？

没戏。隐式溶剂模型的参数是冲着标况下的实验溶解自由能来拟合的。

(5)怎么得到溶解焓？

没戏。原因同上。要怪就怪隐式溶剂模型的提出者没有冲着溶解焓来拟合参数。溶解熵什么的同样也得不到。

(6)有没有必要考虑气相和溶剂下平动和转动熵的差异？

气相下和溶剂下体系的平动、转动特征会有所不同，因为溶剂下分子没法整体自由运动。有一些文章鼓吹必须考虑额外的校正项体现平动和转动熵受溶剂环境的影响，不同的人提出了不同的校正方法。但这是其实是纯属多余的，在J. Phys. Chem. A, 122, 1392 (2018)中通过细致的研究，有信服力地证明考虑这些额外校正不仅没对结果有什么改进，反倒往往更糟糕。

(7)关于基于分子表面静电势直接预测溶解自由能

分子的溶解自由能与体系表面静电势有非常密切的关系。通过Multiwfn程序（http://sobereva.com/multiwfn）计算的分子范德华表面上静电势的统计特征，代入到前人拟合的经验公式里，对于有机分子可以直接方便地计算出精度满意的水中的溶解自由能，原理和操作步骤详见《使用Multiwfn预测晶体密度、蒸发焓、沸点、溶解自由能等性质》（http://sobereva.com/337）。

(8)关于单原子离子的溶解自由能的计算

对于不是很小的离子，直接按照上述做法靠SMD模型去算溶解自由能虽然误差比中性体系大得多，但往往还算能接受。如果要算比如质子、Li离子、Mg离子等的溶解自由能，那是绝对不能只靠SMD模型来算的（即直接靠溶剂模型下与气相下的单点能求差来得到）。因为这样的体系电荷密度非常高，与溶剂分子结合非常强烈，比如质子在水中会形成(H3O)+，显然这样的强相互作用不可能靠溶剂模型定性合理地体现。质子在水中的溶解自由能已经有现成的准确值，完全不需要自己算，见Tissandier等人的J. Phys. Chem. A, 102, 7787 (1998)，以及后来Cramer等人对其数据的验证文章J. Phys. Chem. B, 110, 16066 (2006)，数值是-265.9 kcal/mol（前后都是1M浓度）。如果想准确地自行计算，可以利用显式+隐式溶剂模型，参考比如Dixon等人的J. Phys. Chem. A, 105, 11534 (2001)，对于其它单原子离子在各种溶剂下的计算也可以效仿这种方式，其中构造离子+显式溶剂的团簇模型这一步可以利用molclus程序（http://www.keinsci.com/research/molclus.html）做团簇搜索。对于计算质子在你要用的溶剂中的溶解自由能时，如果你嫌考虑一堆显式溶剂分子的方式费事，为得到定性正确的结果，最最最起码也得考虑一个溶剂分子，这种做法算质子在各种溶剂中的溶解自由能的例子可以参看比如J. Mol. Struct. (THEOCHEM), 952, 25 (2010)和Comput. Theor. Chem., 1077, 11 (2016)。如果你要算过渡金属离子的溶解自由能，则必须考虑溶剂分子与这个金属离子的实际配位结构，将之作为它在实际溶剂中的状态。

(9)关于uESE溶剂模型

如果你研究的是净电荷为+1或-1的多原子离子，又懒得用显式溶剂模型，那么计算溶解自由能强烈建议用uESE溶剂模型，比用SMD准确得多，详见《比SMD算溶解自由能更好的溶剂模型uESE的使用》（http://sobereva.com/593）。
     
     

3 溶剂环境下溶质自由能的计算方法

溶剂环境下溶质的自由能就是溶质在气相下的自由能加上溶解自由能。为了避免歧义，这里说得更明确、严谨一些：

溶质在溶剂环境下标况（298.15K、1M）时的自由能 = 溶质在1atm气相下的自由能 + 直接通过隐式溶剂模型计算的溶解自由能（始末都是1M） + 1.89kcal/mol（标况下1atm→1M自由能变）

有几点要注意：
(1)默认设定下，通过freq关键词或者热力学组合方法得到的气相下的自由能就是标况（298.15K、1atm）状态下的。
(2)用什么方式、什么级别算气相下的自由能，跟用什么级别算溶解自由能没有丝毫关系！这两部分都分别尽量算准了，溶剂下溶质的自由能才能算准。
(3)对绝大多数中性且局部不显离子性的分子，气相和溶剂下结构差异甚微，此时结构优化及之后的任务用气相下优化的结构就够了。但如果溶剂效应对当前溶质结构影响可能较大，或者想追求稳妥，结构优化、振动分析获得热力学校正量都应当在隐式溶剂模型下进行，并且之后算溶解自由能和高精度单点能的时候也用这个结构。

怎么计算标况下气相下自由能，初学者总是没完没了犯错，故有必要再提一下。计算方式分为以下三种情况，结果准确度依次提升。注意这里B3LYP/6-31G*仅用来泛指中低计算级别，最适合用什么中低计算级别视具体体系而定，参看《谈谈量子化学中基组的选择》（http://sobereva.com/336）和《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）：
(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-D3(BJ)/def2-TZVP级别，关键词是B2PLYPD3/def2TZVP，或更好的比如CCSD(T)/cc-pVTZ）。之后将二者相加。
(3)讲究的人：先用B3LYP/6-31G*优化体系，然后计算高精度单点。然后将Shermo程序的settings.ini里的E=后面填上高精度单点能的数值，sclZPE后面填上B3LYP/6-31G*级别的ZPE校正因子，ilowfreq设2使用Grimme的准RRHO模型，然后启动Shermo并读取输出的自由能即可。对于柔性体系，通过Shermo用这种准RRHO模型算自由能比用“粗人”的做法准确得多！因为Gaussian用的RRHO模型不能恰当考虑低频模式对熵的贡献。不熟悉Shermo程序的话见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）。若不懂频率校正因子相关知识的话，仔细阅读《谈谈谐振频率校正因子》（http://sobereva.com/221）。

如果算的体系本身很小（热力学组合方法能算得动），且溶剂效应对结构影响很小时，吐血推荐直接用热力学组合方法算气相下自由能，不仅超省事而且精度高。比如关键词就写个G4，其它什么都不写，输出文件最后给出的G4 Free Energy=就是要取的气相下自由能的值了。比较值得使用的热力学组合方法的精度&耗时关系是：W1>CBS-APNO>=G4>G4(MP2)>CBS-QB3>G3(MP2B3)。但即便最便宜的G3(MP2B3)，在比如30多核的服务器下计算30个原子以上还是极难算得动的。虽然还有更便宜的热力学组合方法CBS-4M，但精度在现在来看已经过时了，不建议再用。

另外，一种粗略但比较省事的计算溶剂下焓或自由能的方法如下：
(1)用opt freq优化结构同时获得焓或自由能的热校正量（获得焓的热校正量时建议用ZPE校正因子）。溶剂模型是否需要用看情况
(2)在(1)的结构下结合溶剂模型用高级别计算单点
然后将1和2的结果相加，再加上1.89 kcal/mol即可。此种做法步骤少，但相当于是在高级别下计算了溶解自由能，会比在溶剂模型专门参数化的级别下算的溶解自由能精度低。

4 溶剂环境下溶质自由能的计算实例

上面的文字已经把环境下溶质的自由能的计算方式说得很透彻、详细了，但总是有一些初学者还是翻来覆去问怎么算，还各种误读、歪曲上面文字所传递出的信息。于是笔者干脆在这里给出一个具体例子，是环氧乙烷在乙醇中的自由能计算，其中气相自由能的计算用的是上文中的“粗人”的计算流程（对柔性体系，强烈建议用“讲究的人”的做法）。这个例子很有代表性，把这个例子看明白了，算其它体系也肯定不会遇到什么障碍。

首先下载文件包：http://sobereva.com/attach/327/file.rar。其中有4个Gaussian输入文件，按顺序将1.gjf、2.gjf、3.gjf、4.gjf都计算完之后，就有了获得环氧乙烷在乙醇中自由能所需的全部数据了。输出的.out文件在文件包里也提供了。

这四个文件的内容挨个说一下：
(1)关键词是B3LYP/6-311G* em=GD3BJ opt freq scrf(SMD,solvent=ethanol)，即对应SMD溶剂模型表现乙醇环境的情况下用B3LYP-D3(BJ)/6-311G*对体系进行优化和振动分析。这个级别不昂贵，而用于opt freq目的足够给出可靠的结果。
(2)关键词是B2PLYPD3/def2TZVP geom=allcheck，即使用较好的B2PLYP-D3(BJ)/def2-TZVP级别基于(1)的结构计算真空下较高精度电子能量
(3)关键词是M052X/6-31G* geom=allcheck，即基于(1)的结构用M05-2X/6-31G*在真空下计算电子能量
(4)关键词是M052X/6-31G* scrf(SMD,solvent=ethanol) geom=allcheck，即基于(1)的结构用M05-2X/6-31G*在SMD溶剂模型表现乙醇环境的情况下计算电子能量

我们先计算环氧乙烷在气相下、标况下（298.15 K、1 atm）的自由能。在1.out中搜索Gibbs可找到下面这行
 Thermal correction to Gibbs Free Energy=         0.033857
即自由能的热校正量为0.033857 Hartree。然后在2.out中读取电子能量，结果是-153.7254582 Hartree，因此此体系的气相标况自由能为-153.725458+0.033857=-153.691601 Hartree。大家读双杂化泛函算的电子能量的时候千万别读错地方，不懂的话看《谈谈该从Gaussian输出文件中的什么地方读电子能量》（http://sobereva.com/488）。
注：虽然1.gjf是带着溶剂模型算的，但对于中性且没有局部显离子性的体系，带不带隐式溶剂模型对结构和自由能的热校正量影响甚微，因此上面的Thermal correction to Gibbs Free Energy可以直接当成是气相的自由能热校正量。之所以这里刻意带了溶剂模型，是因为这一节试图给读者一套普适性的关键词、可以应用于各种体系。对于整体或局部显离子性体系，不带溶剂模型优化出的结构与溶剂环境下的实际状态可能定性不符而无法接受，因此这套模板关键词索性就总是带着溶剂模型。

接下来计算环氧乙烷在乙醇中的溶解自由能。从3.out中读取真空下的M05-2X/6-31G*级别的电子能量，为-153.758808 Hartree。然后在4.out中读取溶剂模型下的M05-2X/6-31G*级别的电子能量，为-153.765311 Hartree，因此前后都是1M标准态的溶解自由能为-153.765311-(-153.758808)=-0.006503 Hartree，折合-4.08 kcal/mol。

最终，1 M标准态的环氧乙烷在乙醇中的自由能就是1 atm标准态气相下的自由能加上溶剂模型算的溶解自由能，再加上1.89 kcal/mol，即-153.691601-0.006503+1.89/627.51=-153.695092 Hartree。

下面是一些零碎的讨论，想多了解一些建议看。

大家算其它体系的时候，可以将上面这四个文件作为模板，直接把自己的体系套进去。上面的计算用的计算级别比较有普适性，至少对于不含第四周期以后元素的主族体系肯定是妥当的，计算弱相互作用体系如分子团簇也完全没问题。但注意Gaussian算大体系的双杂化泛函巨慢，还特别耗硬盘，可以改用ORCA来算，省时得多，这里专门提到了《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）。如果只能用Gaussian，当体系超过50原子的话，建议用M06-2X代替B2PLYP-D3(BJ)，会便宜得多，而精度通常下降不多（对主族体系而言）。

如果你想计算其它温度下的溶质在溶剂下的自由能，只需在1.gjf里加上比如temperature=320关键词即可，之后按上面步骤，最终算出来的溶剂下溶质的自由能就对应320K的。但这里相当于假定溶解自由能受温度的影响可忽略不计。

对于环氧乙烷这个例子，1.gjf里其实不需要加D3色散校正，因为色散校正对这个体系的opt freq任务的结果影响几乎为0。但考虑到这些文件会被一些读者当做模板，所以加了D3校正来提升B3LYP的普适性。而B2PLYP是必须加D3的，即便算的不是弱相互作用，对热力学量计算方面也有不可忽视的改进。更多相关介绍看《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）和《谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的》（http://sobereva.com/557）。

此例在1.gjf中可以不加溶剂模型，因为环氧乙烷是中性、没有局部显离子性的体系，因此溶剂模型对opt freq的影响微乎其微。但为了使得1.gjf的关键词更有普适性，也避免之后某些审稿人吐槽为什么优化不在溶剂模型下进行，所以就带了溶剂模型，反正也不会令耗时增加太多。

我很建议在opt freq过程中用Gaussian默认的IEFPCM而非SMD，也即把1.gjf里scrf中的SMD去掉。因为SMD对于某些体系，特别是柔性体系，容易造成优化难收敛，以及出现小虚频，这在《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（http://sobereva.com/278）中专门说了。几何优化用的溶剂模型和计算溶解自由能使用不同的隐式溶剂模型在原理上没有任何冲突之处。

上面在计算时候为了图省事，没有考虑频率校正因子，也没有利用Shermo程序基于更好的准RRHO模型计算热力学量，如果追求更准确的话这两点皆应当考虑，参见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552），也即使用前述的“讲究的人”的计算流程来完成。顺带一提，Shermo可以特别方便地给你不同温度下的热力学数据。