補足３：じゅりあ「そのままで…いいよ？」ぱいそん「コンプレックス、あるんだ…」

tight-binding5

元記事はこちら
http://buhin-blog.blogspot.jp/2018/03/blog-post_29.html

In [1]:
#5. グラフェンと呼ばれるタイトバインディング模型：アームチェア版
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt

#アームチェア型グラフェンの結晶構造
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-3,4)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-4,4)]
x0= [2*j*1.5*a for j in range(-2,3)]
x1= [(2*j+1)*1.5*a for j in range(-2,2)]

c="red"
lw=0.5
ls="--"

for y in y0:
    for x in x0:
        plt.plot([x,x+a               ],[y,y  ],color=c, lw=lw, ls=ls)
        plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
        plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)

for y in y1:
    for x in x1:
        plt.plot([x,x+a               ],[y,y  ],color=c, lw=lw, ls=ls)
        plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
        plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
        
plt.gca().set_aspect("equal")
plt.xlim([-3,3])
plt.ylim([-2.5,3])
Out[1]:
(-2.5, 3)
In [2]:
#2種類の原子の色分け
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-2,3)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-2,3)]
x0= [2*j*1.5*a for j in range(-1,3)]
x1= [(2*j+1)*1.5*a for j in range(-1,2)]

c="red"
b="blue"
lw=0.5
ls="--"

for y in y0:
    for x in x0:
        plt.plot([x,x+a               ],[y,y  ],color=c, lw=lw)
        plt.plot([x],[y],color=b,marker="o")
        plt.plot([x+a],[y],color=c,marker="o")
        plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
        plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)

for y in y1:
    for x in x1:
        plt.plot([x,x+a               ],[y,y  ],color=c, lw=lw)
        plt.plot([x],[y],color=b,marker="o")
        plt.plot([x+a],[y],color=c, marker="o")
        plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
        plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
        
plt.gca().set_aspect("equal")
plt.xlim([-1.8,3])
plt.ylim([-2.2,2])
Out[2]:
(-2.2, 2)
In [3]:
#ユニットセルのプロット
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-2,3)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-2,3)]
x0= [2*j*1.5*a for j in range(-1,3)]
x1= [(2*j+1)*1.5*a for j in range(-1,2)]

c="red"
b="blue"
lw=0.5
ls="--"

x=0
y=0

plt.plot([x,x+a               ],[y,y  ],color=c, lw=lw)
plt.plot([x],[y],color=b,marker="o")
plt.plot([x+a],[y],color=c,marker="o")

plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw)

plt.plot([x-a/2],[y-np.sqrt(3)*a/2],color=c,marker="o")
plt.plot([x-a/2],[y+np.sqrt(3)*a/2],color=c,marker="o")

plt.plot([x+a,x-a/2+2*a               ],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x+a,x-a/2+2*a               ],[y,y+np.sqrt(3)*a/2],color=c, lw=lw)

plt.plot([x-a/2+2*a],[y-np.sqrt(3)*a/2],color=b,marker="o")
plt.plot([x-a/2+2*a],[y+np.sqrt(3)*a/2],color=b,marker="o")
        
plt.gca().set_aspect("equal")
plt.xlim([-0.5,1])
plt.ylim([-1,1])
Out[3]:
(-1, 1)
In [4]:
#アームチェア型グラフェンのハミルトニアン
def calc_HGrapheneA(Nx,Ny,μ):
    N=Nx*Ny*4
    mat_Htb=np.zeros([N,N])
    mat_Htb += (-μ)*np.eye(N)
    t=1.0
    
    for ix in range(1,Nx+1):
        for iy in range(1,Ny+1):
            for dx in range(-1,2):
                for dy in range(-1,2):
                    jx=ix+dx
                    jx+=-Nx if jx>Nx else 0
                    jx+=Nx if jx<1 else 0
                    jy=iy+dy
                    jy+=-Ny if jy>Ny else 0
                    jy+=Ny if jy<1 else 0
                    
                    for a in range(1,5):
                        for b in range(1,5):
                            hop=False
                            ii=((iy-1)*Nx+ix-1)*4+a
                            if a==1:
                                if dx==-1 and dy==0 and b==4:
                                    hop = True
                                elif dx==0 and dy==1 and b==2:
                                    hop=True
                                elif dx==0 and dy==0 and b==2:
                                    hop=True
                            elif a==2:
                                if dx==0 and dy==0 and b==1:
                                    hop = True
                                elif dx==0 and dy==0 and b==3:
                                    hop=True
                                elif dx==0 and dy==-1 and b==1:
                                    hop=True
                            elif a==3:
                                if dx==0 and dy==0 and b==2:
                                    hop = True
                                elif dx==0 and dy==0 and b==4:
                                    hop=True
                                elif dx==0 and dy==-1 and b==4:
                                    hop=True
                            elif a==4:
                                if dx==0 and dy==0 and b==3:
                                    hop = True
                                elif dx==1 and dy==0 and b==1:
                                    hop=True
                                elif dx==0 and dy==1 and b==3:
                                    hop=True
                                    
                            if hop and 1<=jx<=Nx and 1<=jy<=Ny:
                                jj=((jy-1)*Nx+jx-1)*4+b
                                
                                mat_Htb[ii-1,jj-1]=t
                                
    return mat_Htb
