DMol3基本原理和參數(shù)設置

1、DMol3 基本原理和參數(shù)設置SCF(DFT)制備材料表征材料各種性能參數(shù)制備工藝結構能量電子波函數(shù)能帶，態(tài)密度，電荷密度一階導數(shù)應力，彈性力常數(shù)，體模量.二階導數(shù)聲子頻率，散射譜H = E力場怎么使SCF收斂？如何使SCF的結果準確？Walter KohnE To U Exc )()()(,2rrrvkikikieff)(),.(1rErrENHohenberg-Kohn theoremKohn-Sham equations Exact only for ground state Needs approximation to ExcKohn-Sham theorem生成KS勢求解KS方程得到

2、新nout(r)和之前的n(r)比較收斂與否 ?nstart(r) 輸出結果)()()(,2rrrnvkikikieffrrrrrr)(rnveffr rrfrniiiirrr電荷密度依據(jù)電子自旋方向的不同，一分為二：自旋密度定義為:總的磁矩定義為:存在自旋極化的體系中，交換相關勢也相應的變?yōu)椴煌孕较螂姾擅芏鹊暮瘮?shù):PS： 由于在自洽計算中需要分別考慮電荷密度和自旋密度的收斂，因此計算時間會延長，收斂難度會增加。)()()(rnrnrnrrr)(),(rnrnvvxcxcrr)()()(rnrnrrrrrdrnrnM3)()(rrDmol3的基本原理和參數(shù)設置)(rcjjiji適合于分子、

3、團簇、分子篩、分子晶體、聚合物等“開放類結構”Tips：對于空體積較大的晶體，使用DMol3的效率要高于CastepRcut周期性和非周期性體系lmlmnljYrRr),()()(Radial portion atomic DFT eqs. numericallyAngular Portion 在DMol3模塊中，電子密度實際上由各個原子軌道的平方和來確定： 在這里，電子密度實際上是由所有占據(jù)的分子軌道i.來決定。分子軌道可能由上自旋電子(Alpha電子)和下自旋電子(Beta電子)占據(jù)。當Alpha電子和Beta電子的數(shù)目相等的時候，我們可以用單一的分子軌道i.來進行表述，這類體系稱為閉殼層

4、體系(Closed-shell)，在DMol3中不需要選中Spin restricted前面的選項。當Alpha電子和Beta電子的數(shù)目不相等的時候，我們將會使用不同的i.來表述Alpha電子和Beta電子，這類體系稱之為開殼層體系(Opened-shell)或者自旋極化。在DMol3計算的時候需要選中Spin restricted的選項，并指定自旋數(shù)目。 在開殼層體系中，會有兩個不同的電子密度：一個是Alpha電子的電子密度，一個是Beta電子的電子密度。它們的和就是整個體系的總電荷，它們的差就是自旋密度。betaalphatotalbetaalphaspin其中，動能項的方程為：注意，動能

5、項實際上是一個常數(shù)項，在第一次計算完成后，該數(shù)值基本上可以確定，后繼計算中，則可以忽略這一步驟。勢能項處理： 當我們確定了電子密度之后，傳統(tǒng)的薛定諤方程會從對電子波函數(shù)的處理轉換為對電子密度進行處理的函數(shù)。 需要注意的是，在勢能項中的電子-電子相互作用，指的是兩個電子間的相互作用。但是，在整個體系中，還有三電子、四電子之間的相互作用，這一部分的內容從數(shù)學上是沒有辦法得到精確解的，在密度泛函理論中，將這一部分的內容歸入了Exc這一項。 對這一項的處理，才是密度泛函理論處理的核心。針對不同的體系，有LDA和GGA兩種處理方法。 LDA (Local density approximation)局域

6、密度近似方法假定在原子尺度電子密度變化非常緩慢，也就是說，在整個分子區(qū)域內，整個體系表現(xiàn)為連續(xù)的電子氣狀態(tài)。 那么，整個電子交換-相關能就可以表示為對整個電子氣的積分。 在DMol3模塊中，常用的兩種LDA方法是VWN和PWC： VWN：最常用的LSD (Local spin density)相關勢函數(shù)。用來擬和電子氣的精確數(shù)值結果。 PWC：近期發(fā)展PWC泛函是在對VWN泛函的某些錯誤校正后的結果，是DMol3模塊的默認泛函。 LSD方法可以精確預測共價體系共價體系的結構預測、頻率計算和相關能量。但是，鍵能往往會高估。LDA方法不能用于處理弱健體系不能用于處理弱健體系，如氫鍵。 LDA的這些

7、缺陷，可以使用更大展開的Exc處理來校正。也稱之為梯度校正方法。 GGA (General gradient-corrected)也成為NLSD(Non-local spin density)方法，近二十年來的計算工作表明，使用梯度校正交換-相關能Exc, d ()可以很好的描述分子體系的熱力學性質。 需要注意的是，GGA方法實際上是一種經(jīng)驗性的描述函數(shù)。對于不同的研究體系，計算所使用的泛函的精確度實際上是不一樣的。如果要求計算的結果準確可靠，則需要對相關的函數(shù)進行查閱，或者從相關文獻中查找所使用的泛函。 P91, BP, BLYP, BOP：也稱為廣義梯度近似方法。一般是Becke交換函數(shù)(

8、B88)與Perdew-Wang相關函數(shù) (BP)或者Lee-Yang-Parr相關函數(shù)(BLYP)組合使用。 PBE：PBE (Perdew, Burke和Enzerhof) 泛函具有較強的物理背景，主要用于固體計算，可靠的數(shù)值計算性能，在DFT計算中經(jīng)常被使用。交換項與Becke相類似，相關項與Pedew-Wang函數(shù)相接近。 RPBE：在PBE的基礎上修改得到，對熱力學計算結果較為可靠。 HCTH： 對無機物和氫鍵體系的熱力學計算結果較為可靠。 VWN-BP：COSMO或者COSMO-RS計算推薦使用。 盡管NLSD計算比起LSD計算有很大的改善，但是，在計算反應能壘的時候，得到的能壘數(shù)

9、據(jù)和實驗結果相比仍任會有所低估。PW91Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992).BOPTsuneda, T.; Suzumura, T.; Hirao, K. J. Chem. Phys., 110, 10664 (1999).BPBecke, A. D. J. Chem. Phys., 88, 2547 (1988).Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992).VWN-BPVosko, S. J.; Wilk, L.; Nusair, M. Can. J.

