DFT方法和HF方法的比较.doc

DFT方法和HF方法的比较计算中非常常用的方法。DFT方法考虑了电子相关，而且算起来也快，所以很多人喜欢用，其中又以B3LYP最为常用，所以清楚每种方法的优劣和适用范围有利于我们对问题的理解。 HF方法忽略的大部分的电子相关。相反，在很多时候DFT方法常常是过多的考虑了电子相关，这会使得过渡态的能量偏低，造成算出来的活化能偏低而且计算氢键的键能也会偏低。 并且在计算有机分子的芳香性也不好，DFT会过多考虑电子离域，导致计算出来的能量偏低 但对于过渡金属、有机生物分子，DFT方法都能很好的处理，这是它比其它方法好的地方。 由于HF方法忽略了电子相关，所以在处理弱键的时候是不好的，所以我们一般也不用它来计算活化能（偏高）和解离能（偏低）并且在体系中如果有孤对电子、共扼体系的话，电子比较松散或者离域的比较厉害，这时候电子相关也很重要，HF方法也不能很好的定量解决。 B3LYP与其它泛函相比，对分子基态得到的特性一般相差不大。所以B3LYP很常用也无可厚非。 但是一些人用TDDFT算激发能也用B3LYP，而且只用B3LYP，很差的结果也拿出来发表，根本不考虑原因，这就有问题了。上个世纪末，很多使用TDDFT算激发能的文章都得到一个相同的结论，就是B3LYP作TDDFT激发能计算的结果是不可靠的：对不同的分子体系，有的时候跟实验值相当接近，有的时候却差得不得了。因此在做TDDFT激发能计算的时候，应该多试几种泛函，特别是没有实验值，或者B3LYP的结果比较差的时候。如果图省事，想找一个通用的泛函作TDDFT计算，推荐用PBE极其几种改进版，Gaussian 98的用户可以用B3PW91。它们比较稳定，虽然不总是最好，但也不会太差。另外，HCTH系列泛函（特别是HCTH407）也值得一试。 B3LYP之所以计算TS能量会偏低,主要在于其交换相关势不够准确,特别是在长程区的渐近行为不够好,也正是如此,b3lyp是不可能准确计算氢键.上次徐昕老师在作报告时就提到了b3lyp在计算氢键和van der wal力结果是poor and bad,他发展的 X3LYP是可以得到较好的结果的.用B3LYP计算弱相互作用体系He2是不成键的,而X3LYP计算是成键的,所以X3LYP在远程势的渐近行为是比较好的,如果用来计算过渡态的话也比b3lyp要准确的多.当然,现在出现了很多的泛函结果都会比b3lyp准确,不过要加入程序还尚待时日,希望大家在计算前要仔细考虑所用的方法合不合适. 量子化学计算方法的演进主要分为(1)分子轨道法（简称MO法）和(2)价键法。以下只介绍MO法。分子轨道法的核心是哈特里福克罗特汉方法，简称HFR方程，它是以三个在分子轨道法发展过程中做出卓越贡献的人命名的方程。1928年DR哈特里（Hartree）提出了一个将N个电子体系中的每一个电子都看成是由其余的N-1个电子所提供的平均势场运动的假设。这样对于体系中的每一个电子都得到了一个单电子方程（表示这个电子运动状态的量子力学方程），称为哈特里方程。使用自洽场迭代方式求解这个方程（自洽场分子轨道法），就可以得到体系的电子结构和性质。哈特里方程未考虑由于电子的自旋而需要遵守泡利原理。1930年，BA福克（Fock）和JC斯莱特（Slater）分别提出了考虑泡利原理的自洽场迭代方程，称为哈特里福克方程。它的单电子轨函数（即分子轨道）取为自旋轨函数（即电子的空间函数与自旋函数的乘积）。泡利原理要求，体系的总电子波函数要满足反对称化要求，即对于体系的任何两个粒子的坐标交换都使总电子波函数改变正负号，而斯莱特行列式波函数正是满足反对称化要求的波函数。将哈特里福克方程用于计算多原子分子，会遇到计算上的困难。CCJ罗特汉（Roothaan）提出将分子轨道向组成分子的原子轨道（AO）展开，这样的分子轨道称为原子轨道的线性组合（简称LCAO）。使用LCAO-MO，原来积分微分形式的哈特里福克方程就变为易于求解的代数方程，称为哈特里福克罗特汉方程，简称HFR方程。RHF方程 闭壳层体系是指体系中所有的电子均按自旋相反的方式配对充满某些壳层（壳层指一个分子能级或能量相同的即简并的两个分子能级）。这种体系的特点，是可用单斯莱特行列式表示多电子波函数（分子的状态），描述这种体系的HFR方程称为限制性的HFR方程，所谓限制性，是要求每一对自旋相反的电子具有相同的空间函数。限制性的HFR方程简称RHF方程。UHF方程 开壳层体系是指体系中有未成对的电子（即有的壳层未充满）。描述开壳层体系的波函数一般应取斯莱特行列式的线性组合，这样计算方案就将很复杂，然而对于开壳层体系的对应极大多重度（所谓多重度，指一个分子因总自旋角动量的不同而具有几个能量相重的状态）的状态（即自旋角动量最大的状态）来说，可以保持波函数的单斯莱特行列式形式（近似方法）以描述这些体系。从头计算法 原则上讲，有了HFR方程（不论是RHF方程或是UHF方程），就可以计算任何多原子体系的电子结构和性质。真正严格的计算称之为从头计算法。在从头计算法里，分子轨道由组成体系的原子的全部原子轨道线性组合而成。对于原子轨道有不同的选法。斯莱特型轨道适于描写电子云的分布，但在计算一些积分时包含对无穷级数的积分，十分麻烦，所以在从头计算法里，常取高斯型函数做为基函数，取一个高斯型函数或数个高斯型函数的线性组合模拟一个原子轨道。已经有了不少进行多原子体系的从头计算法的标准的计算机程序，如JA波普尔研究集体推出计算机程序系列：高斯系列。该系列第一个公开的版本为高斯70（指1970年），以后差不多每年更新一次，功能逐次增多，算法也日趋完善。赝势价轨道从头计算法 从直观的化学观出发，可以想见，在原子形成分子时，仅仅原子的价层电子发生了较大的变形，而内层电子分布则改变较小。为了节省计算时间，而又不失去计算精度，从20世纪70年代以来，开始只考虑原子的价电子，而把内层电子和原子核看成一个凝固的原子实，用一个核模型势来代替内层电子与价电子的从头计算法，称为赝势价轨道从头计算法。不同的方案对于赝势的取法不尽相同，且大都能得到与全电子从头计算法相近的结果，因而大大节省了计算时间，特别对于含有重原子的体系，恰好是全电子从头算难以处理的。赝势价轨道从头计算法将会在过渡金属络合催化的量子化学研究方面发挥重要作用。组态相互作用法 比从头计算法精确度更高的组态相互作用法有两个优点：一是不依赖于试探波函数的形式，就能原则上提非相对论薛丁谔方程的精确解；二是原则上可用于原子或分子体系的任何稳定态。近似计算法 在半经验计算法中，目前最常用的近似是零微分重叠（简称ZDO）；近似程度最高的是全略微分重叠（简称CNDO）近似，是JA波普尔在1965年提出的。此外，在近似计算法中，还有休克尔分子轨道法和推广的分子轨道法。Chem3D使用心得这是一个三维分子结构演示软件，上面有和Gaussian的链接，可以直接从 Chem3D打开Gaussian运算，本身也带有分子力场计算功能，还可以和Mopac 链接 可以打开的文件主要有以下几种： *.gjf,*.fch,*.cub,这几种是Gaussian的输入、formcheck和cube文件 *.dat,*.mop,*.mpc,mopac的输入文件 *.zmt,这个是内坐标格式的分子构型文件，和HyperChem、mopac通用 *.pdb,*.ent,Protein DB格式，和HyperChem、AlChem等等通用，Gaussian 本身也可以用这两种格式来生成内坐标文件 *.mol,MDL的mol格式，和HyperChem稍微有点不兼容，但是如果用AlChem打开，在保存一下（仍然是*.mol），Chem3D和HyperChem都可以打开了。最有用的可能是*.zmt,*.mol,*.gjf,*.fch四种格式 ，比如用Gaussian优化以后，在Ultilities选项里面选formcheck，就可以把check文件*.chk转化成*.fch，而*.fch是可以直接用Chem3D打开的。*.zmt是一内坐标的分子结构文件，我这里有个小软件可以把这种格式的文件转换成*.gjf，是Gaussian的分子内坐标文件，但没有route section等等。Chem3D除了生成Gaussian输入文件以外，还有许多不错的功能，比如可以手动地转动某一个键，可以用鼠标显示某个键长，键角和二面角等地大小，等等，还可以显示所以内坐标。显示也有各种模型方式，比如球棍、stick等等，可以显示原子编号原子名称，但是没有GaussView和HyperChem以及AlChem，WinMopac等等漂亮。Gaussian常见问题分析1 检查是否有初始文件错误在命令行中加入 %kJob L301 or %kJob L302如果通过则一般初始文件ok。常见初级错误：a. 自旋多重度错误b. 变量赋值为整数c. 变量没有赋值或多重赋值d. 键角小于等于0度，大于等于180度e. 分子描述后面没有空行f. 二面角判断错误，造成两个原子距离过近g. 分子描述一行内两次参考同一原子,或参考原子共线2 SCF(自洽场)不收敛则一般是L502错误省却情况做64个cycle迭代(G03缺省128 cycles)a. 修改坐标，使之合理b. 改变初始猜 Guess=Huckel 或其他的，看Guess关键词。c. 增加叠代次数 SCFCYC=N （对小分子作计算时最好不要增加，很可能结构不合理）d. iop（5/13=1）这样忽略不收敛，继续往下做。3 分子对称性改变a. 修改坐标，强制高对称性或放松对称性b. 给出精确的、对称性确定的角度和二面角。 如CH4的角度给到109.47122c. 放松对称性判据 Symm=loosed. 不做对称性检查 iop(2/16=1) （最好加这个选项） iop(2/16=2) 则保持新的对称性来计算4 Opt时收敛的问题a. 修改坐标，使之合理b. 增加叠代次数optcyc=N5 优化过渡态，若势能面太平缓，则不好找到。iop(1/8=10) 默认30（下一个结构和该结构的差别0.3Å），可改成10。如果每一步都要用到小的步长，应该加opt(notrustupdate)6 在CI（组态）方法中如QCISD(T),CCSD(T),CID方法中，省却最大循环50，若出错（L913错误）解决方法：#P QCISD(maxcyc=N) 注：N5127 优化过渡态opt=TS (给出过渡态) opt=qst2 (给出反应物和产物) opt=qst3 (给出反应物和产物和过渡态)a. 用G03时的出错 opt=ts 必须加FC (force constant)写法：opt=(TS, calcFc)or opt=(TS,calchffc)计算HF力常数，对QCISD，CCSD等方法用；or opt=(TS,modRedundant) (最好写这个)b. 如果计算采用QCISD计算（不好计算FC）则写为QCISD opt=(TS, calcHFFC) (用HF计算FC)8. 无法写大的Scratch文件RWFa. 劈裂RWF文件 %rwf=loc1,size1,loc2,size2,.,locN,-1b. 改变计算方法 MP2=Direct可以少占硬盘空间c. 限制最大硬盘 maxdisk=N GB，*MB,有些系统写2GB会出错，可以写2000MB9. FOPT出错 原因是变量数与分子自由度数不相等。 可用POPT 或直接用OPT10. 优化过渡态只能做一个STEP 原因是负本征数目不对 添加 iop（1/11）=1或者noeigentest量化小节（转）关于自旋多重度定义: 多重度2S+1, S=n*1/2,n为单电子数。所以，关键是单电子的数目是多少。当有偶数个电子时，例如 O2，共有16个电子，那么单电子数目可能是0，即8个alpha和8个beta电子配对，对应单重态，但是也可能是有9个电子和7个电子，那么能成对的是7对，还剩2个没有配对，于是n=2,对应的是多重度3。同理还可以有多重度5，7，9， .一般而言,是多重度低的能量低,最稳定,所以,一般来说,偶数电子的体系多重度就是1。但是也有例外，例如O2就是一个大家都知道的例子，它的基态是三重态，其单重态反而是激发态。所以对于未知的体系，还是算几个保险一点，看哪个能量更低。所以，总结一下，就是电子数目是偶数，未成对电子数目n=0,2,4,6,.自旋多重度是1,3,5,7,.电子数目是奇数，未成对电子数目n=1,3,5,7,.自旋多重度是2,4,6,8,.多数情况是多重度低的能量低，有时（特别是有“磁”性的时候，例如顺磁的O2，以及Fe啊什么的），可能会高多重度的能量低，所以需要都算算，看哪个能量更低。关于赝势：简单来说，赝势就是不计算内层电子，而是把内层电子的贡献用一个势来描述，放在哈密顿里面。适用于重元素。赝势基组，实际上包括 赝势和基组两个部分，内层电子采用赝势，即effective core potential (ECP)，外层价电子采用一般的基组。 比如：LanL2DZ: D95V on first row, Los Alamos ECP plus DZ on Na-Bi.就是对第一行原子是 D95V (这个是非赝势基组)，对Na-Bi是使用一个叫做Los Alamos的有效核势加上一个DZ基组。所以Lanl2dz就是对前面的原子全电子基组，对后面的原子是赝势基组。 （再次说明，量化里面，C，O那一行，算周期表的第一行）使用赝势的3个原因：1。没有相应的全电子基组。2。减少计算量。3。赝势可以包含重金属相对论效应的修正。在高斯中，lanl2dz基组，在手册中可以查到其定义为：LanL2DZ: D95V on first row , Los Alamos ECP plus DZ on Na-Bi，也就是说，对于C，O等元素来说（量化中，H大概算第0周期，C，O才是第一周期），Lanl2dz实际上还是全电子基组，而对于Na以后才是对内层电子用Los Alamos ECP赝势，外层电子用DZ 基组。使用赝势的输入文件：1.所有原子使用lanl2dz-#HF/lanl2dz optlanl2dz for all atoms0 2O 0.0 0.0 0.0C 0.0 0.0 1.2Cu 0.0 0.0 3.2-2.所有原子使用lanl2dz的另一种输入方法。-#HF/genecp optlanl2dz for all atoms0 2O 0.0 0.0 0.0C 0.0 0.0 1.2Cu 0.0 0.0 3.2C O 0lanl2dz*Cu 0lanl2dz /定义价电子的基组, C O 0 是碳，氧，零，其中零用作终止符号。*Cu 0lanl2dz /定义内层电子的赝势-3.混和基组，即有的使用全电子，有的使用Lanl2dz。格式同2。-#HF/genecp optlanl2dz for Cu, 6-31G(d) for C and O0 2O 0.0 0.0 0.0C 0.0 0.0 1.2Cu 0.0 0.0 3.2C O 06-31G(d) /另一种全电子基组*Cu 0lanl2dz /定义价电子的基组*Cu 0lanl2dz /定义内层电子的赝势-开壳层和闭壳层闭壳层计算就是对于多重度是1的体系，此时和的电子数目相同，可以把和配对，成对的和使用同一个轨道，一个轨道上填充2个电子。开壳层计算就是对和电子分别计算，一个轨道上只填充1个电子，一般来说，多重度是1时，开壳层计算和闭壳层计算会给出相同的结果。限制开壳层计算是对多重度大于1的体系，此时和的电子数目不同，设有m个和n个电子，mn，那么让前n个轨道上每个填充一个和一个，剩下的m-n个电子再填充m-n个轨道。即前n个轨道是闭的（每个轨道2个电子），后m-n个轨道是开的（每个轨道1个电子）在高斯中，以HF为例，闭、开、限制开壳层计算分别是RHF，UHF，ROHF。如果只写HF，则按下面的方式取默认方法：对多重度是1的体系，默认为RHF，对多重度大于1的体系，默认是UHF。关于收敛问题 (L502, L508, L9999)对于一个优化计算，它的过程是先做一个SCF计算，得到这个构型下的能量，然后优化构型，再做SCF，然后再优化构型。因此，会有两种不收敛的情况：一是在某一步的SCF不收敛(L502错误)，或者构型优化没有找到最后结果（L9999错误）。预备知识:计算时保存chk文件,可以在后续计算中使用guess=read读初始猜测.对于SCF不收敛，通常有以下的解决方法：1. 使用小基组，或低级算法计算，得到scf收敛的波函数，用guess=read读初始波函数。2. 使用scf=qc，这个计算会慢，而且需要用stable关键字来测试结果是否波函数稳定。如果这个还不收敛，会提示L508错误。3. 改变键长，一般是缩小一点，有时会有用。4. 计算相同体系的其他电子态，比如相应的阴离子、阳离子体系或单重态体系，得到的收敛波函数作为初始猜测进行计算。5. 待补充.对于优化不收敛,即L9999错误，实际上是在规定的步数内没有完成优化，即还没有找到极小值点。（或者对于过渡态优化，还没有找到过渡态）这有几种可能性：1. 看一下能量的收敛的情况，可能正在单调减小，眼看有收敛的趋势，这样的情况下，只要加大循环的步数（opt(maxcycle=200)），可能就可以解决问题了。2. 加大循环步数还不能解决的（循环步数有人说超过200再不收敛，再加也不会有用了，这虽然不一定绝对正确，但200步应该也差不多了），有两种可能。一是查看能量，发现能量在振荡了，且变化已经很小了，这时可能重新算一下，或者构型稍微变一下，继续优化，就可以得到收敛的结果（当然也有麻烦的，看运气和经验了）；二是构型变化太大，和你预计的差别过大，这很可能是你的初始构型太差了，优化不知道到哪里去了，这时最好检查一下初始构型，再从头优化。3. 对于L9999快达到收敛时,考虑减小优化步长有时对于能量振荡的情况也是有用的，opt(maxstep=1).(flyingheart )一个建议是，对于大体系，难收敛体系，先用小基组，低精度算法优化一下，以得到较好的初始构型，再用高精度的计算接着算。如果前面的方法保留了chk文件，重新计算时需要使用 geom=allcheck 读入构型（就不必麻烦地写构型了）， guess=read（读入初始波函数，可以加快第一步SCF收敛）。关于对称性：（原贴 26/cgi-bin/topic.cgi?forum=3&topic=1736）刚刚在zixia上看到zhuoliu的大作。结合自己以前知道的东西，总结一下。1。gaussian中输入什么对称性，一般优化的结果仍然还是那个对称性，比如CO2，如果初始两个CO键长输入不是完全相等（比如一个1.214,一个1.215），那么程序就会判断为Cv 对称，那么优化结果虽然键长几乎相等，但仍然认为是Cv ，这个从振动频率或者分子轨道对称性上可以看出来。-我们知道，CO2实际上是直线的两边对称的构型，其对称性应该是Dv 。因此，为了得到高的对称性，必须输入的时候，精确地输入数值，比如sqrt(2)，就要保留很多的小数点，180.0角，就不能写成179.9。2。有时计算过程中对称性会变化，比如做过渡态的时候，这时需要用 IOP(2/16=3)，否则计算会出错退出。3。比如用直角坐标输入一个正三角形构型，其对称性应该是D3h，但是如果输入的小数点后面的数字不够多，那么常常得到的是C2v或其它。为了消除输入文件中坐标的有效位数的影响，得到较高的对称性，可以降低对称性判断的严格性。一般可以用symm=loose，这等价于 IOP (2/17=4, 2/18=3)。还可以减小这4和3这两个数值，使得更加loose，但不能过小，否则会出错。symm=loose只是在第一步判断输入构型的对称性时用到。此外，也可以用gaussview来调整设置初始构型的对称性。4。如果要降低对称性，那么可以用symm(PG=C3v)等等来做。使判断出来的对称性为C3v的一个子群。即由PG来限制最高对称性。附用到的IOP的详细解释。IOp(2/16)action taken if the point group changes during an optimization. IOP(2/16=0)Abort the job. IOP(2/16=1)Keep going. IOP(2/16=2)Keep going and leave symmetry on, using the old symmetry. IOP(2/16=3)Keep going and leave symmetry on, using the new symmetryIOp(17)Tolerance for distance comparisons in symmetry determination.0 Default (determined in the symmetry package, currently 1.