我們研究群所使用基於密度泛函理論的計算軟體CASTEP來研究材料物性在許多系統的計算上相當成功，但是對於過度金屬氧化物這有些局域化d電子態的材料其倍頻係數計算有不精確的情形，所以引入了LDA+U方法來修正，我們研究群使用的 old CASTEP 程式已存在了LDA+U方法，並且也給出合理的結果，我們想要更進一步讓 New CASTEP程式也有LDA+U方法。在嘗試建構New CASTEP LDA+U程式之後，我們也成功的在New CASTEP 上建立了我們的LDA+U程式，並進行了一些測試。測試後我們的LDA+U程式能夠正常運作，得到的結果也能夠符合理論數據，雖然我們的LDA+U程式收斂性仍有待加強，並且仍有一些未完成的部分(不能進行能帶結構計算)，但整體來說仍是一個不錯的參考。

The DFT bases code CASTEP was used in our group for many successful studies of materirals properties of various system.However for transition metal Oxides there localized d-electron casue the calculate SHG coefficient to be not accurate.We used LDA+U to correct the problem which was demonstrated using old CASTEP code.
Now we would like to implement LDA+U into New CASTEP.We succeed implement LDA+U code into New CASTEP, and make some test.Our code can normally complete LDA+U calculate,and the result conform with other’s paper.Although our LDA+U code not realy complete(can not be used on band structure calculate) ,but the whole was still a good reference.

目錄
第一章	簡介  -----------------------------------------------------------------------------  1
1.1 動機   ------------------------------------------------------------------------------  1
1.2 LDA+U 理論  ---------------------------------------------------------------------  1
1.3 一個簡化的分析  ----------------------------------------------------------------- 10

第二章	LDA+U 在 CASTEP v4.2 下的探討  ------------------------------------  13
2.1 CASTEP 簡介  -------------------------------------------------------------------  13
2.2 CASTEP 程式架構  -------------------------------------------------------------  13
2.2.1電子結構收斂計算electronic_minimisation副程式的分析  --------  15
2.2.2 electronic_minimisation副程式底下所使用的副程式的分析  ------  20
2.2.3能帶結構計算分析  ---------------------------------------------------------  31
2.3 LDA+U 程式實作  ---------------------------------------------------------------  32

第三章	LDA+U程式實測  -------------------------------------------------------------  41
  3.1 NiO LDA+U計算  ----------------------------------------------------------------  41
  3.2 ZnO LDA+U計算  ----------------------------------------------------------------  48
  3.3 FeO LDA+U計算  ----------------------------------------------------------------  52
  3.4原子軌域波函數的切割與切割修飾  ------------------------------------------  57
  3.5 Ultrasoft 贗勢下LDA+U計算  ------------------------------------------------  67

第四章	與其他 LDA+U 方法的比較  -----------------------------------------------  71
4.1英國CASTEP發展群的LDA+U方法  ---------------------------------------  71
4.2英國CASTEP發展群的LDA+U程式分析  ---------------------------------  73
4.2.1 electronic_minimisation副程式的分析  ----------------------------------  73
4.2.2 electronic_minimisation副程式所使用的副程式分析  ----------------  80
4.2.3 LDA+U功能所使用的副程式分析  --------------------------------------  86
4.3與英國CASTEP發展群LDA+U方法的比較  ------------------------------  94

第五章 結論  ------------------------------------------------------------------------------  96

參考文獻  ----------------------------------------------------------------------------------  98

附錄A  -------------------------------------------------------------------------------------  99
附錄B  -------------------------------------------------------------------------------------  100

圖目錄