10、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).PBEPerdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett., 77, 3865 (1996).RPBEHammer, B.; Hansen, L. B.; Norskov, J. K. Phys. Rev. B, 59, 7413 (1999).BLYPBecke, A. D. J

11、. Chem. Phys., 88, 2547 (1988). Lee, C.; Yang, W.Parr, R. G. Phys. Rev. B, 37, 786 (1988).HCTHBoese, A. D.; Handy, N. C. J. Chem. Phys., 114, 5497 (2001). 當我們確定了如何處理電子的交換-相關能作用后，我們就可以將初始的薛定諤方程進行轉換，轉換后的Kohn-Sham方程就可以用下式來進行描述： 這里的Et是整個分子的總能量。而如果要精確的確定整個分子的總能量，就必須通過改變體系的電子密度，并由此得到不同的能量數(shù)值，當體系的能量不再發(fā)生變化的時

12、候，就可以認為整個體系趨于穩(wěn)定，最終得到的能量，就是整個體系的最終能量。 獲得最終能量的一系列計算工作我們稱之為SCF (Self-Consistent field自洽反應場)計算，整個SCF迭代的計算流程，如下圖所示：1. 選擇初始的Ci。 2. 根據(jù)Ci構建初始的分子軌道i。3. 構建電子密度。4. 使用電子密度，計算電子間的勢能項，并考察Exc。5. 計算整個體系的哈密頓能量。6. 解自洽方程以獲得新的一組Ci系數(shù)。7. 構建新的分子軌道i和電子密度 8. 如果new = old，那么通過新的薛定諤方程計算總能量，并結束。9. 如果new old，那么返回步驟4。 對于有機分子，一般10

