nano_exit

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

XMCDのSum ruleをpythonで処理

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

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

手順としては

  1. pre edgeを差し引く。
  2. 引いたものをpost edgeで割る。
  3. XMCDを(再度)作る。
  4. ピークの位置を求める。
  5. step関数を作る。
  6. step関数をXASから引いて積分する。
  7.  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 )
ax1.set_xlabel( "Energy [eV]")

ax2.plot( energy, norm0, "r", label="norm0" )
ax2.plot( energy, post )
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。
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関数を作る。

c = 0.67
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()

本当は係数 cをフィッティングによってきちんと決めたいが、ここでは簡単のため手で適当に合わせた値 c=0.67を採用する。
f:id:koideforest:20181010055019p:plain
それらしく階段関数が置けている。

階段関数を差し引いた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 - subract )
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] )
print( int_xmcd[ 483 ] )  # ~ 738 eV
# 18.914...
# - 4.3069...

f:id:koideforest:20181010060456p:plain
ちょっと後ろの方でフニャッとしているが、こんなもんだろう。
 L_3 L_2の境目は、力尽きてしまったので、目で見て大体の位置の値を取って来てしまっているが、730~740 eVの間で二階微分が最大のところを持ってくるのが良さそうに思う。
ひとまずはこれで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

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