XMCDのSum ruleをpythonで処理

とあるスクールで、講師が事前に用意したOriginのマクロを使ってXMCDのSum ruleから磁気モーメントを求めるという講習があったが、Originを使いたいと全く思わなかったため、一人でpython スクリプトを黙々と作っていた。

ローレンツ関数やら誤差関数やらで適当に作った架空のXASとXMCD。
f:id:koideforest:20181010044729p:plain
これを処理してSum ruleを適用する。

手順としては

pre edgeを差し引く。
引いたものをpost edgeで割る。
XMCDを（再度）作る。
ピークの位置を求める。
step関数を作る。
step関数をXASから引いて積分する。
$L_3$ と $L_2$ の境目を決めてXMCDを積分する。

という感じだろう。

pre edge、post edgeを以下の関数を使って、求める。

import numpy as np

def pre_edge( x, y, start, end ):
    pre_x = []
    pre_y = []
    for ix, xx in enumerate( x ):
        if start < xx:
            pre_x.append( xx )
            pre_y.append( y[ ix ] )
        if end < xx:
            break
    pre_x = np.array( pre_x )
    pre_y = np.array( pre_y )
    A = np.c_[ pre_x, np.ones( len( pre_x ) ) ]
    a, b = np.linalg.lstsq( A, pre_y, rcond=None )[0]
    return a * x + b

def post_edge( x, y, start, end ):
    post_x = []
    post_y = []
    for ix, xx in enumerate( x ):
        if start < xx:
            post_x.append( xx )
            post_y.append( y[ ix ] )
        if end < xx:
            break
    post_x = np.array( post_x )
    post_y = np.array( post_y )
    A = np.c_[ post_x, np.ones( len( post_x ) ) ]
    a, b = np.linalg.lstsq( A, post_y, rcond=None )[0]
    return a * x + b

要は、各エネルギー領域を指定してあげて、その中のデータ点を用いて最小二乗法で直線を定めるというもの。

今の場合だと、

pre edge: 680 ~ 700 eV
pre edge: 760 ~ 800 eV

くらいで設定した。
とりあえずプラスhelicityのXASにのみ適用してみる。

# energy: array for energy points
# spec_p: array for XAS spectrum with plus helicity

pre   = pre_edge( energy, spec_p, 680, 700 )
norm0 = spec_p - pre

post = post_edge( energy, norm0, 760, 800 )
norm_p = norm0 / post

fig = plt.figure( figsize=( 8, 2.5 ) )
ax1 = fig.add_subplot( 1, 3, 1 )
ax2 = fig.add_subplot( 1, 3, 2 )
ax3 = fig.add_subplot( 1, 3, 3 )

ax1.plot( energy, spec_p, "r", label="spec_p" )
ax1.plot( energy, pre, "green" )
ax1.set_xlabel( "Energy [eV]")

ax2.plot( energy, norm0, "r", label="norm0" )
ax2.plot( energy, post, "purple" )
ax2.set_xlabel( "Energy [eV]")

ax3.plot( energy, norm_p, "r", label="norm_p" )
ax3.set_xlabel( "Energy [eV]")

plt.tight_layout()
plt.show()

f:id:koideforest:20181010051755p:plain
うまく引けている感じがする。

規格化後のXASとXMCD。

xas = ( norm_p + norm_m ) / 2.
xmcd = norm_p - norm_m

f:id:koideforest:20181010053707p:plain

ピークの位置を、次の関数で求める。
最初のピークは710 ~ 730 eV、次のピークは730 ~ 750 eVの範囲にあるので、そのように設定。

def white_line_position( x, y, x0, x1 ):
    p = max([ yy for xx, yy in zip( x, y ) if x0 < xx < x1 ])
    for ix, xx in enumerate( x ):
        if x0 < xx < x1:
            if abs( y[ ix ] - p ) < 1e-7:
                break
    return ix, xx

_, e1 = white_line_position( energy, xas, 710, 730 )
_, e2 = white_line_position( energy, xas, 730, 750 )
print( e1, e2 )
# 720.36..., 740.54...

d状態以外の寄与を差し引くためのstep関数を次のように用意。

from scipy.special import erf

def step_func( x, x0, d ):
    return ( erf( ( x - x0 ) / d ) + 1. ) / 2.

def subtract_func( x, x01, x02, d, c ):
    return c * step_func( x, x01, d ) + ( c / 2. ) * step_func( x, x02, d)

