はじめに

座標データ計算用のプログラムをC言語で作成した。逆格子点の座標を計算してファイルに記録する。回折を起こす逆格子点にはflagを立てる。
Tcl/TkによるGUIのフロントエンドを合わせて作成した。gnuplotを利用してグラフを表示する。
計算方法を知りたい人のためにプログラムソースも公開する（洗練されていないのは許してほしい）。

X線回折理論を学ぶとき，逆格子とEwald球と回折の関係は解ってしまえば気持ち良いほどすっきりしているが，最初は解り辛い。
学生さんがイメージできないと訴えるので説明のために絵を描こうとすると，見やすく描くのが難しい。
ええい，コンピューターで3Dで描かせればぐりぐり回せて解りやすかろう，そんなソフトは誰かが作っているだろう，と最初に検索した時には検索の仕方が悪かったようで適当なソフトが見つからなかった。
まあ，無ければ作るか，というのが作成動機となります。

このプログラムのテストバージョンを作った後で落ち着いて検索してみたら，そんなソフトはやはり存在した。例えばXRayViewである。似たような洗練されたソフトはいくつかあった。微妙にやりたいことと違っていたが。逆格子点の座標データが必要でなければそちらを使うのもよいだろう。

必要とする環境

Mac，あるいはWindowsまたはLinux（Unix系ならたぶん何でも）の動くPC

座標の計算だけでよい場合

WindowsとMacは各環境用に用意した計算プログラムのバイナリのみでできる。特にインストールしておかなければならないものはない。
自分の環境に合わせてCのソースをコンパイルする方がよいが，その場合はもちろんCのコンパイラが必要。
Linuxは自分でコンパイルしましょう。

GUIフロントエンドを使用する場合

