あなたと、フィットしたい…したくない？


【2020/3/29追記】Git-hub https://github.com/Buhin-Mara/Fittings

【イントロ】
　Publish or Perish.
　研究は、論文を書くまでが研究です。
　つまり、実験(計算)してデータを取得するだけではなく、それが何を意味するか考えるまでが研究です。
　データの意味を考える上で役に立つのがデータのプロット、フィッティングによる理論との比較です。
　そこで、今回はPythonをつかった様々な物理量のプロットとフィッティングを行ってみました。
なお、温度や周波数依存性を観たかったので、パラメータは適当です。手を抜いてこそ・・・

【結果】

まずは練習として線形関数のフィッティングです。コレは簡単。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#線形関数のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#x軸の配列
list_linear_x = range(0,20,2)
 
array_error = np.random.normal(size=len(list_linear_x))
 
array_x = np.array(list_linear_x)
 
#y軸の配列
array_y = array_x + array_error
 
#フィッティング
def linear_fit(x,a,b):
    return a*x+b
 
param, cov=curve_fit(linear_fit,array_x,array_y)
 
array_y_fit=param[0]*array_x + param[1]
 
#フィッティング結果のプロット
plt.plot(array_x, array_y,marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="Fitting")
plt.savefig('linear')
plt.legend()
plt.show()

図1、線形関数のフィッティング

練習をもう一つ、非線形関数のフィッティングです。コレで雰囲気がつかめてきました。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#非線形関数のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#x軸の配列
list_linear_x = range(0,20,2)
 
array_x = np.array(list_linear_x)
 
#y軸の配列
list_y=[]
for num in array_x:
    list_y.append( 1.0*np.exp(num/(1.0+num))+np.random.rand()/10)
 
array_y=np.array(list_y)
 
#非線形関数のフィッティング
def nonlinear_fit(x,a,b):
    return b * np.exp(x / (a+x)  )
 
param, cov=curve_fit(nonlinear_fit,array_x,array_y)
 
list_y_fit=[]
for num in array_x:
    list_y_fit.append(param[1]*np.exp(num/(param[0]+num)))
 
array_y_fit=np.array(list_y_fit)
 
#フィッティング結果のプロット
plt.plot(array_x, array_y,marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="Fitting")
plt.legend()
plt.savefig('nonlinear')
plt.show()

図2、非線形関数のフィッティング

ここからは物理量のフィッティングです。
まずは、磁化率。 キュリー則として知られる、温度に反比例する振る舞いです。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
##磁化率のキュリー則のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#温度Tの配列
list_linear_x = range(1,100,1)
 
array_x = np.array(list_linear_x)
 
#磁化率キュリー則の温度依存性
list_y=[]
for num in array_x:
    list_y.append( 1.0 + 1/num + np.random.rand()/10)
 
array_y=np.array(list_y)
 
#キュリー則のフィッティング
def Curie_fit(x,a,b):
    return a + b/x 
 
param, cov=curve_fit(Curie_fit,array_x,array_y)
 
list_y_fit=[]
for num in array_x:
    list_y_fit.append(param[0] + param[1] /num)
 
array_y_fit=np.array(list_y_fit)
 
#キュリー則のプロット
plt.plot(array_x, array_y, marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="fitting")
plt.title("Curie law")
plt.xlabel('temperature (K)')
plt.ylabel('magnetization(emu)')
plt.legend()
plt.savefig('Curie')
plt.show()

図3、キュリー則のフィッティング

次は、電気抵抗率。
電子間相関に起因するフェルミ液体的な電気抵抗は温度の二乗に比例します。
非フェルミ液体では色々な温度依存性が生じるのも面白いですね。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#電子相関によるフェルミ液体の電気抵抗のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#温度の配列
list_linear_x = range(1,50,2)
 
array_x = np.array(list_linear_x)
 
#電気抵抗の配列
list_y=[]
for num in array_x:
    list_y.append( 1.0 + 0.1*num*num + np.random.rand()/10)
 
array_y=np.array(list_y)
 
#フェルミ液体的電気抵抗のフィッティング
def Fermiliquid_fit(x,a,b):
    return a + b * x * x 
 
param, cov=curve_fit(Fermiliquid_fit,array_x,array_y)
 
list_y_fit=[]
for num in array_x:
    list_y_fit.append(param[0] + param[1] * num * num)
 
