密度行列繰り込み群法(DMRG)で1次元多体系の基底状態を求めよう

エレクトロ二クスにつかう材料として磁性体(=磁石)はとても重要な材料として注目されてきた。磁石などの性質は一般に相互作用によって発現するが、相互作用のあるシステム(多体系)は計算量が膨大で現実的な規模での計算には近似が必要となる。
一番簡単なのは平均場近似である：

しかし、この近似精度ではとらえられない現象もあり、そんなときに使う方法の１つが密度行列繰り込み群法である。これは１次元のモデルにしか使えないという欠点もあるが、化学系ではポリマーのシミュレーションに使われ、また最近話題のテンソルネットワークの基礎になっている手法であるという点でも重要だと思い、ここで紹介したい。

モデル

横磁場イジングモデルを例に考える($J>0$)。

H = -J \sum_i \sigma^z_i \sigma^z_{i+1} -h\sum_i \sigma^x_i 

アイデア

適当な$j$に対して

H_L = -J\sum_{i=0}^{j-1} \sigma^z_i \sigma^z_{i+1} -h\sum_{i=0}^{j-1} \sigma^x_i\\
H_R = -J\sum_{i=j+1}^{N} \sigma^z_i \sigma^z_{i+1} -h\sum_{i=j+1}^{N} \sigma^x_i\\
H_C = -J \sigma^z_j \sigma^z_{j+1} \\
H = H_L + H_R + H_C

として$H_L,H_R$の固有状態$(E^L_i, |L_i\rangle), (E^R_i, |R_i\rangle) $を考える(どうやって求めたかはひとまず置いておく)。
このとき、

\langle L_{i_l} R_{i_r} | H | L_{j_l} R_{j_r} \rangle = E^L_{i_l} \delta_{i_l,j_l} + E^R_{i_r} \delta_{i_r,j_r} + \langle L_{i_l} R_{i_r} | H_C | L_{j_l} R_{j_r} \rangle

を(i_l,i_r)について辞書順に並べて成分とした行列が$H$となる。それを対角化すれば全体の固有状態が分かる。
ただ、$| L_{i_l} R_{i_r} \rangle$を全部考えると行列のサイズが膨大になってしまうので、$| L_{i_l} \rangle, |R_{i_r} \rangle$それぞれのエネルギーが低い$m$個だけをとってきて行列を作れば、$H$の行列サイズは$m^2 \times m^2$まで小さくできる。$Mが$十分小さければ、$H$の対角化の計算量はそこまで大きくならない。

これで長さ$N$の近似固有状態が分かったので、これを$H_L,H_R$としてもう一度今のプロセスを行えば、長さ$2N$のハミルトニアンの近似固有状態がわかる。$N=1$(単一サイト)から初めてこのプロセスを繰り返せば長さ$N=2^n$の固有状態を求めることができる。

簡単な例

簡単のために$h=0$(ただのイジングモデル)で考える。$M=2$とする。

n=0

$H=0$なので$|\uparrow\rangle,|\downarrow\rangle$をとればよい。

n=1

$|LR\rangle = ( |\uparrow\uparrow\rangle, |\uparrow\downarrow\rangle, |\downarrow\uparrow\rangle, |\downarrow\downarrow\rangle )$となるので、

H = \begin{pmatrix}
 -J & & & \\
& J & & \\
& & J & \\
& & & -J
\end{pmatrix}

となる。そこで$E=-J$である$|\uparrow\uparrow\rangle, |\downarrow \downarrow\rangle$をとる。

n=2

$|LR\rangle = ( |\uparrow\uparrow\uparrow\uparrow\rangle, |\uparrow\uparrow\downarrow\downarrow\rangle, |\downarrow\downarrow\uparrow\uparrow\rangle, |\downarrow \downarrow\downarrow\downarrow\rangle)$だから、

H = -2J + \begin{pmatrix}
 -J & & & \\
& J & & \\
& & J & \\
& & & -J
\end{pmatrix}

となって$E=-3J$の$|\uparrow\uparrow\uparrow\uparrow\rangle,|\downarrow \downarrow\downarrow\downarrow\rangle$をとる。

あるn

帰納的に長さ$2^{n-1}$で$|\uparrow\cdots\uparrow\rangle,|\downarrow \cdots\downarrow\rangle$が$E=-(2^{n-1}-1)J$で最安定となる。

