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

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

2、和之前的n(r)比較收斂與否 ?nstart(r) 輸出結果)()()(,2rrrnvkikikieffrrrrrrjej=+)(rnveffrDFT和磁學性能電荷密度依據電子自旋方向的不同，一分為二：自旋密度定義為:總的磁矩定義為:存在自旋極化的體系中，交換相關勢也相應的變?yōu)椴煌孕较螂姾擅芏鹊暮瘮?shù):PS： 由于在自洽計算中需要分別考慮電荷密度和自旋密度的收斂，因此計算時間會延長，收斂難度會增加。Dmol3的根本原理和參數(shù)設置DMol3: 原子軌道線形組合法(LCAO)適合于分子、團簇、分子篩、分子晶體、聚合物等“開放類構造Tips：對于空體積較大的晶體，使用DMol3的效率要高于Cast

3、epRcut周期性和非周期性體系Radial portion atomic DFT eqs. numericallyAngular Portion 在DMol3模塊中，電子密度實際上由各個原子軌道的平方和來確定： 在這里，電子密度實際上是由所有占據的分子軌道i.來決定。分子軌道可能由上自旋電子(Alpha電子)和下自旋電子(Beta電子)占據。當Alpha電子和Beta電子的數(shù)目相等的時候，我們可以用單一的分子軌道i.來進展表述，這類體系稱為閉殼層體系(Closed-shell)，在DMol3中不需要選中Spin restricted前面的選項。當Alpha電子和Beta電子的數(shù)目不相等的時候

4、，我們將會使用不同的i.來表述Alpha電子和Beta電子，這類體系稱之為開殼層體系(Opened-shell)或者自旋極化。在DMol3計算的時候需要選中Spin restricted的選項，并指定自旋數(shù)目。 在開殼層體系中，會有兩個不同的電子密度：一個是Alpha電子的電子密度，一個是Beta電子的電子密度。它們的和就是整個體系的總電荷，它們的差就是自旋密度。其中，動能項的方程為：注意，動能項實際上是一個常數(shù)項，在第一次計算完成后，該數(shù)值根本上可以確定，后繼計算中，那么可以忽略這一步驟。勢能項處理： 當我們確定了電子密度之后，傳統(tǒng)的薛定諤方程會從對電子波函數(shù)的處理轉換為對電子密度進展處理的

5、函數(shù)。 需要注意的是，在勢能項中的電子-電子相互作用，指的是兩個電子間的相互作用。但是，在整個體系中，還有三電子、四電子之間的相互作用，這一局部的內容從數(shù)學上是沒有方法得到準確解的，在密度泛函理論中，將這一局部的內容歸入了Exc這一項。 對這一項的處理，才是密度泛函理論處理的核心。針對不同的體系，有LDA和GGA兩種處理方法。 LDA (Local density approximation)局域密度近似方法假定在原子尺度電子密度變化非常緩慢，也就是說，在整個分子區(qū)域內，整個體系表現(xiàn)為連續(xù)的電子氣狀態(tài)。 那么，整個電子交換-相關能就可以表示為對整個電子氣的積分。 在DMol3模塊中，常用的兩種

6、LDA方法是VWN和PWC： VWN：最常用的LSD (Local spin density)相關勢函數(shù)。用來擬和電子氣的準確數(shù)值結果。 PWC：近期開展PWC泛函是在對VWN泛函的某些錯誤校正后的結果，是DMol3模塊的默認泛函。 LSD方法可以準確預測共價體系的構造預測、頻率計算和相關能量。但是，鍵能往往會高估。LDA方法不能用于處理弱健體系，如氫鍵。 LDA的這些缺陷，可以使用更大展開的Exc處理來校正。也稱之為梯度校正方法。 GGA (General gradient-corrected)也成為NLSD(Non-local spin density)方法，近二十年來的計算工作說明，使用

7、梯度校正交換-相關能Exc, d ()可以很好的描述分子體系的熱力學性質。 需要注意的是，GGA方法實際上是一種經歷性的描述函數(shù)。對于不同的研究體系，計算所使用的泛函的準確度實際上是不一樣的。如果要求計算的結果準確可靠，那么需要對相關的函數(shù)進展查閱，或者從相關文獻中查找所使用的泛函。 P91, BP, BLYP, BOP：也稱為廣義梯度近似方法。一般是Becke交換函數(shù)(B88)與Perdew-Wang相關函數(shù) (BP)或者Lee-Yang-Parr相關函數(shù)(BLYP)組合使用。 PBE：PBE (Perdew, Burke和Enzerhof) 泛函具有較強的物理背景，主要用于固體計算，可靠的

