多軌道系のRPA感受率を計算しよう

\def\k{\mathbf{k}}
\newcommand{\kq}{\mathbf{k}+\mathbf{q}}

のもう少し複雑な例をやってみます。

１体模型

で紹介されていたRPA計算を再現してみます。

まず、相互作用なしのハミルトニアンは消滅演算子を

\mathbf{c}(\k) = (c_{d_{zx}\uparrow},c_{d_{zx}\downarrow},c_{d_{yz}\uparrow},c_{d_{yz}\downarrow})^T

とすれば、

H_0(\k) = \mathbf{c}^\dagger(\k) \left[ \begin{pmatrix}
-2t\cos{k_x} & 4t'\sin{k_x}\sin{k_y} \\
4t'\sin{k_x}\sin{k_y} & -2t\cos{k_y}
\end{pmatrix}\otimes\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} \right] c(\k)

です。波数は

k_x,k_y = \frac{2m\pi}{L} \ (m = 0,1,\cdots,L-1)

です。

多軌道模型のクーロン相互作用

原子軌道が複数あるとき、同一原子上の電子管にはどの原子軌道をとるかに応じて種類の電子管相互作用が生じます。副格子、原子軌道、スピンの自由度を$r,l,s$と表すことにすれば、異方性がなければ

\begin{align}
H_{int} &= \frac{U}{2} \sum_{r,l,\sigma} c^\dagger_{r,l,\sigma} c^\dagger_{r,l,-\sigma} c_{r,l,-\sigma} c_{r,l,\sigma} \\

