Research Notes

T. Matsumura

結晶構造因子の計算

Mathematica routines

X線回折や中性子回折実験の基本である結晶構造因子の計算プログラム.たいした内容ではないが,いろいろな意味で重要.

LaRu2Al10 –– Cmcm (No. 63) ––

原子位置座標

(* Wyckoff Positions of Atoms *)
yLa = 0.1256;
xAl1 = 0.23; yAl1 = 0.35;
xAl2 = 0.349; yAl2 = 0.137;
yAl3 = 0.149; zAl3 = 0.601;
yAl4 = 0.38; zAl4 = 0.043;
xAl5 = 0.227;
wcpLa = {{0, yLa, 1/4}, {0, -yLa, 3/4}, {1/2, 1/2 + yLa, 1/4}, {1/2, 1/2 - yLa, 3/4}};
wcpRu = {{1/4, 1/4, 0}, {3/4, 3/4, 1/2}, {3/4, 1/4, 1/2}, {1/4, 3/4, 0}, {1/2 + 1/4, 1/2 + 1/4, 0}, {1/2 + 3/4, 1/2 + 3/4, 1/2}, {1/2 + 3/4, 1/2 + 1/4, 1/2}, {1/2 + 1/4, 1/2 + 3/4, 0}};
wcpAl1 = {{xAl1, yAl1, 1/4}, {-xAl1, -yAl1, 3/4}, {-xAl1, yAl1, 1/4}, {xAl1, -yAl1, 3/4}, {1/2 + xAl1, 1/2 + yAl1, 1/4}, {1/2 - xAl1, 1/2 - yAl1, 3/4}, {1/2 - xAl1, 1/2 + yAl1, 1/4}, {1/2 + xAl1, 1/2 - yAl1, 3/4}};
wcpAl2 = {{xAl2, yAl2, 1/4}, {-xAl2, -yAl2, 3/4}, {-xAl2, yAl2, 1/4}, {xAl2, -yAl2, 3/4}, {1/2 + xAl2, 1/2 + yAl2, 1/4}, {1/2 - xAl2, 1/2 - yAl2, 3/4}, {1/2 - xAl2, 1/2 + yAl2, 1/4}, {1/2 + xAl2, 1/2 - yAl2, 3/4}};
wcpAl3 = {{0, yAl3, zAl3}, {0, -yAl3, zAl3 + 1/2}, {0, yAl3, -zAl3 + 1/2}, {0, -yAl3, -zAl3}, {1/2, 1/2 + yAl3, zAl3}, {1/2, 1/2 - yAl3, zAl3 + 1/2}, {1/2, 1/2 + yAl3, -zAl3 + 1/2}, {1/2, 1/2 - yAl3, -zAl3}};
wcpAl4 = {{0, yAl4, zAl4}, {0, -yAl4, zAl4 + 1/2}, {0, yAl4, -zAl4 + 1/2}, {0, -yAl4, -zAl4}, {1/2, 1/2 + yAl4, zAl4}, {1/2, 1/2 - yAl4, zAl4 + 1/2}, {1/2, 1/2 + yAl4, -zAl4 + 1/2}, {1/2, 1/2 - yAl4, -zAl4}};
wcpAl5 = {{xAl5, 0, 0}, {-xAl5, 0, 1/2}, {-xAl5, 0, 0}, {xAl5, 0, 1/2}, {1/2 + xAl5, 1/2, 0}, {1/2 - xAl5, 1/2, 1/2}, {1/2 - xAl5, 1/2, 0}, {1/2 + xAl5, 1/2, 1/2}};

格子定数と単位格子ベクトル

(* lattice parameters *)
aorth = 9.1442;
borth = 10.1873;
corth = 9.11611;

(* unit vectors of the lattice in the xyz-coordinate *)
a1 = {1, 0, 0}*aorth;
a2 = {0, 1, 0}*borth;
a3 = {0, 0, 1}*corth;

逆格子ベクトル

(* reciprocal lattice vectors *)
b1 = Cross[a2, a3]/(a1.Cross[a2, a3])*2*Pi;
b2 = Cross[a3, a1]/(a1.Cross[a2, a3])*2*Pi;
b3 = Cross[a1, a2]/(a1.Cross[a2, a3])*2*Pi;