• gnuplot
  グラフの作成についてはgnuplotに全て任せるようになっているので，WindowsとMacではgnuplotのインストールと環境変数にgnuplotへのパスの登録が必要である。Macの場合はgnuplotを使用するために予めXQuartz[https://www.xquartz.org]のインストールも必要である。
  Linuxではgnuplotがインストールされてなかったら入れる。普通は入れてるよね？
• Tcl/Tk
  フロントエンド部分（Tcl/Tkで記述）は，MacとLinuxでは（通常はTcl/Tkが入っているので）特に何もしなくても動くはずである。
  Windowsで配布のStarpack（Tcl/Tkの実行環境＋実行スクリプト）はそれだけで動くように作成してあるのでTcl/Tkをインストールしなくても動く。
  Starpackが動かない場合はTcl/Tkをインストールするか，Tcl/Tkの実行環境を提供するtclkitを利用して実行スクリプトを動かす。
• 1600x900 (HD+，WXGA++) 以上の画素数を持つモニター

自分でコンパイルする場合　　CのコンパイラはWindowsの場合は選択肢が多いので好みのものを使おう。Macの場合はCommand Line Toolsをインストールすれば（フルセットの開発環境であるXcodeまでは必要ない）よいであろう。
Tcl/Tk or tclkit　　WindowsでStarpackが動かなかった場合，またはMacとLinuxで既に入っているものより新しいバージョンを使用したい場合
Tcl/TkはTcl Developer Xchangeの各バージョンのTcl/TkかActiveStateが配布しているActiveTclをインストールすればよい。Tcl/Tkを使用したいがインストールしたくない場合はtclkitを利用する。

使用法

作業用のディレクトリを準備し，そこに計算プログラムのバイナリ「calcrl.exe」を置く。GUIを使用する場合は，Starpackの「ViewRL.exe」またはTclスクリプトの「ViewRL.tcl」も同じ階層に置く。

コマンドライン入力の表記について

Cのソースファイルから自分でバイナリを作成する場合

コマンドラインでの命令

慣れればこちらが早いし，バッチで少しずつ角度を変えた図を一度に作成する等ができるようになる。

ターミナルやコマンドプロンプトで作業ディレクトリへ移動し，計算の条件を引数で与えると「diffractionRLP.dat」というファイルを書き出す。この中には計算に使った逆格子点の指数と座標，回折を起こす逆格子点かどうかのflag，平板写真に含まれるかどうかのflagと平板写真中の回折点の座標のデータが入っている。データファイルの書式は「詳細」を参照のこと。この「diffractionRLP.dat」にある座標データをgnuplotなど，グラフ作成ソフトで表示させれば回折が起こる様子が見やすい。

Defaultでは格子定数のα，β，γが90.0°，散乱角2θの上限が45°，Mosaicity(Full width)が1.0°になっている。このままで良い場合は入力を省略できる。 逆に格子定数のa，b，c，結晶の方位（定義は下方を参照）φ，χ，ψ，X線の波長λは必ず入力する。

case 1: 単位胞が直交格子で散乱角2θが45°までのとき

引数として格子定数（a,b,cのみ:Å），結晶の方位（定義は下記参照:°），X線の波長（Å）の7個を与える。 例えば格子定数a=4.5, b=7.0, c=12.0, α=β=γ=90.0で結晶の方位が0.0，10.0，20.0，波長0.71の場合，

$ ./calcrl.exe  4.5 7.0 12.0  0.0 10.0 20.0  0.71

とする（Mac，Linuxのbashやtcshの場合）。Windowsのコマンドプロンプトでは「./」を除く。

case 2: 単位胞が直交格子ではなく，散乱角2θの上限が45°までのとき

引数として格子定数，結晶の方位（定義は下記参照），X線の波長の10個を与える。例えば格子定数a=4.5, b=7.0, c=12.0, α=γ=90.0, β=85.0で結晶の方位が0.0，10.0，20.0，波長0.71の場合，

$ ./calcrl.exe  4.5 7.0 12.0  90.0 85.0 90.0  0.0 10.0 20.0  0.71

とする。

case 3: 散乱角2θの上限が45°以外のとき

引数として格子定数，結晶の方位（定義は詳細を参照），X線の波長，散乱角の上限を指定する引数2個の合計12個を与える。 散乱角の上限は1つめの引数が2つめの引数の種類を示す。0の場合は2つめの引数が2θ（°）で，その入力範囲は0.0から180.0となり，入力が範囲を越えている場合は180.0に置き換えられる。1つめの引数が1の場合，2つめの引数は（1/波長）を単位とした逆格子空間の球の半径を表し，その内部の逆格子点にflagを立てる。その入力範囲は0.0から2.5（理論的に回折が起こりうるのは限界球である半径2.0の内部）としている。入力が範囲を越えている場合は2.0に置き換えられる。 例えば格子定数a=4.5, b=7.0, c=12.0, α=β=γ=90.0で結晶の方位が0.0，10.0，20.0，波長0.71，散乱角2θが60.0までの様子を見たい場合，

$ ./calcrl.exe  4.5 7.0 12.0  90.0 90.0 90.0  0.0 10.0 20.0  0.71  0 60.0

とする。

case 4: Mosaicityを変えたいとき

引数として格子定数，結晶の方位（定義は詳細を参照），X線の波長，逆格子点の指定，Mosaicityを指定する引数2個の合計14個を与える。 Mosaicity指定の1つめの引数は1とする。2つめの引数はFull Width（°）で，その入力範囲は0.0から10.0となり，入力が範囲を越えている場合は1.0に置き換えられる。1つめの引数が1以外の場合もMosaicityは1.0に置き換えられる。 例えば格子定数a=4.5, b=7.0, c=12.0, α=β=γ=90.0で結晶の方位が0.0，10.0，20.0，波長0.71，散乱角2θが60.0まで，Mosaicityが1.5の場合，

$ ./calcrl.exe  4.5 7.0 12.0  90.0 90.0 90.0  0.0 10.0 20.0  0.71  0 60.0  1 1.5

とする。

gnuplotでの表示

Ewald球と逆格子点と回折を起こす逆格子点と回折光の方向を手っ取り早く見るには，gnuplotを起動してデータをグラフ表示すれば良い。
計算の際に波長をλ=0.71（モリブデンのKα線）にしていたとする。

> set parametric  
> lambda=0.71  

と入力，（lambdaは実際に計算で使用した波長にする）した後，

> splot (cos(u)*cos(v)-1.0)/lambda, 1.0*sin(u)*cos(v)/lambda, 1.0*sin(v)/lambda title "Ewald Shpere",\
> "diffractionRLP.dat" using 4:5:6 title "Reciprocal lattice point",\
> "diffractionRLP.dat" using ($11>0 ? -$7 : 1/0):(0):(0):($4+$7):5:6 with vector title "Diffraction beam",\
> "diffractionRLP.dat" using ($10>0 ? $4 : 1/0):5:6 title "RLPoint to cause diffraction"

と入力（"\"は「次の行に続く」という意味なので上から３行の最後には必ずタイプする）するとEwald球，逆格子点，回折のwavevector，回折を起こす逆格子点（消滅則は考慮していない）が表示される。

回折イメージ（消滅則は考慮していない）は

> plot "diffractionRLP.dat" using ($11>0 ? $8 : 1/0):9

と入力すれば表示される。gnuplotの操作方法に詳しければ，さらに凝った図も作成できる。

TkによるGUIフロントエンド

Tcl/Tkがインストールしてあれば，ターミナルやコマンドプロンプトで作業ディレクトリに移り，

$ wish ViewRL.tcl

と入力（wishではなくwish8.5やwish8.6等の場合もある）すればGUIフロントエンドが起動する。
WindowsのStarpackはアイコンのダブルクリックで起動させる。MacのStarpackはダブルクリックではなくターミナルから呼び出す。tclkitの場合，tclkitの実行ファイルの名前が「tclkit.exe」であり，作業ディレクトリに置いてあるとすると，コマンドプロンプトで

$ tclkit.exe ViewRL.tcl

とすれば起動する。
使用方法は各パラメータを設定して座標計算をするなら「Calculate coordinates」ボタンを，グラフを表示するなら「Calc. & View Graph」（座標計算をしてから表示）ボタンか「View Graph」（既にあるデータを表示）ボタンをクリックすれば良い。「Exit」ボタンで終了である。
「Save png & dat」チェックボックスをチェックした状態で「Calc. & View Graph」または「View Graph」をクリックすると，その時の座標データファイルのコピー（「diffractionRLP%04d.dat」）とPNG形式のグラフ（「output%04d.png」）の2つが保存される。この保存されるファイルはフォーマットが「%04d」の共通の連番（「0000」から始まり「9999」まで）が付いている。

詳細

座標

グローバル座標系は

• X線の入射は水平方向で，X軸の負から正の方向とする
• Y軸は水平方向で，X線とは垂直にする
• Z軸は鉛直方向に上が正となる方向にする

と取っている。

グローバル座標の原点と逆格子の原点を一致させる。結晶はX線の入射方向に1/λだけ戻った(-1/λ, 0, 0)に位置すると考える。 回折イメージのみは別で，水平方向をx，垂直方向をyとし，X線の入射方向からスクリーンを眺めている形でビームセンターが原点となるように座標をとる。

結晶の方位はa軸をX軸の正の方向，b軸をXY平面と平行でYの正の方向，c軸がZの正の方向に存在するように置いたところを基準(0, 0, 0)とし，Z軸回りにφ，次にX軸回りにχ，最後にZ軸回りにψ回転したとする（Z-X-Z extrinsic rotations (φ,χ,ψ)）。この結晶の方位φ(gamma), χ(beta), ψ(alpha)は剛体の回転でよく使われるオイラー角（alpha, beta, gamma)とはφ=gamma, χ=beta, ψ=alphaの関係である（z-x'-z'' intrinsic rotations （alpha, beta, gamma)）。