8、數(shù)值計算性能，在DFT計算中經常被使用。交換項與Becke相類似，相關項與Pedew-Wang函數(shù)相接近。 RPBE：在PBE的根底上修改得到，對熱力學計算結果較為可靠。 HCTH： 對無機物和氫鍵體系的熱力學計算結果較為可靠。 VWN-BP：COSMO或者COSMO-RS計算推薦使用。 盡管NLSD計算比起LSD計算有很大的改善，但是，在計算反響能壘的時候，得到的能壘數(shù)據和實驗結果相比仍任會有所低估。Generalized Gradient Approximation (GGA)PW91Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992

9、).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. Phys., 58, 1200-1211 (1980). Becke, A. D. J. Chem. Phys., 88, 2547 (1988)

10、.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. Chem. Phys., 88, 2547 (1988). Lee, C.; Yang, W.Parr, R. G. Phys. Rev. B

11、, 37, 786 (1988).HCTHBoese, A. D.; Handy, N. C. J. Chem. Phys., 114, 5497 (2001). 當我們確定了如何處理電子的交換-相關能作用后，我們就可以將初始的薛定諤方程進展轉換，轉換后的Kohn-Sham方程就可以用下式來進展描述： 這里的Et是整個分子的總能量。而如果要準確確實定整個分子的總能量，就必須通過改變體系的電子密度，并由此得到不同的能量數(shù)值，當體系的能量不再發(fā)生變化的時候，就可以認為整個體系趨于穩(wěn)定，最終得到的能量，就是整個體系的最終能量。 獲得最終能量的一系列計算工作我們稱之為SCF (Self-Consist

12、ent field自洽反響場)計算，整個SCF迭代的計算流程，如以下圖所示：選擇初始的Ci。 根據Ci構建初始的分子軌道i。構建電子密度。使用電子密度，計算電子間的勢能項，并考察Exc。計算整個體系的哈密頓能量。解自洽方程以獲得新的一組Ci系數(shù)。構建新的分子軌道i和電子密度 如果new = old，那么通過新的薛定諤方程計算總能量，并完畢。如果new old，那么返回步驟4。 對于有機分子，一般10部左右可以保證整個new old 10-6，整個體系能量收斂。但是對于金屬體系，一般需要更多的迭代來確保體系收斂。DMol3中的SCF迭代過程DMol3: 參數(shù)控制SetupTask: Energy

13、, geometry optimization, TS-search, Overall quality controlXC-Option: LDA, GGASpin optionCharged systems (add or remove electrons)DMol3:參數(shù)控制Task： Energy：直接進展SCF迭代計算，迭代的結果就是體系的總能量。 Geometry Optimization：求解體系的勢能與坐標的一階導數(shù)，當該數(shù)值為0時，各個原子受力為零，此時的構造為能量最低構造。 Dynamics：考慮為各個原子分配速率，模擬體系在外界溫度場作用下的動力學行為； TS Search

14、、TS Optimization、TS Conformation： 對反響體系的計算工作。TS Search是尋找反響過程中的過渡態(tài)構造，并給出反響的能壘和反響熱；TS Optimization是對過渡態(tài)構造進展進一步處理，尋找真實的一階鞍點；TS Conformation那么是對搜索到的反響過渡態(tài)進展進一步確認，搜索可能的中間體和過渡態(tài)。構造優(yōu)化 在DMol3計算中，采用的是BFGS計算方法。BFGS計算是進展構造優(yōu)化的標準算法。此方法適合于獲得能量最低的構造優(yōu)化算法。Tips：單擊More后，可以看到Optimize Cell選項，只有當整個晶胞的晶胞參數(shù)都需要優(yōu)化時，才可以選擇此選項。孤