13、部左右可以保證整個new old 10-6，整個體系能量收斂。但是對于金屬體系，一般需要更多的迭代來確保體系收斂。DMol3中的中的SCF迭代過程迭代過程Task: Energy, geometry optimization, TS-search, Overall quality controlXC-Option: LDA, GGASpin optionCharged systems (add or remove electrons)Task： Energy：直接進行SCF迭代計算，迭代的結果就是體系的總能量。 Geometry Optimization：求解體系的勢能與坐標的一階導數(shù)，當該數(shù)

14、值為0時，各個原子受力為零，此時的結構為能量最低結構。 Dynamics：考慮為各個原子分配速率，模擬體系在外界溫度場作用下的動力學行為； TS Search、TS Optimization、TS Conformation： 對反應體系的計算工作。TS Search是尋找反應過程中的過渡態(tài)結構，并給出反應的能壘和反應熱；TS Optimization是對過渡態(tài)結構進行進一步處理，尋找真實的一階鞍點；TS Conformation則是對搜索到的反應過渡態(tài)進行進一步確認，搜索可能的中間體和過渡態(tài)。結構優(yōu)化結構優(yōu)化 在DMol3計算中，采用的是BFGS計算方法。BFGS計算是進行結構優(yōu)化的標準算法。

15、此方法適合于獲得能量最低的結構優(yōu)化算法。Tips：單擊More后，可以看到Optimize Cell選項，只有當整個晶胞的晶胞參數(shù)都需要優(yōu)化時，才可以選擇此選項。孤立體系無此選項。 在結構優(yōu)化中，可以勾中Use starting Hessian來指定各個原子的受力方向。Integration: 各個基函數(shù)的精度SCF tolerance：SCF迭代誤差K-points ：布里淵區(qū)采樣Core treatment (all electron, PP,)Basis set：基組選擇Orbital cutoff quality：各個原子軌道的半徑設置More All Electron 不對內核電子進

16、行特殊處理，所有的電子都被包含在計算體系中進行處理。Effective Core Potentials (ECP) 使用單個有效勢替代內核電子，以縮減計算成本。ECPs會在內核處理中引入相對論校正。All Electron Relativistic 處理體系中的所有電子，并在內核電子的處理中引入相對論效應。此方法更為精確，但是計算成本最長。 DFT Semi-core Pseudopots (DSPP) 使用單個有效勢替代內核電子，以縮減計算成本。DSPPS會在內核處理中引入相對論校正。 注意：ECP和DSPP都是對21號以后的重元素進行處理，如果體系內包含C、H、O、Zn，那么程序將只會考慮

17、Zn的相對論效應，其它元素則使用全電子進行處理。其中，DSPP特別針對DMol3模塊開發(fā)，而ECP則來源于Hartree-Fock勢。 Core treatment參數(shù)控制那些原子軌道上的電子需要進行處理。默認的設置是對于所有的電子進行處理。對于重原子而言，內核電子的速率接近光速，就必須考慮到相對論效應。 Min Min 一個占位原子軌道之使用一個數(shù)值軌道基組。速度更快，精度最低。DNDN 雙數(shù)值軌道基組。在Min基組的基礎上加入了第二組價電子原子軌道。DNDDND 雙數(shù)值軌道基組+d軌道極化函數(shù)。與DN基組類似，但是為所有非氫原子加入了d軌道函數(shù)極化。為DMol3模塊默認基組，在可接受的計算

18、時內，確保精度在可接受的范圍內。也是DMol3計算的最低可信精度基組。DNPDNP 雙數(shù)值軌道基組+軌道極化函數(shù)。與DND基組類似，但是為所有氫原子加入了p軌道函數(shù)極化。 精度更高，計算時間更長。對于氫鍵計算更為重要。TNP TNP 三三數(shù)值軌道基組+軌道極化函數(shù)。 與DNP基組類似，為所有原子加入極化函數(shù)。目前只能夠使用與H到Cl的元素，且不包括He和Ne。 精度更高，計算時間更長。 從理論上講，體系的電子密度和相關函數(shù)都是針對整個空間來進行處理的。由于電子密度在遠離原子核的時候迅速衰減，因此，可以限制積分的范圍。這一設置能夠大大減少計算的時間，但不會影響到計算的精度。注意，全局半徑要按照體