入力パラメータ

   a             double:格子定数/Å
   b             double:格子定数/Å
   c             double:格子定数/Å
   alpha         double:格子定数/°
   beta          double:格子定数/°
   gamma         double:格子定数/°
   phi           double:結晶の方位を示す回転角/°
   chi           double:結晶の方位を示す回転角/°
   psi           double:結晶の方位を示す回転角/°
   lambda        double:X線の波長/Å
   flag_reso     int:逆格子の範囲の種類
   resolim       double:逆格子の範囲/°(flag=0)，逆格子の範囲/(1/lambda) (flag=1)
   flag_mosaic   int:mosaicityの入力flag
   mosaic        double:mosaicity/°

の14個あるが，条件によっては減らすことができる。

case 1: 単位胞が直交格子で逆格子の範囲として散乱角2θの上限が45°までのとき

$ ./calcrl.exe (a) (b) (c) (φ) (χ) (ψ) (λ)

引数が7個の場合は，直交格子で逆格子の範囲として散乱角2θが45°までであると判断する。

case 2: 単位胞が直交格子ではなく，逆格子の範囲として散乱角2θの上限が45°までのとき

$ ./calcrl.exe (a) (b) (c) (α) (β) (γ) (φ) (χ) (ψ) (λ)

引数が10個の場合は，逆格子の範囲として散乱角2θが45°までであると判断する。