15、立體系無此選項。 在構造優(yōu)化中，可以勾中Use starting Hessian來指定各個原子的受力方向。DMOL: 參數(shù)設置Electronic optionsIntegration: 各個基函數(shù)的精度SCF tolerance：SCF迭代誤差K-points ：布里淵區(qū)采樣Core treatment (all electron, PP,)Basis set：基組選擇Orbital cutoff quality：各個原子軌道的半徑設置More Core Treatment?All Electron 不對內核電子進展特殊處理，所有的電子都被包含在計算體系中進展處理。Effective Cor

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

17、l3模塊開發(fā)，而ECP那么來源于Hartree-Fock勢。 Core treatment參數(shù)控制那些原子軌道上的電子需要進展處理。默認的設置是對于所有的電子進展處理。對于重原子而言，內核電子的速率接近光速，就必須考慮到相對論效應。 基組設置Min 一個占位原子軌道之使用一個數(shù)值軌道基組。速度更快，精度最低。DN 雙數(shù)值軌道基組。在Min基組的根底上參加了第二組價電子原子軌道。DND 雙數(shù)值軌道基組+d軌道極化函數(shù)。與DN基組類似，但是為所有非氫原子參加了d軌道函數(shù)極化。為DMol3模塊默認基組，在可承受的計算時內，確保精度在可承受的范圍內。也是DMol3計算的最低可信精度基組。DNP 雙數(shù)值

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

19、次能量之間的差值，以電子密度的差異為評判標準Smearing: 選中后，意味著將熱拖尾效應應用到軌道占據上，來加速收斂速度(SCF迭代失敗后，第一項需要調的參數(shù))Density mixing (SCF迭代失敗后，第二項需要調的參數(shù))DMOL: 參數(shù)設置SCF 選項Smearing - Fermi (thermal) occupancySmearing參數(shù)允許電子在所有軌道中按照指定的能量差E 進展拖尾。類似于物理上的熱占位現(xiàn)象。此方法能夠通過允許軌道馳豫而大大加速SCF迭代的收斂速度。會導致虛軌道與占據軌道進展混合，因此，會有一些軌道出現(xiàn)分數(shù)占位。Smearing 熱占位模式的效果初始構造 紅

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

21、) ： 也稱為Pulay mixing方法，由Peter Pulay開展得到，是一種外推方法。用于Harter-Fock SCF迭代方法的加速并穩(wěn)定體系的收斂過程中。 中選中時，能夠加速整個體系的SCF迭代速度。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 pra

