计算化学公社

标题: 谈谈隐式溶剂模型下溶解自由能和体系自由能的计算 [打印本页]
作者
Author: sobereva    时间: 2016-5-25 02:48
标题: 谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
On the calculation of solvation free energy and system free energy via implicit solvent model

文/Sobereva @北京科音
First release: 2016-May-24   Last update: 2024-Jun-20


量子化学中，怎么在隐式溶剂模型下计算体系的溶解自由能和体系在溶剂下的自由能是频繁被问及的问题，其实计算方式非常简单，但很多人还是没有搞清楚或认识上有误区，这里就说说如何正确地计算。首先在第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）。

(10)关于SMD18溶剂模型
隐式溶剂模型算溶解自由能的准确度和构造溶质孔洞用的原子半径关系极大，SMD的作者后来发现SMD原文对Br和I用的半径不合适，导致对涉及Br和I的反应的自由能计算误差很大，因此对这俩元素提出了新的半径来显著改进结果，这称为SMD18模型。此模型的细节以及在Gaussian中的实现见《显著改进SMD溶剂模型描述Br和I元素精度的溶剂模型SMD18的介绍》（http://sobereva.com/613）。



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可以特别方便地给你不同温度下的热力学数据。


作者
Author: tobeant    时间: 2016-5-25 10:21
拜读！
作者
Author: 陌上芬菲    时间: 2016-5-25 20:03
本帖最后由 陌上芬菲 于 2016-5-25 20:25 编辑

# opt freq M052X/6-31G* scrf(SMD,solvent=acetone) 结果的 Sum of electronic and thermal Free Energies 是不是 “巨懒无敌懒人” 的溶质在溶剂中的自由能？？
作者
Author: sobereva    时间: 2016-5-26 00:34
陌上芬菲 发表于 2016-5-25 20:03
# opt freq M052X/6-31G* scrf(SMD,solvent=acetone) 结果的 Sum of electronic and thermal Free Energies ...

是
作者
Author: lastzealot    时间: 2016-5-28 18:58
拜读了，看来我也是巨懒无敌懒的类型:)
作者
Author: 鸢翔flybird    时间: 2016-6-1 22:42
本帖最后由 鸢翔flybird 于 2016-6-2 00:35 编辑

sobereva老师，本帖的内容适合算带电荷的原子分子吗？比如我试着算了一下质子的气相自由能和在水中的溶解自由能，前者与文献值一致，后者差别太大——在M062X/6-31G(d)级别下分别计算气相和溶剂模型（smd）下的单点能差，并加1.89kcal/mol校正，得到质子在水中溶解自由能为-133.14kcal/mol，而文献中的值为-265.9kcal/mol，大约是真实值的一半。差在哪呢？
作者
Author: sobereva    时间: 2016-6-2 01:43
鸢翔flybird 发表于 2016-6-1 22:42
sobereva老师，本帖的内容适合算带电荷的原子分子吗？比如我试着算了一下质子的气相自由能和在水中的溶解自 ...


你不能直接拿隐式溶剂模型算质子，这种 电荷/体积 比非常小的体系用隐式溶剂模型误差极大。这个问题怎么算看此文：Computational and Theoretical Chemistry 1077 (2016) 11–17，必须得考虑起码一个溶剂分子。

多原子的离子用隐式溶剂模型还凑合，虽然误差会比中性大不少。
作者
Author: 鸢翔flybird    时间: 2016-6-2 01:47
sobereva 发表于 2016-6-2 01:43
你不能直接拿隐式溶剂模型算质子，这种 电荷/体积 比非常小的体系用隐式溶剂模型误差极大。这个问题怎 ...

谢谢Sob老师！
作者
Author: cgchen    时间: 2016-6-14 16:41
Sob 看到这个帖子 我顿时想起了当初专门发邮件请教这块问题 大约是14年的事情了   解释的很清楚 关键问题说的很明白
博士论文致谢唯一遗漏的就是您  
今天补上 谢谢sob老师！



作者
Author: 鸢翔flybird    时间: 2016-6-21 20:36
本帖最后由 鸢翔flybird 于 2016-6-21 21:14 编辑

