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

1、DFT 方法和 HF 方法的比较计算中非常常用的方法。DFT 方法考虑了电子相关，而且算起来也快，所以很多人喜欢用，其中又以 B3LYP 最为常用，所以清楚每种方法的优劣和适用范围有利于我们对问题的理解。 HF 方法忽略的大部分的电子相关。相反，在很多时候 DFT 方法常常是过多的考虑了电子相关，这会使得过渡态的能量偏低，造成算出来的活化能偏低而且计算氢键的键能也会偏低。 并且在计算有机分子的芳香性也不好，DFT 会过多考虑电子离域，导致计算出来的能量偏低 但对于过渡金属、有机生物分子，DFT 方法都能很好的处理，这是它比其它方法好的地方。 由于 HF 方法忽略了电子相关，所以在处理弱键的时候

2、是不好的，所以我们一般也不用它来计算活化能（偏高）和解离能（偏低）并且在体系中如果有孤对电子、共扼体系的话，电子比较松散或者离域的比较厉害，这时候电子相关也很重要，HF 方法也不能很好的定量解决。 B3LYP 与其它泛函相比，对分子基态得到的特性一般相差不大。所以B3LYP 很常用也无可厚非。 但是一些人用 TDDFT 算激发能也用 B3LYP，而且只用 B3LYP，很差的结果也拿出来发表，根本不考虑原因，这就有问题了。上个世纪末，很多使用TDDFT 算激发能的文章都得到一个相同的结论，就是 B3LYP 作 TDDFT 激发能计算的结果是不可靠的：对不同的分子体系，有的时候跟实验值相当接近，有

3、的时候却差得不得了。因此在做 TDDFT 激发能计算的时候，应该多试几种泛函，特别是没有实验值，或者 B3LYP 的结果比较差的时候。如果图省事，想找一个通用的泛函作 TDDFT 计算，推荐用 PBE 极其几种改进版，Gaussian 98的用户可以用 B3PW91。它们比较稳定，虽然不总是最好，但也不会太差。另外，HCTH 系列泛函（特别是 HCTH407）也值得一试。 B3LYP 之所以计算 TS 能量会偏低,主要在于其交换相关势不够准确 ,特别是在长程区的渐近行为不够好,也正是如此,b3lyp 是不可能准确计算氢键 .上次徐昕老师在作报告时就提到了 b3lyp 在计算氢键和 van de

4、r wal 力结果是 poor and bad,他发展的 X3LYP 是可以得到较好的结果的.用 B3LYP 计算弱相互作用体系 He2是不成键的,而 X3LYP 计算是成键的,所以 X3LYP 在远程势的渐近行为是比较好的,如果用来计算过渡态的话也比 b3lyp 要准确的多 .当然,现在出现了很多的泛函结果都会比 b3lyp 准确,不过要加入程序还尚待时日 ,希望大家在计算前要仔细考虑所用的方法合不合适. 量子化学计算方法的演进主要分为(1)分子轨道法（简称 MO 法）和(2)价键法。以下只介绍 MO 法。分子轨道法的核心是哈特里福克罗特汉方法，简称 HFR 方程，它是以三个在分子轨道法发展

5、过程中做出卓越贡献的人命名的方程。1928 年DR哈特里（Hartree）提出了一个将 N 个电子体系中的每一个电子都看成是由其余的 N-1 个电子所提供的平均势场运动的假设。这样对于体系中的每一个电子都得到了一个单电子方程（表示这个电子运动状态的量子力学方程） ，称为哈特里方程。使用自洽场迭代方式求解这个方程（自洽场分子轨道法） ，就可以得到体系的电子结构和性质。哈特里方程未考虑由于电子的自旋而需要遵守泡利原理。1930 年，BA福克（Fock）和 JC斯莱特（Slater）分别提出了考虑泡利原理的自洽场迭代方程，称为哈特里福克方程。它的单电子轨函数（即分子轨道）取为自旋轨函数（即电子的空间

6、函数与自旋函数的乘积） 。泡利原理要求，体系的总电子波函数要满足反对称化要求，即对于体系的任何两个粒子的坐标交换都使总电子波函数改变正负号，而斯莱特行列式波函数正是满足反对称化要求的波函数。将哈特里福克方程用于计算多原子分子，会遇到计算上的困难。CCJ罗特汉（Roothaan）提出将分子轨道向组成分子的原子轨道（AO）展开，这样的分子轨道称为原子轨道的线性组合（简称 LCAO） 。使用 LCAO-MO，原来积分微分形式的哈特里福克方程就变为易于求解的代数方程，称为哈特里福克罗特汉方程，简称 HFR 方程。RHF 方程 闭壳层体系是指体系中所有的电子均按自旋相反的方式配对充满某些壳层（壳层指一个