case 3: 逆格子の範囲として散乱角2θの上限が45°以外のとき

$ ./calcrl.exe (a) (b) (c) (α) (β) (γ) (φ) (χ) (ψ) (λ) (flag) (resolim)

逆格子点の範囲は，flag=0の場合はresolimは2θ（°）で与える。その値の範囲は0.0から180.0となり，入力が範囲を越えている場合は180.0に置き換えられる。flag=1の場合，resolimは（1/波長）を単位とした逆格子空間の球の半径で与え，その入力範囲は0.0から2.5としている。入力が範囲を越えている場合は2.0（限界球に当たる）に置き換えられる。

例えば格子定数a=4.5, b=7.0, c=12.0, α=β=γ=90.0，結晶の方位が0.0，10.0，20.0，波長0.71，逆格子は限界球(半径=2.0/λ)の内部すべての様子を見たい場合，

$ ./calcrl.exe 4.5 7.0 12.0 90.0 90.0 90.0 0.0 10.0 20.0 0.71 1 2.0

という入力になる。これは

$ ./calcrl.exe 4.5 7.0 12.0 90.0 90.0 90.0 0.0 10.0 20.0 0.71 0 180.0

すなわち2θ=180.0°と同等である。

case 4: Mosaicityを変えたいとき （DefaultでFull widthが1.0°に設定）

Mosaicityは逆格子点の広がりをFull Width（°）で与える。実際の振動写真を考えると，結晶のMosaicityと振動角を足したものをここで与える値にするのが妥当であろう。 Mosaicity指定の1つめの引数に1を置くことでMosaicityの入力がされると判断される。2つめはFull Width（°）で，その入力範囲は0.0から10.0となり，入力が範囲を越えている場合は1.0に置き換えられる。1つめの引数が1以外の場合もMosaicityは1.0に置き換えられる。

$ ./calcrl.exe (a) (b) (c) (α) (β) (γ) (φ) (χ) (ψ) (λ) (flag) (resolim) 1 (mosaic)

という入力になる。
例えば格子定数a=4.5, b=7.0, c=12.0, α=β=γ=90.0で結晶の方位が0.0，10.0，20.0，波長0.71，散乱角2θが60.0まで，Mosaicityが1.5の場合，

$ ./calcrl.exe 4.5 7.0 12.0 90.0 90.0 90.0 0.0 10.0 20.0 0.71 0 60.0 1 1.5

とする。

出力ファイルの書式

「diffractionRLP.dat」はファイル内で

 h k l x y z scale px py flagD flagP phi chi psi