$|LR\rangle = ( |\uparrow\cdots\uparrow\uparrow\cdots\uparrow\rangle, |\uparrow\cdots\uparrow\downarrow\cdots\downarrow\rangle, |\downarrow\cdots\downarrow\uparrow\cdots\uparrow\rangle, |\downarrow\cdots\downarrow\downarrow\cdots\downarrow\rangle)$だから、

H = -(2^n-2)J + \begin{pmatrix}
 -J & & & \\
& J & & \\
& & J & \\
& & & -J
\end{pmatrix}

となって$E=-(2^n-1)J$の$|\uparrow\cdots\uparrow\uparrow\cdots\uparrow\rangle,|\downarrow\cdots \downarrow\downarrow\cdots\downarrow\rangle$をとる。

確かに強磁性状態が最安定となることが確かめられる。

密度行列

先は波動関数$|\psi\rangle$だけで考えていたが、今度は密度行列$\rho$を利用してハミルトニアンの行列サイズ圧縮を行う。
密度行列はいくつかの波動関数$|\psi_1\rangle,\cdot,|\psi_n\rangle$への射影演算子を適当な割合で重み付き平均したものであり、

\rho = \sum_i p_i |\psi_i \rangle\langle \psi_i | \\
\sum_i p_i = 1 \quad (p_i \geq 0)

で表されて、波動関数の拡張概念になっている。
特別な場合としてある$i$について$p_i=1$のときは

\rho = |\psi_i \rangle\langle \psi_i |

となって、波動関数と1-1対応するようになる。このような$\rho$を純粋状体の密度行列、そうでない場合は混合状態の密度行列という。
例えば有限温度の系に対するボルツマン分布

\rho(T) = \frac{1}{Z}\sum_i |i \rangle\langle i | \exp\left(-\frac{E_i}{k_B T} \right)

は混合状態の密度演算子である。

また、物理量$A$の平均値は

\langle A \rangle = Tr[A\rho]

で計算される。

シュミット分解

まず,
波動関数について左右のブロック全体のシステムに対する波動関数

|\Psi\rangle = \sum_{ij} \Psi_{ij} |i \rangle_{L} \otimes |j \rangle_{R}

を考える。

行列$\Psi$を特異値分解する：$\Psi = U^T \sqrt{W} V$。
ここで$ \sqrt{W} = Diag[w_1 \cdots w_{N_{schmidt}}, 0,0,\cdots ]$, $w_1 \geq w_2\geq \cdots \geq w_{N_{schmidt}}$

すると、

|w_\mu \rangle_L = \sum_i U_{\mu i} |i \rangle_{L}\\ 
|w_\mu \rangle_R = \sum_j V_{\mu j} |j \rangle_{R}

|\Psi \rangle = \sum_{ij} \sum_{\mu} U^T_{i\mu} \sqrt{w_\mu} V_{\mu j} |i \rangle_{L} \otimes |j \rangle_{R} \\
= \sum_{\mu} \sqrt{w_\mu} |w_\mu \rangle_L \otimes |w_\mu \rangle_R

となる。これに対応する密度行列$\rho = |\Psi \rangle\langle \Psi |$で左右について部分トレースをとり、

\rho_L = Tr_R[ |\Psi \rangle\langle \Psi | ] = \sum_{\mu} w_\mu |w_\mu \rangle_L \langle w_\mu|_L\\
\rho_R = Tr_L[ |\Psi \rangle\langle \Psi | ] = \sum_{\mu} w_\mu |w_\mu \rangle_R \langle w_\mu|_R

として、左右のブロックに対する密度行列を取得できる。
ちなみに$w_\mu$は

\rho = \sum_{iji'j'} \Psi_{ij}\Psi^*_{i'j'} |i \rangle_{L} |j \rangle_{R} \langle i' |_{L} \langle j |_{R}

の係数行列$(\Psi_{ij}\Psi^*_{i'j'})_{(i,j),(i'j')}$の対角化でも得られる。

DMRG

図のようなシステムを考える。このシステムでは長さ$l$の左右のブロックの間で２つのサイトが結合している。このシステムのハミルトニアンを以下のように表す。

H^{2l+2} = H^l_L + H_{L\circ} + H_{\circ\circ} + H_{\circ R} + H^l_R

$\circ$はブロックの間に挟んだ２つのサイトを表しており、全体として左から$l+1,l+2$番目のサイトを表している。