请教sob老师：
带正或负电荷的有机分子(100个原子以内，含氨基、羟基)，气相下自由能的计算是否有必要考虑弱相互作用，即根据教学贴《谈谈谐振频率校正因子》(http://sobereva.com/221) 由原来的B3LYP/6-31G*改为M062X/6-31G*级别下优化结构并振动分析？

另外看到教学帖《DFT-D色散校正的使用》(http://sobereva.com/210) 和《乱谈DFT-D》(http://sobereva.com/83) 中建议B3LYP-D3(BJ)或M062X-D3计算弱相互作用，而G09 D.01版D3(BJ)计算频率时色散校正的贡献是不完全准确的，这是否意味着对于G09 D.01版计算气相自由能，得用M062X-D3优化结构并振动分析(基组6-31G*)？但此级别下计算比B3LYP/6-31G*慢好多，能否用B3LYP-D3(BJ)优化结构，B3LYP-D3振动分析(基组6-31G*)？(另一个与此类似的问题是，若已在DFT-D3下优化了结构，那么计算热力学校正量即振动分析时是否必须也要加上DFT-D3校正？)

若最终目的是高精度计算气相自由能，可否计算高精度单点能时用D3(BJ)校正，而振动分析这块用零阻尼D3校正？又：可否计算高精度单点能加DFT-D3校正，而振动分析无色散校正？或者，您能给出此处关于DFT-D3校正的较优建议？
作者
Author: sobereva    时间: 2016-6-21 21:38
鸢翔flybird 发表于 2016-6-21 20:36
请教sob老师：
带正或负电荷的有机分子(100个原子以内，含氨基、羟基)，气相下自由能的计算是否有必要考虑 ...


优化用了D3，振动分析也必须用，否则势能面根本不一样，结果无意义
用B3LYP-D3(BJ)优化结构，而B3LYP-D3做振动分析完全说不通，依然不是相同的级别
B3LYP本身能不能用得看具体体系，有个结构图最能说清楚。如果有弱相互作用但都是静电主导的，B3LYP做opt freq并无问题。
没必要优化和单点用不同阻尼，这样显得略诡异，也没明显好处。必须用B3LYP-D3的话，想得到好点结果，用B3LYP-D3结合6-31G*优化，算单点时候用大点基组诸如TZVP。
作者
Author: 鸢翔flybird    时间: 2016-6-21 22:26
本帖最后由 鸢翔flybird 于 2016-6-21 22:34 编辑
sobereva 发表于 2016-6-21 21:38
优化用了D3，振动分析也必须用，否则势能面根本不一样，结果无意义
用B3LYP-D3(BJ)优化结构，而B3LYP- ...

谢谢sob老师指教。不过我没太理解老师的意思，我的体系弱相互作用主要就是氨基、羟基可能与水溶剂形成氢键，所以根据您的建议，(A) 可以用B3LYP不加D3 opt freq，然后双杂化泛函大基组级别下加D3(BJ)校正算高精度单点能，由此得比较精确的气相自由能，是这个意思吗？
(B) 还是B3LYP不加D3 opt freq，然后双杂化泛函大基组级别下不加D3校正算高精度单点能？
(C) 还是B3LYP加D3 opt freq，然后双杂化泛函大基组级别下加D3校正算高精度单点能？
(D)还是B3LYP加D3(BJ) opt freq，然后双杂化泛函大基组级别下加D3(BJ)校正算高精度单点能？（考虑到"G09 D.01版D3(BJ)计算频率时色散校正的贡献是不完全准确的"）
如果答案是A、C、D中的某一个，为了得到在水溶剂中的自由能，后面溶解自由能的计算时单点能的计算是否也要加上相应D校正？
作者
Author: sobereva    时间: 2016-6-21 23:49
鸢翔flybird 发表于 2016-6-21 22:26
谢谢sob老师指教。不过我没太理解老师的意思，我的体系弱相互作用主要就是氨基、羟基可能与水溶剂形成氢 ...


B3LYP-D3 opt freq，然后双杂化泛函大基组级别下加D3(BJ)校正算高精度单点能。
如果发现加了D3几何优化不好收敛，就去掉，对当前体系结构优化影响不会很大。

溶解自由能计算不要加D3，就按此文用M05-2X/6-31G*计算
谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
http://sobereva.com/327
作者
Author: 鸢翔flybird    时间: 2016-6-22 14:24
sobereva 发表于 2016-6-21 23:49
B3LYP-D3 opt freq，然后双杂化泛函大基组级别下加D3(BJ)校正算高精度单点能。
如果发现加了D3几何优 ...

感谢sob老师耐心解答。请问B3LYP-D3 freq (基组为6-31G*)对应的谐振频率校正因子是否能用B3LYP/6-31G*对应的数值？
作者
Author: sobereva    时间: 2016-6-22 14:27
鸢翔flybird 发表于 2016-6-22 14:24
感谢sob老师耐心解答。请问B3LYP-D3 freq (基组为6-31G*)对应的谐振频率校正因子是否能用B3LYP/6-31G*对 ...

可以
作者
Author: dreamyeye    时间: 2016-6-23 22:23
是否可以认为，只要溶剂化自由能是负的，就说明溶质溶剂化导致了体系的稳定呢？还是溶剂化自由能必须要达到一定的大小才能这么说？
作者
Author: sobereva    时间: 2016-6-24 03:49
dreamyeye 发表于 2016-6-23 22:23
是否可以认为，只要溶剂化自由能是负的，就说明溶质溶剂化导致了体系的稳定呢？还是溶剂化自由能必须要达到 ...

只要是负的就说明溶剂能稳定化体系，即把溶质分子放入到溶剂这个过程中整个体系自由能会下降。
作者
Author: dreamyeye    时间: 2016-6-24 06:19
sobereva 发表于 2016-6-24 03:49
只要是负的就说明溶剂能稳定化体系，即把溶质分子放入到溶剂这个过程中整个体系自由能会下降。

谢谢sob老师
作者
Author: 三石草祭    时间: 2016-11-14 10:31
我竟然一直在用巨懒无比之人的算法计算
作者
Author: wugaxp    时间: 2016-11-25 23:55
多谢sobereva的详细介绍，受益良多。
文中提到，从隐式溶剂模型无法得到不同温度下的溶解自由能。那是否有别的不这么直接的方法可以得到这个值呢？

例如，是否可以用这个模型：溶质分子+若干个显式溶剂分子覆盖住溶质，加上隐式溶剂模型，然后在某个温度下跑一段时间的MD，或者更严格地就是从若干个初始态跑MD，取能量的平均作为这个温度下的自由能？对需要的显式溶剂分子个数也可以做个收敛测试
作者
Author: sobereva    时间: 2016-11-26 11:41
wugaxp 发表于 2016-11-25 23:55
多谢sobereva的详细介绍，受益良多。
文中提到，从隐式溶剂模型无法得到不同温度下的溶解自由能。那是否有 ...


准确计算溶解自由能一般用TI、EFP方法靠MD来做
作者
Author: wugaxp    时间: 2016-11-26 23:59
本帖最后由 wugaxp 于 2016-11-27 00:20 编辑
sobereva 发表于 2016-11-26 11:41
准确计算溶解自由能一般用TI、EFP方法靠MD来做

MD的最大的问题是是否存在准确可靠的力场，换个元素都会需要重新拟合和验证。有的时候遇到棘手的系统，例如包含f电子的体系，或者混合了多种成键类型（当然粗糙的做法就是直接用现有的力场分块处理，但是可信度成问题），力场模型本身都会是个大问题。况且TI的计算量也不见得比cluster模型下的量化计算来得省劲。所以我才想是否可以有比较合理且快速的方法去做这件事，或者经过讨论能找到这样的方法。
事实上我担心的是熵的贡献。对于单个分子，也许量化里面的做法还算合理，毕竟分子内的构型基本不会有大的变化，所以误差抵消了许多。但是如果用cluster模型，也许这部分的贡献就要重新考虑。或者有别的方法可以绕开这个的话就更好了，例如直接先用隐式溶剂模型计算不同温度下溶剂分子自身的溶解能量，辅以气态情况下的熵贡献。

EFP似乎是个合理的选择，因为溶剂分子一般都还是有经验势的。但是QM和MM区之间的连接精度不知如何。Sobereva对EFP有了解吗？
作者
Author: sobereva    时间: 2016-11-27 00:38
wugaxp 发表于 2016-11-26 23:59
MD的最大的问题是是否存在准确可靠的力场，换个元素都会需要重新拟合和验证。有的时候遇到棘手的系统，例 ...


前面打错了，EFP我是指FEP，不是GAMESS-US里的EFP
FEP和TI基本没什么区别，都是基于MD做的，用这个不是为了省计算量去的
没有合适的力场显然没法用TI或FEP。我说用TI/FEP + MD是对于有机小分子体系，这种方式算自由能极为成熟。对于配合物，糙点就直接用SMD，准确点就用杂化溶剂模型并考虑多种溶剂层排布。
作者
Author: wugaxp    时间: 2016-11-27 22:34
sobereva 发表于 2016-11-27 00:38
前面打错了，EFP我是指FEP，不是GAMESS-US里的EFP
FEP和TI基本没什么区别，都是基于MD做的，用这个不 ...

如果是FEP那我就明白了。
力场总是个问题，而且TI和FEP的计算量都太大，如果需要快速筛选一大批分子的情况下很不实际。以前用过TI，和Flory–Huggins这种快速模型比起来，定型差别不大，计算量则完全不可比，大太多了，准确度则太依赖于力场了。
另外如果分子表面不是球形，不知道怎么处理溶剂层排布？sobereva了解这方面的工作吗？我之前说用abinitio-MD，就是想避开这个。另外"杂化"溶剂模型是什么？极化？

我去看了一下EFP，感觉如果可以和TI/FEP之类的连起来，也许能给出比较精确的预言。基本思想在那里呢，只要能把QM和MM区的耦合强度（vdW,charge-transfer,Electrostatic）连续变化，就能用FEP的思想了。问题是似乎EFP的限制挺多，计算量也不小。
作者
Author: sobereva    时间: 2016-11-28 01:12
wugaxp 发表于 2016-11-27 22:34
如果是FEP那我就明白了。
力场总是个问题，而且TI和FEP的计算量都太大，如果需要快速筛选一大批分子的情 ...

杂化溶剂模型就是显式溶剂模型和隐式溶剂模型一起用，显式溶剂只考虑第一配位层
没有要求必须是球形，初始结构可以靠经典动力学模拟产生再量化优化，也可以用一些团簇结构生成工具来产生

如果大批量快速筛选，还是直接用SMD为宜，对于中性体系，尤其是有机体系，误差并不大，计算量也很小。
作者
Author: wugaxp    时间: 2016-11-29 11:19
本帖最后由 wugaxp 于 2016-11-29 11:20 编辑
sobereva 发表于 2016-11-28 01:12
杂化溶剂模型就是显式溶剂模型和隐式溶剂模型一起用，显式溶剂只考虑第一配位层
没有要求必须是球形，初 ...

谢谢sobereva的解释。

不过还是回到那个问题了，不是常温的话SMD就麻烦了。
而且上面也说了，力场不一定总是有的，所以经典分子动力学不一定可以用。我最近遇到的体系就是含有f电子的，而且不只是要考虑常温。

还有，你说的中性体系SMD误差不大，但是如果是有比较强的偶极矩呢？误差能有多大？
作者
Author: sobereva    时间: 2016-11-29 20:30
wugaxp 发表于 2016-11-29 11:19
谢谢sobereva的解释。

不过还是回到那个问题了，不是常温的话SMD就麻烦了。

偶极矩大没问题，不如说偶极矩大反倒是好事
主要是SMD算带电体系误差大
作者
Author: 1130240115    时间: 2016-11-29 21:17
请问老师，我计算50个原子的体系（主要为C,H,O)吸附阴离子的吸附能，考虑溶剂化效应，优化使用B3LYP-D3/6-311G**，单点使用B3LYP-D3/6-311+G**（均考虑溶剂化效应），溶解自由能使用您推荐的M05-2X/6-31G*计算，这样可以吗？谢谢老师。
作者
Author: sobereva    时间: 2016-11-29 23:17
1130240115 发表于 2016-11-29 21:17
请问老师，我计算50个原子的体系（主要为C,H,O)吸附阴离子的吸附能，考虑溶剂化效应，优化使用B3LYP-D3/6-3 ...