22、ctice, only the most 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, DMOL: 參數(shù)設置K-points專門針對周期性體系需要針對不同布里淵區(qū)的K點路徑進展設置(如果不確定，可以

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

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

25、電子的電子密度差。ElectrostaticsElectrostatic potentialElectrostatic moments: dipole, quadrupole, octupole, and hexadecapoleNuclear electric field gradientsFukui Functions - 預測反響性FFs 提供了親電或者親核攻擊的相關信息那個區(qū)域反響更易發(fā)生？ FF可以繪制出3D圖，以顯示反響區(qū)域。 這些類似構造，那個更有反響性？濃縮的FFs顯示方式，可以以原子為中心的等勢面，定量比照不同區(qū)域的反響性能。親核攻擊親電攻擊Fukui index for e

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

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

28、通過頻率計算后，可以將得到的最低能量振動得到的零點能校正參加到之前計算的能量數(shù)值來獲得較為可靠的能量數(shù)值。而體系的熵、焓數(shù)值，那么可以分解為由分子平動、振動和轉動帶來的能量變化。Transition State Search TS Search TS Optimization TS ConformationReaction CoordinateEnergyWith catalystZ1Z2PGRGA*P*Ea,0Ea,1Ea,2DHRWithout catalyst過渡態(tài)計算方法一覽 局域法: (EF, ScanPES) 從單一構造出發(fā)計算勢能面TS 計算較為準確 插值方法: (LST/QST

29、, Ridge, NEB)需要給定初始反響物和產物構造 計算單個構造的能量和能量梯度一般需要作更進一步的驗證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極大值化計算，并使用完全共軛梯度法對構造進展后繼處理。 Complet

30、e LST/QST：執(zhí)行單個LST極大值化計算，并重復使用共軛梯度法優(yōu)化和QST極大值方法進展后繼處理，直到確定過渡態(tài)構造。 QST / Optimization：從QST路徑出發(fā)，重復使用共軛梯度法優(yōu)化和QST極大值方法進展后繼處理，直到確定過渡態(tài)構造。Reaction PathEnergy (Ha)LST Plus CG Peak Refinement:Ring opening of CaprolactoneLSTMaximizationConjugate Gradient MinimizationConjugate gradientrefinement of saddlepoint is

31、 important!0.6eVReactantProductHalgren Libscomb過渡態(tài)構造的頻率分析真實的過渡態(tài)構造應該是一個鞍點，反映在頻率計算上，應該只有一個虛頻率虛頻率的振動模式，應當是在反響物和產物之間發(fā)生震蕩，由此保證反映過程的一致性。 當過渡態(tài)構造中存在多個虛頻，有時候就需要使用TS Optimization方法將多余的虛頻消除，以獲得真正的一階鞍點。 在使用TS Optimization的時候，首先需要使用Tools | Vibrational Analysis對整個體系的Hessian矩陣進展計算處理。 在設置參數(shù)的時候，需要將Vibrational Analys

32、is的窗口激活，并選中最小的正頻率，之后開場DMol3計算。 當體系發(fā)生震蕩的時候，可以將Max. Step size的數(shù)值調小，以確保收斂。假定反響路徑AM1M2BTS 1TS 2A B involves 2 unexpected minima假定 QST 路徑點能否在反響物、產物和TS 2之間找到其他能量最小的構造？AM1M2BTS 1TS 2過渡態(tài)驗證的方法傳統(tǒng)方法:由TS 出發(fā)執(zhí)行最陡下降法計算簡單的最陡下降法優(yōu)化阻尼MD計算 要求大量的步驟二階導數(shù)法 要求進展大量的Hessian矩陣計算在第一個找到的最小點處完畢NEB方法.搜索完整的反響路徑需要計算的步數(shù)少不需要計算Hessian矩

33、陣依賴于離域內坐標優(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工作。TS 計算的完整流程執(zhí)行TS Search計算工作在TS Search得到的xtd路徑上，進展TS Conformaiton工

34、作完成后會給出IRP路徑，可能有一些局域最小點對新找到的最小點構造進展優(yōu)化如果必須，在新得到的最小點構造之間執(zhí)行TS Search計算工作CASTEP的根本原理和參數(shù)設置CASTEP平面波贗勢法(Planewave-Pseudopotential)CASTEP將電子波函數(shù)通過平面波函數(shù)展開贗勢pseudo wave functionpseudopotential有效的減少平面波數(shù)目只適用于周期性體系CASTEP的參數(shù)設置Setup計算任務(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)體系電荷CASTEP的參數(shù)設置Electronic動能截斷值(Energy cutoff)自洽計算收斂精度(SCF tolerance)K點設置(k- point set)贗勢(Pseudopotential) Ultrasoft, Normconserving, On the fly贗勢的表達，實空間 o

36、r 倒空間Core hole的選擇 MoreCASTEP的參數(shù)設置Electronic-MoreBasis 基組的相關參數(shù)設置 動能截斷值的設定 FFT格子設定 基組的測試(Finite basis correction ) SCF 自洽計算的相關參數(shù)設置k-points Potential 選擇贗勢CASTEP的參數(shù)設置Electronic-MoreBasis 基組的相關參數(shù)設置SCF 自洽計算的相關參數(shù)設置 收斂精度的具體設定 自洽計算允許的步數(shù) Density mixing 軌道占據數(shù) 初始自旋值的優(yōu)化k-points Potential 選擇贗勢Density Mixing Densi

37、ty Mixing參數(shù)控制體系中，如何根據特征方程來如何構造新的電子密度。在整個體系中，通過參加阻尼振蕩來確保整個體系的平滑收斂。對于一個最簡單的阻尼方法，遵循如下方程： A就是阻尼系數(shù)。in是當前前一步SCF迭代計算得到的電子密度，根據這一項來構建特征矩陣。out是沒有考慮阻尼系數(shù)，根據新的分子軌道系數(shù)構建得到的電子密度。而new那么是真正用于下一步迭代的電子密度。DIIS(direct inversion in an iterative subspace) ： 也稱為Pulay mixing方法，由Peter Pulay開展得到，是一種外推方法。用于Harter-Fock SCF迭代方法的