実装して学ぶデータ処理システムその1 - XDS の座標系

XDS の主要なアルゴリズムを手を動かして実装することを通じて、X 線回折のデータ処理について理解を深めよう。サンプルデータは Thaumatin / Diamond Light Source I04 user training を用いることにする。

まず、XDS の座標系を学ぶ。nXDS は回転軸の概念がなくなった分、単純に見えて、計算がややこしくなっているので、通常の XDS を対象にする。リファレンスは Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 だ。今回は、2.1 節を実装する。

回折を記述するのに必要なパラメータは、XPARM.XDS に全て含まれているので、これを R でパースする。

# Utility functions 便利関数

outer_prod <- function(x, y) {
  c(x[2] * y[3] - x[3] * y[2], 
    x[3] * y[1] - x[1] * y[3], 
    x[1] * y[2] - x[2] * y[1])
}

inner_prod <- function(x, y) sum(x * y)
vec_norm <- function(x) sqrt(sum(x * x))

# Rodrigues rotation formula に基づく
# https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula 
rotate_vec <- function(axis, phi, vec) {
  ip <- inner_prod(axis, vec)
  op <- outer_prod(axis, vec)
  axis * ip + (vec - axis * ip) * cos(phi) + op * sin(phi)
}

deg2rad <- function(d) d / 180.0 * pi
rad2deg <- function(r) r / pi * 180.0

# 数字として読み込む。最初の行はバージョン情報なので飛ばす

XPARM <- scan("~/Desktop/Datasets/DIALS-thaumatin/xds/XPARM.XDS", skip=1)

#XPARM.XDS    VERSION March 1, 2015

# ROTATION PARAMETERS 回転角度と回転軸 (m2)
#1        0.0000    0.1500  0.999998 -0.001616 -0.001308

IMAGE0 <- XPARM[1]
PHI0 <- XPARM[2]
DELTA_PHI <- XPARM[3]
M2 <- XPARM[4:6]

# BEAM PARAMETERS X 線の波長と入射方向 (S0)
#0.976250       0.001299       0.001023       1.024326

LAMBDA <- XPARM[7]
S0 <- XPARM[8:10]

# Form the goniostat system. m2 と S0 から goniostat 座標系が作られる
M1 <- outer_prod(M2, S0)
M1 <- M1 / vec_norm(M1)
M3 <- outer_prod(M1, M2)

G <- cbind(M1, M2, M3)

# CRYSTAL PARAMETERS 結晶の情報。空間群は含まれない。
#1     57.8654     57.8673    150.1588  90.000  90.013  89.940
#15.074243      53.742535      15.261243
#-19.928106      -9.533589      53.484589
#135.412857     -49.825455      41.574448

B1 <- XPARM[18:20]
B2 <- XPARM[21:23]
B3 <- XPARM[24:26]

B <- cbind(B1, B2, B3)
B_STAR <- solve(B) # 逆格子の基底も求めておく

# 12..17 番目にある格子定数と一致することを確認しておく
vec_norm(B1)
1 / (vec_norm(B_STAR[1,]))
vec_norm(B2)
1 / (vec_norm(B_STAR[2,]))
vec_norm(B3)
1 / (vec_norm(B_STAR[3,]))

# DETECTOR PARAMETERS 検出器の大きさとピクセルサイズ
#1      2463      2527    0.172000    0.172000

NX <- XPARM[28]
NY <- XPARM[29]
QX <- XPARM[30] * 1E7 # mm to A
QY <- XPARM[31] * 1E7

# 検出器の原点
# 1225.349976    1193.469971     265.269989

X0 <- XPARM[32] # ORGX px
Y0 <- XPARM[33] # ORGY px
F <- XPARM[34] * 1E7 # mm to A

# 検出器系の基底ベクトル
#1.000000       0.000000       0.000000
#0.000000       1.000000       0.000000
#0.000000       0.000000       1.000000

# Form the detector system

D1 <- XPARM[35:37] # DIRECTION_OF_DETECTOR_X-AXIS
D2 <- XPARM[38:40] # DIRECTION_OF_DETECTOR_Y-AXIS
D3 <- XPARM[41:43] # D3 <- outer_prod(D1, D2) なのでやや冗長

D <- cbind(D1, D2, D3)

# SEGMENT PARAMETERS マルチセグメント検出器については今回はパス
#1         1      2463         1      2527
#0.00    0.00    0.00  1.00000  0.00000  0.00000  0.00000  1.00000  0.00000

ここから分かるように、XPARAM.XDS の項目の順序もよく考えられているのだ。

ところで、なぜ、laboratory system, goniostat system, detector system(, segment system) という複数の座標系を用意する必要があるのか。詳しくは精密化のところで明らかになるが、パラメータの依存性を整理するためである。例えば、回折がどの方向( $\mathbf{S}$ )に起こるかは、回転後の結晶の基底ベクトルの向きと、入射ビームの波数ベクトルだけで決まり、検出器とは無関係である。なにしろ検出器が存在しなくても、回折は同じように起こるのだから。そのため、 $\mathbf{S}$ ベクトルを計算する部分と、それが検出器のどこで観測されるか―を別々に扱ったほうが都合が良いのである。

それでは、この座標系を使ってみよう。INTEGRATE.HKL から適当に一行抜き出して、それを再現することにする。

# INTEGRATE.HKL から適当な一行を取ってきた
#  -31 4 -35 6.284E+00 2.354E+00 554.2 435.1 23.8 0.48608 100.00000 60 5 554.6 434.2 23.5 39.83 0.10 -131.51 33.29 8.73 1

hkl <- c(-31, 4, -35)
xcalc <- 554.2
ycalc <- 435.1 
zcalc <- 23.8

xy2S <- function(x, y) {
# 検出器系座標から、laboratory system へ
  s <- ((x - X0) * QX * D1 + (y - Y0) * QY * D2 + F * D3)
  s / vec_norm(s) / LAMBDA
}

hkl2P0_STAR <- function(h, k, l) {
# hkl から laboratory system へ
  B_STAR[1, ] * h + B_STAR[2, ] * k + B_STAR[3, ] * l
}

# 一致することを確認する
p0_star <- hkl2P0_STAR(hkl[1], hkl[2], hkl[3]) # 回転前の逆空間座標

# zcalc に基づいて回転。zcalc の単位は DELTA_PHI 度であることに注意
p_star <- rotate_vec(M2, deg2rad(PHI0 + zcalc * DELTA_PHI), p0_star)

S_calc <- p_star + S0
S_calc_from_XY <- xy2S(xcalc, ycalc)

# 確認。5E-4 くらいの精度があるようだ
all.equal(1 / vec_norm(S0), 1 / vec_norm(S_calc))
all.equal(S_calc, S_calc_from_XY)

なお、上の計算では実験室系の座標でベクトルの成分を表示して計算した。しかし、goniostat system で計算したほうが話が単純になる。なぜなら、回転軸 $\mathbf{m_2}$ は goniostat system の第2基底なので、その周りの回転によって第1成分と第3成分だけが変化する(= $\mathbf{m_1, m_3}$ 平面上での回転となる)からである。これが 2.2 節の前半部分である。

次回、ZCALC の求め方を説明する。# 図を作るのが面倒……。コードとしては Kabsch transform までできているので、ゆっくり進めます