が並ぶ。ここで(x, y, z)はミラー指数h, k, l の逆格子点の座標である。scaleは1/lambdaである。(px, py)はカメラ平面上の座標である。この座標の計算はtan（2θ）を使っており，原点から距離が1.0の円は散乱角2θが45°，距離がsqrt(3)の円は2θが60°となる所を示す。flagDは回折を起こすかどうかのflagでflagD=1で回折を起こす点，0で起こさない点であることを示す。flagPは回折イメージに入るかどうかのflagで，厳密には0°=< 散乱角2θ < 90°のときに1，それ以外のときに0となる。phi，chi，psiは結晶の方位として入力された値である。
消滅則を考慮したい場合は，h，k，lの情報があるので例えばawkで強度が0になる回折点のflagDを0に書き換えれば良い。

gnuplotでの表示

gnuplotを起動して

> set parametric  
> lambda=0.71  

と入力（例としてλ=0.71:モリブデンのΚα線の波長としたが，実際に計算で入力した波長にする）した後，

> splot (cos(u)*cos(v)-1.0)/lambda, 1.0*sin(u)*cos(v)/lambda, 1.0*sin(v)/lambda title "Ewald Sphere", \
> "diffractionRLP.dat" using 4:5:6 title "Reciprocal lattice point" w p ps 0.2, \
> "diffractionRLP.dat" using ($11>0 ? $4 : 1/0):5:6 title "RLPoint to cause diffraction", \
> "diffractionRLP.dat" using (-$7):(0.0):(0.0):(2.0 + $7):(-1.0 * (2.0 + $7) * $8):((2.0 + $7) * $9) with vector title "Wavevector of diffraction", \
> "diffractionRLP.dat" using ($11>0 ? (2.0) : 1/0):(-1.0 * (2.0 + $7) * $8):((2.0 + $7) * $9) title "Diffraction Point"

でX=2.0平面に平板カメラがあった場合の回折イメージ（但し消滅則は考慮していない）を重ねた図ができる。

バッチによるデータ処理

MacとLinuxではデフォルトのシェルがbashなので，以下のようなスクリプトを作成して少しずつ結晶の向きを変えた計算とグラフの作成を一気にさせることができる。WindowsではBATファイルで同様のことができるはずであるが…筆者はDOSに疎いのでご勘弁ください（筆者であればCygwinかMsysを入れてこのスクリプトを走らせる）。

シェルスクリプト command.sh

#!/bin/bash

# 格子定数
C_A=20.0
C_B=40.0
C_C=60.0
ALPHA=90.0
BETA=110.0
GAMMA=90.0
# 結晶の方位 psiはrotdegreeとして変化させる
PHI=20.0
CHI=40.0
rotdegree=0
# X線の波長lambdaを適切に変更する
LAMBDA=1.00
# 散乱ベクトルの最大値
MSV=0.5
# Mosaicity
MOSAIC=0.9

# 最初は上のfor文でうまく動くことを確認する
#for i in `seq 0 1`;do  # for test
for i in `seq 0 359`;do
 rotdegree=`echo $i \* 0.5 |bc`
 COMMAND1=$(printf "./calcrl.exe %.1f %.1f %.1f  %.1f %.1f %.1f  %.1f %.1f %.1f  %f  1 %f 1 %f" $C_A $C_B $C_C $ALPHA $BETA $GAMMA $PHI $CHI $rotdegree $LAMBDA $MSV $MOSAIC)
 COMMAND2=$(printf "cp diffractionRLP.dat diffractionRLP%04d.dat" $i)
 COMMAND3=$(printf "sed s/OUTFILE/out%04d/" $i)
 COMMAND4=$(printf "sed s/IN_LAMBDA/%f/" $lambda)
 $COMMAND1
 $COMMAND2
 cat template.gplt | $COMMAND3 | $COMMAND4 >CMMND.gplt
 gnuplot <CMMND.gplt
done

gnuplotスクリプト template.gplt
command.shにより出力するファイル名が置換されたものがgnuplotに読み込まれる。

# pngcairoで出力する環境がなければ「pngcairo」の部分を「png」に変える。
set terminal pngcairo enhanced size 1500,1200
set output "OUTFILE.png"
lambda=IN_LAMBDA

set border lw 3
set zeroaxis
set parametric
set ticslevel 0
set isosamples 40,40