可以
作者
Author: 1130240115    时间: 2016-11-29 23:38
谢谢老师回复。
作者
Author: paradise    时间: 2016-12-19 21:25
请问老师金属存在的情况下，应该是m052x/genecp .
此时的金属应该用高精度的么？其他元素还是6-31g*?
作者
Author: sobereva    时间: 2016-12-20 01:32
paradise 发表于 2016-12-19 21:25
请问老师金属存在的情况下，应该是m052x/genecp .
此时的金属应该用高精度的么？其他元素还是6-31g*?

能用6-31G*就用6-31G*，用不了的时候用赝势基组
作者
Author: MrD    时间: 2017-5-9 13:21
本帖最后由 MrD 于 2017-5-9 14:08 编辑

请教社长两个问题：

1.  文中提到"如果算的体系本身很小（热力学组合方法能算得动），且溶剂效应对结构影响很小时，吐血推荐直接用热力学组合方法算气相下自由能，不仅超省事而且精度高。比如关键词就写个G4，其它什么都不写，输出文件最后给出的G4 Free Energy=就是要取的气相下自由能的值了。"

请问，可以B3LYP/6-31G*(或其它中等级别)优化结构，然后G4计算单点能，得出气相下的自由能么 ？
因为这样就免去了计算气相高精度单点能，以及计算频率得出自由能校正，这两项麻烦。
但不了解G4方法，自由能校正是用频率计算得出的，G4是不需要计算频率的么？


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

请问，如果对于反应 A + B -- C 求自由能变，那么反应物天然就比产物多了1.89kcal/mol的自由能，是这样么？因为在求ABC每一项的液相自由能时候，都涉及到加1.89.

作者
Author: sobereva    时间: 2017-5-9 19:57
MrD 发表于 2017-5-9 13:21
请教社长两个问题：

1.  文中提到"如果算的体系本身很小（热力学组合方法能算得动），且溶剂效应对结构 ...

1 不能。G4根本就不是像一般理论方法那样算单点能的，本身已经包含了几何优化过程。至少要读一篇热力学组合方法的文章再去用G4，例如Gn系列的综述：DOI: 10.1002/wcms.59

2 是
作者
Author: MrD    时间: 2017-5-9 21:37
sobereva 发表于 2017-5-9 19:57
1 不能。G4根本就不是像一般理论方法那样算单点能的，本身已经包含了几何优化过程。至少要读一篇热力学组 ...

