DMol3基本原理和参数设置

1、DMol3 基本原理和参数设置 SCF(DFT) 制备材料表征材料 各种性能参数 制备工艺 结构 能量 电子波函数能带，态密度，电荷密度 一阶导数 应力，弹性力常数， 体模量. 二阶导数 声子频率，散射谱 H = E 力场 怎么使SCF收敛？ 如何使SCF的结果准确？ Walter Kohn E To U Exc )()()( , 2 rrrv kikiki eff )(),.( 1 rErrE N Hohenberg-Kohn theorem Kohn-Sham equations Exact only for ground state Needs approximation to Exc

2、Kohn-Sham theorem 生成KS势 求解KS方程 得到新nout(r) 和之前的n(r)比较 收敛与否 ? nstart(r) 输出结果 )()()( , 2 rrrnv kikiki eff rrr rrr )(rnveff r rrfrn i i ii rrr 电荷密度依据电子自旋方向的不同，一分为二： 自旋密度定义为: 总的磁矩定义为: 存在自旋极化的体系中，交换相关势也相应的变为不同自旋方向电荷密度 的函数: PS： 由于在自洽计算中需要分别考虑电荷密度和自旋密度的收敛，因此 计算时间会延长，收敛难度会增加。 )()()(rnrnrn rrr )(),(rnrnvv xcx

3、c rr )()()(rnrnr rrr rdrnrnM 3 )()( rr Dmol3的基本原理和参数设置 )( rc j j iji 适合于分子、团簇、分子筛、分子晶体、聚合物等“开放类结构” Tips：对于空体积较大的晶体，使用DMol3的效率要高于Castep Rcut 周期性和非 周期性体系 lm lmnlj YrRr),()()( Radial portion atomic DFT eqs. numerically Angular Portion 在DMol3模块中，电子密度实际上由各个原子轨道的平方和来确定： 在这里，电子密度实际上是由所有占据的分子轨道i.来决定。分子轨道可 能

4、由上自旋电子(Alpha电子)和下自旋电子(Beta电子)占据。当Alpha电子和 Beta电子的数目相等的时候，我们可以用单一的分子轨道i.来进行表述，这类 体系称为闭壳层体系(Closed-shell)，在DMol3中不需要选中Spin restricted前面 的选项。当Alpha电子和Beta电子的数目不相等的时候，我们将会使用不同的 i.来表述Alpha电子和Beta电子，这类体系称之为开壳层体系(Opened-shell)或 者自旋极化。在DMol3计算的时候需要选中Spin restricted的选项，并指定自旋 数目。 在开壳层体系中，会有两个不同的电子密度：一个是Alpha电

5、子的电子密 度，一个是Beta电子的电子密度。它们的和就是整个体系的总电荷，它们的差 就是自旋密度。 betaalphatotal betaalphaspin 其中，动能项的方程为： 注意，动能项实际上是一个常数项，在第一次计算完成后，该数值基本上可 以确定，后继计算中，则可以忽略这一步骤。 势能项处理： 当我们确定了电子密度之后，传统的薛定谔方程会从对电子波函数的处理转 换为对电子密度进行处理的函数。 需要注意的是，在势能项中的电子-电子相互作用，指的是两个电子间的相互作用。 但是，在整个体系中，还有三电子、四电子之间的相互作用，这一部分的内容从数学上 是没有办法得到精确解的，在密度泛函理论

6、中，将这一部分的内容归入了Exc这一项。 对这一项的处理，才是密度泛函理论处理的核心。针对不同的体系，有LDA和 GGA两种处理方法。 LDA (Local density approximation)局域密度近似方法假定在原子尺度电子密度变 化非常缓慢，也就是说，在整个分子区域内，整个体系表现为连续的电子气状态。 那么，整个电子交换-相关能就可以表示为对整个电子气的积分。 在DMol3模块中，常用的两种LDA方法是VWN和PWC： VWN：最常用的LSD (Local spin density)相关势函数。用来拟和电子气的精确数 值结果。 PWC：近期发展PWC泛函是在对VWN泛函的某些错误

7、校正后的结果，是DMol3模 块的默认泛函。 LSD方法可以精确预测共价体系共价体系的结构预测、频率计算和相关能量。但是，键能 往往会高估。LDA方法不能用于处理弱健体系不能用于处理弱健体系，如氢键。 LDA的这些缺陷，可以使用更大展开的Exc处理来校正。也称之为梯度校正方法。 GGA (General gradient-corrected)也成为NLSD(Non-local spin density)方法，近 二十年来的计算工作表明，使用梯度校正交换-相关能Exc, d ()可以很好的描述分子 体系的热力学性质。 需要注意的是，GGA方法实际上是一种经验性的描述函数。对于不同的研究体系， 计

8、算所使用的泛函的精确度实际上是不一样的。如果要求计算的结果准确可靠，则需要 对相关的函数进行查阅，或者从相关文献中查找所使用的泛函。 P91, BP, BLYP, BOP：也称为广义梯度近似方法。一般是Becke交换函数(B88)与 Perdew-Wang相关函数 (BP)或者Lee-Yang-Parr相关函数(BLYP)组合使用。 PBE：PBE (Perdew, Burke和Enzerhof) 泛函具有较强的物理背景，主要用于固体 计算，可靠的数值计算性能，在DFT计算中经常被使用。交换项与Becke相类似，相关 项与Pedew-Wang函数相接近。 RPBE：在PBE的基础上修改得到，对

9、热力学计算结果较为可靠。 HCTH： 对无机物和氢键体系的热力学计算结果较为可靠。 VWN-BP：COSMO或者COSMO-RS计算推荐使用。 尽管NLSD计算比起LSD计算有很大的改善，但是，在计算反应能垒的时候，得到 的能垒数据和实验结果相比仍任会有所低估。 PW91 Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992). BOP Tsuneda, T.; Suzumura, T.; Hirao, K. J. Chem. Phys., 110, 10664 (1999). BP Becke, A. D. J. Chem. Phys.,

