他のプログラムも智慧とデータが拓くエレクトロニクス新材料開発拠点　公開・非公開プログラム情報で公開しています

2024年度 応用物理学会 結晶工学分科会
結晶工学スクール 関連資料

東京工業大学 科学技術創成研究院　神谷利夫

★注意★　以下の資料、pythonプログラムは教育目的で作成したもので、
間違いやバグ、計算精度の不足、プログラムの最適化などの保証は全くありません。
自己責任でお使いください。

1. pythonについて、必要なモジュール (numpy, scipyなど)をインストールしてください
2. 下記の 「レベル」の★の数は、プログラムの難易度を表しています
3. *.py のリンクをクリックすると、pythonプログラムのダウンロードができます。
   うまくいかない場合は、リンクをマウスで右クリックし、ファイルに保存してください。
4. (プログラムコード) あるいは (プログラムコード・実行結果) をクリックすると、
   web　browser上でプログラムコードを表示します。
   コメントが色付けされているので、プログラムを理解する場合はこちらをご覧ください。 

• 神谷講義スライド    Kamiya-Slides2024.pdf (2024/7/30 14:33更新)
• 講義資料本文     　20240729Text.pdf (2024/7/29 17:09更新)
  　　　　　　　　付録  20210805Appendix.pdf

• 2024年5月12日　 チュートリアル：実空間像から理解するバンド理論 (↓と同じページ)
• 2024年5月17日　 チュートリアル： 第一原理計算で何がわかるか ( 資料・録画)
• 2023年7月28日　 チュートリアル：どのように第一原理バンド計算の条件を決めるか (資料・録画)
• 2024年7月18日　 チュートリアル： 欠陥の第一原理計算 (資料・録画)

関連資料： 神谷利夫、「機能性化合物の設計」、金属2020年10月号
　　　848ページ以降に、周期構造の透過率スペクトルとバンド構造の関係についての説明があります

結晶学関係・一般

• pythonについて

結晶関係

• VESTAによる見やすい結晶構造図の描き方 (VESTAについては http://jp-minerals.org/vesta/jp/)
  　複雑な結晶構造: VESTA-DrawComplexCrystal.pdf
  　結晶の表面構造: VESTA-DrawSurfaceStructure.pdf
  　　ミラー指数の表記法： MillerIndex.pdf
  　　結晶格子とベクトル演算 －非直交格子の法線ベクトルの計算法－: Crystal-VectorAnalysis.pdf
   
• 結晶構造関係プログラム
  
• レベル★★★　ピークフィッティングプログラム　peakfit-scipy-minimize.py  (プログラムコード・実行結果)
  説明： 単一のGauss関数で、共役勾配法を用いてピークフィッティングを行う。
  　　　　フィッティングさせるデータは ycal()関数で生成し、目的関数 (残差二乗和) CalS2()を最小化する。 
  
  詳細

  使用しているアルゴリズム： 非線形最適化 (共役勾配法, scipy.minimize())
  実行方法： python peakfit-scipy-minimize.pｙ  
  資料： Optimization.pdf
• レベル★★　エピタキシャル薄膜、超格子薄膜の対称反射XRDパターン slxrd.py (プログラムコード・事項結果)
  説明： 結晶構造因子を定数として与え、cellsリスト変数で与えられる単位格子構造と、
           layersリスト変数で定義される層構造から、対象配置2θ/θ回折パターンを計算する。
  　　　　画面出力は、単位格子からの回折強度、cells[i]で定義されている単位格子の繰り返し単位からの回折強度、
           layersで定義される層構造からの回折強度を出力。 
  
  詳細

  使用しているアルゴリズム： なし
  実行方法： python slxrd.py

半導体統計関係

• レベル★★　金属の電子密度の計算　 N-integration-metal.py (プログラムコード・実行結果)
  説明：　金属の電子密度を、状態密度関数 N(E) と Fermi-Dirac分布関数 f(E)の積 N(E)f(E) の数値積分によって求める時間を比較する。
  　　　　　数値積分関数として、integrate.quad()とintegrate.romberg()を比較している。
  
  詳細

  使用しているアルゴリズム： 数値積分 (多項式適合法 scipy.integrate.quad(), Romberg積分 scipy.integrate.romberg()) 
  Usage: python N-integration-metal.py T EF
  実行例: python N-integration-metal.py 300 5.0
           　300K、EF=5.0eVの場合に、integrate.quad関数で相対精度 10^-8 を指定して300回積分したときの計算時間を比較する。
  　　　　　積分区間 (1): E = 0 ~ E_F + 6k_BT
  　　　　　積分区間 (2): E = 0 ~ E_F - 6k_BT
  　　　　　積分区間 (3): E = E_F - 6k_BT ~ E_F + 6k_BT
  　　　　　　　積分区間 (1) ＝ 積分区間 (1) ＋ 積分区間 (2)
  資料： EF-Metal.pdf、NumericalIntegration.pdf
• レベル★★★　金属のE_Fの温度依存性　 EF-T-metal.py (プログラムコード・実行結果)   
  説明：　金属のフェルミエネルギーの温度依存性を数値積分とNewton法を使って求める。
  　　　　初期値として 0 K の E_F から始め、温度を上げながら、前回温度のE_F(T) を次の初期値として使うことで、
  　　　　Newton法でも安定して計算できる。
  
  詳細

  使用しているアルゴリズム： 数値積分 (多項式適合法 scipy.integrate.quad())
  　　　　　　　　　　　　　　　　　　方程式の数値解法（Newton法 scipy.optmize.newton())
  実行方法: python EF-T-metal.py
  資料： EF-Metal.pdf、 NumericalIntegration.pdf、 Optimization.pdf
  