array_y_fit=np.array(list_y_fit)
 
#フィッティング結果のプロット
plt.plot(array_x, array_y, marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="Fitting")
plt.title("Fermi-liquid")
plt.xlabel('temperature (K)')
plt.ylabel('resistivity(mohm cm)')
plt.legend()
plt.savefig('Fermi-liquid')
plt.show()
図4、フェルミ液体的電気抵抗のフィッティング

次は、同じく電気抵抗率ですが、フォノン散乱に起因する成分です。
こちらは、Bloch-Gruneisen関数とよばれる関数形を示します。
低温では温度の5乗に比例するため、上記の電子間相関が抵抗の主要因となります。 積分関数でフィッティングする際はScipyのOptimize.leastsq()をうまく使うのが大切のようです。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#フォノンによる電気抵抗の計算
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
from scipy import optimize
from scipy import integrate
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#温度の配列
list_linear_x = range(1,300,1)
 
array_x = np.array(list_linear_x)
 
def Bloch(x):
    return (x**5)/((np.exp(x)-1)*(1-np.exp(-1*x)))
 
#Bloch-Gruneisen関数の配列
list_y=[]
for num in array_x:
    p=integrate.quad(Bloch,0,300/num)
    list_y.append(((num/300)**5)*p[0] +np.random.rand()/100)
 
array_y=np.array(list_y)
 
#ここから最小自乗によるフィッティング
def Bloch_fit(parameter,xs,ys):
    a=parameter[0]
    b=parameter[1]
     
    fxs=[(b*(x/a)**5)*integrate.quad(Bloch,0,a/x)[0] for x in xs]
    residual=ys-fxs 
    return residual
 
#フィッティングの初期値
parameter0=[300,1]
result=optimize.leastsq(Bloch_fit,parameter0,args=(array_x,array_y))
print(result)
a_fit=result[0][0]
b_fit=result[0][1]
 
#フィッティングパラメータを使用した関数の配列
list_y_fit=[]
for num in array_x:
    p=integrate.quad(Bloch,0,a_fit/num)
    list_y_fit.append((b_fit*(num/a_fit)**5)*p[0] )
 
array_y_fit=np.array(list_y_fit)
 
#フィッティング結果のプロット
plt.plot(array_x, array_y, marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="Fitting")
plt.title("Bloch-Gruneisen")
plt.xlabel('Temperature(K)')
plt.ylabel('resistivity(mOhm cm)')
plt.legend()
plt.show()
図5、フォノンによる電気抵抗のフィッティング

Fermi-Dirac関数のフィッティングです。
あるエネルギーの状態は全て埋まっていて、 それ以上は空っぽのような占有状態を教えてくれます。
色々な実験に現れる関数ですが、電子の占有をみる光電子分光などでよく見かけますね。

タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#Fermi-Dirac関数のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#エネルギーの配列
list_linear_x = range(1,100,1)
 
array_x = np.array(list_linear_x)
 
#FD関数の配列
list_y=[]
for num in array_x:
    list_y.append( 1/(np.exp(2*(num-50))+1) + np.random.rand()/60)
 
array_y=np.array(list_y)
 
#FD関数のフィッティング
def Fermi_Dirac_fit(x,a,b):
    return  1/(np.exp(a*(x-b))+1) 
 
param, cov=curve_fit(Fermi_Dirac_fit,array_x,array_y)
 
list_y_fit=[]
for num in array_x:
    list_y_fit.append( 1/(np.exp(param[0]*(num-param[1]))+1))
 
array_y_fit=np.array(list_y_fit)
 
#フィッティング結果のプロット
plt.plot(array_x, array_y, marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="",label="Fitting")
plt.title("Fermi-Dirac")
plt.xlabel('Energy (meV)')
plt.ylabel('f(E)(a. u. )')
plt.legend()
plt.savefig('Fermi-Dirac')
plt.show()

図6、Fermi-Dirac関数のフィッティング

光学伝導度を表すDrude-Lorentzモデルによるフィッティングです。
Drude項が自由電子による応答、Lorentz項がフォノンやバンド間遷移を表しています。 急峻なピークを示すLorentz項は初期値を適切に選んであげないと収束しづらいです。
実際には散乱率のせいでピークはブロードになっていることが多いです。
他にも、ディラック電子による応答や、散乱率に周波数依存性があると仮定する拡張Drudeモデルなどもあり、 色々な情報を引き出すことができます。
なんちゃってバルク敏感の某分光と違って、こちらは本物のバルク敏感なプローブです。

タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#光学伝導のフィッティング
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#周波数の配列
list_linear_x = np.arange(1,200,0.2)
 
array_x = np.array(list_linear_x)
 
#Drude-Lorentz関数の配列
list_y=[]
for num in array_x:
    list_y.append( 1.0 + 1/(1+num**2) \
    + num**2/(5*num**2+((100*100-num**2)**2))\
    + np.random.rand()/60)
 
array_y=np.array(list_y)
 
 
#Drude-Lorentz関数のフィッティング
def Drude_Lorentz_fit(x,a,b):
    return a + 1/(1+x**2) + x**2/(5*x**2+(b**2-x**2)**2) 
 
#初期値の設定,ローレンツピークの位置を指定しないと上手くいかない
a_initial=1
b_initial=100
param, cov=curve_fit(Drude_Lorentz_fit,array_x,array_y,[a_initial,b_initial])
 
list_y_fit=[]
for num in array_x:
    list_y_fit.append( param[0] + 1/(1.0+num**2) \
    + num**2/(5*num**2+(param[1]**2-num**2)**2))
 
array_y_fit=np.array(list_y_fit)
 
#光学伝導度のプロット
plt.plot(array_x, array_y, marker="o",label="Numerical calculation")
plt.plot(array_x, array_y_fit, marker="", label="Fitting")
plt.title("Drude-Lorentz")
plt.xlabel('frequency (cm^-1)')
plt.ylabel('conductivity(Ohm^-1 cm^-1)')
plt.legend()
plt.savefig('Drude')
plt.show()
図7、Drude-Lorentzモデルのフィッティング

最後は超伝導転移に伴う比熱のフィッティングです。
比熱からは超伝導ギャップの大きさや、ギャップのノード(節、ある特定のk方向でギャップがゼロになる部分)の有無、 常伝導状態の電子比熱成分からは電子間相関の強さ(いわゆる重い電子系)を推定できたりと得られる情報がたくさんあります。
ただし、実用的にはフォノン成分の取り除き方が難しいのが難点と思います。 (Tcの低い物質では磁場で超伝導を潰してフォノン成分を見積もり、それを引けば良いです。Tcが高い物質だと、 臨界磁場も大きく完全な常伝導状態のフォノン成分を推定するのが難しいため、フィッティングの任意性が大きくなってしまいます。)
その場合、ゼロ温度成分の磁場依存性からギャップ構造を推定することもあります。
いわゆるVolovik効果などですね。
タイトル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
##比熱の温度依存性の計算
 
##フィッティングに使うもの
from scipy.optimize import curve_fit
from scipy import optimize
from scipy import integrate
import numpy as np
 
## 図示のために使うもの
import matplotlib.pyplot as plt
 
#Temperature
list_linear_T = range(1,100,1)
array_T = np.array(list_linear_T)
 
#gap^2の温度依存性の数値計算
list_gap=[]
for num in array_T:
    def f(d):
        def gap(x):
           return np.tanh(1/num*np.sqrt(x**2+d**2))/np.sqrt(x**2+d**2)
 
             
        fxs=integrate.quad(gap,0,500)
        residual=3-fxs[0] 
        return residual
     
    p=optimize.fsolve(f,0)
    list_gap.append(p[0]**2)
 
gap2=np.array(list_gap)
 
dg2=np.gradient(gap2)
 
#比熱の温度依存性の計算
list_C=[]
for num in array_T:
    def f(E):
        return (np.exp(np.sqrt(E**2+gap2[num-1])/num))/(np.exp(np.sqrt(E**2+gap2[num-1])/num)+1)**2\
        *(E**2+gap2[num-1]-num/2*dg2[num-1])
     
    p=integrate.quad(f,0,1000)
    list_C.append(2/num/num/num*p[0])
 
C_gap=np.array(list_C)
 
#比熱の温度依存性のプロット       
plt.plot(array_T, C_gap, marker="", label="Numerical calculation")
plt.title("Electronic specific heat")
plt.xlabel('Temperature(K)')
plt.ylabel('C_el(J/K)')
plt.legend()
plt.savefig('C_el(T)')
plt.show()

