import numpy as np import matplotlib.pyplot as plt import xrayutilities as xu import random # RMC 的にフィッチングするため乱数を使う # GaN, AlN の結晶構造をCIFから読む GaN = xu.materials.Crystal.fromCIF("CIFs/GaN-Hex.cif") AlN = xu.materials.Crystal.fromCIF("CIFs/AlN-Hex.cif") # 今作った結晶(GaN, AlN)は、結晶構造の情報しか持ってないので # 積層構造にした時の歪の計算ができるように、弾性に関わる情報を追加する # ## どちらも Hexagonal GaN.symm_class_name = 'Hexagonal' AlN.symm_class_name = 'Hexagonal' # ## GaN, AlN の弾性定数 Cij の値を、名前のついた値として準備 C11_GaN, C12_GaN, C13_GaN, C33_GaN, C44_GaN = 390, 145, 106, 398, 105 C11_AlN, C12_AlN, C13_AlN, C33_AlN, C44_AlN = 396, 137, 108, 373, 116 # ## xu にサポートしてもらって 2階テンソルの形に組み上げる cij_array_GaN = xu.materials.material.HexagonalElasticTensor( C11_GaN, C12_GaN, C13_GaN, C33_GaN, C44_GaN ) cij_array_AlN = xu.materials.material.HexagonalElasticTensor( C11_AlN, C12_AlN, C13_AlN, C33_AlN, C44_AlN ) # ## さらに4階テンソルの形にもしておく ij to ijkl (プログラムでは '2' を 'to' の意味でよく使う) cijkl_tensor_GaN = xu.materials.material.Cij2Cijkl(cij_array_GaN) cijkl_tensor_AlN = xu.materials.material.Cij2Cijkl(cij_array_AlN) # ## Crystal オブジェクト(GaN, AlN)に4階テンソルを属性として設定 setattr(GaN, 'cijkl', cijkl_tensor_GaN) setattr(AlN, 'cijkl', cijkl_tensor_AlN) # ## 2階テンソルも elastic 属性に代入 GaN.elastic = cij_array_GaN AlN.elastic = cij_array_AlN # さらに混晶も準備する Alx = 0.2 # Al組成 x = 0.2 Target = xu.materials.material.Alloy( GaN, AlN, Alx ); # AlN 0.2 + GaN (1-0.2) # ここまでで、弾性定数の情報も持った GaN, AlN という結晶が準備できた。 # 次は、この物質のを積み上げていくためにその部品として Layer を作り、 # それを積み上げて積層構造にする # ## まずは 基板(厚さ無限のGaN Layer)、AlN層を準備 sub = xu.simpack.Layer( GaN, float('inf') ) # 基板 (GaN) lTarget = xu.simpack.Layer( Target, 1500, relaxation=0.0 ) # 膜 (Al0.2GaN) # relaxation=1.0 => 0.0 にすると回折ピークは広角側に動く # ## 準備した Layer を積み上げて AlN層 / GaN 基板という積層構造を構築 Target = xu.simpack.PseudomorphicStack001( "AlN/GaN", sub+lTarget ) # ω-2θ 強度計算 ## 波長、エネルギー、計算する角度範囲指定 wavelength = xu.wavelength('CuKa1') # 波長は Cu Ka1 の波長(1.5405...A) energy = 12390 / wavelength # エネルギーは波長から計算 twotheta = np.linspace(32,38,200) # 横軸になる角度範囲のベクトルを作っておく # 動的モデルの設定とシミュレーション ## リファレンスに基板(GaN)だけの計算もする Target_Model = xu.simpack.DynamicalModel( Target, energy=energy, resolution_width=0.0001 ) Int_Target = Target_Model.simulate( twotheta/2, hkl=(0, 0, 2) ) # この結果をフィッティング対象と考えて AlN 組成x と、緩和率 relaxation をパラメータとして # フィッティングをかけてみる Alx = 0.0 # 初期値 Alx = 0 (フィッティング対象は 0.2) relax = 1.0 # 初期値 完全緩和を仮定(フィッティング対象は 0.0) AlGaN_Try = xu.materials.material.Alloy(GaN, AlN, Alx); # Alx を仮定した結晶生成 lAlGaN_Try = xu.simpack.Layer( AlGaN_Try, 1500, relaxation=relax ) # Alx と緩和を仮定した layer 生成 Try = xu.simpack.PseudomorphicStack001( "Try/GaN", sub+lAlGaN_Try ) Try_Model = xu.simpack.DynamicalModel( Try, energy=energy, resolution_width=0.0001 ) Int_Try = Try_Model.simulate( twotheta/2, hkl=(0, 0, 2) ) Int_Try_Now = Int_Try # ループに入る前の初期値(比較対象) Res = np.sum( ( np.log(Int_Target) - np.log(Int_Try) )**2 ) # 仮の残差 log を取ると縦軸logのグラフで見える微細な差が影響する # Res = np.sum( ( Int_Target - Int_Try )**2 ) # 仮の残差 log を取らないと、細かなところは無視される Limit = 0.01 # log をとると残差は大きい # Limit = 0.0001 # log をとらないと残差は小さい if ( Res > 5 ) : RRange = 0.5 else : RRange = Res / 10; print ( "Res", Res ) # 初期のプロット plt.figure() plt.plot( twotheta/2, Int_Target ) plt.plot( twotheta/2, Int_Try_Now ) plt.title( 'Dynamical Simulation of AlGaN/GaN and GaN Sub') plt.yscale('log') plt.xlabel(r'$\omega\ (deg)$') plt.ylabel('Intensity (arb. u.)') plt.show() total_count = 0 # 全トライの数 count = 0 # 最後にパラメータが更新されてからのトライの数 relax_New = relax # ループ内で生成しない(ランダムウォークの対象にしない)場合でもこの変数が存在するようにする Alx_New = Alx # ループ内で生成しない(ランダムウォークの対象にしない)場合でもこの変数が存在するようにする while ( ( Res > Limit ) and ( count < 5000 ) ) : # 残差が十分小さくなるまで実行, 但し残差が 1000 回更新されなければ終了 # Alx_New = Alx + random.uniform( -0.3, 0.3 ) # 組成の現在の値 + [-0.3, 0.3]の一様乱数 Alx_New = Alx + random.uniform( -RRange, RRange ) # 組成の現在の値 + +/-RRange の一様乱数 if ( Alx_New < 0.0 ) : Alx_New = 0.0 # 下限(0.0)を越えた時の補正 if ( Alx_New > 1.0 ) : Alx_New = 1.0 # 上限(1.0)を越えた時の補正 # relax_New = relax + random.uniform( -0.1, 0.1 ) # relaxの現在の値 + [-0.1, 0.1]の一様乱数 relax_New = relax + random.uniform( -RRange, RRange ) # relaxの現在の値 + +/-RRange の一様乱数 if ( relax_New < 0.0 ) : relax_New = 0.0 # 下限(0.0)を越えた時の補正 if ( relax_New > 1.0 ) : relax_New = 1.0 # 上限(1.0)を越えた時の補正 AlGaN_Try_New = xu.materials.material.Alloy(GaN, AlN, Alx_New); # Alx_New を仮定した結晶生成 lAlGaN_Try_New = xu.simpack.Layer( AlGaN_Try_New, 1500, relaxation=relax_New ) # Alx_New と緩和を仮定した layer 生成 Try_New = xu.simpack.PseudomorphicStack001( "Try_New/GaN", sub+lAlGaN_Try_New ) Try_Model_New = xu.simpack.DynamicalModel( Try_New, energy=energy, resolution_width=0.0001 ) Int_Try_New = Try_Model_New.simulate( twotheta/2, hkl=(0, 0, 2) ) Res_New = np.sum( ( np.log(Int_Target) - np.log(Int_Try_New) )**2 ) # 残差 log を取ると縦軸logのグラフで見える微細な差が影響する # Res_New = np.sum( ( Int_Target - Int_Try_New )**2 ) # 残差 log を取らないと、細かなところは無視される total_count += 1 count += 1 print( total_count, count, ":", Res, Alx, relax ) # 新しいパラメータで生成された回折スペクトルのほうがよく合っていれば、各パラメータを新しいパラメータに置き換える if ( Res_New < Res ) : print( Res, Alx, relax, " : ", Res_New, Alx_New, relax_New ) Alx = Alx_New relax = relax_New Int_Try_Now = Int_Try_New Res = Res_New if ( Res > 5 ) : RRange = 0.5 else : RRange = Res / 10; count = 0 # 中間結果のプロット plt.figure() plt.plot( twotheta/2, Int_Target ) plt.plot( twotheta/2, Int_Try_Now ) plt.plot( twotheta/2, Int_Try_New ) plt.title( 'Dynamical Simulation of AlGaN/GaN and GaN Sub') plt.yscale('log') plt.xlabel(r'$\omega\ (deg)$') plt.ylabel('Intensity (arb. u.)') plt.show() # 結果のプロット print( "Final : ", Res, Alx, relax, ) plt.figure() plt.plot( twotheta/2, Int_Target ) plt.plot( twotheta/2, Int_Try_Now ) plt.title( 'Dynamical Simulation of AlGaN/GaN and GaN Sub') plt.yscale('log') plt.xlabel(r'$\omega\ (deg)$') plt.ylabel('Intensity (arb. u.)') plt.show()