左のブロックのエネルギーの低い$M^L$個の固有状態$|m_L \rangle_L$、右のブロックのエネルギーの低い$M^R$個の固有状態$|m_R \rangle_R$をとり、単一サイトの固有状態は$|\sigma \rangle$とする。
このシステム全体の基底状態を最もよく近似する波動関数が

|\Psi\rangle = \sum_{m_L,\sigma_L}\sum_{m_R,\sigma_R} \Psi_{(m_L,\sigma_L),(m_R,\sigma_R)} |m_L \rangle_{L}|\sigma_L \rangle \otimes |m_R \rangle_{R}|\sigma_R \rangle \\
= \sum_{ij} \Psi_{ij} |i \rangle_{L} \otimes |j \rangle_{R}

であるとする。ここで$i=(m,\sigma)$である。これは$H^{(2l+2)}$を対角化して、最低エネルギー状態を求めれば分かる。
そして右側について部分トレースをとった

\rho^L=Tr_R[|\Psi \rangle \langle \Psi| ] \\
\Leftrightarrow \rho^L_{ii'}= \sum_j \Psi_{ij}\Psi_{i'j}

が左側のブロックを１サイト分長くしたときの近似基底状態の密度行列になる。これを対角化して大きい方から$M^L$個の固有状態をとる：

w^L_1 \geq w^L_2 \geq \cdots \geq w^L_{M^L}

そしてこれを並べた変換行列

T_L = \begin{pmatrix}
_L\langle 1| w^L_1 \rangle & \cdots & \langle 1_L| w^L_{M^L} \rangle \\
\vdots & \ddots & \vdots \\
_L\langle N^L| w^L_1 \rangle & \cdots & _L\langle (N^L)| w^L_{M^L} \rangle
\end{pmatrix}

を用いて、ハミルトニアンを

H^{l+1}_L = T_L^\dagger (H^l_L + H_{L\circ}) T_L

と、サイズ$M^L \times M^L$に圧縮したものに置き換える。

同様にして

\rho^R=Tr_L[|\Psi \rangle \langle \Psi| ] \\
\Leftrightarrow \rho^R_{jj'}= \sum_i \Psi_{ij}\Psi_{ij'}

を対角化して

w^R_1 \geq w^R_2 \geq \cdots \geq w^R_{M^R}

をとり、

H^{l+1}_R = T_R^\dagger (H_{\circ R} + H^l_R) T_R

とする。これを繰り返してより長いシステムの波動関数を計算していく。

ハミルトニアンの表し方

システムサイズを大きくしていく際に$H_{L\circ},H_{\circ\circ},H_{\circ R}$をどうすればよいか説明しよう。
ここでは簡単のために外場項を落として説明する。
更新の際には

• A:$|[\sigma]_l \rangle = | \sigma_1 \sigma_2 \cdots \sigma_l \rangle$を基底にして表したオブジェクト
• B:$| w_i^{(l)} \rangle$を基底にしたオブジェクト
• C:$| w_i^{(l)} \sigma \rangle = | w_i^{(l)} \rangle \otimes |\sigma \rangle$を基底にしたオブジェクト

の３通りの表し方を互いに変換しながら表す必要がある(右側のブロックも同様)。

基底A,Bの変換は行列

V_{L/R}^{(l)} = ( \langle w^{L/R}_{m^{(l)}} | \sigma_1 \sigma_2 \cdots \sigma_l \rangle )

で行う。また、基底B,Cの変換は前節の$T_{L/R}$で行う。

HLoの表し方

$H_L$はサイズ圧縮のために$|w_i^L\rangle$で表されている。

追加する部分に必要な右端のスピンは基底$|[\sigma]_l \rangle$を用いれば、簡単に表すことができて

s_z^{(l)} = |[\sigma]_l \rangle \langle [\sigma]_l | s_z^{(l)} |[\sigma]_l \rangle \langle [\sigma]_l | \\
= |[\sigma]_{l-1} \rangle \langle [\sigma]_{l-1} | \otimes |\sigma_l \rangle \langle \sigma_l | s_z^{(l)} |\sigma_l \rangle  \langle \sigma_l |

となる。ここで

|[\sigma]_{l-1} \rangle \langle [\sigma]_{l-1} | = I_{2^{l-1} \times 2^{l-1}} \quad (単位行列) \\