圖1.1  U參數形式簡易說明圖 ----------------------------------------------------------  2
圖2.1  CASTEP程式架構圖 -----------------------------------------------------------  14
圖2.3  electronic_minimisation副程式說明 -----------------------------------------  15
圖2.2  CASTEP電子結構收斂計算程式流程圖 -----------------------------------  16
圖2.4  electronic_initialise副程式說明 -----------------------------------------------  17
圖2.5  electronic_prepare_H副程式說明 --------------------------------------------  17
圖2.6  electronic_apply_H副程式說明 ----------------------------------------------  17
圖2.7  electronic_minimisation進行自洽計算，並呼叫electronic_iter_diag_H --------------------------------------------------------------------------------------------------- 18
圖2.8  electronic_minimisation呼叫electronic_find_occupancie計算新的波函數的佔據數，並呼叫electronic_apply_H_energy來計算新的波函數總能 ------------ 18
圖2.9  electronic_minimisation呼叫electronic_prepare_H、 electronic_apply_H
建立新的Hamiltonian與新的梯度方向 ----------------------------------------------  19
圖2.10  自洽場計算後，如果收斂且Density Mix，則復原混合過的電荷密度 -------------------------------------------------------------------------------------------------  19
圖2.11  electronic_minimisation副程式指定Etot = total_energy ----------------  20
圖2.12  electronic_apply_H副程式的總類 ------------------------------------------  21
圖2.13  electronic_apply_H_wv又再細分為實空間波函數與導空間波函數 --- 22
圖2.14  electronic_apply_H_recip_wv副程式計算 的交換相干能(exchange potential)部分 ------------------------------------------------------------------------------  22
圖2.15  electronic_apply_H_recip_wv副程式將各項的本徵值加起來而成為Kohn-Sham波函數的本徵值 -----------------------------------------------------------  23
圖2.16  electronic_apply_H_recip_wv副程式將各項的能量加起來當作總能 -------------------------------------------------------------------------------------------------  23
圖2.17  electronic_iter_diag_H副程式的說明 -------------------------------------  23
圖2.18  electronic_iter_diag_H副程式，Density Mixing的流程 --------------  24
圖2.19  electronic_find_eigenstate副程式將波函數wvfn的某個能帶(band)取出儲存至bnd，並歸一化bnd -----------------------------------------------------------------  25