• レベル★★★★　半導体の統計物性の温度依存性　 EF-T-semiconductor.py (プログラムコード・実行結果)
  説明：　半導体のフェルミエネルギー、自由電子濃度、自由正孔濃度、イオン化アクセプター濃度、イオン化ドナー濃度の
  　　　　　温度依存性を二分法を使って求める。非縮退近似を用い、ドナー・アクセプター準位の縮重は考慮していない。
  　　　　　Newton法では E_Fの初期値がかなり真値に近くないと発散するので、
  　　　　　価電子帯上端と伝導帯下端エネルギーを初期として二分法を使う。
  
  詳細

  使用しているアルゴリズム： 方程式の数値解法 (二分法 main()関数中に組み込み) 
  Usage: python EF-T-semiconductor.py EA NA ED ND Ec Nv Nc
  実行例: python EF-T-semiconductor.py 0.05 1.0e15 0.95 1.0e16 1.0 1.2e19 2.1e18
  　　　　　Ec = 1.0 eV (= バンドギャップ)、E_A = 0.05 eV, N_A = 10¹⁵ cm^-3, E_D = 0.95 eV, N_D = 10¹⁶ cm^-3　　　 N_c = 1.2x10¹⁹ cm^-3, N_v = 2.1x10¹⁸ cm^-3 (バンドギャップを1.12 eVとすれば、ほぼ Si の物性値) 
  資料： EF-Semiconductor.pdf
  
• レベル★★★★　⁴HeのBose凝縮における化学ポテンシャルの温度依存性　 bose_condensation.py (プログラムコード・実行結果)
  説明：極低温の⁴HeのBose凝縮における化学ポテンシャルを二分法を使って求める。
  
  詳細

  使用しているアルゴリズム: Γ関数 (Gamma(sigma))
  　　　　　　　　　　　　　　　　　　数値積分 (多項式適合法 scipy.integrate.quad())
  　　　　　　　　　　　　　　　　　　方程式の数値解法 (二分法 CalEF() 関数中に組み込み)
  Usage: python bose_condensation.py [Fs|mu] Tmin Tmax Tstep alphamin alphamax 
  実行例１: python bose_condensation.py Fs 3 4.5 0.2 0 0.5 0.002 
  　　　　3.0～4.5 Kの範囲で0.2K毎に、Fs(σ, α)をα = -μ/k_B/T = 0～0.5, 0.02ステップで計算しプロット 
  実行例２: python bose_condensation.py mu 2.5 4.5 0.01 
  　　　　2.5～4.5 Kの範囲で0.01K毎に µ, N’, n0 などを計算してプロット。近似式と比較。
  資料： BoseCondensation.pdf
  参考資料： 阿部 龍蔵 著、熱統計力学 (裳華房、1995)、他の文献は上記資料を参照
  

半導体物性関係

• レベル★　高ドープ半導体のバンドフィリング効果 (Burstein-Mossシフト) BMshift.py　(プログラムコード・実行結果)
  説明： 自由電子モデルによる、高ドープ半導体の光学バンドギャップ変化
  使用しているアルゴリズム： なし
  実行方法： python BMshift.py
   

