====== xrayUtilities の使用例 ======

python のライブラリとして提供されている[[https://xrayutilities.sourceforge.io/|xrayUtilies]]の使用例を幾つか紹介しています。
応用物理学会 結晶工学分科会が企画する[[https://annex.jsap.or.jp/kessho/event/seminar260309/|結晶工学セミナー「PythonとAIで進める結晶工学 ～結晶成長から評価技術まで～」(2026.3.9)]]での講演内容に対応して準備したものです。

主な目的は 「1.使用例」 を示すことです。\\
python の使用やプログラムを書くことについてはセミナーでは、より詳細な説明がありますのでそちらに従って下さい
([[http://d2mate.mdxes.iir.isct.ac.jp/Lecture/jsap-crystal/seminar2025/|セミナーが推奨する神谷先生のwebページ]])。

後日、セミナーの受講とは別に xrayUtilities を試してみたいという時には「[[prepareXrayUtils|xrayUtilities の準備]]」を参照してみて下さい。


===== - 使用例 =====

==== - 結晶構造の定義 ====
  - 結晶構造の定義と表示 ++ ここをクリックすると表示 |
<code python example-01.py>
#!/usr/bin/env python3
#
# 結晶構造の設定と表示
# https://xrayutilities.sourceforge.io/examples.html より
#

import matplotlib.pyplot as plt
import xrayutilities as xu

# xrayUtilities が持つ materials.elements の定義を使って原子を定義する
# (これは原子のX線光学的な性質を保持している)
In = xu.materials.elements.In
P = xu.materials.elements.P

# この In 原子、P 原子 を使って
# 閃亜鉛鉱と六方晶の InP を作ってみる。

## ここでは使わないが閃亜鉛鉱型の InP については弾性定数も設定してみるのでその準備
elastictensor = xu.materials.CubicElasticTensor(10.11e+10, 5.61e+10, 4.56e+10)

# 閃亜鉛鉱の InP の定義
InP = xu.materials.Crystal(
    "InP", xu.materials.SGLattice(216, 5.8687, atoms=[In, P],
                                  pos=['4a', '4c']),
    elastictensor)   # ここでは不要だが、弾性定数も追加で設定(設定しなくても良い)

# 六方晶の InP の定義
InPWZ = xu.materials.Crystal(
    "InP(WZ)", xu.materials.SGLattice(186, 4.1423, 6.8013,
                                      atoms=[In, P], pos=[('2b', 0),
                                                          ('2b', 3/8.)]))

# 以下は、定義した InP の 3D 表示
f = plt.figure()
InP.show_unitcell(fig=f, subplot=121)
plt.title('InP zincblende')
InPWZ.show_unitcell(fig=f, subplot=122)
plt.title('InP wurtzite')

plt.show()
</code>

{{:tabuchi:example-01-res-01.png?600|}}

++
  - 結晶構造の定義のバリエーション ++ ここをクリックすると表示 |
<code python example-01-01.py>
import matplotlib.pyplot as plt
import xrayutilities as xu

# 1. xu はある程度の物質は知っている。xu の materials のライブラリから取り出せる。
GaAs_1 = xu.materials.GaAs
InP_1  = xu.materials.InP

# 2. xu が持っている原子の情報を使って結晶を組み立てる
In = xu.materials.elements.In
P = xu.materials.elements.P
elastictensor = xu.materials.CubicElasticTensor( 10.11e+10, 5.61e+10, 4.56e+10 )
InP_2 = xu.materials.Crystal( "InP",
                             xu.materials.SGLattice( 216, 5.8687,        # 結晶の対称性を指定する群の番号、格子定数
                                                     atoms=[In, P],      # 指定した対称性で実際に配列する原子
                                                     pos=['4a', '4c']),  # 原子の位置
                                                     elastictensor )     # 弾性定数

# 3. cif ファイルを使う、cif ファイルから情報を読む
#   対称性の指定とか考えなくていいので楽。但し、弾性定数は cif ファイルにはないので必要なら後で指定する
# GaAs_2 = xu.materials.Crystal.fromCIF( "GaAs.cif" )
# InAs_1 = xu.materials.Crystal.fromCIF( "InAs.cif" )
# InP_3 = xu.materials.Crystal.fromCIF( "InP.cif" )
# SiC4H_1 = xu.materials.Crystal.fromCIF( "SiC-4H.cif" )
# GaN_1 = xu.materials.Crystal.fromCIF( "GaN-Hex.cif" )
# AlN_1 = xu.materials.Crystal.fromCIF( "AlN-Hex.cif" )


# 試しに 3D 表示
f = plt.figure()
InP_1.show_unitcell(fig=f)      # InP_1 を他のものに変えて試して見て下さい

plt.show()

</code>
cif ファイルを使った例は、このプログラムを実行する場所に cif ファイルを置いておかないとエラーになるので
コメントアウトしてあります (先頭に #　がある)。
試したい場合には、何か自分でダウンロードしてきて、そのファイルを読み込むように書き換えてください。

cif ファイルを入手したい場合には「結晶名、物質名、化学式等 cif」のように検索するとすぐ見つかります。
[[https://legacy.materialsproject.org/|Materials Project の web ページ]]はよく利用されます(登録が必要かも)。
例えば[[https://legacy.materialsproject.org/materials/mp-2534/|「GaAs cif」で検索してでてきた Materials Project のページ]]を例として挙げます。

\\ \\ [[https://next-gen.materialsproject.org/about/terms|Materials Project のファイルは CC-BY-4.0]] だとわかったので
ここで例に上げた cif ファイルをここからもダウンロードできるようにします。\\ 
1. {{ :tabuchi:gaas.cif | GaAs.cif }}　2. {{ :tabuchi:inas.cif | InAs.cif }}　3. {{ :tabuchi:inp.cif | InP.cif (閃亜鉛鉱)}}　4. {{ :tabuchi:inp-wz.cif | InP-wz.cif (六方)}}　5. {{ :tabuchi:sic-4h.cif | SiC-4h.cif (4H)}}　6. {{ :tabuchi:gan-hex.cif | GaN.cif }}　7. {{ :tabuchi:aln-hex.cif | AlN.cif }}\\
これらは全て、上記 Materials Project の web ページから入手したものです。
++

==== - 構造因子の計算 ====
  - 構造因子の計算と表示 ++ ここをクリックすると表示 |
<code python example-02.py>
# xrayUtilities を使うことの宣言。今回はグラフ表示はしないので matplotlib は呼ばない
import xrayutilities as xu

# xu が持ってる結晶の情報を借りる。
# 別の方法で作った結晶についても試してみよう
InP  = xu.materials.InP
GaAs = xu.materials.GaAs

# 構造因子の大きさは
#   1. G=dk となる G、あるいは同じことだが dk に依存する。
#   2. dk = k'-k が同じなら (G=dkになるなら)良いので、本当は k の大きさには依存しない
#   3. 言い換えると、X線の波長、エネルギーには依存しない。
#   4. はずなのだが、原子散乱因子には電子密度分布のフーリエ変換の項(主たる項)以外に異常分散項が存在し
#      これはエネルギー依存性があるため、厳密な議論には X線のエネルギー、波長を指定する必要がある。
#   5. 異常分散の起源は、原子に束縛された電子が電磁波によって揺さぶられる効果
energy= 8048 # eV

#  構造因子を計算する逆格子点のリスト (好きな逆格子点を好きなだけ並べればいい) 
hkllist = [[0,0,0], [1, 1, 1], [2, 2, 2], [3, 3, 3], [0, 0, 0], [1, 0, 0], [2, 0, 0], [4, 0, 0]]


# これから InP の構造因子の計算結果を表示するよ、という表示
print( f"Structual Factors of InP at {energy}[eV]" )

for hkl in hkllist:                               # hkllistから取り出した逆格子点を順番に hkl として
  qvec = InP.Q( hkl )                           # その hkl にたどり着くのに Q すなわち dk を求める
  F = InP.StructureFactor( qvec, energy )       # InP について、その Q の時の構造因子を計算する
  print( f"|F|{hkl} = {abs(F):8.3f}" )          # 結果表示


# これから GaAs の構造因子の計算結果を表示するよ、という表示
print( f"Structual Factors of GaAs at {energy}[eV]" )

for hkl in hkllist:
  qvec = GaAs.Q( hkl )
  F = GaAs.StructureFactor(qvec, energy)
  print(f"|F|{hkl} = {abs(F):8.3f}")
</code>

このコードを実行すると、以下のような結果が画面に表示される。

<code>
Structual Factors of InP at 8048[eV]
|F|[0, 0, 0] =  258.476
|F|[1, 1, 1] =  182.826
|F|[2, 2, 2] =  103.439
|F|[3, 3, 3] =  112.501
|F|[0, 0, 0] =  258.476
|F|[1, 0, 0] =    0.000
|F|[2, 0, 0] =  121.886
|F|[4, 0, 0] =  163.484
Structual Factors of GaAs at 8048[eV]
|F|[0, 0, 0] =  247.148
|F|[1, 1, 1] =  148.498
|F|[2, 2, 2] =    6.527
|F|[3, 3, 3] =   92.194
|F|[0, 0, 0] =  247.148
|F|[1, 0, 0] =    0.000
|F|[2, 0, 0] =    7.106
|F|[4, 0, 0] =  154.631
</code>
++
  - 構造因子のエネルギー依存性の計算と表示 ++ ここをクリックすると表示 |
<code python example-02-2.py>
#
# 構造因子の計算 (エネルギー依存)
#

import xrayutilities as xu
import numpy
import matplotlib.pyplot as plt

GaAs = xu.materials.GaAs

energy= numpy.linspace(500, 20000, 5000) # 500 - 20000 eV
F = GaAs.StructureFactorForEnergy(GaAs.Q(1, 1, 1), energy)

plt.figure(); plt.clf()
plt.plot(energy, F.real, '-k', label='Re(F)')
plt.plot(energy, F.imag, '-r', label='Imag(F)')
plt.xlabel("Energy (eV)"); plt.ylabel("F"); plt.legend()

plt.show()
</code>
{{:tabuchi:ex-02-2.png?600|}}
++




==== - 回折パタンの計算と逆格子マップの表示 ====
  - 回折パタンの計算と逆格子マップの表示 ++ ここをクリックすると表示 |
<code python example-03-1.py>
# 逆格子空間内でのブラッグ点の表示
# X線のエネルギーが何も書かれていないが、デフォルトで Cu Ka1 になってるはず。

import xrayutilities as xu
import matplotlib.pyplot as plt

GaAs = xu.materials.GaAs
hGaAs = xu.HXRD( GaAs.Q( 1, 0, 0 ), GaAs.Q( 0, 0, 1 ) )              # GaAs の x=(1,0,0), y=(0,0,1) 断面のプロットと
ax, h = xu.materials.show_reciprocal_space_plane( GaAs, hGaAs )

InP  = xu.materials.InP
hInP = xu.HXRD( InP.Q( 1, 0, 0 ), InP.Q( 0, 0, 1 ) )                 # InP の x=(1,0,0), y=(0,0,1) 断面のプロットを重ねている
ax, h = xu.materials.show_reciprocal_space_plane( InP, hInP, ax=ax )

plt.show()     # これがないと show_reciprocal_space_plane で作った「図」が表示されない
</code>
{{:tabuchi:ex-03-1.png?600|}}
++
==== - 回折スペクトルの計算 ====
  - AlN/GaN(Sub)の $\omega-2\theta$スペクトルの計算と表示 ++ ここをクリックすると表示 |
<code python example-04-1.py>
# いつものおまじない。numpy, matplotlib, xrayutilities を使用可能にする
import numpy as np
import matplotlib.pyplot as plt
import xrayutilities as xu

#### ここからが本体 ####

# GaN, AlN の結晶構造をCIFから読む
GaN = xu.materials.Crystal.fromCIF("GaN.cif")  # 所持している GaN の cif ファイルの名前と場所 "場所/名前.cif"
AlN = xu.materials.Crystal.fromCIF("AlN.cif")  # 無ければここにのせている GaN, 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

# ここまでで、弾性定数の情報も持った GaN, AlN という結晶が準備できた。
# 次は、この物質のを積み上げていくためにその部品として Layer を作り、
# それを積み上げて積層構造にする
#
## まずは 基板(厚さ無限のGaN Layer)、AlN層を準備
sub = xu.simpack.Layer( GaN, float('inf') )            # 基板 (GaN)
lAlN = xu.simpack.Layer( AlN, 1500, relaxation=1.0 )   # 膜 (AlN)
#
## 準備した Layer を積み上げて AlN層 / GaN 基板という積層構造を構築
Epi_Sub = xu.simpack.PseudomorphicStack001( "AlN/GaN", sub+lAlN )

# ω-2θ 強度計算
## 波長、エネルギー、計算する角度範囲指定
wavelength = xu.wavelength('CuKa1')           # 波長は Cu Ka1 の波長(1.5405...A)
energy     = 12390 / wavelength               # エネルギーは波長から計算
twotheta   = np.linspace(32,38,200)           # 横軸になる角度範囲のベクトルを作っておく

# 動的モデルの設定とシミュレーション
Epi_Sub_Model     = xu.simpack.DynamicalModel( Epi_Sub, energy=energy, resolution_width=0.0001 )
Int_Epi_Sub_Model = Epi_Sub_Model.simulate( twotheta/2, hkl=(0, 0, 2) )

# 結果のプロット
plt.figure()
plt.plot( twotheta/2, Int_Epi_Sub_Model )
plt.title( 'Dynamical Simulation of AlN(1500nm)/GaN Sub')
plt.yscale('log')
plt.xlabel(r'$\omega\ (deg)$')
plt.ylabel('Intensity (arb. u.)')
plt.show()
</code>
\\ \\ [[https://next-gen.materialsproject.org/about/terms|Materials Project のファイルは CC-BY-4.0]] だとわかったので
ここで例に上げた cif ファイルをここからもダウンロードできるようにします。\\ 
6. {{ :tabuchi:gan-hex.cif | GaN.cif }}　7. {{ :tabuchi:aln-hex.cif | AlN.cif }}\\
これらは全て、上記 Materials Project の web ページから入手したものです。
{{:tabuchi:ex-04-1.png?600|}}\\
++
  - AlGaN/GaN(Sub)の$\omega-2\theta$スペクトルの計算と表示 ++ ここをクリックすると表示 |
<code python example-04-2.py>
##### プログラムはほとんど example-03-1.py と同じ #####
#### ここからしばらく一緒です ####
import numpy as np
import matplotlib.pyplot as plt
import xrayutilities as xu

# 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

#### ここまで一緒です。 ####

# さらに混晶も準備する
Alx = 0.2                                                # Al組成 x = 0.2
AlGaN = xu.materials.material.Alloy(GaN, AlN, Alx);      # AlN 0.2 + GaN (1-0.2)

## まずは 基板(厚さ無限のGaN Layer)、AlGaN層を準備
sub    = xu.simpack.Layer( GaN, float('inf') )            # 基板 (GaN)
lAlGaN = xu.simpack.Layer( AlGaN, 1500, relaxation=1.0 )  # 膜 (Al0.2GaN)
                                                          # relaxation=1.0 => 0.0 にすると回折ピークは広角側に動く
#
## 準備した Layer を積み上げて AlGaN層 / GaN 基板という積層構造を構築
Epi_Sub = xu.simpack.PseudomorphicStack001( "AlGaN/GaN", sub+lAlGaN )

#### この先はまた一緒になります。(コメントや画面表示用の文字列とかは別にして)

# ω-2θ 強度計算
## 波長、エネルギー、計算する角度範囲指定
wavelength = xu.wavelength('CuKa1')           # 波長は Cu Ka1 の波長(1.5405...A)
energy     = 12390 / wavelength               # エネルギーは波長から計算
twotheta   = np.linspace(32,38,200)           # 横軸になる角度範囲のベクトルを作っておく

# 動的モデルの設定とシミュレーション
## リファレンスに基板(GaN)だけの計算もする
Sub_Model         = xu.simpack.DynamicalModel( sub,     energy=energy, resolution_width=0.0001 )
Int_Sub_Model     = Sub_Model.simulate( twotheta/2, hkl=(0, 0, 2) )
## こちらは Epi/Sub = AlGaN/GaN の計算
Epi_Sub_Model     = xu.simpack.DynamicalModel( Epi_Sub, energy=energy, resolution_width=0.0001 )
Int_Epi_Sub_Model = Epi_Sub_Model.simulate( twotheta/2, hkl=(0, 0, 2) )

# 結果のプロット
plt.figure()
plt.plot( twotheta/2, Int_Epi_Sub_Model )
plt.plot( twotheta/2, Int_Sub_Model )
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()
</code>
cif ファイルは一つ前の例と一緒に置いてあります。\\
\\
次の様に緩和(relaxation)を1にした場合(上記の例そのままの場合)。おそらく完全緩和。
   lAlGaN = xu.simpack.Layer( AlGaN, 1500, relaxation=1.0 )  # 膜 (Al0.2GaN)
{{:tabuchi:ex-04-2.png?600|}}\\
\\
次の様に緩和(relaxation)を0にした場合(上記の例から次の1行だけ変更)。
   lAlGaN = xu.simpack.Layer( AlGaN, 1500, relaxation=0.0 )  # 膜 (Al0.2GaN)
{{:tabuchi:ex-04-2-2.png?600|}}\\
\\
緩和(relaxation)が1の場合と0の場合を見比べるため、下記の様に追加変更して両方を計算して表示。
  lAlGaN0 = xu.simpack.Layer( AlGaN, 1500, relaxation=0.0 )  # 膜 (Al0.2GaN)
  lAlGaN1 = xu.simpack.Layer( AlGaN, 1500, relaxation=1.0 )  # 膜 (Al0.2GaN)

  Epi0_Sub = xu.simpack.PseudomorphicStack001( "AlGaN/GaN", sub+lAlGaN0 )
  Epi1_Sub = xu.simpack.PseudomorphicStack001( "AlGaN/GaN", sub+lAlGaN1 )

  plt.plot( twotheta/2, Int_Epi0_Sub_Model )
  plt.plot( twotheta/2, Int_Epi1_Sub_Model )
{{:tabuchi:ex-04-2-3.png?600|}}\\
緩和していないと、(GaNの方が格子定数が長くて、AlNの方が短いので、AlNは横方向に引っ張られて)AlNのc軸方向の面間隔が短くなるので　$2d_{\rm c}\sin\theta = \lambda \Rightarrow \sin\theta = \lambda/2d_{\rm c}$ の関係から回折角 $\theta$ は大きくなるはずだが、実際そうなっていることがわかる。
++
  - 仮に生成したAlGaN/GaN(Sub)の$\omega-2\theta$スペクトルに対するなんちゃってフィッティング ++ ここをクリックすると表示 |
\\ この例では、真ん中あたりに出てくる
  # この結果をフィッティング対象と考えて AlN 組成x と、緩和率 relaxation をパラメータとして
  # フィッティングをかけてみる
というコメント行までは、一つ前の例と同じコードです(一部変数名が違う)。\\
ここでは、Al 組成 $x=0.2$、緩和の度合い $relaxation = 0.0$ でスペクトルをまず計算し、これをフィッティング対象にしています。
このモデルに関連した変数は Target という名前が付いています。\\
これに対して、Al組成と、緩和の度合いを未知パラメータとして RMC でフィッティングを試みています。\\
「あってるか?」と試しに計算して受け入れられたもの(暫定的に合ってると判断されたもの)に関連する変数には Try_Now という名前が付いています。
またTry_Now の値と優劣を比較する、別のパラメータやそれで計算されるものに関連する変数には Try_New という名前が付いています。\\
実際に走らせると、グラフ表示の時点で止まります。グラフ表示のウインドウで 'q' を入力すると先に進みます。\\
残差が Limit (プログラム内で設定)より小さくなるか、パラメータが連続 5000回更新されなければ終了します。\\
\\
スペクトルの形状に与える組成と緩和の影響は強い相関があります。
Al組成が大きくなればピークは広角側に移動しますが、緩和が進むとピークは低角側に移動します。
このため、元のパラメータと比較して、
「Al組成が大きくて緩和が進んだ(1.0に近い)」パラメータと「Al組成が小さくてあまり緩和していない(0.0に近い)」パラメータで計算すると
どちらも同じ様な場所にピークが現れます。フィッティングの観点から言うと、どちらも同じ様に「良い」のでフィッティングが進みにくくなります。
更に細かくみるとAl 組成が違えばスペクトルの形は違うので、最終的には正しい Al 組成と緩和度合いに収束するはずですが、
このプログラムでそこに到れる保証は全くありません(フィッティングのプログラムとして、収束性を上げる工夫を全くしていないので)。
<code python example-04-3.py>
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()
</code>
++
  - AlGaN/GaN(Sub)のロッキングカーブスペクトルの計算と表示 ++ ここをクリックすると表示 |
後半に出てくる、
  #### ここまでは先にやった omega-2theta のシミュレーションの時と同一 ####
と書かれた行までは、2つ目の例と同じコードです。
<code python example-04-4.py>
import numpy as np
import matplotlib.pyplot as plt
import xrayutilities as xu

# 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
AlGaN = 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)
lGaN   = xu.simpack.Layer( GaN,   1500, relaxation=1.0 )  # 膜 (GaN)       # ここでは使わない
lAlN   = xu.simpack.Layer( AlN,   1500, relaxation=1.0 )  # 膜 (AlN)       # ここでは使わない
lAlGaN = xu.simpack.Layer( AlGaN, 1500, relaxation=0.0 )  # 膜 (Al0.2GaN)  # relaxation=1.0 => 0.0 にすると回折ピークは広角側に動く
#
## 準備した Layer を積み上げて AlN層 / GaN 基板という積層構造を構築
Epi_Sub = xu.simpack.PseudomorphicStack001( "AlN/GaN", sub+lAlGaN )

# ω-2θ 強度計算
## 波長、エネルギー、計算する角度範囲指定
wavelength = xu.wavelength('CuKa1')           # 波長は Cu Ka1 の波長(1.5405...A)
energy     = 12390 / wavelength               # エネルギーは波長から計算
two_theta  = np.linspace( 32, 38, 200 )       # 横軸になる角度範囲のベクトルを作っておく

#### ここまでは先にやった omega-2theta のシミュレーションの時と同一 ####

# 動的モデルの設定とシミュレーション
## リファレンスに基板(GaN)だけの計算もする
sub_model         = xu.simpack.DynamicalModel( sub,     energy=energy, resolution_width=0.0001 )
int_w2th_sub      = sub_model.simulate( two_theta/2, hkl=(0, 0, 2) )  # 基板だけの時の omega-2theta スペクトル
peak_index_sub    = np.argmax( int_w2th_sub )                         # そのピーク位置(2theta)の指数
peak_sub          = two_theta[ peak_index_sub ]                       # そのピーク位置(2theta)の角度
omega_sub         = np.linspace( peak_sub/2-0.5, peak_sub/2+0.5, 400 )  # 今わかった角度+/-0.5 omega の範囲を対称に
int_RC_sub        = sub_model.simulate( omega_sub, hkl=(0,0,2), geometry="omega" ) # ロッキングカーブ方向の強度計算

## こちらは Epi/Sub = AlGaN/GaN の計算
epi_sub_model      = xu.simpack.DynamicalModel( Epi_Sub, energy=energy, resolution_width=0.0001 )
int_w2th_epi_sub   = epi_sub_model.simulate( two_theta/2, hkl=(0, 0, 2) ) # omega-2theta スペクトル
peak_index_epi_sub = np.argmax( int_w2th_epi_sub )                        # そのピーク index
peak_epi_sub       = two_theta[ peak_index_epi_sub ]                      # ピーク角度
omega_epi_sub      = np.linspace( peak_epi_sub/2-0.5, peak_epi_sub/2+0.5, 400 )
int_RC_epi_sub     = epi_sub_model.simulate( omega_epi_sub, hkl=(0,0,2), geometry="omega" )

# 結果のプロット
plt.figure()
plt.plot( omega_sub     - peak_sub/2,     int_RC_sub )
plt.plot( omega_epi_sub - peak_epi_sub/2, int_RC_epi_sub )
plt.title( 'Dynamical Simulation of Rocking Curve of AlGaN/GaN and GaN Sub')
plt.yscale('log')
plt.xlabel(r'$\omega\ (deg)$')
plt.ylabel('Intensity (arb. u.)')
plt.show()
</code>
{{:tabuchi:ex-04-4.png?600|}}\\
青線が基板(GaN)のロッキングカーブで、オレンジ線がAlGaN/GaNのロッキングカーブ。
2シータはそれぞれについてシミュレーションで計算された $\omega-2\theta$ スペクトルのピーク位置に合わせて計算されてます
(単純に「一番値が大きいところ」なので、AlGaN/GaN も基板のピークを捕まえてるかも)。
いずれにしても、AlGaN/GaN のロッキングカーブは、基板(GaN)のスペクトルを含んでしまってます。
これは、$2\theta$を制限しない計算になっているためで、現実の測定で言うと、検出器前スリットを前回にして測定している状況に相当します。
現実でも、この様な設定で測定を行うとロッキングカーブを見ているのか$\omega-2\theta$を見ているのかよくわからなくなりますが、
この計算はそれを再現していることになります。
++
------
[[講演資料|田渕の講演資料の一部]]\\
[[start|田渕のページのルート]]\\
当 web ページとその下のページに関するお問い合わせ等ございましたら、[[連絡先|連絡先]]にご連絡をお願いします。 \\
今日: {{counter|today}} / 昨日: {{counter|yesterday}} / 総計: {{counter|total}}
~~NOCACHE~~