図8、超伝導ギャップの温度依存性

図9、比熱の電子成分の温度依存性

【まとめ】
いろいろな物理量をプロット、フィッティングしました。
Pythonのモジュールをうまく組み合わせれば、大体の関数をフィッティングできるようです。ただし、でてきたパラメータが正しい値かは慎重な考察が必要です。
正しいかどうかの判断ができるのは、今の所人間だけですから…
データを上手に解析して、あなただけの自然の真理に迫りましょう！

【参考文献】
参考にしたWeb記事
・Pythonで非線形関数モデリング
・測定データを積分関数でフィッティングしたいです。
・Scilabで関数フィッティング: 金属の電気抵抗
・importコマンドの使用例(1) 超伝導ギャップと比熱
・Introduction to Superconductivity Chapter 3 The BCS Theory
参考にした論文
・光学伝導度：D.N.Basov et al., "Electrodynamics of correlated electron materials", Rev. Mod. Phys., 83, 471(2011)
・比熱：G. Mu et al., "Low temperature specific heat of the hole-doped Ba0.6K0.4Fe2As2 single crystals", PRB 79, 174501 (2009).

コメント

コメントを投稿

このブログの人気の投稿

LLM技術を利用した常圧高温超伝導体の発見～AIリテラシーを高めるぞ！！！～

イメージ
【イントロダクション】 　超伝導は、固体中の電子の協奏的な相互作用から生じる、創発的な巨視的量子現象である。超伝導はゼロ抵抗とマイスナー効果を特徴とし、基礎研究的な観点だけでなく、超伝導電線、超伝導電磁石、そして最近では量子コンピュータ[1]への応用など工学的な注目も高い。 　超伝導の社会的な応用への課題は、その転移温度である。これまで、高圧下では水素化物において250K付近での超伝導転移が観測[2]されているが、常圧では銅酸化物の135Kが最高である[3]。超伝導体のさらなる活用には、冷却が不要でかつ常圧で利用できる高い転移温度をもつ超伝導体の発見が望まれている。 　本研究では、近年進展の著しいLLM技術を活用し、常圧で超伝導転移が203Kで生じる物質、硫酸化ウラン(USO)の発見を報告する。 【結果】 　本研究では、常圧で高い転移温度を示す銅酸化物に含まれる酸素O、高圧で更に高い転移温度を示す硫化水素H3Sに含まれる硫黄S[4]、そして量子効果に伴う非従来型超伝導の発現が期待されている二テルル化ウランUTe2[5]に含まれるUに着目し、U-S-Oの3元系の物質探索を行った。 　図1に示すように、合成された単結晶試料のX線構造解析を行ったところ、正方晶の結晶構造であることが明らかになった。 図1、X線構造回折。 　得られた回折データをリートベルト法により解析したところ、図2のようにU, S, Oの各原子がランダムに結晶格子点を占有していることが明らかになった。組成比はほぼ１：１：１であることから、本物質がUSOの組成であることを確認した。 図2、USOの結晶構造の模式図。赤がU、青がS、白がOを表す。 　図3には、電気抵抗の温度依存性を示している。電気抵抗は高温では温度に線形な振る舞いを示している。抵抗は203Kで超伝導転移を生じ、測定精度の範囲でゼロ抵抗を示すことを確認した。 図3、電気抵抗の温度依存性。T=203Kの赤い破線は転移温度を表す。 　図4には、磁化率の温度依存性を示している。磁化率は高温では磁気転移は確認できず、常磁性体的な振る舞いを示した。電気抵抗と同じく、203K付近から磁化率の低下が生じ、160K以下でほぼ体積分率100%の反磁性を示した。 図4、磁化率の温度依存性。T=203Kの赤い破線は転移温度を表す。 【議論】 　電気抵抗の温度依存性に...
続きを読む

拙著はリジェクト また…流れるでござるよ～逆オッカムの剃刀を使って論文を魅力的に～