デバイス関係

• レベル★★★　TFTの飽和移動度 TFTiv.py 入力データ: TransferCurve.csv (プログラムコード・実行結果)
  説明： TFTのIds-Vgs-Vdsデータから飽和移動度を求める
  使用しているアルゴリズム： 線形最小二乗法 (numpy.polyfit())
  実行方法: TransferCurve.csv をカレントフォルダーに入れて、python TFTiv.py
  資料： TFT.pdf
   
• レベル★★★★　線形最小二乗法の誤差とTFTの飽和移動度の誤差
  　　　　　　　　　　　　　　　　TFTiv_errorcheck.py 入力データ: TransferCurve.csv (プログラムコード・実行結果)
  説明： TFTのIds-Vgs-Vdsデータから飽和移動度を求める。
  　　　　線形最小二乗法の誤差から、フィッティング範囲と飽和移動度／閾値電圧の誤差を比較する
  使用しているアルゴリズム： 誤差評価付き一次多項式線形最小二乗法 (lsq1())
  実行方法: TransferCurve.csv をカレントフォルダーに入れて、python TFTiv_errorcheck.py
  資料： TFT-LSQError.pdf
  

バンド理論

• レベル★　二次元波動関数の描画　wavefunction2D.py (プログラムコード・実行結果)
  説明： 自由電子モデル、箱型ポテンシャル量子ドット、水素原子用モデルの２次元波動関数を描画
  
  詳細

  使用しているアルゴリズム： なし
  Usage1: python wavefunction2D.py pw kx0 ky0 kz0 kx ky kz
    自由電子。
  　波数ベクトル(量子数) (kx0, ky0, kz0) (任意実数) をもち、
  　Bloch波数 (kx, ky, kz) の波動関数 (平面波) を描画
  Usage2: python wavefunction2D.py qbox nx ny nz kx ky kz
  　無限大ポテンシャルに閉じ込められた量子ドット。 
  　量子数 (nx, ny, nz) (= 1, 2, 3…)をもち、 
  　Bloch波数 (kx, ky, kz) の 波動関数 (正弦・余弦関数) を描画 
  Usage3: python wavefunction2D.py H n l m kx ky kz
  　水素原子様モデル。
  　量子数 (n, l, m) (n = 1, 2, …, l = 0, 1, …, n-1, m = -l, -l+1, …., l-1, l) をもち、
  　Bloch波数 (kx, ky, kz) の波動関数 を描画。見やすいように、動径関数のサイズは適当に変えている
  実行方法： python wavefunction2D.py pw 0.1 0.0 0.0 0.0 0.0 0.0
  　波数ベクトル (0.1, 0, 0),　Bloch波数 (0, 0, 0) の平面波波動関数を描画
  実行方法： python wavefunction2D.py pw 0.1 0.0 0.0 0.5 0.0 0.0
  　波数ベクトル (0.1, 0, 0),　Bloch波数 (1/2, 0, 0) の平面波波動関数を描画
  実行方法： python wavefunction2D.py pw 0.1 0.0 0.0 1.0 0.0 0.0
  　波数ベクトル (0.1, 0, 0),　Bloch波数 (1, 0, 0) の平面波波動関数を描画
  実行方法： python wavefunction2D.py qdot 1 1 1 0.0 0.0 0.0
  　量子数 (1, 1, 1) の量子ドットの波動関数を描画
  実行方法： python wavefunction2D.py H 1 0 0 0.0 0.0 0.0
  　H 1s (n = 1, l = 0, m = 0)の水素用原子の波動関数を描画
  実行方法： python wavefunction2D.py H 2 1 1 0.0 0.0 0.0
  　H 2p_x (n = 2, l = 1, m = 1)の水素用原子の波動関数を描画
  
• レベル★★　三次元自由電子バンド free_electron_band.py (プログラムコード・実行結果)
  説明： 自由電子モデル (ゼロポテンシャル)による三次元バンド構造
  
  詳細

  使用しているアルゴリズム： なし
  実行方法： python free_electron_band.py
  資料： pw.pdf
  