元々の誤差関数は $\pm1$ の値を持つため、step_funcでは $[0,1$ ]になるように調整した。
強度比は軌道モーメントがなければ $L_3/L_2 = 2$ であるため、大きさを決める係数 $c$ は一つで良い。（この関数は、（遷移金属原子における）d状態以外の寄与のためのモデルであり、それらは一般に軌道モーメントは小さい。）

求めたピークの位置を使って、step関数を作る。
step関数のfitting領域は760 eV以降のpost edgeに設定。

def index_energy( energy, e0 ):
    for i, e in enumerate( energy ):
        if e >= e0:
            if abs( energy[ i ] - e0 ) <= abs( energy[ i - 1 ] - e0 ):
                return i
            else:
                return i - 1

i_post = index_energy( energy, 760 )

from scipy.optimize import curve_fit

c, _ = curve_fit(
    lambda e, c: subtract_func( e, e1, e2, de, c ),
    energy[ i_post : ],
    xas[ i_post : ]
    )
subtract = subtract_func( energy, e1, e2, de, c )

plt.plot( energy, subtract )
plt.plot( energy, xas )
plt.xlabel( "Energy [eV]" )
plt.ylabel( "Intensity [arb.unit]" )
plt.show()

f:id:koideforest:20181021055450p:plain
それらしく階段関数が置けている。
curve_fitの使い方は以下のサイトを参照。
pythonでフィッティングをする - おっぱいそん！

階段関数を差し引いたXAS、そしてXMCDを次のように積分する。
ここではSimpson法を用いる。

from scipy.integrate import simps

def integrate( x, y ):
    result = [ 0. ]
    for ix in range( len( x ) )[ 1: ]:
        result.append( simps( y[ :ix ], x[ :ix ] ) )
    return np.array( result )

int_xas = integrate( energy, xas - subtract )
int_xmcd = integrate( energy, xmcd )
plt.plot( energy, int_xas )
plt.plot( energy, int_xmcd )
plt.xlabel( "Energy [eV]" )
plt.show()

print( int_xas[-1] )
# 18.914...

f:id:koideforest:20181010060456p:plain
ちょっと後ろの方でフニャッとしているが、こんなもんだろう。

$L_3$ と $L_2$ の境目は、XASの一回微分がゼロになるところを持って来た。

def get_index( xs, x0 ):
    for i, x in enumerate( xs ):
        if x >= x0:
            if abs( xs[ i ] - x0 ) <= abs( xs[ i - 1 ] - x0 ):
                return i
            else:
                return i - 1

def index_zero( x, y, x_start, x_end ):
    iy_start = get_index( x, x_start )
    iy_end   = get_index( x, x_end )
    y_zero       = min( abs( y[ iy_start : iy_end ] ) )
    for iy, yy in enumerate( y ):
        if not iy_start <= iy <= iy_end: continue
        if abs( yy ) - y_zero <= 1e-7:
            return iy

i_zero = index_zero( energy[ 1 : ], np.diff( xas, n=1 ), 730, 735 )
plt.plot( energy, xas )
plt.plot( energy[ im ], xas[ im ], "o" )
plt.show()

print( int_xmcd[ i_zero ] )
# -4.5744...

f:id:koideforest:20181021064816p:plain
差分による近似的に微分した配列が得られるnp.diffの使い方は、以下のサイトを参照。
要素の差分、足し合わせを計算するNumPyのdiff関数とcumsum関数の使い方 - DeepAge

ひとまずはこれでSum ruleに必要なXAS及びXMCDの積分値を求めることが出来た。

ちなみに、適当に作ったスペクトルは以下のスクリプトを使用。

def lorentzian_norm( x, x0, gamma ):
    return gamma**2 / ( ( x - x0 )**2 + gamma**2 ) 

energy = np.linspace( 680, 800, 1000 )
de = energy[1] - energy[0]

peak1 = 720
peak2 = 740

back  = ( 0.01 * ( energy - 900 ) )**2
base  = subtract_func( energy, peak1, peak2, 10 * de, 1 )

peak_p  = 2 * lorentzian_norm( energy, peak1, 2 )
peak_p += 1 * lorentzian_norm( energy, peak2, 2 )
spec_p = base + peak_p + back

peak_m  = 3.0 * lorentzian_norm( energy, peak1, 2 )
peak_m += 0.5 * lorentzian_norm( energy, peak2, 2 )
spec_m = base + peak_m + back

ローレンツ関数の大きさ等は以下を参照。
ローレンツ関数とは？

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

XMCDのSum ruleをpythonで処理