実空間での原子位置

(* atomic positions in the xyz-coordinate *)
Lasite = Table[wcpLa[[i, 1]]*a1 + wcpLa[[i, 2]]*a2 + wcpLa[[i, 3]]*a3, {i, 1, Length[wcpLa]}];
Rusite = Table[wcpRu[[i, 1]]*a1 + wcpRu[[i, 2]]*a2 + wcpRu[[i, 3]]*a3, {i, 1, Length[wcpRu]}];
Al1site = Table[wcpAl1[[i, 1]]*a1 + wcpAl1[[i, 2]]*a2 + wcpAl1[[i, 3]]*a3, {i,1, Length[wcpAl1]}];
Al2site = Table[wcpAl2[[i, 1]]*a1 + wcpAl2[[i, 2]]*a2 + wcpAl2[[i, 3]]*a3, {i,1, Length[wcpAl2]}];
Al3site = Table[wcpAl3[[i, 1]]*a1 + wcpAl3[[i, 2]]*a2 + wcpAl3[[i, 3]]*a3, {i,1, Length[wcpAl3]}];
Al4site = Table[wcpAl4[[i, 1]]*a1 + wcpAl4[[i, 2]]*a2 + wcpAl4[[i, 3]]*a3, {i,1, Length[wcpAl4]}];
Al5site = Table[wcpAl5[[i, 1]]*a1 + wcpAl5[[i, 2]]*a2 + wcpAl5[[i, 3]]*a3, {i,1, Length[wcpAl5]}];
Alsite = Join[Al1site, Al2site, Al3site, Al4site, Al5site];

原子散乱因子

(データベースからダウンロードしてきたものを読み込む.具体的な数値を出すために使う.)
(* atomic scattering factors *)
dataLa = ReadList["Calculation/RT2Al10/atfacLa", Number, RecordLists -> True];
dataRu = ReadList["Calculation/RT2Al10/atfacRu", Number, RecordLists -> True];
dataAl = ReadList["Calculation/RT2Al10/atfacAl", Number, RecordLists -> True];
f0La = Interpolation[dataLa]
f0Ru = Interpolation[dataRu]
f0Al = Interpolation[dataAl]
Plot[{f0La[x], f0Ru[x], f0Al[x]}, {x, 0, 1.4}]

(* An example to draw a good looking figure *)
fig = Plot[{f0La[x], f0Ru[x], f0Al[x]}, {x, 0, 1.4}, PlotRange -> {{0, 1.5}, {0, All}}, PlotStyle -> Thick, Frame -> True, FrameStyle -> Thickness[0.002], BaseStyle -> Large, FrameLabel -> {"sin\[Theta]/\[Lambda]", Subscript[Style["f", Italic, 22], 0]}, AspectRatio -> 0.75, ImageSize -> 400];figf0 = Show[{Graphics[fig], Graphics[Text[Style[La, 20], {0.2, 52}]],Graphics[Text[Style[Ru, 20], {0.2, 30}]], Graphics[Text[Style[Al, 20], {0.2, 13}]]}]
Export["figf0.png", figf0]

X線のエネルギー,波長,波数

(* X-ray energy, wavelength, wavenumber *)
engy = 7.246;
lambda = 12.398/engy;
kxray = 2*Pi/lambda;

反射指数,逆格子ベクトル,散乱角2θ,結晶構造因子