In [5]:
#アームチェア型グラフェンのエネルギー固有値
Nx=20
Ny=40
μ=0

mat_H=calc_HGrapheneA(Nx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
plt.hist(energy,bins=100,label="y1",ec="black")
Out[5]:
(array([  27.,   24.,   28.,   32.,   28.,   24.,   36.,   14.,   38.,
          34.,   26.,   30.,   34.,   26.,   32.,   36.,   24.,   34.,
          34.,   36.,   44.,   32.,   30.,   38.,   40.,   48.,   32.,
          42.,   46.,   44.,   60.,   50.,   50.,  105.,   52.,   56.,
          30.,   32.,   32.,   26.,   24.,   12.,   24.,   16.,   12.,
          10.,    2.,   10.,    2.,    2.,    2.,    2.,   10.,    2.,
          10.,   12.,   16.,   24.,   12.,   24.,   26.,   32.,   32.,
          30.,   56.,   52.,  105.,   50.,   50.,   60.,   44.,   46.,
          42.,   32.,   48.,   40.,   38.,   30.,   32.,   44.,   36.,
          34.,   34.,   24.,   36.,   32.,   26.,   34.,   30.,   26.,
          34.,   38.,   14.,   36.,   24.,   28.,   32.,   28.,   24.,   27.]),
 array([-3.  , -2.94, -2.88, -2.82, -2.76, -2.7 , -2.64, -2.58, -2.52,
        -2.46, -2.4 , -2.34, -2.28, -2.22, -2.16, -2.1 , -2.04, -1.98,
        -1.92, -1.86, -1.8 , -1.74, -1.68, -1.62, -1.56, -1.5 , -1.44,
        -1.38, -1.32, -1.26, -1.2 , -1.14, -1.08, -1.02, -0.96, -0.9 ,
        -0.84, -0.78, -0.72, -0.66, -0.6 , -0.54, -0.48, -0.42, -0.36,
        -0.3 , -0.24, -0.18, -0.12, -0.06,  0.  ,  0.06,  0.12,  0.18,
         0.24,  0.3 ,  0.36,  0.42,  0.48,  0.54,  0.6 ,  0.66,  0.72,
         0.78,  0.84,  0.9 ,  0.96,  1.02,  1.08,  1.14,  1.2 ,  1.26,
         1.32,  1.38,  1.44,  1.5 ,  1.56,  1.62,  1.68,  1.74,  1.8 ,
         1.86,  1.92,  1.98,  2.04,  2.1 ,  2.16,  2.22,  2.28,  2.34,
         2.4 ,  2.46,  2.52,  2.58,  2.64,  2.7 ,  2.76,  2.82,  2.88,
         2.94,  3.  ]),
 <a list of 100 Patch objects>)
In [6]:
#片側フーリエ変換したハミルトニアン
def calc_HGrapheneAkx(kx,Ny,μ):
    N=Ny*4
    mat_Htb=np.zeros([N,N],dtype=np.complex )#複素数の行列
    mat_Htb += (-μ)*np.eye(N)
    t0=1.0
    
    #for ix in range(1,Nx+1):
    for iy in range(1,Ny+1):
        #for dx in range(-1,2):
            for dy in range(-1,2):
                #jx=ix+dx
                #jx+=-Nx if jx>Nx else 0
                #jx+=Nx if jx<1 else 0
                jy=iy+dy
                jy+=-Ny if jy>Ny else 0
                jy+=Ny if jy<1 else 0
                t=0.0
                
                for a in range(1,5):
                    for b in range(1,5):
                        hop=False
                        t=0.0
                        ii=(iy-1)*4+a
                        if a==1:
                            if dy==0 and b==4:
                                hop = True
                                t+=t0*np.exp(-1j*kx)
                            elif dy==1 and b==2:
                                hop=True
                                t+=t0
                            elif dy==0 and b==2:
                                hop=True
                                t+=t0
                        elif a==2:
                            if dy==0 and b==1:
                                hop = True
                                t+=t0
                            elif dy==0 and b==3:
                                hop=True
                                t+=t0
                            elif dy==-1 and b==1:
                                hop=True
                                t+=t0
                        elif a==3:
                            if dy==0 and b==2:
                                hop = True
                                t+=t0
                            elif dy==0 and b==4:
                                hop=True
                                t+=t0
                            elif dy==-1 and b==4:
                                hop=True
                                t+=t0
                        elif a==4:
                            if dy==0 and b==3:
                                hop = True
                                t+=t0
                            elif dy==0 and b==1:
                                hop=True
                                t+=t0*np.exp(1j*kx)
                            elif dy==1 and b==3:
                                hop=True
                                t+=t0
                                    
                        if hop and 1<=jy<=Ny:
                            jj=(jy-1)*4+b
                                
                            mat_Htb[ii-1,jj-1]=t
                                
    return mat_Htb
In [7]:
#kxを横軸としてエネルギー固有値をプロット
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*4])
cnt=0

