import numpy as np import matplotlib.pyplot as plt import xrayutilities as xu import random # RMC 的にフィッチングするため乱数を使う # GaN, AlN の結晶構造をCIFから読む GaN = xu.materials.Crystal.fromCIF("GaN.cif") AlN = xu.materials.Crystal.fromCIF("AlN.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 # さらに混晶も準備する Target_Alx = 0.2 # フィッティング対象の Al組成 x = 0.2 Target_relax = 0.0 # フィッティング対象の 緩和 Target_thick = 1500 # フィッティング対象の 厚さ Target = xu.materials.material.Alloy( GaN, AlN, Target_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, Target_thick, relaxation=Target_relax ) # ## 準備した 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 をパラメータとして # フィッティングをかけてみる Try_Alx = 0.0 # 初期値 Alx = 0 (フィッティング対象は 0.2) Try_relax = 1.0 # 初期値 完全緩和を仮定(フィッティング対象は 0.0) Try_thick = 1500 # 膜厚の初期値 AlGaN_Try = xu.materials.material.Alloy(GaN, AlN, Try_Alx); # Alx を仮定した結晶生成 lAlGaN_Try = xu.simpack.Layer( AlGaN_Try, Try_thick, relaxation=Try_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.text( 0.60, 0.90, 'Start!', transform=plt.gca().transAxes ) 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 # 最後にパラメータが更新されてからのトライの数 Try_Alx_New = Try_Alx # ループ内で生成しない(ランダムウォークの対象にしない)場合でもこの変数が存在するようにする Try_relax_New = Try_relax # ループ内で生成しない(ランダムウォークの対象にしない)場合でもこの変数が存在するようにする Try_thick_New = Try_thick # ループ内で生成しない(ランダムウォークの対象にしない)場合でもこの変数が存在するようにする xRange = 1.0 # 組成の変化のスケール rRange = 1.0 # 緩和の変化のスケール tRange = 10 # 厚さの変化のスケール while ( ( Res > Limit ) and ( count < 5000 ) ) : # 残差が十分小さくなるまで実行, 但し残差が 5000 回更新されなければ終了 Try_Alx_New = Try_Alx + random.uniform( -RRange*xRange, RRange*xRange ) # 組成の現在の値 + +/-RRange の一様乱数 if ( Try_Alx_New < 0.0 ) : Try_Alx_New = 0.0 # 下限(0.0)を越えた時の補正 if ( Try_Alx_New > 1.0 ) : Try_Alx_New = 1.0 # 上限(1.0)を越えた時の補正 Try_relax_New = Try_relax + random.uniform( -RRange*rRange, RRange*rRange ) # relaxの現在の値 + +/-RRange の一様乱数 if ( Try_relax_New < 0.0 ) : Try_relax_New = 0.0 # 下限(0.0)を越えた時の補正 if ( Try_relax_New > 1.0 ) : Try_relax_New = 1.0 # 上限(1.0)を越えた時の補正 Try_thick_New = Try_thick + random.uniform( -RRange*tRange, RRange*tRange ) # 膜厚は少し大きく変化させる if ( Try_thick_New < 0.0 ) : Try_thick_New = 0; AlGaN_Try_New = xu.materials.material.Alloy(GaN, AlN, Try_Alx_New); # Alx_New を仮定した結晶生成 lAlGaN_Try_New = xu.simpack.Layer( AlGaN_Try_New, Try_thick_New, relaxation=Try_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( f"{total_count} {count} : R= {Res:.5f} x= {Try_Alx:.2f} relax= {Try_relax:.2f} T= {Try_thick:.1f}" ) # 新しいパラメータで生成された回折スペクトルのほうがよく合っていれば、各パラメータを新しいパラメータに置き換える if ( Res_New < Res ) : print( f" Old: x= {Try_Alx:.2f} relax= {Try_relax:.2f} T= {Try_thick:.2f} : R= {Res:.5f}" ) print( f" New: x= {Try_Alx_New:.2f} relax= {Try_relax_New:.2f} T= {Try_thick_New:.2f} : R= {Res_New:.5f}" ) print( f" Ans: x= {Target_Alx:.2f} relax= {Target_relax:.2f} T={Target_thick:.2f}" ) Try_Alx = Try_Alx_New Try_relax = Try_relax_New Try_thick = Try_thick_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.text( 0.55, 0.90, f'Result: Count= {total_count} R= {Res:.5f}', transform=plt.gca().transAxes ) plt.text( 0.55, 0.85, f' x= {Try_Alx:.2f} (x= {Target_Alx:.2f})', transform=plt.gca().transAxes ) plt.text( 0.55, 0.80, f' r= {Try_relax:.2f} (r= {Target_relax:.2f})', transform=plt.gca().transAxes ) plt.text( 0.55, 0.75, f' t= {Try_thick:.2f} (t= {Target_thick:.2f})', transform=plt.gca().transAxes ) 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( f"Final : R= {Res:.5f} x= {Try_Alx:.2f} relax= {Try_relax:.2f} T= {Try_thick:.2f}" ) plt.figure() plt.plot( twotheta/2, Int_Target ) plt.plot( twotheta/2, Int_Try_Now ) plt.text( 0.55, 0.90, f'Result: Count= {total_count} R= {Res:.5f}', transform=plt.gca().transAxes ) plt.text( 0.55, 0.85, f' x= {Try_Alx:.2f} (x= {Target_Alx:.2f})', transform=plt.gca().transAxes ) plt.text( 0.55, 0.80, f' r= {Try_relax:.2f} (r= {Target_relax:.2f})', transform=plt.gca().transAxes ) plt.text( 0.55, 0.75, f' t= {Try_thick:.2f} (t= {Target_thick:.2f})', transform=plt.gca().transAxes ) 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()