構造因子の数値も大事だが,どの原子からの散乱がどのように構造因子に入ってくるかをしっかり見るためにも,文字で表しておくことは重要.特に,秩序化を起こす磁性原子はひとくくりにせず,サイトごとに番号をつけて構造因子を書き下し,符号や係数を見ておくのは重要.
{h, k, l} = {6, 2, 0};
kappa = b1*h + b2*k + b3*l;
tau = Sqrt[Sum[kappa[[i]]^2, {i, 1, 3}]];
dd = 2*Pi/tau;
tth = ArcSin[tau/2/kxray]*180/Pi*2;
stfac = Expand[Chop[N[
Sum[fLa*Exp[I*(kappa.Lasite[[i]])], {i, 1,Length[Lasite]}]
+Sum[fRu*Exp[I*(kappa.Rusite[[i]])], {i, 1,Length[Rusite]}]
+Sum[fAl*Exp[I*(kappa.Alsite[[i]])], {i, 1,Length[Alsite]}]
]]] ;
stfacLa = Expand[Chop[N[
fLa1*Exp[I*(kappa.Lasite[[1]])]+fLa2*Exp[I*(kappa.Lasite[[2]])]
+fLa3*Exp[I*(kappa.Lasite[[3]])]+fLa4*Exp[I*(kappa.Lasite[[4]])]
]]];
stfacval = stfac /. {fLa->f0La[tau/4/Pi],fRu->f0Ru[tau/4/Pi],fAl->f0Al[tau/4/Pi]};
Print["kappa=", kappa]
Print["tau=", tau, Tab, "tau/4Pi=", tau/4/Pi, " d=", dd]
Print["two-theta=", tth]
Print["F=", stfac, Tab, Abs[stfacval]^2]
Print["F_La=", stfacLa]

〔注〕最後に構造因子を数値で出すときに,原子散乱因子f0を代入しているが,もちろんここに異常散乱因子を追加で入れれば,より正確な計算になる.磁気散乱や共鳴散乱の原子散乱因子を入れれば,それらの散乱による構造因子が得られる.

YbAl3C3 –– P63/mmc (No. 194) ––

原子位置座標

(* Wyckoff Positions of Atoms *)
zAl4f = 0.13;
zC4f = 0.59;
wcpYb2a = {{0, 0, 0}, {0, 0, 1/2}};
wcpAl2d = {{2/3, 1/3, 1/4}, {1/3, 2/3, 3/4}};
wcpAl4f = {{1/3, 2/3, zAl4f}, {2/3, 1/3, -zAl4f}, {2/3, 1/3, 1/2 + zAl4f}, {1/3, 2/3, 1/2 - zAl4f}};
wcpC2c = {{1/3, 2/3, 1/4}, {2/3, 1/3, 3/4}};
wcpC4f = {{1/3, 2/3, zC4f}, {2/3, 1/3, -zC4f}, {2/3, 1/3, 1/2 + zC4f}, {1/3, 2/3, 1/2 - zC4f}};

格子定数と単位格子ベクトル

(* lattice parameters *)
ahex = 3.39804;
bhex = 3.39804;
chex = 17.1243;
(* unit vectors of the lattice in the xyz-coordinate *)
a1 = {1, 0, 0}*ahex;
a2 = {-1/2, Sqrt[3]/2, 0}*ahex;
a3 = {0, 0, 1}*chex;

逆格子ベクトル

(* reciprocal lattice vectors *)
b1 = Cross[a2, a3]/(a1.Cross[a2, a3])*2*Pi;
b2 = Cross[a3, a1]/(a1.Cross[a2, a3])*2*Pi;
b3 = Cross[a1, a2]/(a1.Cross[a2, a3])*2*Pi;

実空間での原子位置

(* atomic positions in the xyz-coordinate *)
Yb2asite = Table[wcpYb2a[[i, 1]]*a1 + wcpYb2a[[i, 2]]*a2 + wcpYb2a[[i, 3]]*a3, {i, 1, Length[wcpYb2a]}];
Al2dsite = Table[wcpAl2d[[i, 1]]*a1 + wcpAl2d[[i, 2]]*a2 + wcpAl2d[[i, 3]]*a3, {i, 1, Length[wcpAl2d]}];
Al4fsite = Table[wcpAl4f[[i, 1]]*a1 + wcpAl4f[[i, 2]]*a2 + wcpAl4f[[i, 3]]*a3, {i, 1, Length[wcpAl4f]}];
C2csite = Table[wcpC2c[[i, 1]]*a1 + wcpC2c[[i, 2]]*a2 + wcpC2c[[i, 3]]*a3, {i,1, Length[wcpC2c]}];
C4fsite = Table[wcpC4f[[i, 1]]*a1 + wcpC4f[[i, 2]]*a2 + wcpC4f[[i, 3]]*a3, {i,1, Length[wcpC4f]}];

Hexagonalなので,実空間での単位格子ベクトルの書き方が違う以外は,LaRu2Al10の計算と同じである.