感谢社长的回答和提供的文献。这就去学习。
作者
Author: 研道    时间: 2018-5-3 15:09
感谢分享
作者
Author: ZHANGZY    时间: 2018-8-5 23:12
sobereva：
“对于混合溶剂，应该将几种溶剂分子的eps和epsinf按照体积比例进行混合来定义自定义溶剂。”

有没有相关文献，这么做？我找了半天，没有发现已发表的SCI文章，有这方面的说明，能帮忙提供一下相关文献吗？
作者
Author: sobereva    时间: 2018-8-6 05:51
ZHANGZY 发表于 2018-8-5 23:12
sobereva：
“对于混合溶剂，应该将几种溶剂分子的eps和epsinf按照体积比例进行混合来定义自定义溶剂。”
...

我一时半会儿在硬盘里找不到现成的，这么做原理上很正当，不会有问题
作者
Author: ZHANGZY    时间: 2018-8-6 06:43
我遇到一个审稿人意见：
(4) How do you define the DMSO:wáter 1:1 mixture in Gaussian?. G09 code does not include this solvent combination by default. As the g09 manual remarks: "We list the <epsilon> values here for convenience, but be aware it is only one of many internal parameters used to define solvents. Thus, simply changing the <epsilon> value will not define a new solvent properly." (http://wild.life.nctu.edu.tw/~jsyu/compchem/g09/g09ur/k_scrf.htm).

我是按照sobereva老师的帖子，定义了EpsInf 和 Eps，没有定义radii。我觉得没什么问题，但是要回复审稿人意见，需要举例子，所以要找已发表的SCI文献。
作者
Author: sobereva    时间: 2018-8-6 08:19
ZHANGZY 发表于 2018-8-6 06:43
我遇到一个审稿人意见：
(4) How do you define the DMSO:wáter 1:1 mixture in Gaussian?. G09 code doe ...


根本不需要定义radii

如果你只考虑溶剂模型的极性部分，定义eps（有时也要定义epsinf）就足够了，不依赖于别的参数；如果你还要同时考虑非极性部分的参数，没法严格定义，因为算这部分还需要溶剂的一大堆其它参数，一般都不可能找全。看本帖1.5节
作者
Author: xinren    时间: 2018-11-29 21:34
最近在计算激发态的电荷传递，有考虑到溶剂的影响，主要是质子性溶剂与极性非质子性溶剂。在计算非质子性溶剂随极性的增加讨论影响时，变化趋势与实验相反，但是在极性较强的溶剂中添加一个溶剂分子就可以得到实验的趋势。想请问一下什么情况下需要隐性溶剂化模型与溶剂分子一起使用？@ sobereva 。谢谢老师。
作者
Author: sobereva    时间: 2018-11-29 23:19
xinren 发表于 2018-11-29 21:34
最近在计算激发态的电荷传递，有考虑到溶剂的影响，主要是质子性溶剂与极性非质子性溶剂。在计算非质子性溶 ...

溶剂分子与溶质存在较强的、隐式溶剂模型没法定性正确表现的作用时就需要加入显式溶剂分子，你当前的情况正是如此
作者
Author: xinren    时间: 2018-11-30 16:54
谢谢老师!
作者
Author: sslc1985    时间: 2019-4-17 15:10
Sob老师，隐式模型下计算溶解自由能，能算3个大气压的吗，我看文献有人做298k，3atm的，我记得您培训的时候说只能做动力学算压力，文献是咋做到的呢，他没做动力学。
作者
Author: sobereva    时间: 2019-4-17 17:44
sslc1985 发表于 2019-4-17 15:10
Sob老师，隐式模型下计算溶解自由能，能算3个大气压的吗，我看文献有人做298k，3atm的，我记得您培训的时候 ...

溶解自由能只能得到标况的
文章说法不实，要么你理解错了作者的意思
作者
Author: Somnus荣小荣    时间: 2019-5-7 23:25
求助！老师，我在文献上看到溶解自由能的定义是溶液和气相中溶质自由能的差异，而溶剂环境下溶质的自由能是溶质在气相下的自由能加上溶解自由能，那这样的话，溶剂环境下溶质的自由能不就是等于溶液中溶质的自由能+1.89kcal/mol（标况下1atm->1M自由能变）吗？这样的话直接使用例如，M062X/6-31G*-SMD溶剂模型,计算得到的G值+1.89kcal/mol（标况下1atm->1M自由能变）是不是就是溶剂环境下溶质自由能。为何我们用溶剂模型下计算的单点能减去气相下计算的单点能得到溶解自由能呢？是因为溶剂模型本身对溶解自由能的定义是溶剂模型下计算的单点能减去气相下计算的单点能吗？
作者
Author: sobereva    时间: 2019-5-8 09:34
Somnus荣小荣 发表于 2019-5-7 23:25
求助！老师，我在文献上看到溶解自由能的定义是溶液和气相中溶质自由能的差异，而溶剂环境下溶质的自由能是 ...

文中专门说了

例如，计算某分子在乙醇中的溶解自由能，就用比如# M052X/6-31G* scrf(SMD,solvent=ethanol)计算得到的单点能，减去# M052X/6-31G*计算得到的单点能。这种计算溶解自由能的方式是绝对合理的，那些专门搞溶剂模型的人计算溶解自由能也是用这种方式计算，并将得到的结果与实验结果相对比来拟合参数。如果对这点还心存疑虑的话，一定仔细去看J. Phys. Chem. A, 114, 13442 (2010)，包括其中补充材料部分，文章专门说了这个问题。
作者
Author: 韩同学    时间: 2019-5-8 09:35
sob老师，你好，想跟您请教一个问题，我想用高斯加溶剂：乙酸乙酯，这个溶剂不是高斯内置的，也没有Abraham提供全部计算参数，但还是希望加上 氢键酸度 和 氢键碱度 ，这种情况怎么办呢？可以用跟乙酸乙酯结构差不多的溶剂分子的参数代替莫

要是凑合一下，怎么看是不是凑成功了，一般来说，是加上 氢键酸度 和 氢键碱度好，还是干脆就不加了呐？谢谢您！
作者
Author: Somnus荣小荣    时间: 2019-5-8 10:16
sobereva 发表于 2019-5-8 09:34
文中专门说了

好的，赶紧仔细看文章，多谢sob 老师啦
作者
Author: Somnus荣小荣    时间: 2019-5-8 16:32
求助！谢谢老师推荐的文献看懂了溶质在溶剂环境下自由能的计算方法。一下是在计算的过程中遇到的一个疑问，求助老师！
使用M062X/6-31G*进行溶剂以及气相单点计算，以下简写M062X
(1)溶解自由能=M062X(SMD)单点-M062X(气相g)单点
(2)溶质在溶剂环境下自由能=溶解自由能+M062X(g)自由能
(3)M062X(g)自由能=M062X(g)单点+M062X-G校正项
将(1),(3)带入(2)中得到溶质在溶剂环境下自由能=M062X(SMD)单点-M062X(气相g)单点+M062X(g)单点+M062X-G校正项。
结果得到溶质在溶剂环境下自由能=M062X(SMD)单点+M062X-G校正项。

如果使用更高级别如CCSD(T)进行更高级别的单点校正。
(1)溶解自由能=M062X(SMD)单点-CCSD(T)(气相g)单点
(2)溶质在溶剂环境下自由能=溶解自由能+CCSD(T)(g)自由能
(3)CCSD(T)(g)自由能=CCSD(T)(g)单点+M062X-G校正项
(4)将(1),(3)带入(2)中得到溶质在溶剂环境下自由能=M062X(SMD)单点-CCSD(T)(g)单点+CCSD(T)(g)单点+M062X-G校正项
结果得到溶质在溶剂环境下自由能=M062X(SMD)单点+M062X-G校正项。

这样的话，按照将每一步都算准确的原理，无论使用多高级别的单点校正似乎都是没有用，而最终的结果只与G校正项的计算级别有关？
还是说高级别的单点校正不需要加在(1)中，只需加在(3)中？
作者
Author: sobereva    时间: 2019-5-9 09:15
韩同学 发表于 2019-5-8 09:35
sob老师，你好，想跟您请教一个问题，我想用高斯加溶剂：乙酸乙酯，这个溶剂不是高斯内置的，也没有Abraham ...

不可能不加，不设置那些就不是SMD计算需要的完整的溶剂参数，根本没法算
作者
Author: sobereva    时间: 2019-5-9 09:18
Somnus荣小荣 发表于 2019-5-8 16:32
求助！谢谢老师推荐的文献看懂了溶质在溶剂环境下自由能的计算方法。一下是在计算的过程中遇到的一个疑问， ...

即便你用CCSD(T)算单点，依然是
溶解自由能=M062X(SMD)单点-M062X(气相g)单点

建议用SMD文中测试过的M05-2X算溶解自由能，M06-2X在SMD原文中未经测试
作者
Author: Somnus荣小荣    时间: 2019-5-9 10:52
sobereva 发表于 2019-5-9 09:18
即便你用CCSD(T)算单点，依然是
溶解自由能=M062X(SMD)单点-M062X(气相g)单点

好嘞，明白，多谢sob老师！
作者
Author: Somnus荣小荣    时间: 2019-5-10 11:40
求助！抱歉老师，还是有一个关于溶解自由能的疑问，在J. Phys. Chem. A, 114, 13442 (2010)文章，它给出的计算溶解自由能的公式是：△Gsolv=（Esoln+Gnes）-Egas （3），这里边的Gnes是非静电贡献的总和，而在SMD溶剂化模型中这一项是指的哪个值呢？，是SCF 得到的单点值包括了吗？我认为是scf 得到的单点值就包括了非静电贡献了，所以Esoln+Gnes=scf能量。像cpcm这些没有考虑非静电贡献，Gnes就是0,所以Esoln就是SCF能量，但是我老师不相信我的说法，让我找到具体的文献说明Gnes确实被包含进了scf 得到的能量中。并且文献还提到很多人误认为（Esoln+Gnes）是total free energy in solution，那这样做肯定会引入额外的误差，但是还是没有说明（Esoln+Gnes）究竟是哪个值，此外J. Phys. Chem. B 2009, 113, 6378–6396，SMD的文章也只是说到将非静电贡献加入到计算溶解自由能中，还是有点不清楚Gnes去哪了，恳请老师帮忙解答
作者
Author: sobereva    时间: 2019-5-10 16:31
Somnus荣小荣 发表于 2019-5-10 11:40
求助！抱歉老师，还是有一个关于溶解自由能的疑问，在J. Phys. Chem. A, 114, 13442 (2010)文章，它给出的 ...

非极性部分已经包含在了SMD下计算的单点里，这一项在输出文件里也单独给出了，而且提示得极其明确

你们老师根本不懂，让他仔细去看SMD原文去
作者
Author: Somnus荣小荣    时间: 2019-5-10 19:30
sobereva 发表于 2019-5-10 16:31
非极性部分已经包含在了SMD下计算的单点里，这一项在输出文件里也单独给出了，而且提示得极其明确

你 ...

嗯嗯，明白了，多谢sob老师！
作者
Author: Somnus荣小荣    时间: 2019-5-13 10:35
求助！老师我最近算液相反应的自由能变，对反应物以及产物气相优化的时候发现，气相无法优化到稳定的complex（反应物是自由基和离子，优化的过程中两种反应物的距离越来越远，达到十几埃了，而且能量是一直下降一直不收敛），但是液相可以得到稳定的complex，这样的情况是不是就得用液相的优化结构算气相的单点了。
作者
Author: sobereva    时间: 2019-5-13 13:09
Somnus荣小荣 发表于 2019-5-13 10:35
求助！老师我最近算液相反应的自由能变，对反应物以及产物气相优化的时候发现，气相无法优化到稳定的comple ...

是
凡是液相的反应，建议优化和单点一律在溶剂模型下进行
作者
Author: sobereva    时间: 2019-5-13 14:00
今日对文中大量细节做了补充和更新
作者
Author: FANYI    时间: 2019-6-12 11:22
老师，请教个问题：
我在用隐式溶剂模型计算溶剂模型下的单点能（m052x，scrf(SMD,solvent=...)）时，不收敛；在输入文件里加关键词guess=read，去读优化的chk文件，还是报错，不收敛。报错内容如下：

>>>>>>>>>> Convergence criterion not met.
Error on total polarization charges =  0.04078
SCF Done:  E(UM052X) =  -3563.93254300     A.U. after  129 cycles
            NFock=128  Conv=0.43D-03     -V/T= 2.0194
<Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 1.0000 <S**2>= 2.0159 S= 1.0053
<L.S>= 0.000000000000E+00
KE= 3.496260347002D+03 PE=-1.309128207973D+04 EE= 3.722008145278D+03
SMD-CDS (non-electrostatic) energy       (kcal/mol) =      -5.14
(included in total energy above)
Annihilation of the first spin contaminant:
S**2 before annihilation     2.0159,   after     2.0001
Convergence failure -- run terminated.
Error termination via Lnk1e in /home/test/gaussian//g09/l502.exe at Wed Jun 12 10:59:06 2019.
Job cpu time:       0 days  9 hours 35 minutes 32.9 seconds.
File lengths (MBytes):  RWF=    315 Int=      0 D2E=      0 Chk=     20 Scr=      1

请老师指点，谢谢！
作者
Author: sobereva    时间: 2019-6-12 19:56
FANYI 发表于 2019-6-12 11:22
老师，请教个问题：
我在用隐式溶剂模型计算溶剂模型下的单点能（m052x，scrf(SMD,solvent=...)）时，不收 ...

解决SCF不收敛问题的方法
http://sobereva.com/61
作者
Author: FANYI    时间: 2019-6-13 08:26
sobereva 发表于 2019-6-12 19:56
解决SCF不收敛问题的方法
http://sobereva.com/61

谢谢老师！看了博文，放宽收敛标准，问题解决啦。
作者
Author: xiaos    时间: 2019-7-23 14:59
本帖最后由 xiaos 于 2019-7-23 15:04 编辑
sobereva 发表于 2019-6-12 19:56
解决SCF不收敛问题的方法
http://sobereva.com/61

sob老师您好。我的体系为174个原子的有机阴离子体系（含四个烷基磺酸根侧链），分子内存在明显的pi-pi相互作用。计算级别：B3LYP-D3BJ/6-311+G**//M062X-D3/6-311+G**（想对这两个泛函进行筛选，确定后续体系的计算级别）。在进行基态优化时，发现SMD隐式溶剂模型（溶剂为水）下优化的得到的构型（与晶体结构基本一致）与气相优化的构象存在明显差异。

1、为减小计算量，进行基态和激发态结构优化时，基组能否由6-311+G**改变为6-311G**？
2、以下溶解自由能的计算水平是否合理？基组是否应该使用6-311+G**与正常计算任务（基态几何优化，激发态计算，激发态几何优化等）保持一致? opt和freq能否分开计算（opt freq计算较慢）？
气相：#p M062X/6-31G* opt freq int=ultrafine scale=0.9580 empiricaldispersion=gd3
溶剂相：#p M062X/6-31G* opt freq int=ultrafine scale=0.9580  scrf=(SMD,solvent=warter) empiricaldispersion=gd3
3、优化得到的气相分子构型与晶体结构相差较大时，溶解自由能的计算是否合理？


作者
Author: sobereva    时间: 2019-7-23 18:52
xiaos 发表于 2019-7-23 14:59
sob老师您好。我的体系为174个原子的有机阴离子体系（含四个烷基磺酸根侧链），分子内存在明显的pi-pi相 ...

B3LYP-D3BJ/6-311+G**//M062X-D3/6-311+G**这种组合完全莫名其妙。//后头的是优化级别，哪能优化用的级别比单点还贵？反过来还差不多。

1 可以
2 既然你的那个级别是算溶解自由能，拿优化后的结构在气相和溶剂模型下各算个单点求差就完了，干嘛还写opt freq？
opt freq比opt和freq单独算通常花的时间更少
3 带电体系，气相和极性溶剂环境下优化出的结果有显著差异十分正常。不要多想
作者
Author: xiaos    时间: 2019-7-23 19:13
sobereva 发表于 2019-7-23 18:52
B3LYP-D3BJ/6-311+G**//M062X-D3/6-311+G**这种组合完全莫名其妙。//后头的是优化级别，哪能优化用的级别 ...

好的好的   谢谢sob老师的解答
作者
Author: linqiaosong    时间: 2019-8-28 11:44
sob老师，我想请教一个问题：在用不起显式溶剂模型的时候，能不能直接用隐式模型下的自由能代替气相反应TST的自由能，用来近似计算溶液体系反应的速率常数呢？
作者
Author: sobereva    时间: 2019-8-28 23:14
linqiaosong 发表于 2019-8-28 11:44
sob老师，我想请教一个问题：在用不起显式溶剂模型的时候，能不能直接用隐式模型下的自由能代替气相反应TST ...

可以。这是非常常规的做法
作者
Author: linqiaosong    时间: 2019-8-29 11:15
sobereva 发表于 2019-8-28 23:14
可以。这是非常常规的做法

非常感谢sob老师的回答
作者
Author: wr52139    时间: 2019-8-30 03:12
Sob老师好！关于您这篇帖子里：“溶质在溶剂环境下标况（298.15K、1M）时的自由能 = 溶质在1atm气相下的自由能 + 直接通过隐式溶剂模型计算的溶解自由能（始末都是1M） + 1.89kcal/mol（标况下1atm→1M自由能变）”有两个问题有点不懂：
1. 我要计算一个反应的能垒和反应晗 A + B --> C + D，我之前的做法是，在溶剂模型下用M062x/6-311G(d,p)优化得到Thermal Correction to Free Energy, 然后用M062x/def2TZVP计算单点（就是您说的一般人的做法），再将两者相加得到其自由能。 再用TS的自由能减去A和B的自由能。但看了您的帖子之后，发现我没有加上“溶质在1atm气相下的自由能” 和 “1.89kcal/mol”，我想问，计算能垒的时候，是否必须要加上“溶质在1atm气相下的自由能”？是否要在气相条件下再做一次结构优化和计算单点能得到它的自由能？

2. 双分子反应，如果没有考虑到1.89kcal/mol的校正，那算出来的能垒就会高1.89 kcal/mol。如果有两个或三个水分子参与了反应，那考虑这个校正，就会是(TS+1.89)-(A+1.89)-(B+1.89)-N(水分子个数)*1.89？
作者
Author: sobereva    时间: 2019-8-30 04:10
wr52139 发表于 2019-8-30 03:12
Sob老师好！关于您这篇帖子里：“溶质在溶剂环境下标况（298.15K、1M）时的自由能 = 溶质在1atm气相下的自 ...

1 你的问题很奇怪。“溶质在1atm气相下的自由能”不可能加上，否则等于电子能量算了两遍
仔细、反复阅读文章，里面的描述很充分且没有歧义
2 几分子反应和反应机理计算中有几个水完全不是一码事
作者
Author: sobereva    时间: 2019-8-30 05:48
鉴于总有初学者搞不明白计算流程，干脆给此文增加了第4节，是一个详细的计算例子，写得没法更清楚了。
作者
Author: wr52139    时间: 2019-8-30 07:07
sobereva 发表于 2019-8-30 05:48
鉴于总有初学者搞不明白计算流程，干脆给此文增加了第4节，是一个详细的计算例子，写得没法更清楚了。

感谢Sob老师的耐心解答，这次彻底明白了！
作者
Author: ch199412    时间: 2019-10-23 09:51
本帖最后由 ch199412 于 2019-10-23 20:34 编辑

sob老师，为什么我仿照你的例子中2.out那个，刚运行就出现2070错误呢？
作者
Author: ch199412    时间: 2019-10-23 21:58
ch199412 发表于 2019-10-23 09:51
sob老师，为什么我仿照你的例子中2.out那个，刚运行就出现2070错误呢？

老师，我已解决，看到l李老师的问题总结了，感谢平台
作者
Author: 国风寒洋    时间: 2019-10-30 16:55
感谢sob老师
作者
Author: sobereva    时间: 2019-11-26 22:40
鉴于时不时有人问怎么算质子的溶解自由能，很多人犯了严重误区，今日在此文第二节添加了“(8)关于单原子离子的溶解自由能的计算”
作者
Author: linima    时间: 2019-11-27 13:35
sob老师 ， 我把您的算环氧乙烷的实例1号文件在自己电脑上跑了一下出现这个图一的情况，我用的高斯09 ，另外我主要是想跑我的CO2自由基负离子，想算这个自由基负离子在乙腈环境下的自由能：

#P B3LYP/6-311G* em=GD3BJ opt freq scrf(SMD,solvent=acetonitrile)

Title Card Required

-1 2

这和结果跑下来还是图一的情况。
作者
Author: sobereva    时间: 2019-11-28 02:24
linima 发表于 2019-11-27 13:35
sob老师 ， 我把您的算环氧乙烷的实例1号文件在自己电脑上跑了一下出现这个图一的情况，我用的高斯09 ，另 ...

提问Gaussian问题的常识：绝对不要问l几几几报错怎么回事，更不要问2070提示怎么处理，应当从输出文件末尾逐渐往上找对应报错的确切的提示语句并贴出来，关键词也最好一并贴出来。说不清楚的时候直接上传输出文件。
作者
Author: linima    时间: 2019-11-28 15:01
本帖最后由 linima 于 2019-11-28 15:02 编辑
sobereva 发表于 2019-11-28 02:24
提问Gaussian问题的常识：绝对不要问l几几几报错怎么回事，更不要问2070提示怎么处理，应当从输出文件末 ...

em=GD3BJ  sob老师这个短语下面有个黑点
%chk=C:\mol.chk
------------------------------------------------------------
#P B3LYP/6-311G* em=GD3BJ opt freq scrf(SMD,solvent=ethanol)
------------------------------------------------------------
QPErr --- A syntax error was detected in the input line.
#P B3LYP/6-311G* em=GD3BJ opt freq scrf(
                                 '
Last state="GCL"
TCursr= 1177 LCursr=   19
Error termination via Lnk1e in d:\gaussian\l1.exe at Thu Nov 28 14:58:50 2019.
Job cpu time:  0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      1 Int=      0 D2E=      0 Chk=      1 Scr=      1
作者
Author: linima    时间: 2019-11-28 16:53
sobereva 发表于 2019-11-28 02:24
提问Gaussian问题的常识：绝对不要问l几几几报错怎么回事，更不要问2070提示怎么处理，应当从输出文件末 ...

我用热力学组合方法的时候也是出现这种错误，在W1下面标了一个点

%chk=C:\Users\Administrator\Desktop\O2-gas.chk
-----------------------
#p W1 geom=connectivity
-----------------------
QPErr --- An ambiguous keyword was detected.
#p W1 geom=connectivity
        '
Last state="GCL"
TCursr=  430 LCursr=    3
Error termination via Lnk1e in d:\gaussian\l1.exe at Thu Nov 28 14:22:58 2019.
Job cpu time:  0 days  0 hours  0 minutes  0.0 seconds.
File lengths (MBytes):  RWF=      1 Int=      0 D2E=      0 Chk=      1 Scr=      1
作者
Author: sobereva    时间: 2019-11-29 08:29
linima 发表于 2019-11-28 16:53
我用热力学组合方法的时候也是出现这种错误，在W1下面标了一个点

%chk=C:%users\Administrator\Deskto ...

没法直接写W1
看手册就知道，要么W1U要么W1BD
作者
Author: sobereva    时间: 2019-11-29 08:30
linima 发表于 2019-11-28 15:01
em=GD3BJ  sob老师这个短语下面有个黑点
%chk=C:\mol.chk
----------------------------------------- ...

G09 D.01之前的版本不支持DFT-D3，看
DFT-D色散校正的使用
http://sobereva.com/210
作者
Author: naonao5205    时间: 2019-12-14 18:46
sob老师你好，请问在文中提到：
“但是，在优化和振动分析过程中使用SMD时，极个别情况下可能造成本来是高对称性的结构出现虚频、优化收敛变得困难，而默认的IEFPCM则没这个问题。因此使用SMD优化和振动分析时碰到这种情况，可以尝试改用默认的IEFPCM再试。”
请问在实际过程中，在一条反应路径中，使用IEFPCM优化和震动分析的结果可以和SMD模型的结果进行比较吗？（因为SMD难以得到优化的结果）
作者
Author: sobereva    时间: 2019-12-15 10:45
naonao5205 发表于 2019-12-14 18:46
sob老师你好，请问在文中提到：
“但是，在优化和振动分析过程中使用SMD时，极个别情况下可能造成本来是高 ...

某一类型任务在不同体系间横向对比时，要用IEFPCM就都用，要么都用SMD。否则容易被审稿人批
作者
Author: naonao5205    时间: 2019-12-16 13:47
本帖最后由 naonao5205 于 2019-12-21 11:06 编辑
sobereva 发表于 2019-12-15 10:45
某一类型任务在不同体系间横向对比时，要用IEFPCM就都用，要么都用SMD。否则容易被审稿人批

谢谢sob老师解答。
使用IEFPCM模型优化好的结构作为smd模型初始结构，opt(gdiis,maxstep=5,notrust) 后能正常收敛了。（猜测是分子太柔了给的初始构型不够好）修改：
老师要求下换成09（fine） 没遇到SMD无法震荡错误了
作者
Author: leo1992    时间: 2019-12-20 23:14
原文：在SMD模型下计算溶解自由能，一定要用M05-2X/6-31G*  对于体系带有负电荷，M05-2X/6-31G*是否需要改成M05-2X/6-31+G*
作者
Author: sobereva    时间: 2019-12-22 00:34
leo1992 发表于 2019-12-20 23:14
原文：在SMD模型下计算溶解自由能，一定要用M05-2X/6-31G*  对于体系带有负电荷，M05-2X/6-31G*是否需要改 ...

不需要
作者
Author: a815648905    时间: 2019-12-26 21:56
老师，想请教一下，如果想计算反应势垒，按照您上面的方法需要把每个分子都纳入吗？比如反应物和生成物有气体的，是不是不用计算溶解自由能，谢谢老师了。
作者
Author: sobereva    时间: 2020-1-3 12:07
a815648905 发表于 2019-12-26 21:56
老师，想请教一下，如果想计算反应势垒，按照您上面的方法需要把每个分子都纳入吗？比如反应物和生成物有气 ...

单独的气体状态的分子不用考虑溶解自由能。但如果是已经进入了溶液里面才参加的反应（诸如计算溶液中的反应物复合物），需要把复合物当做溶液中的体系考虑溶解自由能。
作者
Author: 九曜    时间: 2020-2-11 16:53
社长，实验里的反应是在CDCl3溶剂里进行的，NMR检测的产物比例。在计算该反应机理的时候溶剂里没有CDCl3，请问写SCRF(SMD,solvent=chloroform)可以么。
作者
Author: sobereva    时间: 2020-2-13 17:39
九曜 发表于 2020-2-11 16:53
社长，实验里的反应是在CDCl3溶剂里进行的，NMR检测的产物比例。在计算该反应机理的时候溶剂里没有CDCl3， ...

溶剂模型不区分同位素
就这么写
作者
Author: 金城羽Alex    时间: 2020-5-10 21:03
sobereva老师您好，我按照您说的方法算溶质在溶剂中吉布斯自由能（以溶剂化中结构算thermal correction和一溶剂化优化的结构为基础算高精度单点等等）
但是我的导师不认同这种方法，认准了Ho文章J. Phys. Chem. B 2016, 120, 1319−1329中提出的要用气相优化结构。也听不进去我解释的，按照您说溶剂化优化结构，也听不进去我给他看您在前面说的“次年Cramer和Truhlar在JPCB,115,14556(2011)中打脸”的文章。
有没有其他的文献能证明溶剂化优化结构是比气相中优化结构更合适的文献让我跟他讲道理吗？
我今年毕业，老板全盘否定我算的吉布斯自由能让我全部重算，但是我并不觉得他的方法合理。他让我用已经用溶剂化优化的结构做M052X/6-31g(d)的气相优化结构，用这个结构做高精度单点算电子能。
诚挚地等待您的回复
作者
Author: sobereva    时间: 2020-5-11 14:05
金城羽Alex 发表于 2020-5-10 21:03
sobereva老师您好，我按照您说的方法算溶质在溶剂中吉布斯自由能（以溶剂化中结构算thermal correction和一 ...

有Truhlar这一篇就足够顶事了，Truhlar+Cramer的名气比JPCB那篇的作者大不知多少
作者
Author: 金城羽Alex    时间: 2020-5-12 10:08
sobereva 发表于 2020-5-11 14:05
有Truhlar这一篇就足够顶事了，Truhlar+Cramer的名气比JPCB那篇的作者大不知多少

Sob老师，我的导师太犟了，他不认为Truhlar这一篇说明溶剂化优化结构更好。现在逼我全部重新算气相下优化结构，不按照他说的做我就后果自负，但是19号就要交论文和答辩材料了，根本算不完。
我没有办法说服他溶剂化优化更好（跟何况我整条反应路径中还有两个带+1电的结构）
作者
Author: 金城羽Alex    时间: 2020-5-12 10:20
sobereva 发表于 2020-5-11 14:05
有Truhlar这一篇就足够顶事了，Truhlar+Cramer的名气比JPCB那篇的作者大不知多少

并且我自己找到了Ho等人在2016年发表的J. Phys. Chem. B 2016, 120, 1319−1329，在这篇文献中他们仍旧采取的是气相中优化的结构。我导师抓住这一点认为如果Ho等人错了，那这篇就不会延续用气相结构（Truhlar那一篇是2011年发表的）
作者
Author: sobereva    时间: 2020-5-13 06:52
金城羽Alex 发表于 2020-5-12 10:08
Sob老师，我的导师太犟了，他不认为Truhlar这一篇说明溶剂化优化结构更好。现在逼我全部重新算气相下优化 ...

就后果自负呗
在气相下搞带净电荷的体系，到时候极有可能挨审稿人批，把Ho那篇搬出来也不顶事
作者
Author: slope    时间: 2020-5-13 11:14
sob老师您好，我在计算某反应的溶剂化自由能时，里面涉及氢氧根离子。氢氧根离子的自由能是通过计算液态的H2O和H2得到，但是H2O和H2又不参加反应，所以我没有进行溶剂化处理。想请教下sob老师，在这里计算H2O和H2的自由能时，是否需要进行溶剂化处理吗？谢谢老师。
作者
Author: leo1992    时间: 2020-5-13 11:34
sob老师，如果反应在水溶液下进行，计算反应的吉布斯自由能变。那么计算高精度单点能是否也需要加上水溶剂校正代替气相中的高精度单点能。有阴离子参与的反应采用M06-2X/ma-TZVPP级别是否可以。
作者
Author: sobereva    时间: 2020-5-14 16:00
slope 发表于 2020-5-13 11:14
sob老师您好，我在计算某反应的溶剂化自由能时，里面涉及氢氧根离子。氢氧根离子的自由能是通过计算液态的H ...

我不晓得氢氧根离子的自由能是如何“通过计算液态的H2O和H2得到”，原理上并没法直接实现，除非借助额外数据。
作者
Author: sobereva    时间: 2020-5-14 16:07
leo1992 发表于 2020-5-13 11:34
sob老师，如果反应在水溶液下进行，计算反应的吉布斯自由能变。那么计算高精度单点能是否也需要加上水溶剂 ...

“水溶剂校正”指代不明。
溶剂下的自由能怎么算在主贴里已经给了具体的例子了
阴离子的电子能量用这个级别没问题



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3