&+ \frac{U'}{2} \sum_{r}\sum_{l_1 \neq l_2}\sum_{\sigma_1,\sigma_2} c^\dagger_{r,l_1,\sigma_1} c^\dagger_{r,l_2,\sigma_2} c_{r,l_2,\sigma_2} c_{r,l_2,\sigma_1} \\

&+ \frac{J}{2} \sum_{r}\sum_{l_1 \neq l_2}\sum_{\sigma_1,\sigma_2} c^\dagger_{r,l_1,\sigma_1} c^\dagger_{r,l_2,\sigma_2} c_{r,l_1,\sigma_2} c_{r,l_2,\sigma_1} \\

&+ \frac{J'}{2} \sum_{r}\sum_{l_1 \neq l_2}\sum_{\sigma} c^\dagger_{r,l_1,\sigma} c^\dagger_{r,l_1,-\sigma} c_{r,l_2,-\sigma} c_{r,l_2,\sigma}
\end{align}

となります。必要に応じて異なるサイトの電子間の相互作用

H'_{int} = \frac{1}{2} \sum_{\langle r_1,r_2 \rangle} V(r_1 - r_2) \sum_{l_1,l_2,\sigma_1,\sigma_2} c^\dagger_{r_1,l_1,\sigma_1} c^\dagger_{r_2,l_2,\sigma_2} c_{r_2,l_2,\sigma_2} c_{r_1,l_1,\sigma_2}

（$\langle r_1,r_2 \rangle$は原子の結合があるサイトのペア）も考える場合もあります。

\begin{align}
J &= 0.1U \\
U' &= U - 2J \\
J' &= J
\end{align}

とするのが一般的なようです。

相関関数の表示

$(i_1,i_2)$を辞書順に並べて１つの添え字とみなしします。すると0次の相関関数は

で計算したように

\begin{align}
\Phi^0_{(i_1,i_2)(i_3,i_4)}(\mathbf{q},i\omega_l)
= \frac{1}{L^d}\sum_\k\sum_{j_1,j_2} 
& U^*_{i_1j_1}(\k)U_{i_3j_1}(\k)U^*_{i_4j_2}(\mathbf{k+q})U_{i_2j_2}(\mathbf{k+q}) \\
&\times \frac{ f[E_{j_2}(\mathbf{k+q})] - f[E_{j_1}(\k)] }{ E_{j_1}(\k) - E_{j_2}(\mathbf{k+q})  + i\omega_l }
\end{align}

です。

ここで

\mathscr{U}_{(i_1,i_2),(j_1,j_2)}(\k,\mathbf{q}) = U^*_{(i_1,j_1)}(\k) U_{(i_2,j_2)}(\kq)

と置くと、

\mathscr{U}^\dagger_{(j_1,j_2),(i_3,i_4)}(\k,\mathbf{q}) = U_{(i_3,j_1)}(\k) U^*_{(i_4,j_2)}(\kq)

だから、

\mathcal{F}(\k,\mathbf{q},\omega_l) =  Diag\left[ \left\{ \frac{ f[E_{j_2}(\mathbf{k+q})] - f[E_{j_1}(\k)] }{ E_{j_1}(\k) - E_{j_2}(\mathbf{k+q})  + i\omega_l } \right\}_{(j_1,j_2)} \right]

とすると、

\begin{align}
\Phi^0_{(i_1,i_2)(i_3,i_4)}(\mathbf{q},i\omega_l)
&= \frac{1}{L^d}\sum_\k\sum_{j_1,j_2,j_3,j_4} 
\mathscr{U}_{(i_1,i_2),(j_1,j_2)}(\k,\mathbf{q})
\mathcal{F}_{(j_1,j_2)(j_3,j_4)}(\k,\mathbf{q},\omega_l) 
\mathscr{U}^\dagger_{(j_3,j_4),(i_3,i_4)}(\k,\mathbf{q}) \\
&= \frac{1}{L^d} \sum_\k \mathscr{U}(\k,\mathbf{q}) \mathcal{F}(\k,\mathbf{q},\omega_l) \mathscr{U}^\dagger(\k,\mathbf{q})
\end{align}

と書けます。

そしてRPA感受率は

によると

\begin{align}
\Phi^{RPA}(\mathbf{q},i\omega_l) &= [1 + \Phi^0(\mathbf{q},i\omega_l)V]^{-1}\Phi^0(\mathbf{q},i\omega_l)\\
\chi_{BA} &= \sum_{i_1,\cdots,i_4} B_{i_1i_2} A^*_{i_3,i_4} \Phi^{RPA}_{(i_1,i_2)(i_3,i_4)}(\mathbf{q},i\omega_l)
\end{align}

です。$V$は上記の記事の手順で相互作用の行列からつくった行列です。
また物理量を表す行列を

\vec{A}_{(i,j)} = A_{ij}

とベクトルとして表すと

\chi_{BA}(\mathbf{q},\omega+i\delta) = \vec{B} \Phi^{RPA}(\mathbf{q},\omega+i\delta) \left( \vec{A} \right)^*

となります。

化学ポテンシャルのチューニング

温度$T$での１サイト当たりの粒子数$\nu$は化学ポテンシャル$\mu$に対して

\nu(\mu;T) = \frac{1}{L^d}\sum_{\k}\sum_{n} f(E_n(\k);T)

で計算できます。

よって

\begin{align}
\mu_{low} &\leq E_{\min} \\
\mu_{sup} &\geq E_{\max} 
\end{align}

から初めて

\begin{align}
\nu_{sup} &= \nu(\mu_{sup};T) \\
\nu_{low} &= \nu(\mu_{low};T) \\
\nu_{mid} &= \frac{\nu_{sup} + \nu_{low}}{2} \\

&\begin{cases}
\mu_{sup} \gets \mu_{mid} \ ( \nu_{mid} \geq \nu_{taeget}  ) \\
\mu_{low} \gets \mu_{mid} \ ( \nu_{mid} < \nu_{taeget}  )
\end{cases}
\end{align} 

と更新して二分探索で化学ポテンシャルを決めます。

実装

RPA感受率を計算をするコンポーネントは以下のようになります。

相互作用なしの相関関数の計算

を元にMPI並列化をしています。

bare.py

import numpy as np
import numpy.linalg as LA
import itertools
from mpi4py import MPI

comm = MPI.COMM_WORLD
procs = comm.Get_size()
rank = comm.Get_rank()

# MPI並列化での配列の分割
class assign:
    
    def __init__(self,procs,size):
        self.size = size # 配列のサイズ
        self.procs = procs # 使用コア数
        self.num = [ size//procs + (1 if i<size%procs else 0) for i in range(procs) ]
        self.id = [0] + list(itertools.accumulate( self.num))

# フェルミ分布関数
class Fermi(object):
    
    def __init__(self,chemicalPotential):
        self.chem = chemicalPotential
    
    def priventInf(self,x):
        if x > 709: #prevent overflow
            return 709
        elif x < -745: #prevent underflow
            return -745
        return x

    def boltzmann_factor(self,energy, kbt):
        inv_temp = 1/kbt
        a = -self.priventInf(inv_temp*(energy))
        return np.exp( a )
    
    def fermi(self,kbt,energy):
        E2 = energy - self.chem
        return 1/(self.boltzmann_factor(-E2, kbt)+1)
    
    def diffFermi(self,kbt,energy):
        E2 = energy - self.chem
        bz = self.boltzmann_factor(-E2, kbt) if E2>0 else self.boltzmann_factor(E2, kbt)
        return -bz/(kbt*(1+bz)**2)

# 化学ポテンシャルの調整
class FermiLevelTune(Fermi):
    
    def __init__(self,size,dim,hamiltonianFunc,threshold):
        super().__init__(0)
        self.threshold = threshold # 収束判定の閾値
        self.hamiltonian = hamiltonianFunc # ハミルトニアンの関数; 引数はkのみになるようにして下さい
        self.size = size # メッシュサイズ
        self.volume = size**dim # 総サイト数　dimは次元
        mesh = [ np.array(k_tuple)/size for k_tuple in itertools.product(list(range(size)),repeat=dim) ]
        self.energy_arr = np.sort(np.array([LA.eigvalsh(hamiltonianFunc(k)).tolist() for k in mesh]).flatten())[::-1]
        self.Emax = np.max(self.energy_arr)
        self.Emin = np.min(self.energy_arr)
        self.Elim = max(abs(self.Emax),abs(self.Emin))
                    
    def NumParticle(self,chem,kbt):
        self.chem = chem
        return np.sum(np.vectorize(lambda E: self.fermi(kbt,E))(self.energy_arr))/self.volume
    
    def tune(self,filling,kbt):
        Esup = 2*self.Elim
        Elow = -2*self.Elim
        while(Esup - Elow > self.threshold):
            Emid = 0.5*(Esup+Elow)
            fillTmp = self.NumParticle(Emid,kbt)
            if(fillTmp<filling): Elow = Emid
            else: Esup = Emid
        return 0.5*(Esup+Elow)

# 相関関数行列の添え字
class IDcorr(object):
    def __init__(self,Nr,Nl):
        # Nr: 副格子の数, Nl: 1原子当たりの原子軌道の数
        self.Nr = Nr; self.Nl = Nl; self.Ns = 2
        self.N = Nr*Nl*self.Ns
        self.N2 = self.N*self.N

    def IdBand(self,r,l,s):
        return r*self.Nl*self.Ns + l*self.Ns + s
    
    def Id2state(self,i1,i2):
        return i1*self.N + i2
    
    def IdInt(self,r1,l1,s1,r2,l2,s2):
        return  self.Id2state(self.IdBand(r1,l1,s1),self.IdBand(r2,l2,s2))

# 相関関数(相互作用なし)の計算
class BareSusceptability(IDcorr,Fermi,assign):
    
    def __init__(self,procs,rank,
                 chemicalPotential,delta,
                 size,dim,Nr,Nl,hamiltonianFunc):
        IDcorr.__init__(self,Nr,Nl)
        Fermi.__init__(self,chemicalPotential)
        self.size = size; self.dim = dim; self.volume = size**dim
        self.delta = delta
        self.hamiltonian = hamiltonianFunc
        self.k_list = [ np.array([kx,ky])/self.size for kx,ky 
             in itertools.product(list(range(self.size)),repeat=self.dim)]
        assign.__init__(self,procs,len(self.k_list))
        self.k_listPerCore = self.k_list[self.id[rank]:self.id[rank+1]]

    def LindhardElement(self,Ek,Ekq,kbt,omega):
        if(abs(omega)<1e-6 and abs(Ek-Ekq)<1e-3):
            if(Ekq>=Ek):
                return -Fermi.diffFermi(self,kbt,Ekq)
            if(Ekq<Ek):
                return Fermi.diffFermi(self,kbt,Ek)
        numerator = Fermi.fermi(self,kbt,Ekq) - Fermi.fermi(self,kbt,Ek)
        denominator = ( Ekq - Ek + omega - 1j*self.delta )
        denominator /= ( ( Ekq - Ek + omega )**2 + self.delta**2 )
        return -numerator*denominator 
        
    def Lindhard(self,k,q,kbt,omega):
        Ek,Uk = LA.eigh(self.hamiltonian(k)); Ekq,Ukq = LA.eigh(self.hamiltonian(k+q))
        UU = np.kron(np.conj(Uk),Ukq)
        Fvec = np.zeros(self.N2).astype("complex")
        for i,j in itertools.product(list(range(self.N)),repeat=2):
            Fvec[IDcorr.Id2state(self,i,j)] = self.LindhardElement(Ek[i],Ekq[j],kbt,omega)
        return UU.dot(np.diag(Fvec)).dot(np.conj(UU.T))
    
    def CorrPerCORE(self,q,kbt,omega):
        return np.array([ self.Lindhard(k,q,kbt,omega)
                for k in self.k_listPerCore 
            ]).sum(axis=0)/self.volume
        
    def Corr(self,q_list,kbt,omega):
        sendbuf = np.array([ self.CorrPerCORE(q,kbt,omega) for q in q_list ])
        Corr_arr = np.empty(sendbuf.shape, dtype="complex128")
        comm.Allreduce(sendbuf, Corr_arr, MPI.SUM)
        return Corr_arr

### 実行部分

# ハミルトニアン
def hamiltonian(k):
    t1 = 1.; t2 = 0.1*t1; p = 2*np.pi
    return np.kron(np.array([
        [-2*t1*np.cos(p*k[0]), 4*t2*np.sin(p*k[0])*np.sin(p*k[1])],
        [4*t2*np.sin(p*k[0])*np.sin(p*k[1]), -2*t1*np.cos(p*k[1])]
    ]),np.identity(2))

# 感受率の波数点
k_list = [ 0.5*np.array([k,0]) for k in np.linspace(0,1,31)[:-1]] + [
    0.5*np.array([1,k]) for k in np.linspace(0,1,31)[:-1]
    ] + [ 0.5*np.array([k,k]) for k in np.linspace(1,0,31)]

# パラメーター
size = 2**11
kbt = 0.02
filling = 8/3
delta = 1000/(size**2)
SaveFile = "CorrBareData"

# 実行部分
chem = FermiLevelTune(2**7,2,hamiltonian,1e-5).tune(filling,kbt)
bare = BareSusceptability(procs,rank,chem,delta,size,2,1,2,hamiltonian)
Corr_bare_arr = bare.Corr(k_list,kbt,0)
# 計算結果を出力
if(rank == 0): np.save(SaveFile,Corr_bare_arr)

RPA感受率への変換

先のコードで計算した相互作用なしの相関関数のデータCorrBareData.npyを読み込んで、これをRPA感受率に変換します。

bareToRPA.py

import numpy as np
import numpy.linalg as LA
import itertools
import pandas as pd

Corr_bare_arr = np.load("CorrBareData.npy")

class IDcorr(object):
    def __init__(self,Nr,Nl):
        self.Nr = Nr; self.Nl = Nl; self.Ns = 2
        self.N = Nr*Nl*self.Ns
        self.N2 = self.N**2

    def IdBand(self,r,l,s):
        return r*self.Nl*self.Ns + l*self.Ns + s
    
    def Id2state(self,i1,i2):
        return i1*self.N + i2
    
    def IdInt(self,r1,l1,s1,r2,l2,s2):
        return  self.Id2state(self.IdBand(r1,l1,s1),self.IdBand(r2,l2,s2))

# RPAの計算で相互作用の効果を与える行列の構成
class INTt2g(IDcorr):
    
    def __init__(self,Nr,Nl):
        super().__init__(Nr,Nl)
    
    def interactionT2g(self,U,U2,J,J2):
        MatInt = np.zeros((self.N2,self.N2)).astype("complex")
        Nl_list = list(range(self.Nl)); Ns_list = list(range(self.Ns))
        for r in range(self.Nr):
            for l,s1 in itertools.product(Nl_list,Ns_list):
                s2 = (s1+1)%2
                MatInt[self.IdInt(r,l,s1,r,l,s2),
                       self.IdInt(r,l,s2,r,l,s1)] = 0.5*U
            for l1,l2 in itertools.permutations(Nl_list,2):
                for s1,s2 in itertools.product(Ns_list,repeat=2):
                    MatInt[self.IdInt(r,l1,s1,r,l2,s2),
                           self.IdInt(r,l2,s2,r,l1,s1)] = 0.5*U2
                    MatInt[self.IdInt(r,l1,s1,r,l2,s2),
                           self.IdInt(r,l1,s2,r,l2,s1)] = 0.5*J
                for s1 in Ns_list:
                    s2 = (s1+1)%2
                    MatInt[self.IdInt(r,l1,s1,r,l1,s2),
                           self.IdInt(r,l2,s2,r,l2,s1)] = 0.5*J2
        return MatInt
    
    def IntMatRPA(self,U,U2,J,J2):
        MatRPA = np.zeros((self.N2,self.N2)).astype("complex")
        MatInt = self.interactionT2g(U,U2,J,J2)
        for i1,i2,i3,i4 in itertools.product(list(range(self.N)),repeat=4):
            i = self.Id2state(i1,i2); j = self.Id2state(i3,i4)
            MatRPA[i,j] += MatInt[self.Id2state(i2,i3),self.Id2state(i4,i1)]
            MatRPA[i,j] += MatInt[self.Id2state(i3,i2),self.Id2state(i1,i4)]
            MatRPA[i,j] -= MatInt[self.Id2state(i2,i3),self.Id2state(i1,i4)]
            MatRPA[i,j] -= MatInt[self.Id2state(i3,i2),self.Id2state(i4,i1)]
        return MatRPA

# bareの相関関数をRPA相関関数に変換
class RPA(INTt2g):
    
    def __init__(self,U,Nr,Nl):
        super().__init__(Nr,Nl)
        J = 0.1*U
        U2 = U - 2*J; J2 = J
        self.MatRPA = self.IntMatRPA(U,U2,J,J2)
    
    def BareToRPA(self,Corr):
        return LA.inv((np.identity(self.N2) + Corr.dot(self.MatRPA))).dot(Corr)
    
    def Corr(self,Corr_arr):
        return np.array([ self.BareToRPA(Corr) for Corr in Corr_arr ])

# 感受率の計算
def susceptibility(A,Corr):
    return A.flatten().dot([e-mail](A.flatten()))

def susceptibilityArr(A,CorrArr):
    return np.array([susceptibility(A,C) for C in CorrArr])

# グラフ描画用の波数点
k_list = [ 0.5*np.array([k,0]) for k in np.linspace(0,1,31)[:-1]] + [
    0.5*np.array([1,k]) for k in np.linspace(0,1,31)[:-1]
    ] + [ 0.5*np.array([k,k]) for k in np.linspace(1,0,31)]

klabel_arr = np.array([0] 
                    + list(itertools.accumulate([ LA.norm(k_list[i] - k_list[i-1]) 
                                                for i in range(1,len(k_list)) ])))
klabel_arr /= klabel_arr[-1]

# 感受率を計算する物理量
s0 = np.identity(2).astype("complex"); s1 = np.array([[0,1],[1,0]]).astype("complex")
s2 = np.array([[0,-1j],[1j,0]]); s3 = np.array([[1,0],[0,-1]]).astype("complex")

nc = np.identity(4).astype("complex")
sz = np.kron(s0,s3)
o2 = np.kron(s3,s0)
oxy = np.kron(s1,s0)
lz = np.kron(s2,s0)
tz = np.kron(s3,s3)
txyz = np.kron(s1,s3)

# 実行部分
for U in np.round(np.linspace(0,2,6),1):
    Corr_rap_arr = RPA(U,1,2).Corr(Corr_bare_arr)
    outputFile = "rpa_U"+str(U)+".csv"

    df_rpa = pd.DataFrame({
        "k":klabel_arr,
        "n":susceptibilityArr(nc,Corr_rap_arr).real,
        "sz":susceptibilityArr(sz,Corr_rap_arr).real,
        "o2":susceptibilityArr(o2,Corr_rap_arr).real ,
        "oxy":susceptibilityArr(oxy,Corr_rap_arr).real ,
        "lz":susceptibilityArr(lz,Corr_rap_arr).real ,
        "tz":susceptibilityArr(tz,Corr_rap_arr).real ,
        "txyz":susceptibilityArr(txyz,Corr_rap_arr).real
        }).to_csv(outputFile,index=False)

これによって相互作用の大きさを変えた感受率のデータファイルrpa_Uxxx.csv(xxx=0,0.4,$\cdots$,2.0)が生成されます。

実行結果

例えば、次のように計算結果を出力したファイルを読み込んで図示します。定義の問題で感受率の値が参考文献の2倍になっているので、$0.5$をかけて調整してあります。

plot.py

import pandas as pd

U = 0.0
df = pd.read_csv("rpa_U"+str(U)+".csv")
df.iloc[:,1:] *=0.5

df.plot("k",["sz","o2","oxy","lz"],xlim=[0,1],ylim=[0,1])

パラメーター

$U$を変化させて計算します。他のパラメーターは以下のもので固定しました。

\begin{align}
t &= 1. \\
t' &= 0.1t \\
L &= 2^{11} \\
k_BT &= 0.02t \\
\nu &= \frac{8}{3} \\
\delta &= \frac{1000}{L^2}
\end{align}

相互作用なし(U=0の場合)

数値誤差がまだ残っていますが、参考文献の計算結果と概ね合っています。

相互作用あり

コメント

• 数値誤差が出やすくて参考資料の計算を再現するにはかなり大きなシステムサイズが必要でした。特に$\mathrm{\Gamma,X}$点の周りでその影響が大きかったです。デモ計算の際はクラスタコンピューターで100並列して数時間かかりました。(手元のPCでできる範囲でもある程度の結果は得られます。)デモ計算でも数値誤差が残っていたのでもっと大きいシステムサイズで計算する必要があるかもしれません。

• 今回は分かりやすさ重視でpythonで実装しました。しかし、３次元系や複雑なハミルトニアンを用いて計算したいときはc++やjuliaで実装する必要があると思います。

• $\delta$の値は$E_n(k)$の値の間隔程度でとるとよいとされています。体積(サイト数)$V=L^2$に対して$V\delta = const.$となるようにとるのが通例となっています。

• 定義の問題で僕のコードでは感受率の値が参考文献の2倍になっています(グラフに描画する前に調整しています)。