10、 88, 2547 (1988). Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992). VWN-BP Vosko, S. J.; Wilk, L.; Nusair, M. Can. J. Phys., 58, 1200-1211 (1980). Becke, A. D. J. Chem. Phys., 88, 2547 (1988). Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992). PBE Perdew, J. P.; Burke, K.; Ernzerhof, M. P

11、hys. Rev. Lett., 77, 3865 (1996). RPBE Hammer, B.; Hansen, L. B.; Norskov, J. K. Phys. Rev. B, 59, 7413 (1999). BLYP Becke, A. D. J. Chem. Phys., 88, 2547 (1988). Lee, C.; Yang, W. Parr, R. G. Phys. Rev. B, 37, 786 (1988). HCTH Boese, A. D.; Handy, N. C. J. Chem. Phys., 114, 5497 (2001). 当我们确定了如何处理电

12、子的交换-相关能作用后，我们就可以将初始的 薛定谔方程进行转换，转换后的Kohn-Sham方程就可以用下式来进行描述： 这里的Et是整个分子的总能量。而如果要精确的确定整个分子的总能量， 就必须通过改变体系的电子密度，并由此得到不同的能量数值，当体系的能 量不再发生变化的时候，就可以认为整个体系趋于稳定，最终得到的能量， 就是整个体系的最终能量。 获得最终能量的一系列计算工作我们称之为SCF (Self-Consistent field自洽反应场)计算，整个SCF迭代的计算流程，如下图所示： 1. 选择初始的Ci。 2. 根据Ci构建初始的分子轨道i。 3. 构建电子密度。 4. 使用电子密度

