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


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

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

【結果】

まずは練習として線形関数のフィッティングです。コレは簡単。 
#線形関数のフィッティング

##フィッティングに使うもの
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、線形関数のフィッティング

練習をもう一つ、非線形関数のフィッティングです。コレで雰囲気がつかめてきました。 
#非線形関数のフィッティング

##フィッティングに使うもの
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、非線形関数のフィッティング

ここからは物理量のフィッティングです。
まずは、磁化率。 キュリー則として知られる、温度に反比例する振る舞いです。 
##磁化率のキュリー則のフィッティング

##フィッティングに使うもの
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、キュリー則のフィッティング

次は、電気抵抗率。
電子間相関に起因するフェルミ液体的な電気抵抗は温度の二乗に比例します。
非フェルミ液体では色々な温度依存性が生じるのも面白いですね。 
#電子相関によるフェルミ液体の電気抵抗のフィッティング

##フィッティングに使うもの
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()をうまく使うのが大切のようです。 
#フォノンによる電気抵抗の計算

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

#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モデルなどもあり、 色々な情報を引き出すことができます。
なんちゃってバルク敏感の某分光と違って、こちらは本物のバルク敏感なプローブです。 

#光学伝導のフィッティング

##フィッティングに使うもの
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効果などですね。
##比熱の温度依存性の計算

##フィッティングに使うもの
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).

コメント

コメントを投稿

このブログの人気の投稿

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

イメージ
【イントロ】 自分のアウトプットを、そして自分自身を魅力的に見せること。 世の中を渡っていく上では大切な技術です。 一方で、それがいき過ぎると嘘になってしまうためバランスが大切になります。 そんなバランスがアカデミックの世界で崩れているのではないかという懸念を表するエッセイ論文が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が挙げられて...
続きを読む

ワイが面白いと思った論文はこれだ！～2024年まとめ～

イメージ
 【イントロ】 1年早すぎる。 これが相対性理論てやつですね。 巷ではAIの進化が猛威をふるった1年でした。 今年も2024年にCond-matに報告された論文の中で興味深かった論文をまとめてみました。 毎年面白い論文が多すぎる！ 【方法】 1年間のブログ記事を見返して、面白かった論文を選抜しました。 基準は霊感です。 AI、かかってこいよ！！！ 【本文の前に】 2024年1月から12月の間に、Cond-matには約2.2万本の論文が投稿されたようです（ぶひん調べ） その中から、本ブログにメモした論文は3077本。 約14％の論文が琴線に触れたみたいです。 こいつちょろすぎるな。。。 メモしただけでアブストしかほとんど読んでないですが、まあ毎日ワンピースやコナンの続きがWebで無料アップされてたらチェックしちゃいますよね。。。 【気になった論文たち】 １，一番刺激的だった論文 誘電率とエネルギーギャップの普遍的関係 Universal relation between energy gap and dielectric constant https://arxiv.org/abs/2401.04180 ”量子重み”の提案 Quantum weight https://arxiv.org/abs/2401.13847 構造因子に対する量子幾何による制限 Geometric bound on structure factor https://arxiv.org/abs/2412.02656 昨年の論文の中で個人的に一番興味深かったのは、MITのOnishi&Fuによる一連の研究、特に絶縁体のエネルギーギャップと誘電率の関係を明らかにした研究と、基底状態の量子ゆらぎの大きさの指標となる”Quantum Weight"の概念を提案、そして量子幾何が構造因子を制約することを指摘した３論文です。Onishi&Fuは昨年から、 トポロジカル絶縁体のギャップサイズがチャーン数の大きさによってキャップされる関係を発見 するなど、エネルギーギャップと各種物理量の関係性を明らかにしており、興味深い方向性を切り開いています。 理論で提案されている物理量が光学測定で求まる値なのも良いですね。 物理量の上下限を定める関係性は実用上も重要となると考えられ、今...
続きを読む

2025年3月の気になった論文（完全版）

イメージ
 オギャ！オギャ！３月！ 3月のまとめ ・室温超伝導をみつける10の方法 ・ハミルトニアンとしての深層ニューラルネットワーク ・マイクロソフトのマヨラナQubit批判 ・原子サイト交代磁性の提案 ・ウテテの太極図ボルテックス ・実験データから物理公式を発見するLLMファインマンの開発 ・カオスはトポロジカル超対称性の自発的破れ ・冷却原子ペア相関におけるBCS理論の破れ ・熱電交流ジョセフソン効果 ・Bi2212における擬ギャップと関連する反位相密度波 ・超高速鞍点トロニクス ・BaCo2(AsO4)2 の50年来の謎の解決 ・YBCOの強磁場非弾性X線散乱 ・タイプ３マルチフェロイクスは存在するか？ ・単結晶RuO2は量子振動で見ると非磁性 ・Rh置換Sr2IrO4におけるネマティック隠れた秩序相 ・トポロジカル近藤絶縁体SmB6におけるトポロジカルダイオード効果 ・ウテテの奇パリティ準粒子干渉 ・ニューラルネットワークの自然な量子化 ・完全調和スピンサイクロイドワイル半金属GdAlSi ・マジックアングルグラフェンのV字型超伝導ギャップの観測 ・スピンネマティック相Na2BaNi(PO4)2における四重極子波 ・四層グラフェンのナノSQUID観測によるアイソスピン磁気構造の可視化 ・二層モテテにおける高温無散逸分数化チャーン絶縁体状態の観測 ・7層グラフェンにおける軌道チャーン絶縁体の電場制御 ・FIB加工による歪勾配を利用したSr2RuO4の歪依存性電子構造のARPES ・MnPS3のクソデカ熱ホール効果 ・畳み込みトランスフォーマー波動関数 (CTWF) の提案 ・駆動超流体における超固体的音波モードの観測 ・冷却原子におけるフェルミ加速の観測 ・テルルの非線形カイラル熱電ホール効果 ・BCS理論的Tc公式はBCS的ペア形成メカニズムを意味しない ・量子流体の臨界温度はコヒーレンス長と有効質量に比例する ・Bi2201におけるハーブボゴリュボンの発見 ・1次元銅酸化物のRIXSでみる2スピノン励起は1次元ハバード模型で説明できない ※Amazonアフィにアクセスしてやってください🙇‍♀🙇‍♀🙇‍♀🙇‍♀🙇‍♀🙇‍♀ Kindle本 セール＆キャンペーン https://amzn.to/3Q2hjzT ‐2025/3/31、202...
続きを読む