for kx in vkx:
    cnt+=1
    mat_H=calc_HGrapheneAkx(kx,Ny,μ)
    energy,mat_v=np.linalg.eigh(mat_H)
    for i in range(1,Ny*4+1):
        #print(energy[i-1])
        ep[cnt-1,i-1]=energy[i-1]
        
plt.plot(vkx,ep)
In [8]:
#端状態を出すためにy方向の周期条件を外したハミルトニアン
def calc_HGrapheneAkxw(kx,Ny,μ):
    N=Ny*4
    mat_Htb=np.zeros([N,N],dtype=np.complex )
    mat_Htb += (-μ)*np.eye(N)
    t0=1.0
    
    #for ix in range(1,Nx+1):
    for iy in range(1,Ny+1):
        #for dx in range(-1,2):
            for dy in range(-1,2):
                #jx=ix+dx
                #jx+=-Nx if jx>Nx else 0
                #jx+=Nx if jx<1 else 0
                jy=iy+dy
                #jy+=-Ny if jy>Ny else 0
                #jy+=Ny if jy<1 else 0
                t=0.0
                
                for a in range(1,5):
                    for b in range(1,5):
                        hop=False
                        t=0.0
                        ii=(iy-1)*4+a
                        if a==1:
                            if dy==0 and b==4:
                                hop = True
                                t+=t0*np.exp(-1j*kx)
                            elif dy==1 and b==2:
                                hop=True
                                t+=t0
                            elif dy==0 and b==2:
                                hop=True
                                t+=t0
                        elif a==2:
                            if dy==0 and b==1:
                                hop = True
                                t+=t0
                            elif dy==0 and b==3:
                                hop=True
                                t+=t0
                            elif dy==-1 and b==1:
                                hop=True
                                t+=t0
                        elif a==3:
                            if dy==0 and b==2:
                                hop = True
                                t+=t0
                            elif dy==0 and b==4:
                                hop=True
                                t+=t0
                            elif dy==-1 and b==4:
                                hop=True
                                t+=t0
                        elif a==4:
                            if dy==0 and b==3:
                                hop = True
                                t+=t0
                            elif dy==0 and b==1:
                                hop=True
                                t+=t0*np.exp(1j*kx)
                            elif dy==1 and b==3:
                                hop=True
                                t+=t0
                                    
                        if hop and 1<=jy<=Ny:
                            jj=(jy-1)*4+b
                                
                            mat_Htb[ii-1,jj-1]=t
                                
    return mat_Htb
In [9]:
#バンド構造。端状態は出ない。
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*4])
cnt=0

for kx in vkx:
    cnt+=1
    mat_H=calc_HGrapheneAkx(kx,Ny,μ)
    energy,mat_v=np.linalg.eigh(mat_H)
    for i in range(1,Ny*4+1):
        #print(energy[i-1])
        ep[cnt-1,i-1]=energy[i-1]
        
plt.plot(vkx,ep)

コメント

コメントを投稿

このブログの人気の投稿

交代磁性があるのかないのかどっちなんだい！～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論文、盛り上がるやつですね」 【どんな実験したの？】
続きを読む