13、，计算电子间的势能项，并考察Exc。 5. 计算整个体系的哈密顿能量。 6. 解自洽方程以获得新的一组Ci系数。 7. 构建新的分子轨道i和电子密度 8. 如果new = old，那么通过新的薛定谔方程计算总能量，并结 束。 9. 如果new old，那么返回步骤4。 对于有机分子，一般10部左右可以保证整个new old 10-6，整个体系能量收敛。但是对于金属体系，一般需要更多的迭 代来确保体系收敛。 DMol3中的中的SCF迭代过程迭代过程 Task: Energy, geometry optimization, TS-search, Overall quality control XC

14、-Option: LDA, GGA Spin option Charged systems (add or remove electrons) Task： Energy：直接进行SCF迭代计算，迭代的结果就是体系的总能量。 Geometry Optimization：求解体系的势能与坐标的一阶导数， 当该数值为0时，各个原子受力为零，此时的结构为能量最低结构。 Dynamics：考虑为各个原子分配速率，模拟体系在外界温度场作 用下的动力学行为； TS Search、TS Optimization、TS Conformation： 对反应体系的计算工作。TS Search是寻找反应过程中的过渡态

15、结 构，并给出反应的能垒和反应热；TS Optimization是对过渡态结构 进行进一步处理，寻找真实的一阶鞍点；TS Conformation则是对搜 索到的反应过渡态进行进一步确认，搜索可能的中间体和过渡态。 结构优化结构优化 在DMol3计算中，采用的 是BFGS计算方法。BFGS计算 是进行结构优化的标准算法。 此方法适合于获得能量最低的 结构优化算法。 Tips：单击More后，可以看到 Optimize Cell选项，只有当整个晶 胞的晶胞参数都需要优化时，才 可以选择此选项。孤立体系无此 选项。 在结构优化中，可以勾中 Use starting Hessian来指定各个 原子的

16、受力方向。 Integration: 各个基函数的 精度 SCF tolerance：SCF迭代误 差 K-points ：布里渊区采样 Core treatment (all electron, PP,) Basis set：基组选择 Orbital cutoff quality： 各个原子轨道的半径设置 More All Electron 不对内核电子进行特殊处理，所有的电子都被包含在计 算体系中进行处理。 Effective Core Potentials (ECP) 使用单个有效势替代内核电子，以缩 减计算成本。ECPs会在内核处理中引入相对论校正。 All Electron Rela

17、tivistic 处理体系中的所有电子，并在内核电子的处 理中引入相对论效应。此方法更为精确，但是计算成本最长。 DFT Semi-core Pseudopots (DSPP) 使用单个有效势替代内核电子， 以缩减计算成本。DSPPS会在内核处理中引入相对论校正。 注意：ECP和DSPP都是对21号以后的重元素进行处理，如果体 系内包含C、H、O、Zn，那么程序将只会考虑Zn的相对论效应， 其它元素则使用全电子进行处理。其中，DSPP特别针对DMol3 模块开发，而ECP则来源于Hartree-Fock势。 Core treatment参数控制那些原子轨道上的电子需要进行处 理。默认的设置是对

18、于所有的电子进行处理。对于重原子而言， 内核电子的速率接近光速，就必须考虑到相对论效应。 Min Min 一个占位原子轨道之使用一个数值轨道基组。速度更快，精度 最低。 DNDN 双数值轨道基组。在Min基组的基础上加入了第二组价电子原 子轨道。 DNDDND 双数值轨道基组+d轨道极化函数。与DN基组类似，但是为所有 非氢原子加入了d轨道函数极化。为DMol3模块默认基组，在可接受的 计算时内，确保精度在可接受的范围内。也是DMol3计算的最低可信精 度基组。 DNPDNP 双数值轨道基组+轨道极化函数。与DND基组类似，但是为所有 氢原子加入了p轨道函数极化。 精度更高，计算时间更长。对于