set title "Ewald Sphere and Reciprocal Lattice"

# ここの範囲は適当な値に変える
set xrange [-1.4:0.8]
set yrange [-1.1:1.1]
set zrange [-0.86:0.86]
set view 70, 310, 1, 1
set xlabel "X"
set ylabel "Y"
set zlabel "Z"
splot (cos(u)*cos(v)-1.0)/lambda, sin(u)*cos(v)/lambda, sin(v)/lambda \
title "Ewald Sphere" w l lw 0.4 lc 2, \
"diffractionRLP.dat" using 4:5:6 title "Reciprocal lattice point" w p lc 1, \
"diffractionRLP.dat" using ($10>0 ? (-$7): 1/0):(0.0):(0.0):($4+$7):5:6 with vector \
lw 2 lc 4 title "Wavevector of Diffraction", \
"diffractionRLP.dat" using ($10>0 ? $4 : 1/0):5:6 w p pt 3 lw 2 lc 3 \
title "RLPoint to cause diffraction"

上の例のスクリプトであれば，ターミナルで

$ bash command.sh

とすれば，Z軸回りに0.5°刻みで結晶の向きを0°から179.5°まで変えたときの逆格子のデータ（diffractionRLP0000.datからdiffractionRLP0359.dat）とその3Dグラフ（out0000.pngからout0359.png）が作成される。数分で終了するであろう。
ffmpegがインストールされていれば

$ ffmpeg -r 5 -i out%04d.png -an video.mp4

とすることでパラパラマンガのように纏めたmp4ムービーもできる。
以下にさらに凝ったスクリプトで作成したムービーの例をリンクする。

小さめの格子定数でMosaicityが変わると

• a=4.0, b=6.0, c=8.0, α=γ=90.0, β=110.0, λ=1.0, mosaicity=1.0 で結晶を回したときの様子: 15.4MB
• a=4.0, b=6.0, c=8.0, α=γ=90.0, β=110.0, λ=1.0, mosaicity=4.0 で結晶を回したときの様子: 18.0MB
• a=4.0, b=6.0, c=8.0, α=γ=90.0, β=110.0, λ=1.0, mosaicity=8.0 で結晶を回したときの様子: 20.1MB

大きめの格子定数（タンパク質などの結晶がこのような感じになる）では

• a=12.0, b=18.0, c=24.0, α=γ=90.0, β=110.0, λ=1.0, mosaicity=0.8 で結晶を回したときの様子: 28.6MB

TkによるGUIフロントエンド

Tcl/Tkがインストールしてあり，作業ディレクトリにViewRLtclとcalcrl.exeが置いてあれば，ターミナルかコマンドプロンプトで

$ wish ViewRL.tcl

とすれば起動する（「wish」ではなく「wish8.5」や「wish8.6」等かもしれない）。WindowsのStarpackの場合はアイコンのダブルクリックで起動する。tclkitを利用する場合は作業ディレクトリで

    [tclkitへのフルパス] ViewRL.tcl

とタイプすれば使用できる。
使用方法は各パラメータを設定して座標計算のみをするなら「Calculate coordinates」ボタンを，グラフを表示するなら「Calc. & View Graph」（座標計算をしてから表示）ボタンか「View Graph」（既にあるデータを表示）ボタンをクリックすれば良い。「Exit」ボタンで終了である。
このGUIフロントエンドではgnuplotのtkcanvasへの出力を利用してキャンバスにグラフを描くようにした。オプションとして「Save png & dat」チェックボックスをチェックした状態で「Calc. & View Graph」または「View Graph」をクリックすると，pngファイルも出力できるようにした。この時，PNG形式のグラフ（「output%04d.png」）とその時の座標データファイルのコピー（「diffractionRLP%04d.dat」）が保存される。この保存されるファイルはフォーマットが「%04d」の共通の連番（起動したところで「0000」から始まり，保存されるごとに1増えていき「9999」まで）が付いている。

回折の起こる逆格子点の選択について:理論と計算式