• レベル★★★　平面波基底による一次元バンド計算 pw1d.py (プログラムコード・実行結果)
  説明： 平面波基底と井戸型ポテンシャルによる一次元バンド構造
  
  詳細

  使用しているアルゴリズム： フーリエ変換 (numpy.fft.fft())、エルミート行列の対角化 (numpy.linalg.eig())
  Usage: python pw1d.py mode (args)
      mode: ft: ポテンシャルのフーリエ変換を表示
               band:　バンド構造を表示
               wf: 波動関数を表示   
  Usage1: python pw1d.py (ft     a na pottype bwidth bpot)
  Usage2: python pw1d.py (band a na pottype bwidth bpot nG kmin kmax nk)
  Usage3: python pw1d.py (wf    a na pottype bwidth bpot nG kw iLevel xwmin xwmax nxw)
  　　pottype: ポテンシャル型 現在は rect のみ実装
  実行例１： python pw1d.py ft 5.4064 64 rect 0.5 10.0 
  　　格子定数 5.4064 Åを 64分割 (FFTを使うため、2ⁿ である必要がある）し、フーリエ変換した係数をプロット。
  　　ポテンシャル障壁 0.5 Å幅、高さ10.0 eV。
  実行例２： python pw1d.py band 5.4064 64 rect 0.5 10.0 3 -0.5 0.5 21
  　　格子定数 5.4064 Åを 64分割。ポテンシャル障壁 0.5 Å幅、高さ10.0 eV。
  　　平面波基底を 3つ 用いて、k = [-0.5, 0.5] の範囲を 21 分割してバンド構造をプロット。
  実行例３： python pw1d.py wf 5.4064 64 rect 0.5 10.0 3 0.0 0 0.0 16.2192 101
  　　格子定数 5.4064 Åを 64分割。ポテンシャル障壁 0.5 Å幅、高さ10.0 eV。
  　　平面波基底を 3つ 用いて、k = 0.0 におけるバンドの 0 番目の準位の波動関数を描く。
  　　　　　（注意： 準位の番号は、エネルギー準位順にはなっていない。コンソール出力で確認すること）
  　　波動関数は 0.0 ～ 16.2192  Å の範囲を 101 分割してプロット。
  資料： pw.pdf
  
• レベル★★★　転送行列法による透過率、波動関数の計算 transfer_matrix.py (プログラムコード・実行結果)
  説明： 転送行列法による、一次元多重井戸型ポテンシャルの電子の透過率
  
  詳細

  使用しているアルゴリズム： なし
  Usage1: python transfer_matrix.py (wf pottype nz Ez0)
  Usage2: python transfer_matrix.py (tr pottype nz Ez0 Emin Emax nE)
  実行例１： python transfer_matrix.py wf mqw 201 0.1 
  　　ポテンシャルを多重量子井戸mqwにする。
  　　座標 (z) を 201分割し、電子のエネルギー　0.1 eV における波動関数の係数およびΨ(z)をプロットする。 
  実行例１'： python transfer_matrix.py wf potenetial_defect.xlsx 201 0.1 
  　　ポテンシャルをExcelファイルpotential_defect.xlsx から読み込む。
  　　座標 (z) を 201分割し、電子のエネルギー　0.1 eV における波動関数の係数およびΨ(z)をプロットする。 
  実行例２： python transfer_matrix.py tr mqw 201 0.1 0.01 1.0 1001
  　　ポテンシャルを多重量子井戸mqwにする。
  　　座標 (z) を 201分割し、電子のエネルギー　0.1 eV における透過率と波動関数をプロットする。 
  　　透過率は、0.01 ～ 1.0 eV の範囲を1001分割する。
  資料： TransferMatrix.pdf　potential_defect.xlsx
  関連資料： 神谷利夫、「機能性化合物の設計」、金属2020年10月号
  　　　　848ページ以降に、周期構造の透過率スペクトルとバンド構造の関係についての説明があります
  
• レベル★★★★　Kronig-Penneyモデルによる一次元バンド計算 kronig_penney.py (プログラムコード・実行結果)
  説明： Kronig-Penneyモデルによる一次元バンド構造と波動関数
  
  詳細

  使用しているアルゴリズム： 多値方程式の解法 (直接探索、セカント法)
  Usage: python kronig_penney.py 
  Usage1: python kronig_penney.py (graph a bwidth bpot k Emin Emax nE) 
  Usage2: python kronig_penney.py (band a bwidth bpot nG kmin kmax nk) 
  Usage3: python kronig_penney.py (wf a bwidth bpot kw iLevel xwmin xwmax nxw)
  実行例１： python kronig_penney.py graph 5.4064 0.5 10.0 0.0 0.0 9.5 51 
  　　格子定数 5.4064 Å、ポテンシャル幅 0.5 Å、高さ 10.0 eVで、k = 0.0 についてのKronig-Penney方程式の
  　　残差 Δ を E = 0.0 ～ 9.5 eVの範囲を51分割してプロット。Δ = 0 のEが固有エネルギー。
  実行例２： python kronig_penney.py band 5.4064 0.5 10.0 -0.5 0.5 21
  　　格子定数 5.4064 Å、ポテンシャル幅 0.5 Å、高さ 10.0 eVで、k = [-0.5,0.5] の範囲を 21分割してバンド構造をプロットする。
  実行例３： python kronig_penney.py wf 5.4064 0.5 10.0 0.0 0 0.0 16.2192 101
  　　格子定数 5.4064 Å、ポテンシャル幅 0.5 Å、高さ 10.0 eVで、k = 0.0 における下から 0 番目の
  　　準位の波動関数をプロットする。波動関数は x = 0.0 ～ 16.2192 Å を 101 分割してプロットする。
  資料： Kronig-Penney.pdf、Optimization.pdf

第一原理計算関係

• レベル★★★★　DFTの自己相互作用誤差: HF近似とLDAによる水素原子1s 軌道: H1s-HF-LDA.py (プログラムコード・実行結果)
  説明： DFTの自己相互作用による誤差を、自己相互作用のないHF近似と比較
  
  詳細

  使用しているアルゴリズム： 数値積分 (多項式適合法 scipy.integrate.quad())
                                         変分原理ではSIMPLEX法 (scipy.optimize(), method='nelder-mead')
  Usage: python H1s-HF-LDA.py mode Z ka Ne 
  　　　mode: 次の文字の組み合わせ d: 基本的なグラフの表示　v: kaを変分原理で最適化 e: 軌道準位ではなく全エネルギーで出力
  　　　　　　　　　　　　　　　　　　　　　　k: kaを変化させたグラフをプロット　n: Neを変化させたグラフをプロット
  実行例１： python H1s-HF-LDA.py ng 1.0 1.0 1.0 
  　　ka = 1.0 (HFの H 1s 軌道波動関数) での１s 軌道準位をNeを0~1と変化させてプロット
  実行例２： python H1s-HF-LDA.py nvg 1.0 1.0 1.0 
  　　実行例１に、kaを変分原理で最適化させた結果を追加
  資料： H1s-SelfInteraction.pdf
   
• レベル★　E(k)から有効質量を計算 (数値微分) EffectiveMass.py  入力データ: band.csv (プログラムコード・実行結果)
  説明： バンド構造 E(k)を二階微分し、有効質量を計算する。数値微分で二階微分をとる際に中央差分の公式を使う。
  
  詳細

  実行例： band.csv をカレントフォルダーに入れて、python EffectiveMass.py
  資料： EffectiveMass.pdf
  　　
• レベル★★★★　全状態密度 D(E) から半導体の統計物性を計算  EF-T-DOS.py 入力データ： DOS.dat (プログラムコード・実行結果)
  説明： DFTで計算した全状態密度 D(E) のデータから、フェルミ準位、自由電子密度、イオン化ドナー密度などを計算する。
  　　　　両無限範囲積分と同等のため、Simpson則 などを使うよりも、単純和 (Rieman和) を取る方が数値積分精度が高い。
  
  詳細

  使用しているアルゴリズム： 数値積分 (Rieman和、integrate())
  　　　　　　　　　　　　　　　　　　方程式の解法 (二分法)
  実行例１： DOS.dat をカレントフォルダーに入れて、python EF-T-DOS.py T
  実行例２： DOS.dat をカレントフォルダーに入れて、python EF-T-DOS.py EF
  注意: プログラム中の Vcell を、DOSの単位が states/cm³ になるように設定すると、Neなどの単位がcm^-3になる。
  　　　  VASPのDOSは states/unit cell になっているので、Vcell としては単位格子体積を ų で入力する。
  　　　　　　　　　（プログラム中で ų => cm³の変換を行っている)
              Vcell = 212.309 # A^3
  資料： EF-Semiconductor.pdf