圖2.20  electronic_find_eigenstate副程式求出單一能帶的本徵值，並呼叫wave_kinetic_eigenvalues計算單一個能帶的動能 ---------------------------------  25
圖2.21  electronic_find_eigenstate呼叫electronic_step_direction_bks、electronic_apply_H、electronic_ideal_step_size進行梯度落下法計算 ----------  26
圖2.22  electronic_find_eigenstate副程式更新bnd與H_bnd -------------------  27
圖2.23  electronic_apply_H_recip_bnd副程式的片段 ----------------------------  27
圖2.24  electronic_step_direction_bks副程式的說明 -----------------------------  28
圖2.25  electronic_ideal_step_size副程式的說明 ---------------------------------  28
圖2.26  electronic_apply_H_energy_local與electronic_apply_H_energy_eigen兩種副程式的輸入值 ---------------------------------------------------------------------------  29
圖2.27  electronic_apply_H_energy_eigen副程式的說明 ------------------------  30
圖2.28  electronic_apply_H_energy_eigen副程式算出最後的總能 ------------  30
圖2.29  band_structure_calculate副程式的說明 -----------------------------------  31
圖2.30  band_structure_apply_H副程式的說明 ------------------------------------  31
圖2.31  band_structure_apply_H副程式最後沒有進行本徵值與能量的加總 -------------------------------------------------------------------------------------------------  32
圖2.32  popn_create_lcao_basis副程式的說明 -------------------------------------  33
圖2.33  原子軌域波函數儲存在波函數的類型裡 ---------------------------------  33
圖2.34  popn_create_lcao_basis副程式所使用的參數說明 ----------------------  33
圖2.35  popn_create_lcao_basis副程式所使用參數的輸出與輸入 -------------  34
圖2.36  建立局域化原子軌域波函數的插入點 ------------------------------------  34
圖2.37  ldau_initialise副程式的位置 ------------------------------------------------  35
圖2.38  ldau_initialise副程式呼叫ldau_create_den_matrix副程式計算局域軌域的佔據數，並呼叫ldau_trans_matrix副程式將所需的軌域取出 -------------------  35
圖2.39  ldau_initialise副程式呼叫ldau_diagonal_density_matrix對角化密度矩陣-------------------------------------------------------------------------------------------------  36
圖2.40  在 更新之前呼叫ldau_initialise副程式更新局域軌域佔據數 --- 36
圖2.41  u_apply副程式分成u_apply_wv與u_apply_bnd -----------------------  37
圖2.42  u_apply副程式的參數輸入 -------------------------------------------------  37
圖2.43  u_apply_wv副程式計算  ------------  38
圖2.44  nlpot_apply_S副程式的說明 ------------------------------------------------  38
圖2.45  Kohn-Sham波函數與  內積以求得本徵值 -----------------------  38
圖2.46  在electronic_apply_H副程式中加入u_apply副程式 ------------------  39
圖2.47  ldau_calc_total_u_energy副程式所輸入的參數 --------------------------  39
圖2.48  ldau_calc_total_u_energy副程式計算Hubbard U總能 -----------------  40
圖2.49  在electronic_apply_H副程式加上ldau_calc_total_u_energy副程式 -------------------------------------------------------------------------------------------------  40
圖2.50  在electronic_apply_H_energy_eigen副程式加上ldau_calc_total_u_energy副程式 ---------------------------------------------------------------------------------------  40
圖3.1  NiO 晶體模型 -------------------------------------------------------------------  41
圖3.2  上圖U=0eV時NiO的d軌域PDOS。 ------------------------------------  42
圖3.3  NiO的d軌域PDOS圖，U值由上而下依序為1 eV、2 eV、3 eV ---  43
圖3.4  上圖U = 4 eV時NiO的d軌域PDOS --------------------------------------  44
圖3.5  英國CASTEP發展群LDA+U方法計算出來的NiO d軌域PDOS圖，由上而下U值為1 eV、1.5 eV、2 eV ----------------------------------------------------  46
圖3.6  ZnO 晶體模型 -------------------------------------------------------------------  48
圖3.7  U = 0 eV時ZnO的d軌域PDOS，橫軸單位為電子伏特(eV)以費米能量為原點0 eV ----------------------------------------------------------------------------------  49
圖3.8  U=0.3675 eV時ZnO的d軌域PDOS，橫軸單位為電子伏特(eV)以費米能量為原點0 eV ----------------------------------------------------------------------------  49

圖3.9  U = 0.735 eV時ZnO的d軌域PDOS，橫軸單位為電子伏特(eV)以費米能量為原點0 eV ----------------------------------------------------------------------------  50
圖3.10  U = 0 eV與U = 0.735 eV的d軌域PDOS圖，橫軸單位為電子伏特(eV)以費米能量為原點0 eV -------------------------------------------------------------------  50
圖3.11  paper的ZnO DOS圖，橫軸單位為電子伏特(eV)以最高佔據態為原點0 eV。上半部為ZnO LDA+U的DOS圖，U = 2.0 eV -------------------------------  51
圖3.12  FeO 晶體模型 ------------------------------------------------------------------  52
圖3.13  U=0eV時，FeO d 軌域PDOS圖 ------------------------------------------  53
圖3.14  U=0eV時，FeO d 軌域PDOS圖 ------------------------------------------  53
圖3.15  U=3eV時，FeO d 軌域PDOS圖 ------------------------------------------  54
圖3.16  FeO d軌域PDOS圖minority spin部分局部放大 -----------------------  55
圖3.17  上方的圖為GGA方法FeO PDOS圖，下方的圖為LDA+U方法FeO
PDOS圖 -------------------------------------------------------------------------------------  56
圖3.18   Ni的d原子軌域波函數，切割半徑固定為1 Å (Ångstom)，隨著截面寬度d不斷減少的改變 ----------------------------------------------------------------------  57
圖3.19  切割過原子軌域波函數的NiO d軌域PDOS圖，由上而下U = 1 eV、2eV、3 eV ------------------------------------------------------------------------------------  58
圖3.20  切割半徑  Å (Ångstom)，NiO d軌域PDOS圖，由上而下U = 1 eV、2 eV、3 eV、4 eV、5eV ------------------------------------------------------------  60
圖3.21  切割半徑  Å (Ångstom)，NiO d軌域PDOS圖， U值依序為1 eV、2 eV、3 eV、4 eV、5 eV ---------------------------------------------------------  61,62
圖3.22  切割半徑  Å (Ångstom)，NiO d軌域PDOS圖， U值依序為1 eV、2 eV、3 eV、4 eV、5 eV --------------------------------------------------------------  63,64
圖3.23  使用Ultrasoft贗勢的NiO d軌域PDOS圖，U值依序為1 eV、2 eV 、3 eV、4 eV、5 eV ----------------------------------------------------------------------  67,68