イメージ
【イントロ】 自分のアウトプットを、そして自分自身を魅力的に見せること。 世の中を渡っていく上では大切な技術です。 一方で、それがいき過ぎると嘘になってしまうためバランスが大切になります。 そんなバランスがアカデミックの世界で崩れているのではないかという懸念を表するエッセイ論文がNature Physics誌に発表されました。 　そこで本記事では、アメリカの物性物理学者 Igor Mazin によって発表されたこのエッセイの概要を説明してみたいとおもいます。 　もちろん現実逃避です。 逆刃刀 と逆オッカムの剃刀の違い（ぜんぜん違う） 【方法】 Nature Physics誌に発表された以下論文（は読めないのでそのArxiv版）を頑張って読みました。 2ページくらいの短いエッセイなので気になった方はぜひ読んで見てください。 I.Mazin, Inverse Occam’s razor , Nature Physics volume 18, pages367–368 (2022) （Arxiv版： arXiv:2204.08284 ） 【概要】 　「 オッカムの剃刀 」とは、 「ある事柄を説明するためには、必要以上に多くを仮定するべきでない」 という指針、思考の原理のことです。ガリレオ以来、これまでの科学はこの原理に従って進歩してきたとMazinは考えます。すなわち、 「2つの競合する理論がある場合、シンプルな説明のほうが好ましい」 という考え方です。 　ところが、ここ10数年近くこの原理よりもむしろ「逆オッカムの剃刀」と呼ばれる思考原理が採用されるようになりました。これは 「2つの競合する解釈が存在する場合、よりエキゾチックな解釈のほうが好ましい」 という考え方です。この風潮は新しいものを好む人間の本性というよりも、 「エキゾチックな解釈な方がNatureやScienceといったより高インパクトファクターな雑誌に掲載される」 という理由によるものです。資本主義の結果ですね。 　例えば線形磁気抵抗や異常ホール効果を観測したとき、そのメカニズムを ディラックバンド分散 に割り当てるか、 その他の従来型バンド分散 に割り当てるかという解釈の選択を迫られたとき、多くの研究者はよりエキゾチックな前者の解釈を採用するのです。 　また別の例として、Sr2RuO4が挙げられて...
続きを読む

交代磁性があるのかないのかどっちなんだい！～RuO2の場合～

イメージ
 【イントロ】 博士「はぁ～交代じゃ、交代じゃ」 学生「どうしたんですか？キャバクラですか？」 博士「違う違う違う違う今日は違う！交代は交代でも、交代磁性のことを調べておったんじゃ」 学生「交代磁性？」 博士「なんじゃ知らんのか？交代磁性とは、強磁性と反強磁性に次ぐ第三の磁性として最近発見された、新しい磁気状態のことなんじゃ。」 学生「知ってますよ、 Šmejkal（流暢な発音）が提案 してるやつでしょ？強磁性体のようにバンド構造にスピン分裂が生じるけれど、反強磁性体のようにNetの磁性はゼロとなっている、副格子間の反強磁性状態と解釈できる磁性体のことですよね。」 博士「そうじゃ、Šmejkal（流暢な発音）達の提案が大きなきっかけとなった研究分野じゃな。なんじゃ、知っとるじゃないか」 学生「当然です。中学入試程度です。交代磁性って、 MnTe や CrSb 、それに RuO2 で確立された概念じゃないんですか？」 よく見ると(b)がない（ 引用元 ） 博士「 どこ中出身 じゃ。。。うむ、それらの物質でARPESが行われ、交代磁性の特徴の1つであるバンド構造の巨大なスピン分裂が観測されて Nature や PRB に報告されるなど研究が加熱しておる。」 学生「ならいいんじゃないですか？流行物質どんとこい。ビッグウェーブに乗って論文を量産するタイミングですよ。h-index稼ぎましょう」 博士「しかし、そのうちの１つ、RuO2が実は交代磁性体ではなく、非磁性体なのではないかという報告がなされ議論を呼んでおるのじゃ。」 学生「どういうことですか？？？RuO2は 異常ホール効果 が見えたり、 REX（共鳴弾性X線散乱） や 中性子散乱実験 で反強磁性秩序がすでに観測されたんじゃないんですか？？？陰謀ですか？？？？」 博士「それでは今回は、独ヴュルツブルク大学のPhilipp Keßler博士らによる論文、 Philipp Keßler, et al. , Absence of magnetic order in RuO2: insights from μSR spectroscopy and neutron diffraction , arXiv:2405.10820 を読んでいくかのお」 学生「Absence論文、盛り上がるやつですね」 【どんな実験し...
続きを読む