19、氢键 计算更为重要。 TNP TNP 三三数值轨道基组+轨道极化函数。 与DNP基组类似，为所有原 子加入极化函数。目前只能够使用与H到Cl的元素，且不包括He和Ne。 精度更高，计算时间更长。 从理论上讲，体系的电子密度和相关函数都是针对整个空间来进行处 理的。由于电子密度在远离原子核的时候迅速衰减，因此，可以限制积 分的范围。这一设置能够大大减少计算的时间，但不会影响到计算的精 度。注意，全局半径要按照体系内具有最大半径的元素标准进行设置。 SCF tolerance：SCF迭代中， 两次能量之间的差值，以电子密 度的差异为评判标准 Smearing: 选中后，意味着将热 拖尾效应应用到轨

20、道占据上，来 加速收敛速度(SCF迭代失败后， 第一项需要调的参数) Density mixing (SCF迭代失败后， 第二项需要调的参数) 1 1 )( Tk e f B )(f Smearing参数允许电子在所有轨道中按 照指定的能量差E 进行拖尾。类似于 物理上的热占位现象。此方法能够通过 允许轨道驰豫而大大加速SCF迭代的收 敛速度。会导致虚轨道与占据轨道进行 混合，因此，会有一些轨道出现分数占 位。 初始结构 红色为未占据轨道 “Made visible” by smearing 默认值默认值! 最终结构 红色为发生占据 Density Mixing参数控制体系中，如何根据特征方程

21、来如何构造新的 电子密度。在整个体系中，通过加入阻尼振荡来确保整个体系的平滑收敛。 对于一个最简单的阻尼方法，遵循如下方程： A就是阻尼系数。in是当前前一步SCF迭代计算得到的电子密度，根 据这一项来构建特征矩阵。out是没有考虑阻尼系数，根据新的分子轨道 系数构建得到的电子密度。而new则是真正用于下一步迭代的电子密度。 DIIS(direct inversion in an iterative subspace) ： 也称为Pulay mixing方法，由Peter Pulay发展得到，是一种外推方法。 用于Harter-Fock SCF迭代方法的加速并稳定体系的收敛过程中。 当选中时，

22、能够加速整个体系的SCF迭代速度。 inoutinnew A The DIIS method assumes that a good approximation to the final solution pf can be obtained as a linear combination of the previous guess vectors where m is the number of previous vectors (in practice, only the most recent few vectors are used). The coefficients ci are

23、obtained by requiring that the associated residual vector approximates the zero vector in a least-squares sense. Furthermore, the coefficients are required to add to one, 专门针对周期性体系 需要针对不同布里渊区 的K点路径进行设置(如 果不确定，可以使用程 序默认参数) K点路径间隔或者格点 设置(同一系列体系需 要使用相同的格点间隔) 在数学和固体物理学中，第一布里渊布里渊 区区（Brillouin zone）是动量空间动

24、量空间中晶体倒倒 易点阵易点阵的原胞。 在点阵点阵空间中，作某一个阵点与其所 有相邻阵点的垂直平分面，这些平面包围 的空间就是包含前述阵点的第一布里渊区； 亦可等价地定义为：在k空间（即波矢空间 或倒易空间倒易空间）中，从原点出发，不穿越任 何布拉格衍射面布拉格衍射面所能到达的点的集合，就 是第一布里渊区。 FCC晶体的第一布里渊区晶体的第一布里渊区 与计算的全局精度设置 相关 可以按照globally或者 element specific来进 行选择 需要调整不同的Rcut设 置来验证Oribtal Cutoff是否合理 DMOL: 参数设置 Orbital Cutoff Band Stru

25、cture Density of states Electron density Electrostatics Frequency Fukui function Orbitals Population analysis Thermodynamic Properties Electron Density Total density: 整个体系的电子密度 Deformation density: 与单个原子的电子密度相比较，整个分子的电子密 度变化值。计算方法是总电子密度减去单个原子的密度差。 Spin density: -自旋 和 b-自旋电子的电子密度差。 Electrostatics Ele

26、ctrostatic potential Electrostatic moments: dipole, quadrupole, octupole, and hexadecapole Nuclear electric field gradients FFs 提供了亲电亲电或者亲核亲核攻击的相关信息 那个区域反应更易发生？ FF可以绘制出3D图，以显示反应区域。 这些类似结构，那个更有反应性？ 浓缩的FFs显示方式，可以以原子为中心的等势面，定量对比不 同区域的反应性能。 )()( 1 rrrf NNN rrr )()( 1 rrrf NNN rrr 亲核攻击 亲电攻击 Od d- attack

27、Ald d+ attack Orbitals - 预测反应性 对于使用数值轨道的DMol3模块，在计算后能够给出整个体系 的分子轨道信息。对于单个分子而言，这些分子能够帮助我们判断 电子的可能跃迁(需要用到群论知识)；对于反应而言，也可以使用 福井谦一的前线轨道对称守恒原理，判断反应发生的区域。 LUMO HOMO .where reaction takes place! Frequency - 热力学性质计算 分子的振动频率在实验上是通过IR和Raman谱图来进行观测的。这些技 术能够帮助我们观察到分子的结构和化学环境。不同的振动模式对应了不同的 峰。一般而言，键的伸缩具有最高能量的振动频率

28、，而弯曲振动的能量较低， 扭转振动会更低一些。能量最低的那些振动方式，对应了刚性体系的扭动或者 大分子的呼吸效应。 DMol3通过对能量梯度的有限差值分析来获得相应的频率。通过对能量对 位移的一阶导数分析，可以得到整个体系中各个原子的受力状况，而二阶导数 分析则对应了相应的振动模式。 谱图中峰的强度，对应了特定波长的光子被吸收或者Raman散射的几率。 在量子力学计算中，需要通过对跃迁偶极矩跃迁偶极矩的计算得到的波函数进行分析得到 这种几率数值。这种计算得到的强度往往只是相对数值，且具有较大的误差。 经过计算得到的能量，往往与真实能量有一定的误差。这是由于真实体 系往往都具有一定的振动模式。通

29、过频率计算后，可以将得到的最低能量振动 得到的零点能校正零点能校正加入到之前计算的能量数值来获得较为可靠的能量数值。而 体系的熵、焓数值，则可以分解为由分子平动、振动和转动带来的能量变化。 Transition State Search TS Search TS Optimization TS Conformation Reaction Coordinate Energy With catalyst Z1 Z2 PG RG A* P* Ea,0 Ea,1Ea,2 HR Without catalyst 局域法: (ScanPESScanPES) 从单一结构出发 计算势能面 TS 计算较为准确

30、插值方法: () 需要给定初始反应物和产物结构 计算单个结构的能量和能量梯度 一般需要作更进一步的验证 LST line search maximum QST line search maximum A schematic view of the LST and QST approach LST Maximum：执行单个LST极大值化计算。找到反应物和产物之 间的能量最高点结构。 Halgren-Lipscomb：执行单个LST极大值化计算，并使用单线搜索最 小点方法进行后继处理。 LST/Optimization：执行单个LST极大值化计算，并使用完全共轭梯 度法对结构进行后继处理。 Com

31、plete LST/QST：执行单个LST极大值化计算，并重复使用共轭 梯度法优化和QST极大值方法进行后继处理，直到确定过渡态结构。 QST / Optimization：从QST路径出发，重复使用共轭梯度法优化和 QST极大值方法进行后继处理，直到确定过渡态结构。 Reaction Path Energy (Ha) LST Maximization Conjugate Gradient Minimization Conjugate gradient refinement of saddle point is important! 0.6eV Reactant Product Halgren

32、 Libscomb 真实的过渡态结构应该 是一个鞍点，反映在频 率计算上，应该只有一 个虚频率 虚频率的振动模式，应 当是在反应物和产物之 间发生震荡，由此保证 反映过程的一致性。 cm-1-189164474118166169190213228248273284313321 cm-1618657675717770794824838858869911945988999 1033 cm-11348 1363 1369 1379 1451 1456 1467 1471 1480 1892 2936 2952 2968 2976 2983 当过渡态结构中存在多个虚频，有时候就需要使用TS Optim

33、ization方法将多 余的虚频消除，以获得真正的一阶鞍点。 在使用TS Optimization的时候，首先需要使用Tools | Vibrational Analysis对整个体系的 Hessian矩阵进行计算处理。 在设置参数的时候，需要将Vibrational Analysis的窗口激活，并选中最小的正频率，之 后开始DMol3计算。 当体系发生震荡的时候，可以将Max. Step size的数值调小，以确保收敛。 A M1 M2 B TS 1 TS 2 A B involves 2 unexpected minima 能否在反应物、产物和TS 2之间找到其他能量 最小的结构？ A M

34、1 M2 B TS 1 TS 2 传统方法: 由TS 出发 执行最陡下降法计算 简单的最陡下降法优化 阻尼MD计算 要求大量的步骤 二阶导数法 要求进行大量的Hessian矩阵计算 在第一个找到的最小点处结束 NEB方法. 搜索完整的反应路径 需要计算的步数少 不需要计算Hessian矩阵 依赖于离域内坐标优化和QST算法 NEB 思路描述 C B A 收敛到“正确答案”！ 由TS Search计算得到的QST路径出 发 在QST路径上设定“映像结构” C B A QST IRP 考察反应路径 A C 反应路径未知 A B C 1) 将“映像结构”相关联 垂直于路径方向进行结构优化 C B A