7、分子能级或能量相同的即简并的两个分子能级） 。这种体系的特点，是可用单斯莱特行列式表示多电子波函数（分子的状态） ，描述这种体系的 HFR 方程称为限制性的 HFR 方程，所谓限制性，是要求每一对自旋相反的电子具有相同的空间函数。限制性的 HFR 方程简称 RHF 方程。UHF 方程 开壳层体系是指体系中有未成对的电子（即有的壳层未充满） 。描述开壳层体系的波函数一般应取斯莱特行列式的线性组合，这样计算方案就将很复杂，然而对于开壳层体系的对应极大多重度（所谓多重度，指一个分子因总自旋角动量的不同而具有几个能量相重的状态）的状态（即自旋角动量最大的状态）来说，可以保持波函数的单斯莱特行列式形式（

8、近似方法）以描述这些体系。从头计算法 原则上讲，有了 HFR 方程（不论是 RHF 方程或是 UHF 方程） ，就可以计算任何多原子体系的电子结构和性质。真正严格的计算称之为从头计算法。在从头计算法里，分子轨道由组成体系的原子的全部原子轨道线性组合而成。对于原子轨道有不同的选法。斯莱特型轨道适于描写电子云的分布，但在计算一些积分时包含对无穷级数的积分，十分麻烦，所以在从头计算法里，常取高斯型函数做为基函数，取一个高斯型函数或数个高斯型函数的线性组合模拟一个原子轨道。已经有了不少进行多原子体系的从头计算法的标准的计算机程序，如 JA波普尔研究集体推出计算机程序系列：高斯系列。该系列第一个公开的版

9、本为高斯 70（指 1970 年） ，以后差不多每年更新一次，功能逐次增多，算法也日趋完善。赝势价轨道从头计算法 从直观的化学观出发，可以想见，在原子形成分子时，仅仅原子的价层电子发生了较大的变形，而内层电子分布则改变较小。为了节省计算时间，而又不失去计算精度，从 20 世纪 70 年代以来，开始只考虑原子的价电子，而把内层电子和原子核看成一个凝固的原子实，用一个核模型势来代替内层电子与价电子的从头计算法，称为赝势价轨道从头计算法。不同的方案对于赝势的取法不尽相同，且大都能得到与全电子从头计算法相近的结果，因而大大节省了计算时间，特别对于含有重原子的体系，恰好是全电子从头算难以处理的。赝势价轨

10、道从头计算法将会在过渡金属络合催化的量子化学研究方面发挥重要作用。组态相互作用法 比从头计算法精确度更高的组态相互作用法有两个优点：一是不依赖于试探波函数的形式，就能原则上提非相对论薛丁谔方程的精确解；二是原则上可用于原子或分子体系的任何稳定态。近似计算法 在半经验计算法中，目前最常用的近似是零微分重叠（简称ZDO） ；近似程度最高的是全略微分重叠（简称 CNDO）近似，是 JA波普尔在 1965 年提出的。此外，在近似计算法中，还有休克尔分子轨道法和推广的分子轨道法。Chem3D 使用心得这是一个三维分子结构演示软件，上面有和 Gaussian 的链接，可以直接从 Chem3D 打开 Gau

11、ssian 运算，本身也带有分子力场计算功能，还可以和 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 稍微有点不兼容，

12、但是如果用AlChem 打开，在保存一下（仍然是*.mol） ，Chem3D 和 HyperChem 都可以打开了。最有用的可能是*.zmt,*.mol,*.gjf,*.fch 四种格式 ，比如用 Gaussian 优化以后，在 Ultilities 选项里面选 formcheck，就可以把 check 文件*.chk 转化成*.fch，而 *.fch 是可以直接用 Chem3D 打开的。*.zmt 是一内坐标的分子结构文件，我这里有个小软件可以把这种格式的文件转换成*.gjf，是 Gaussian 的分子内坐标文件，但没有 route section 等等。Chem3D 除了生成 Gauss

13、ian 输入文件以外，还有许多不错的功能，比如可以手动地转动某一个键，可以用鼠标显示某个键长，键角和二面角等地大小，等等，还可以显示所以内坐标。显示也有各种模型方式，比如球棍、stick 等等，可以显示原子编号原子名称，但是没有GaussView 和 HyperChem 以及 AlChem，WinMopac 等等漂亮。Gaussian 常见问题分析1 检查是否有初始文件错误在命令行中加入 %kJob L301 or %kJob L302如果通过则一般初始文件 ok。常见初级错误：a. 自旋多重度错误b. 变量赋值为整数c. 变量没有赋值或多重赋值d. 键角小于等于 0 度，大于等于 180 度