19、系內具有最大半徑的元素標準進行設置。SCF tolerance：SCF迭代中，兩次能量之間的差值，以電子密度的差異為評判標準Smearing: 選中后，意味著將熱拖尾效應應用到軌道占據(jù)上，來加速收斂速度(SCF迭代失敗后，第一項需要調的參數(shù))Density mixing (SCF迭代失敗后，第二項需要調的參數(shù))11)(TkefB)(fSmearing參數(shù)允許電子在所有軌道中按照指定的能量差E 進行拖尾。類似于物理上的熱占位現(xiàn)象。此方法能夠通過允許軌道馳豫而大大加速SCF迭代的收斂速度。會導致虛軌道與占據(jù)軌道進行混合，因此，會有一些軌道出現(xiàn)分數(shù)占位。初始結構 紅色為未占據(jù)軌道“Made visi

20、ble” by smearing 默認值默認值!最終結構 紅色為發(fā)生占據(jù) Density Mixing參數(shù)控制體系中，如何根據(jù)特征方程來如何構造新的電子密度。在整個體系中，通過加入阻尼振蕩來確保整個體系的平滑收斂。對于一個最簡單的阻尼方法，遵循如下方程： A就是阻尼系數(shù)。in是當前前一步SCF迭代計算得到的電子密度，根據(jù)這一項來構建特征矩陣。out是沒有考慮阻尼系數(shù)，根據(jù)新的分子軌道系數(shù)構建得到的電子密度。而new則是真正用于下一步迭代的電子密度。DIIS(direct inversion in an iterative subspace) ： 也稱為Pulay mixing方法，由Peter

21、 Pulay發(fā)展得到，是一種外推方法。用于Harter-Fock SCF迭代方法的加速并穩(wěn)定體系的收斂過程中。 當選中時，能夠加速整個體系的SCF迭代速度。inoutinnewAThe 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 m

22、ost recent few vectors are used). The coefficients ci are obtained by requiring that the associated residual vectorapproximates the zero vector in a least-squares sense. Furthermore, the coefficients are required to add to one, 專門針對周期性體系需要針對不同布里淵區(qū)的K點路徑進行設置(如果不確定，可以使用程序默認參數(shù))K點路徑間隔或者格點設置(同一系列體系需要使用相同的

23、格點間隔) 在數(shù)學和固體物理學中，第一布里淵布里淵區(qū)區(qū)（Brillouin zone）是動量空間動量空間中晶體倒倒易點陣易點陣的原胞。 在點陣點陣空間中，作某一個陣點與其所有相鄰陣點的垂直平分面，這些平面包圍的空間就是包含前述陣點的第一布里淵區(qū)；亦可等價地定義為：在k空間（即波矢空間或倒易空間倒易空間）中，從原點出發(fā)，不穿越任何布拉格衍射面布拉格衍射面所能到達的點的集合，就是第一布里淵區(qū)。 FCC晶體的第一布里淵區(qū)晶體的第一布里淵區(qū)與計算的全局精度設置相關可以按照globally或者element specific來進行選擇需要調整不同的Rcut設置來驗證Oribtal Cutoff是否合理D

24、MOL: 參數(shù)設置Orbital CutoffBand StructureDensity of statesElectron densityElectrostaticsFrequencyFukui function OrbitalsPopulation analysis Thermodynamic PropertiesElectron DensityTotal density: 整個體系的電子密度Deformation density: 與單個原子的電子密度相比較，整個分子的電子密度變化值。計算方法是總電子密度減去單個原子的密度差。Spin density: -自旋 和 b-自旋電子的電子密度