35、 2) TS Conformation 在执行TS Conformation的时候，需要注意，一定要将TS Search 中得到的xtd轨迹激活，并在这一窗口下执行TS Conformational工作。 执行TS Search计算工作 在TS Search得到的xtd路径上，进行TS Conformaiton工作 完成后会给出IRP路径，可能有一些局域最小点 对新找到的最小点结构进行优化 如果必须，在新得到的最小点结构之间执行TS Search计算工 作 CASTEP的基本原理和参数设置 CASTEP将电子波函数通过平面 波函数展开 G r)Gk( , ) r ( i Gki eC 赝势 p

36、seudo wave function pseudopotential 有效的减少平面波数目 只适用于周期性体系 计算任务(Task) Energy, Geometry optimization, Elastic constant, TS search, Properties 精度控制(Quality) 交换相关函数(Functional) LDA: CA-PZ GGA: PBE, RPBE, PW91, WC, PBEsol B3LYP 0PBE HF, HF-LDA sX, sX-LDA 自旋极化(Spin polarized) 初始自旋值(Initial spin) 体系电荷 动能截断值

37、(Energy cutoff) 自洽计算收敛精度(SCF tolerance) K点设置(k- point set) 赝势(Pseudopotential) Ultrasoft, Normconserving, On the fly 赝势的表达，实空间 or 倒空间 Core hole的选择 More Basis 基组的相关参数设置 动能截断值的设定 FFT格子设定 基组的测试(Finite basis correction ) SCF 自洽计算的相关参数设置 k-points Potential 选择赝势 Basis 基组的相关参数设置 SCF 自洽计算的相关参数设置 收敛精度的具体设定 自洽计算允许的步数 Density mixing 轨道占据数 初始自旋值的优化 k-points Potential 选择赝势 Density Mixing参数控制体系中，如何根据特征方程来如何构造新的 电子密度。在整个体系中，通过加入阻尼振荡来确保整个体系的平滑收敛。 对于一个最简单的阻尼方法，遵循如下方程： A就是阻尼系数。in是当前前一步SCF迭代计算得到的电子密度，根 据这一项来构建特征矩阵。out是没有考虑阻尼系数，根据新的分子轨道 系数构建得到的电子密度。而new则是真正用于下一步迭代的电子密度。 DIIS(direct inversion i