14、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. 给出精确的、对称性确定的角度和二面角。 如

15、 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

16、错误）解决方法：#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) (用

17、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

18、为单电子数。所以，关键是单电子的数目是多少。当有偶数个电子时，例如 O2，共有 16 个电子，那么单电子数目可能是 0，即8 个 alpha 和 8 个 beta 电子配对，对应单重态，但是也可能是有 9 个 电子和7 个 电子，那么能成对的是 7 对，还剩 2 个 没有配对，于是 n=2,对应的是多重度 3。同理还可以有多重度 5，7，9， .一般而言,是多重度低的能量低,最稳定,所以,一般来说,偶数电子的体系多重度就是 1。但是也有例外，例如 O2就是一个大家都知道的例子，它的基态是三重态，其单重态反而是激发态。所以对于未知的体系，还是算几个保险一点，看哪个能量更低。所以，总结一下，就是电

19、子数目是偶数，未成对电子数目 n=0,2,4,6,.自旋多重度是 1,3,5,7,.电子数目是奇数，未成对电子数目 n=1,3,5,7,.自旋多重度是 2,4,6,8,.多数情况是多重度低的能量低，有时（特别是有“磁”性的时候，例如顺磁的O2，以及 Fe 啊什么的） ，可能会高多重度的能量低，所以需要都算算，看哪个能量更低。关于赝势：简单来说，赝势就是不计算内层电子，而是把内层电子的贡献用一个势来描述，放在哈密顿里面。适用于重元素。赝势基组，实际上包括 赝势和基组两个部分，内层电子采用赝势，即 effective core potential (ECP)，外层价电子采用一般的基组。 比如：La

20、nL2DZ: 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 ro

21、w , 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 optl

22、anl2dz 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) /另一种全电子基组*C

23、u 0lanl2dz /定义价电子的基组*Cu 0lanl2dz /定义内层电子的赝势-开壳层和闭壳层闭壳层计算就是对于多重度是 1 的体系，此时 和 的电子数目相同，可以把 和 配对，成对的 和 使用同一个轨道，一个轨道上填充 2 个电子。开壳层计算就是对 和 电子分别计算，一个轨道上只填充 1 个电子，一般来说，多重度是 1 时，开壳层计算和闭壳层计算会给出相同的结果。限制开壳层计算是对多重度大于 1 的体系，此时 和 的电子数目不同，设有 m 个 和 n 个 电子，mn，那么让前 n 个轨道上每个填充一个 和一个 ，剩下的 m-n 个 电子再填充 m-n 个轨道。即前 n 个轨道是闭的（

24、每个轨道 2 个电子） ，后 m-n 个轨道是开的（每个轨道 1 个电子）在高斯中，以 HF 为例，闭、开、限制开壳层计算分别是RHF，UHF ，ROHF。如果只写 HF，则按下面的方式取默认方法：对多重度是 1 的体系，默认为 RHF，对多重度大于 1 的体系，默认是 UHF。关于收敛问题 (L502, L508, L9999)对于一个优化计算，它的过程是先做一个 SCF 计算，得到这个构型下的能量，然后优化构型，再做 SCF，然后再优化构型。 。 。因此，会有两种不收敛的情况：一是在某一步的 SCF 不收敛 (L502 错误)，或者构型优化没有找到最后结果（L9999 错误） 。预备知识:

25、计算时保存 chk 文件,可以在后续计算中使用 guess=read 读初始猜测.对于 SCF 不收敛，通常有以下的解决方法：1. 使用小基组，或低级算法计算，得到 scf 收敛的波函数，用 guess=read 读初始波函数。2. 使用 scf=qc，这个计算会慢，而且需要用 stable 关键字来测试结果是否波函数稳定。如果这个还不收敛，会提示 L508 错误。3. 改变键长，一般是缩小一点，有时会有用。4. 计算相同体系的其他电子态，比如相应的阴离子、阳离子体系或单重态体系，得到的收敛波函数作为初始猜测进行计算。5. 待补充.对于优化不收敛,即 L9999 错误，实际上是在规定的步数内没