25、差。ElectrostaticsElectrostatic potentialElectrostatic moments: dipole, quadrupole, octupole, and hexadecapoleNuclear electric field gradientsFFs 提供了親電親電或者親核親核攻擊的相關信息那個區(qū)域反應更易發(fā)生？ FF可以繪制出3D圖，以顯示反應區(qū)域。 這些類似結構，那個更有反應性？濃縮的FFs顯示方式，可以以原子為中心的等勢面，定量對比不同區(qū)域的反應性能。 )()(1rrrfNNNrrr )()(1rrrfNNNrrr親核攻擊親電攻擊Od d- attac

26、kAld d+ attackOrbitals - 預測反應性 對于使用數(shù)值軌道的DMol3模塊，在計算后能夠給出整個體系的分子軌道信息。對于單個分子而言，這些分子能夠幫助我們判斷電子的可能躍遷(需要用到群論知識)；對于反應而言，也可以使用福井謙一的前線軌道對稱守恒原理，判斷反應發(fā)生的區(qū)域。LUMOHOMO.where reaction takes place!Frequency - 熱力學性質計算 分子的振動頻率在實驗上是通過IR和Raman譜圖來進行觀測的。這些技術能夠幫助我們觀察到分子的結構和化學環(huán)境。不同的振動模式對應了不同的峰。一般而言，鍵的伸縮具有最高能量的振動頻率，而彎曲振動的能量

27、較低，扭轉振動會更低一些。能量最低的那些振動方式，對應了剛性體系的扭動或者大分子的呼吸效應。 DMol3通過對能量梯度的有限差值分析來獲得相應的頻率。通過對能量對位移的一階導數(shù)分析，可以得到整個體系中各個原子的受力狀況，而二階導數(shù)分析則對應了相應的振動模式。 譜圖中峰的強度，對應了特定波長的光子被吸收或者Raman散射的幾率。在量子力學計算中，需要通過對躍遷偶極矩躍遷偶極矩的計算得到的波函數(shù)進行分析得到這種幾率數(shù)值。這種計算得到的強度往往只是相對數(shù)值，且具有較大的誤差。 經(jīng)過計算得到的能量，往往與真實能量有一定的誤差。這是由于真實體系往往都具有一定的振動模式。通過頻率計算后，可以將得到的最低能

28、量振動得到的零點能校正零點能校正加入到之前計算的能量數(shù)值來獲得較為可靠的能量數(shù)值。而體系的熵、焓數(shù)值，則可以分解為由分子平動、振動和轉動帶來的能量變化。Transition State Search TS Search TS Optimization TS ConformationReaction CoordinateEnergyWith catalystZ1Z2PGRGA*P*Ea,0Ea,1Ea,2HRWithout catalyst 局域法: (ScanPESScanPES) 從單一結構出發(fā)計算勢能面TS 計算較為準確 插值方法: ()需要給定初始反應物和產物結構 計算單個結構的能量和能

29、量梯度一般需要作更進一步的驗證LST line searchmaximumQST line search maximumA schematic view of the LST and QST approach LST Maximum：執(zhí)行單個LST極大值化計算。找到反應物和產物之間的能量最高點結構。 Halgren-Lipscomb：執(zhí)行單個LST極大值化計算，并使用單線搜索最小點方法進行后繼處理。 LST/Optimization：執(zhí)行單個LST極大值化計算，并使用完全共軛梯度法對結構進行后繼處理。 Complete LST/QST：執(zhí)行單個LST極大值化計算，并重復使用共軛梯度法優(yōu)化和QS

30、T極大值方法進行后繼處理，直到確定過渡態(tài)結構。 QST / Optimization：從QST路徑出發(fā)，重復使用共軛梯度法優(yōu)化和QST極大值方法進行后繼處理，直到確定過渡態(tài)結構。Reaction PathEnergy (Ha)LSTMaximizationConjugate Gradient MinimizationConjugate gradientrefinement of saddlepoint is important!0.6eVReactantProductHalgren Libscomb真實的過渡態(tài)結構應該是一個鞍點，反映在頻率計算上，應該只有一個虛頻率虛頻率的振動模式，應當是在反

31、應物和產物之間發(fā)生震蕩，由此保證反映過程的一致性。cm-1-189164474118166169190213228248273284313321cm-1618657675717770794824838858869911945988999 1033cm-11348 1363 1369 1379 1451 1456 1467 1471 1480 1892 2936 2952 2968 2976 2983 當過渡態(tài)結構中存在多個虛頻，有時候就需要使用TS Optimization方法將多余的虛頻消除，以獲得真正的一階鞍點。 在使用TS Optimization的時候，首先需要使用Tools | Vi

32、brational Analysis對整個體系的Hessian矩陣進行計算處理。 在設置參數(shù)的時候，需要將Vibrational Analysis的窗口激活，并選中最小的正頻率，之后開始DMol3計算。 當體系發(fā)生震蕩的時候，可以將Max. Step size的數(shù)值調小，以確保收斂。AM1M2BTS 1TS 2A B involves 2 unexpected minima能否在反應物、產物和TS 2之間找到其他能量最小的結構？AM1M2BTS 1TS 2傳統(tǒng)方法:由TS 出發(fā)執(zhí)行最陡下降法計算簡單的最陡下降法優(yōu)化阻尼MD計算 要求大量的步驟二階導數(shù)法 要求進行大量的Hessian矩陣計算在第

33、一個找到的最小點處結束NEB方法.搜索完整的反應路徑需要計算的步數(shù)少不需要計算Hessian矩陣依賴于離域內坐標優(yōu)化和QST算法NEB 思路描述CBA 收斂到“正確答案”！ 由TS Search計算得到的QST路徑出發(fā) 在QST路徑上設定“映像結構”CBAQSTIRP考察反應路徑 A C反應路徑未知 A B C1) 將“映像結構”相關聯(lián) 垂直于路徑方向進行結構優(yōu)化CBA2)TS Conformation 在執(zhí)行TS Conformation的時候，需要注意，一定要將TS Search中得到的xtd軌跡激活，并在這一窗口下執(zhí)行TS Conformational工作。執(zhí)行TS Search計算工作

