談談隱式溶劑模型下溶解自由能和體系自由能的計算
量子化學中，怎么在隱式溶劑模型下計算體系的溶解自由能和體系在溶劑下的自由能是頻繁被問及的問題，其實計算方式非常簡單，但很多人還是沒有搞清楚或認識上有誤區，這里就說說如何正確地計算。首先在第1節簡要說一些和隱式溶劑模型有關的常識，這幾個常識不知道的話是不可能正確進行計算的。關于溶劑模型更深入、系統的知識以及各種實際應用，筆者會在北京科音初級量子化學培訓班（http://www.keinsci.com/workshop/KEQC_content.html）和基礎量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）進行詳細的講解。
1 幾個與隱式溶劑模型有關的常識
1.1 什么是隱式溶劑模型
簡單來說，隱式溶劑模型就是不具體描述溶質附近的溶劑分子的具體結構和分布，而是把溶劑環境簡單地當成可極化的連續介質來考慮的。這種考慮溶劑效應的好處是可以表現溶劑的平均效應而不需要像顯式溶劑模型那樣需要考慮各種可能的溶劑層分子的排布方式，而且不至于令計算耗時增加很高，因此被廣泛用于量子化學和分子模擬領域。在分子模擬方面的應用筆者以前專門寫過一個博文有興趣可以看看：《談談分子模擬中的隱式溶劑模型與GB模型》（http://www.shanxitv.org/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，而是其它溶劑模型結合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://www.shanxitv.org/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://www.shanxitv.org/multiwfn）計算的分子范德華表面上靜電勢的統計特征，代入到前人擬合的經驗公式里，對于有機分子可以直接方便地計算出精度滿意的水中的溶解自由能，原理和操作步驟詳見《使用Multiwfn預測晶體密度、蒸發焓、沸點、溶解自由能等性質》（http://www.shanxitv.org/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://www.shanxitv.org/593）。
3 溶劑環境下溶質自由能的計算方法
溶劑環境下溶質的自由能就是溶質在氣相下的自由能加上溶解自由能。為了避免歧義，這里說得更明確、嚴謹一些：溶質在溶劑環境下標況（298.15K、1M）時的自由能 = 溶質在1atm氣相下的自由能 + 直接通過隱式溶劑模型計算的溶解自由能（始末都是1M） + 1.89kcal/mol（標況下1atm→1M自由能變）
有幾點要注意：
(1)默認設定下，通過freq關鍵詞或者熱力學組合方法得到的氣相下的自由能就是標況（298.15K、1atm）狀態下的。
(2)用什么方式、什么級別算氣相下的自由能，跟用什么級別算溶解自由能沒有絲毫關系！這兩部分都分別盡量算準了，溶劑下溶質的自由能才能算準。
(3)對絕大多數中性且局部不顯離子性的分子，氣相和溶劑下結構差異甚微，此時結構優化及之后的任務用氣相下優化的結構就夠了。但如果溶劑效應對當前溶質結構影響可能較大，或者想追求穩妥，結構優化、振動分析獲得熱力學校正量都應當在隱式溶劑模型下進行，并且之后算溶解自由能和高精度單點能的時候也用這個結構。
怎么計算標況下氣相下自由能，初學者總是沒完沒了犯錯，故有必要再提一下。計算方式分為以下三種情況，結果準確度依次提升。注意這里B3LYP/6-31G*僅用來泛指中低計算級別，最適合用什么中低計算級別視具體體系而定，參看《談談量子化學中基組的選擇》（http://www.shanxitv.org/336）和《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/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://www.shanxitv.org/552）。若不懂頻率校正因子相關知識的話，仔細閱讀《談談諧振頻率校正因子》（http://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/413）和《談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的》（http://www.shanxitv.org/557）。
此例在1.gjf中可以不加溶劑模型，因為環氧乙烷是中性、沒有局部顯離子性的體系，因此溶劑模型對opt freq的影響微乎其微。但為了使得1.gjf的關鍵詞更有普適性，也避免之后某些審稿人吐槽為什么優化不在溶劑模型下進行，所以就帶了溶劑模型，反正也不會令耗時增加太多。
我很建議在opt freq過程中用Gaussian默認的IEFPCM而非SMD，也即把1.gjf里scrf中的SMD去掉。因為SMD對于某些體系，特別是柔性體系，容易造成優化難收斂，以及出現小虛頻，這在《Gaussian中幾何優化收斂后Freq時出現NO或虛頻的原因和解決方法》（http://www.shanxitv.org/278）中專門說了。幾何優化用的溶劑模型和計算溶解自由能使用不同的隱式溶劑模型在原理上沒有任何沖突之處。
上面在計算時候為了圖省事，沒有考慮頻率校正因子，也沒有利用Shermo程序基于更好的準RRHO模型計算熱力學量，如果追求更準確的話這兩點皆應當考慮，參見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552），也即使用前述的“講究的人”的計算流程來完成。順帶一提，Shermo可以特別方便地給你不同溫度下的熱力學數據。