|\sigma_l \rangle \langle \sigma_l | s_z^{(l)} |\sigma_l \rangle  \langle \sigma_l | =\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix} = \sigma_z \\
s_z^{(l)} = I_{2^{l-1} \times 2^{l-1}} \otimes \sigma_z

である。よって規定を取り換えると

s_z^{(l)} = | w_i^{(l)} \rangle \langle w_i^{(l)} | [\sigma]_l \rangle \langle [\sigma]_l | s_z^{(l)} |[\sigma]_l \rangle \langle [\sigma]_l | w_i^{(l)} \rangle \langle w_i^{(l)} | \\

であり、これを行列で表すと

\Leftrightarrow s_z^{(l)}|_L = V^{(l)}_L ( I_{2^{l-1} \times 2^{l-1}} \otimes \sigma_z ) (V^{(l)}_L)^\dagger

となる。よって右端のスピンを基底$| w_i^{(l)} \rangle$で表せた。

これを使って

\langle w_i^{(l)} \sigma | H^{(l+1)}_L | w_i^{(l)} \sigma \rangle = H^{(l)}_L \otimes I_{2\times 2} -J s_z^{(l)}|_L \otimes \sigma_z \\
\langle w_i^{(l+1)} | H^{(l+1)}_L | w_i^{(l+1)} \rangle = T_L^\dagger [ H^{(l)}_L \otimes I_{2\times 2} -J s_z^{(l)}|_L \otimes \sigma_z ] T_L

また、基底も更新すると

V_L^{(l+1)} =  \langle w^{(l+1)}_i | w^{(l)}_i \sigma \rangle \langle w^{(l)}_i \sigma | \sigma_1 \sigma_2 \cdots \sigma_l \rangle
= T_L^\dagger V_L^{(l)} 

HoRの表し方

同様にして右側のブロックは

 s_z^{(0)}|_R = V^{(l)}_R ( \sigma_z \otimes I_{2^{l-1} \times 2^{l-1}} ) (V^{(l)}_R)^\dagger \\
H^{(l+1)}_R = T_R^\dagger [ I_{2\times 2} \otimes H^{(l)}_R -J \sigma_z \otimes s_z^{(0)}|_R ] T_R \\
V_R^{(l+1)} =T_R^\dagger V_R^{(l)} 

Hの表し方

以上のことから全体のハミルトニアンは $| w_i^{(l)} \sigma \rangle_L \otimes | \sigma w_i^{(l)}\rangle_R$を基底として

H^{2l+2} = [ H^{(l)}_L \otimes I_{2\times 2} -J s_z^{(l)}|_L \otimes \sigma_z ] \otimes I_{2M_R\times 2M_R} + I_{2M_L\times 2M_L} \otimes [ I_{2\times 2} \otimes H^{(l)}_R -J \sigma_z \otimes s_z^{(0)}|_R ] -J I_{M_L\times M_L} \otimes \sigma_z \otimes \sigma_z \otimes I_{M_R\times M_R}

となる。外場項を付け足すと

H^{2l+2} = [ H^{(l)}_L \otimes I_{2\times 2} -J s_z^{(l)}|_L \otimes \sigma_z -h I_{M_L \times M_L } \otimes \sigma_x ] \otimes I_{2M_R \times 2M_R} + I_{2M_L\times 2M_L} \otimes [ I_{2\times 2} \otimes H^{(l)}_R -J \sigma_z \otimes s_z^{(0)}|_R -h \sigma_x \otimes I_{M_R \times M_R } ] -J I_{M_L\times M_L} \otimes \sigma_z \otimes \sigma_z \otimes I_{M_R\times M_R}

となる。

実装

以上のステップをpython3で実装してみたのだが、結果がおかしかった。
一応コードを載せておくが、あくまで参考にとどめておいて使用はしないでほしい。
(ミスを発見したらコメントしていただけると幸いです。)

dmrg.py

import numpy as np
from numpy import linalg as LA

#部分トレース
class pTrace:
    def __init__(self):
        return

    def R(self,M,sizeR):
        N = len(M)
        sizeL = N//sizeR
        Block = np.array([ np.hsplit(m,sizeL) for m in np.vsplit(M,sizeL) ])
        Mp = []
        for i in range(sizeL):
            Mp += [ [ np.trace(m) for m in Block[i] ]  ]
        return np.array(Mp)

    def swap_kron(self,AB,sizeB):
        N = len(AB)
        sizeA = N//sizeB
        Block = np.array([ np.hsplit(m,sizeA) for m in np.vsplit(AB,sizeA) ])
        BA = np.zeros((N,N)).astype("complex")
        for i in range(sizeB):
            for j in range(sizeB):
                I = np.zeros((sizeB,sizeB)); I[i][j]=1.
                BA += np.kron(I,Block[:,:,i,j])
        return BA

    def L(self,M,sizeR):
        N = len(M)
        sizeL = N//sizeR
        return self.R(self.swap_kron(M,sizeR),sizeL)