34、在TS Search得到的xtd路徑上，進行TS Conformaiton工作完成后會給出IRP路徑，可能有一些局域最小點對新找到的最小點結構進行優(yōu)化如果必須，在新得到的最小點結構之間執(zhí)行TS Search計算工作CASTEP的基本原理和參數(shù)設置CASTEP將電子波函數(shù)通過平面波函數(shù)展開Gr)Gk(,) r (iGkieC贗勢pseudo wave functionpseudopotential有效的減少平面波數(shù)目只適用于周期性體系計算任務(Task) Energy, Geometry optimization, Elastic constant, TS search, Properties精

35、度控制(Quality)交換相關函數(shù)(Functional) LDA: CA-PZ GGA: PBE, RPBE, PW91, WC, PBEsol B3LYP 0PBE HF, HF-LDA sX, sX-LDA自旋極化(Spin polarized) 初始自旋值(Initial spin)體系電荷動能截斷值(Energy cutoff)自洽計算收斂精度(SCF tolerance)K點設置(k- point set)贗勢(Pseudopotential) Ultrasoft, Normconserving, On the fly贗勢的表達，實空間 or 倒空間Core hole的選擇 Mo

36、reBasis 基組的相關參數(shù)設置 動能截斷值的設定 FFT格子設定 基組的測試(Finite basis correction ) SCF 自洽計算的相關參數(shù)設置k-points Potential 選擇贗勢Basis 基組的相關參數(shù)設置SCF 自洽計算的相關參數(shù)設置 收斂精度的具體設定 自洽計算允許的步數(shù) Density mixing 軌道占據(jù)數(shù) 初始自旋值的優(yōu)化k-points Potential 選擇贗勢 Density Mixing參數(shù)控制體系中，如何根據(jù)特征方程來如何構造新的電子密度。在整個體系中，通過加入阻尼振蕩來確保整個體系的平滑收斂。對于一個最簡單的阻尼方法，遵循如下方程： A就是阻尼系數(shù)。in是當前前一步SCF迭代計算得到的電子密度，根據(jù)這一項來構建特征矩陣。out是沒有考慮阻尼系數(shù)，根據(jù)新的分子軌道系數(shù)構建得到的電子密度。而new則是真正用于下一步迭代的電子密度。DIIS(direct inversion in an iterative subspace) ： 也稱為P