原子散乱因子のGaussian展開

原子散乱因子のうち、異常分散の影響を受けない項  f_0 は角度に  \sin\theta \over \lambda (sin theta over lambda で stol と略される)の形で依存する。これを、ガウス関数の和に展開して

 f_0(\theta) = \sum_{i = 0}^{N} a_i e^{-B (\frac{\sin\theta}{\lambda})^2} + c

として扱うことが多い。係数は、電子の波動関数を Gaussian-type orbital を基底にして量子化学計算によって求めた後、フーリエ変換すれば自然に得られる(ガウス関数のフーリエ変換 - biochem_fan's note 参照)。

International Tables for Crystallography Volume C (1992) には、4つのガウス関数の和 + 定数項として展開した数表が載っている。これは減衰が速すぎるので、5つのガウス関数の和 + 定数項に拡張したものが Waasmeier & Kirfel (1995) によって提案されたという。一方生体高分子の結晶学では、多くの場合分解能が 1.5-2.5 Å 程度 (stol 0.2-0.3 に対応)なので、展開の項数を減らしても十分な精度が得られる。そこで、cctbx では 分解能に応じて項数の違う数表を利用している(N-gaussian アプローチ)。このあたりの事情は cctbx 公式サイト にある Newsletter "IUCr Computing Commission No. 3, 2004/01" に解説してあるのを読んで知った。実際に R で 計算して確かめてみよう。

# International Tables for Crystallography Volume C (1992) の表
# cctbx/eltbx/xray_scattering/itc1992.h に数値がある
itc_N_a <-c(12.2126, 3.13220, 2.01250, 1.16630)
itc_N_b <- c(0.005700, 9.89330, 28.9975, 0.582600)
itc_N_c <- -11.529

# Waasmeier & Kirfel (1995) の表
# cctbx/eltbx/xray_scattering/wk1995.h に数値がある
wk_N_a <- c(11.893780, 3.277479, 1.858092, 0.858927, 0.912985)
wk_N_b <- c(0.000158, 10.232723, 30.344690, 0.656065, 0.217287)
wk_N_c <- -11.804902

# cctbx の N-gaussian アプローチ
# cctbx/eltbx/xray_scattering/n_gaussian_raw.cpp に数値がある
n_s <- c(6.0, 5.0, 3.0, 1.7, 0.5, 0.17) # 各展開の分解能限界(stol単位)
n_e <- c(0.00174936368848, 0.00378607883957, 0.00812450990726, 0.00974384919663, 0.00712055235301, 0.00877728869515) # 相対誤差

# 1-6項での展開
n_6 <-c(2.77545321643, 15.0644760293, 1.37595750403, 7.17746883597, 1.06289560478, 0.527446769306,
        1.03805703625, 37.9622771317, 0.625821830249, 0.187618748878, 0.120841771171, 0.0471843880812)
n_5 <- c(3.290377358, 10.3009938633, 1.83751629445, 30.4991869731, 1.00843354744, 0.286891277694,
         0.627115494512, 0.765912546206, 0.233020950338, 0.0682199217083)
n_4 <- c(3.16212239104, 9.94082738119,  1.98555889282, 29.2341683858, 1.07984555802, 0.575668823648,
         0.767239655802, 0.151761276156)
n_3 <- c(2.99954939487, 23.2726795155, 2.25583887941, 7.45433091596, 1.7278842283, 0.316224876669)
n_2 <- c(4.0103185728, 19.9718879528, 2.96436307327, 1.75589053867)
n_1 <- c(6.96715024214, 11.4372299305)

stol <- seq(0, 6, by=0.05)
plot(0, 0, type="n", xlim=c(0, 6), ylim=c(0, 8), 
     xlab="sin(theta)/lambda", ylab="Atomic Scattering Factor", main="Atomic Scattering Factor (f0) of N")
lines(y=sapply(stol, function(x){itc_N_c + sum(itc_N_a * exp(-itc_N_b * x * x))}), x=stol, lwd=3)
lines(y=sapply(stol, function(x){wk_N_c + sum(wk_N_a * exp(-wk_N_b * x * x))}), x=stol, col=2, lwd=3)
lines(y=sapply(stol, function(x){sum(n_1[1] * exp(-n_1[2] * x * x))}), x=stol, col=3)
lines(y=sapply(stol, function(x){sum(n_2[(1:2) * 2 - 1] * exp(-n_2[(1:2) * 2] * x * x))}), x=stol, col=4)
lines(y=sapply(stol, function(x){sum(n_3[(1:3) * 2 - 1] * exp(-n_3[(1:3) * 2] * x * x))}), x=stol, col=5)
lines(y=sapply(stol, function(x){sum(n_4[(1:4) * 2 - 1] * exp(-n_4[(1:4) * 2] * x * x))}), x=stol, col=6)
lines(y=sapply(stol, function(x){sum(n_5[(1:5) * 2 - 1] * exp(-n_5[(1:5) * 2] * x * x))}), x=stol, col=7)
lines(y=sapply(stol, function(x){sum(n_6[(1:6) * 2 - 1] * exp(-n_6[(1:6) * 2] * x * x))}), x=stol, col=8)
legend("topright", c("ITC1992", "WK1995", "1-Gaussian", "2-Gaussian", "3-Gaussian", 
                     "4-Gaussian", "5-Gaussian", "6-Gaussian"), col=1:8, lwd=c(3, 3, 1, 1, 1, 1, 1, 1))

プロット結果は以下の通り。

低分子化合物の結晶など超高分解能領域では International Table of Crystallography の4項展開(黒の太線)が早く減衰しすぎることが分かる。一方、生体高分子の分解能領域では、2, 3項もあれば十分そうだ。