圖3.24  MS4.3使用Ultrasoft贗勢的NiO d軌域PDOS圖，U依序為1 eV、2 eV 、3eV、4eV、5eV ------------------------------------------------------------------------  69,70
圖4.1 LDA+U的資訊被宣告為lcao_data的類型(type) ----------------------------  74
圖4.2  lcao_data類型所在的模組 ----------------------------------------------------  74
圖4.3 lcao_data類型(TYPE)所包含的內容 ------------------------------------------  74
圖4.4  呼叫electronic_initialise副程式，並傳入ldau_occm ---------------------  75
圖4.5 呼叫electronic_prepare_H，並傳入ldau_occm ------------------------------  75
圖4.6 呼叫electronic_apply_H副程式，並傳入ldau_occm ----------------------  76
圖4.7 呼叫firstd_calculate_forces，並傳入ldau_occm -----------------------------  76
圖4.8 呼叫electronic_iter_diag_H副程式，並傳入ldau_occm -------------------  76
圖4.9  呼叫electronic_find_occupancies副程式，並傳入ldau_occm -----------  77
圖4.10  呼叫electronic_prepare_H與electronic_apply_H副程式，並傳入ldau_occm -------------------------------------------------------------------------------------- 77
圖4.11 呼叫electronic_scf_footer與electronic_finalise，並傳入ldau_occm ---  78
圖4.12  electronic_minimisation副程式Density Mixing部分的流程圖 -------  79
圖4.13  electronic_initialise副程式的說明 ------------------------------------------  80
圖4.14  electronic_finalise副程式的說明 --------------------------------------------  80
圖4.15  electronic_scf_footer副程式的說明 ----------------------------------------  80
圖4.16 electronic_prepare_H副程式的說明 ------------------------------------------  81
圖4.17  呼叫nlxc_ldau_calculate來計算LDA+U的位能 ------------------------  81
圖4.18  呼叫nlxc_ldau_apply副程式，並將得到的U位勢期望值值加入到
loc_eigenvalues(nb)裡 ---------------------------------------------------------------------  82
圖4.19  呼叫pot_apply副程式，產生局域化位勢本徵值loc_eigenvalues(nb) -------------------------------------------------------------------------------------------------  82
圖4.20  局域化位勢本徵值loc_eigenvalues(nb)加入動能本徵值ek(nb)以及非局
域化交換位勢本徵值nloc_ex_eigenvalues(nb) ---------------------------------------  82
圖4.21   局域化位能本徵值loc_eigenvalues(nb) 加上非局域化遮蔽位勢本徵值
nl_screened_eigenvalues(nb)，而成為Kohn-Sham方程式的本徵值 ------------  83
圖4.22  electronic_apply_H副程式將所有的能量加起來求得總能 ------------  83
圖4.23  呼叫electronic_find_eigenstate，並傳入ldau_occm ---------------------  83
圖4.24  electronic_find_occupancies副程式的說明 -------------------------------  84
圖4.25  Density Mixing部分的內容，可以看到並沒有LDA+U的計算 -----  85
圖4.26  electronic_apply_H_energy副程式將本徵值的總合加上修正項算出最後的總能 ---------------------------------------------------------------------------------------  85
圖4.27  nlxc_ldau_initialise副程式的說明 ------------------------------------------  86
圖4.28  nlxc_ldau_finalise副程式的說明 --------------------------------------------  86
圖4.29  nlxc_ldau_calculate副程式的說明 ------------------------------------------  86
圖4.30  呼叫ldau_calculate_overlap_private計算基底與Kohn-Sham波函數的內積 ---------------------------------------------------------------------------------------------  87
圖4.31  呼叫ldau_calculate_occ_matrix計算局域化原子軌域基底的佔據矩陣 -------------------------------------------------------------------------------------------------  87
圖4.32  nlxc_ldau_calculate副程式計算LDA+U能量的回圈 -----------------  88
圖4.33   nlxc_ldau_energy與nlxc_ldau_correction的變數宣告 ---------------  90
圖4.34  nlxc_ldau_apply副程式的說明 ---------------------------------------------  90
圖4.35  nlxc_ldau_apply副程式計算LDA+U本徵值部分 ----------------------  90
圖4.36  nlxc_ldau_apply副程式計算LDA+U本徵值部分，第三層迴圈的片段 -------------------------------------------------------------------------------------------------  91
圖4.37  nlxc_ldau_apply副程式計算LDA+U本徵值部分，第四層迴圈的片段 -------------------------------------------------------------------------------------------------  91
圖4.38  nlxc_ldau_apply副程式計算LDA+U位能作用到波函數 ------------  92
圖4.39  nlxc_ldau_apply副程式將wvfn_temp儲存至vu_wvfn，並呼叫nlpot_apply_S副程式把S重疊算符作用到vu_wvfn上 ---------------------------  93
圖4.40  nlpot_apply_S副程式的說明 ------------------------------------------------  93