狂気！室温超伝導体は脳内にあった！！！

イメージ
 【イントロ】 博士「これはすごい・・・これはすごいぞ！！！」 学生「どうしたんですか？自分の愚かさのことですか？」 博士「そうではない！ワシは出会ったんじゃ！科学史上最高の論文に！」 学生「どんな論文ですか一体」 博士「なんと 室温超伝導体 が見つかったという話じゃ！」 学生「え！室温超伝導体といえば、室温で電気抵抗がゼロになる物質で、世界を支えるエレクトロニクスに不可避なエネルギーの散逸、発熱の問題を一気に解決し、 世界のエネルギー問題を救うキーテクノロジーになる物質 じゃないですか！」 博士「その通りじゃ！」 学生「気になってきました、早く教えて下さい、その論文！」 博士「うむ…驚くなよ、なんと、 生物の脳の中では超伝導状態が実現しておる のじゃ！！！」 学生「・・・」 博士「ふふふ、あまりの衝撃に言葉もないようじゃな。」 学生「あの」 博士「では、今回は、 P. Mikheenko, Possible Superconductivity in the Brain , Journal of Superconductivity and Novel Magnetism, 2018, No 5, p. 1121-1134 （ Arxiv版 ） を読んでいくぞい！」 学生「帰っていいですか？」 【論文誌がすごい！】 博士「まずは、この論文のなにが凄いかといえば、査読を通って論文誌に掲載されておるのじゃ」 学生「ほんとですよ、どうなってんですか」 博士「そして、掲載された Journal of Superconductivity and Novel Magnetismの2022年度のImpact Factor は1.8なんじゃ。これがどういうことかわかるか？」 学生「そ、それは日本の物理学会が発行するJournal of the Physical Society of Japan、いわゆる JPSJのImpact Factor 、1.7を超えています・・・！」 博士「そうじゃな、つまり論文誌としての影響力としてはJPSJよりも大きな雑誌に掲載された論文ということじゃ」 学生「JPSJもっとがんばってくださいよ」 博士「しかもこの論文は8回も引用されておる」 学生「 みなしPRL じゃないですか・・・」 博士「これだけで如何にこの論文が本物であるか納得してもらえ
続きを読む

まだ最小二乗法で消耗してるの？ ～Passing Bablok法について～

イメージ
【イントロ】 与えられたデータの傾向を見る、つまり回帰分析の一つとして、直線でフィッティングする方法が挙げられます。 ココで問題となるのが、「どうやって直線の傾き(と切片)を求めるか」という点です。 代表的な方法が最小二乗法です。目的の直線を $y=ax$ と仮定します。データ点$y$と求めたい直線との差の二乗を誤差として、すべての点に対する誤差の和を最小化する傾き$a$を求める方法です。 具体的には誤差の和を $E=\sum_{i=1}^{n}(y_i－ax_i)^2$ とすると、$a$について最小値が欲しいので、 $\frac{\partial E}{\partial a}=0$ となる$a$を求めれば良いです。整理すると、 $a=\frac{\sum_{i=1}^{n}(x_i－x_{average})(y_i－y_{average})}{\sum_{i=1}^{n}(x-x_{average})^2}$ となります。 　最小二乗法は直観的にもちょうどよい直線を見つけてくれそうで、最強の方法に見えます。しかし、 いくつかの課題 があります。まず一つが、外れ値に弱いという点です。えっちですね。つまり、ノイズ等の影響でデータ点の集合から大きくハズレた点が存在すると、その点に引きずられて傾き$a$が大きく変化してしまうという点です。 　また、別の問題として、最小二乗法は暗にデータが正規分布に従った誤差を持つことを仮定している点が挙げられます。今回最小化すべき誤差をデータ点と直線の差の二乗和として与えましたが、これはデータ点の周りに分散一定の誤差が正規分布している場合に導かれる最尤推定解と対応しています。詳細は PRML を参照。 さて、それではこうした最小二乗法の呪縛から解き放たれた、 回帰モデルは存在するでしょうか？そうしたモデルとしては 主成分分析 等があります 本記事では、そうした手法の中でもノンパラメトリック回帰モデルの一つである”Passing-Bablok法”について、 原論文 に従った解説とPythonでの実装に取り組みました。 【方法】 Passing-Bablok法はその名の通り、PassingとBablok(PB)により提案された、一次の比例関係を比較するノンパラメトリック回
続きを読む