class DMRG:
    def __init__(self,H0,J,ns):
        self.s0 = np.identity(2).astype("complex")
        self.sz = np.array([[1.,0.],[0.,-1.]]).astype("complex")
        self.J = J
        self.H0 = H0
        self.ns=ns #max size of Hamiltonian

    def I(self,size):
        return np.identity(size).astype("complex")

    def dagger(self,M):
        return np.conj(M).T

    def kron4(self,A,B,C,D):
        AB = np.kron(A,B)
        CD = np.kron(C,D)
        return np.kron(AB,CD)

    def mat_L_end(self,size,M):
        return np.kron(M,self.I(size))

    def mat_R_end(self,size,M):
        return np.kron(self.I(size),M)

    def extend_HL(self,HL,SL):
        HL = self.mat_L_end(2,HL) + self.mat_R_end(len(HL),self.H0)
        HL += -self.J*np.kron(SL,self.sz)
        return HL

    def extend_HR(self,HR,SR):
        HR = self.mat_R_end(2,HR) + self.mat_L_end(len(HR),self.H0)
        HR += -self.J*np.kron(self.sz,SR)
        return HR

    def merge_hamiltonian(self,HL,HR,SL,SR):
        H = self.mat_L_end(len(HR),HL) + self.mat_R_end(len(HL),HR)
        H += -self.J*self.kron4(self.I(len(SL)),self.sz,self.sz,self.I(len(SR)))
        return H

    def update_hamiltonian(self,HL,HR,VL,VR,nL,nR):
        SL = self.mat_R_end(2**(nL-1),self.sz)
        SL = np.conj(VL).dot(SL).dot(VL.T)
        SR = self.mat_L_end(2**(nR-1),self.sz)
        SR = np.conj(VR).dot(SR).dot(VR.T)
        HL = self.extend_HL(HL,SL)
        HR = self.extend_HR(HR,SR)
        H = self.merge_hamiltonian(HL,HR,SL,SR)
        return HL,HR,H

    def update_states(self,density):
        W,P = LA.eigh(density)
        return self.dagger(P)[np.argsort(-W)[0:self.ns]] #id for ns larger W

    def project_op(self,vmin,sizeR):
        trace = pTrace()
        density = np.outer(vmin,np.conj(vmin))
        TL = self.update_states(trace.R(density,sizeR))
        TR = self.update_states(trace.L(density,sizeR))
        return TL,TR

    def DMRG(self,N):
        nL=1; nR=1
        HL=self.H0.copy(); HR = self.H0.copy()
        E1,U1 = LA.eigh(self.H0)
        vmin = U1[np.argmin(E1)]
        ##sets of WF with basis: {|s1s2...sL>} or {|s1s2...sR>}
        TL = U1.copy(); TR = U1.copy()
        VL = U1.copy(); VR = U1.copy()
        while(nL<=N):
            HL,HR,H = self.update_hamiltonian(HL,HR,VL,VR,nL,nR)
            nL+=1; nR+=1
            E,V = LA.eigh(H)
            vmin = V[np.argmin(E)]
            TL,TR = self.project_op(vmin,len(HR))
            HL = TL.dot(HL).dot(self.dagger(TL))
            HR = TR.dot(HR).dot(self.dagger(TR))
            VL = np.kron(VL,self.s0); VR= np.kron(self.s0,VR)
            VLooR = np.kron(VL,VR)
            VL = TL.dot(VL); VR = TR.dot(VR)
        return self.dagger(VLooR)@vmin

J=1.
B=1.
ns=2**5
H0 = -B*np.array([[1.,0.],[0.,-1.]]).astype("complex")
#H0 = -B*np.array([[0.,1.],[1.,0.]]).astype("complex")
dmrg = DMRG(H0,J,ns)

N=2**10
v = dmrg.DMRG(N)

print( np.round(v,3) )

参考文献

• U. Schollwöck, Rev. Mod. Phys., 77, 259 (2005). : https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.77.259