表目錄

表3.1 NiO castep 計算參數 --------------------------------------------------------  41,42
表3.2  ZnO castep 參數設定 -----------------------------------------------------  48,49
表3.3  FeO castep參數設定 ----------------------------------------------------------  52
表3.4  固定截面寬度d，不同切割半徑 的LDA+U計算結果 --------------  65

1.	Vlasdimir I. Anisimov, Jan Zaanen, and Ole K. Andersen , Band theory and Mott insulators: Hubbard U instead of Stoner I, PRB v.44, p.943  (1991)
2.	Vlasdimir I. Anisimov, I. V. Solovyev, and M. A. Korotin , Density-functional theory and NiO photoemission spectra, PRB v.48, p.16929  (1993)
3.	I. V. Solovyev, P. H. Dederichs and V. I. Anisimov, Corrected atomic limit in the local-density approximation and the electronic structure of d impurities in Rb, PRB v.50, p.16861  (1994)
4.	H. Sawada, Y. Morikawa, K. Terakura and N. Hamada, Jahn-Teller distortion and magnetic structure in LaMnO3 , PRB v.56, p.12154  (1997)
5.	David Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, PRB v.41, p.7892  (1990)
6.	Shu-jun Hu, Shi-shen Yan, Ming-wen Zhao, and Liang-mo Mei, First-principles LDA+U calculations of the Co-doped ZnO magnetic semiconductor, PRB v.73, p.245205  (2006)
7.	Matteo Cococcioni* and Stefano de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA+U method, PRB v.71, p.035105  (2005)
8.	Benedict Russell, Micoroscopic mechanisms of iron oxide biomineralisation, St John’s College, Cambridge  (2007)
9.	劉嘉鴻，局域化Kohn-Sham態的計算方法與其在含稀土元素晶體之光學性質研究上的應用，淡江大學