26、有完成优化，即还没有找到极小值点。 （或者对于过渡态优化，还没有找到过渡态）这有几种可能性：1. 看一下能量的收敛的情况，可能正在单调减小，眼看有收敛的趋势，这样的情况下，只要加大循环的步数（opt(maxcycle=200)） ，可能就可以解决问题了。2. 加大循环步数还不能解决的（循环步数有人说超过 200 再不收敛，再加也不会有用了，这虽然不一定绝对正确，但 200 步应该也差不多了） ，有两种可能。一是查看能量，发现能量在振荡了，且变化已经很小了，这时可能重新算一下，或者构型稍微变一下，继续优化，就可以得到收敛的结果（当然也有麻烦的，看运气和经验了） ；二是构型变化太大，和你预计的差别

27、过大，这很可能是你的初始构型太差了，优化不知道到哪里去了，这时最好检查一下初始构型，再从头优化。3. 对于 L9999 快达到收敛时,考虑减小优化步长有时对于能量振荡的情况也是有用的，opt(maxstep=1).(flyingheart )一个建议是，对于大体系，难收敛体系，先用小基组，低精度算法优化一下，以得到较好的初始构型，再用高精度的计算接着算。如果前面的方法保留了chk 文件，重新计算时需要使用 geom=allcheck 读入构型（就不必麻烦地写构型了） ， guess=read（读入初始波函数，可以加快第一步 SCF 收敛） 。关于对称性：（原贴 http:/210.34.15.

28、126/cgi-bin/topic.cgi?forum=3&topic=1736）刚刚在 zixia 上看到 zhuoliu 的大作。结合自己以前知道的东西，总结一下。1。gaussian 中输入什么对称性，一般优化的结果仍然还是那个对称性，比如CO2，如果初始两个 CO 键长输入不是完全相等（比如一个 1.214,一个 1.215） ，那么程序就会判断为 Cv 对称，那么优化结果虽然键长几乎相等，但仍然认为是 Cv ，这个从振动频率或者分子轨道对称性上可以看出来。 -我们知道，CO2 实际上是直线的两边对称的构型，其对称性应该是 Dv 。因此，为了得到高的对称性，必须输入的时候，精确地输入数

29、值，比如 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=l

30、oose 只是在第一步判断输入构型的对称性时用到。此外，也可以用 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 symm

31、etry 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.d-8).N0 10-N.N0 10 -N.N0，修正后的体系能量比原来稍高一点。不考虑 BSSE，我们计算的结合能是用的 Eb=E(AB)E(A

32、0)-E(B0) ，注意，这里 A0 和 B0 分别是 A，B 的各自的优化得到的稳定构型的能量。那么考虑了BSSE 以后，就是(Eb)= Eb+BSSE energy.当然也就是 E6- E(A0)-E(B0).SCF 收敛失败的解决办法1.1. 在 Guess 关键字中使用 Core，Huckel 或 Mix 选项，试验不同的初始猜测；2. 对开壳层体系，尝试收敛到同一分子的闭壳层离子，接下来用作开壳层计算的初始猜测。添加电子可以给出更合理的虚轨道，但是作为普遍的经验规则，阳离子比阴离子更容易收敛。选项 Guess=Read 定义初始猜测从 Gaussian 计算生成的 checkpoin

33、t 文件中读取； 3. 另一个初始猜测方法是首先用小基组进行计算，由前一个波函得到用于大基组计算的初始猜测(Guess=Read 自动进行)； 4. 尝试能级移动(SCF=Vshift)； 5. 如果接近 SCF 但未达到，收敛标准就会放松或者忽略收敛标准。这通常用于不是在初始猜测而是在平衡结构收敛的几何优化。SCF=Sleazy 放松收敛标准，Conver 选项给出更多的控制； 6. 一些程序通过减小积分精度加速 SCF。对于使用弥散函数，长程作用或者低能量激发态的体系，必须使用高积分精度：SCF=NoVarAcc； 7. 尝试改变结构。首先略微减小键长，接下来略微增加键长，接下来再对结构作

34、一点改变； 8. 考虑使用不同的基组； 9. 考虑使用不同理论级别的计算。这并不总是实用的，但除此之外，增加迭代数量总是使得计算时间和使用更高理论级别差不多； 10. 更多的 SCF 迭代( SCF(MaxCycle=N)，其中 N 是迭代数)。这很少有帮助，但值得一试； 11. 关闭 DIIS 外推(SCF=NoDIIS) 。同时进行更多的迭代 ( SCF=(MaxCycle=N) )；12. 使用强制的收敛方法。SCF=QC 通常最佳，但在极少数情况下 SCF=DM 更快。不要忘记给计算额外增加一千个左右的迭代。应当测试这个方法获得的波函，保证它最小，并且正好不是稳定点(使用 Stable 关键字)； 13. 试着改用 DIIS 之外其它方法(SCF=SD 或 SCF=SSD)。