Exact Diagonalization
精确对角化¶
参考材料：
- Anders W. Sandvik 的综述 Computational Studies of Quantum Spin Systems
- 以及他的课程ppt XIV Training Course in the Physics of Strongly Correlated Systems
- Ralf Schneider, et.al. "Computational Many-Particle Physics" 的第 18 章
精确对角化实际上就是把哈密顿量写成矩阵表示，然后再求这个矩阵的最小本征值和本征矢，也就是基态和基态能量。
对于单电子图像的哈密顿量是很容易表示成矩阵形式的，多体哈密顿量的直接表示也不是很难。但是由于多体系统随着尺寸增加，矩阵的维数迅速增大，很容易就会超过计算机甚至超算的内存能力。所以精确对角化方法的核心就是如何利用系统的对称性把哈密顿量的矩阵划分成更小的矩阵。
下面先以海森堡自旋 $\frac{1}{2}$ 链为例说明如何利用对称性划分哈密顿量矩阵。最后再给出一个利用这些对称性的 Fermi-Hubbard 模型例子。其它所有一维电子模型都可以用与这两个例子类似的方法实现。
直接表示哈密顿量的矩阵¶
Heisenberg $S= \frac{1}{2}$ 链 $$ H = J\sum_{i=1}^{N} S_i \cdot S_{i+1} = J \sum_{i=1}^{N} \left( S_i^xS_{i+1}^x+S_i^yS_{i+1}^y+S_i^zS_{i+1}^z \right) $$ $$ H = J \sum_{i=1}^{N} \left[ \frac{1}{2} \left( S_i^+S_{i+1}^-+S_i^-S_{i+1}^+ \right) + S_i^zS_{i+1}^z \right] $$
为了简化起见，下面都取 $J=1$ ，即反铁磁Heisenberg模型。
算符张量积构造哈密顿量¶
一般单粒子的哈密顿量都会由算符的张量积构造得到，如下
import numpy as np
def tensor_product(*args):
if len(args) < 2:
return args[0]
result = args[0]
for i in range(1, len(args)):
result = np.kron(result, args[i])
return result
首先定义每个格点上的算符
Sz = np.array([[1, 0], [0, -1]])
Sp = np.array([[0, 1], [0, 0]])
Sm = np.array([[0, 0], [1, 0]])
II = np.eye(2)
接下来把所有格点的项全部加起来，这里以4个格点为例
H = (
1 / 2 * (tensor_product(Sp, Sm, II, II) + tensor_product(Sm, Sp, II, II))
+ tensor_product(Sz, Sz, II, II)
+ 1 / 2 * (tensor_product(II, Sp, Sm, II) + tensor_product(II, Sm, Sp, II))
+ tensor_product(II, Sz, Sz, II)
+ 1 / 2 * (tensor_product(II, II, Sp, Sm) + tensor_product(II, II, Sm, Sp))
+ tensor_product(II, II, Sz, Sz)
+ 1 / 2 * (tensor_product(Sm, II, II, Sp) + tensor_product(Sp, II, II, Sm))
+ tensor_product(Sz, II, II, Sz)
)
H111 = H
有了哈密顿量矩阵就可以对它进行对角化求基态了，这里调用 Arpack 求最小的本征值
from scipy.sparse.linalg import eigsh
eigvals, eigvecs = eigsh(H, k=1, which="SA")
eigvals
根据基矢构造哈密顿量¶
从上面最一般的过程能看出，构造哈密顿量需要做大量的张量积计算，而且这些张量积中乘的项大多数都是单位矩阵，也就是说哈密顿量矩阵是非常稀疏的，所以我们其实可以用稀疏矩阵来表示哈密顿量，并且直接找到有数字的那些基矢，往哈密顿量里填数字就行了。
在上面的张量积构造过程中我们实际上取的表象是按照格点顺序排列的，即 $$ \begin{array}{lll} |0\rangle &= |\uparrow\uparrow\uparrow\uparrow\rangle &= (0000) \\ |1\rangle &= |\uparrow\uparrow\uparrow\downarrow\rangle &= (0001) \\ |2\rangle &= |\uparrow\uparrow\downarrow\uparrow\rangle &= (0010) \\ |3\rangle &= |\uparrow\uparrow\downarrow\downarrow\rangle &= (0011) \\ \cdots &= \cdots &= \cdots \end{array} $$ 很容易发现基矢的表示自然地与数字的二进制联系起来了，因为每个格点恰好也只有两种可能的基矢，上下自旋分别对应于数字0和1。
这样我们就得到了一个关系 $$ \begin{array}{llll} S^z &\Rightarrow & 0\rightarrow 0 & 1\rightarrow 1 \\ S^- &\Rightarrow & 0\rightarrow \emptyset & 1\rightarrow 0 \\ S^+ &\Rightarrow & 0\rightarrow 1 & 1\rightarrow \emptyset \end{array} $$ 即，$S^z$ 作用无效果，$S^-$ 相当于如果bit是0就翻转为1，如果是1，这个基矢就消失掉。
综合起来，对于最近邻的两个态的四种情况，哈密顿量的作用结果为 $$ \begin{array}{lllll} & (00) & (01) & (10) & (11) \\ S_i^zS_{i+1}^z & (00) & -1(01) & -1(10) & (11) \\ S_i^+S_{i+1}^- &\emptyset & (10) &\emptyset &\emptyset \\ S_i^-S_{i+1}^+ &\emptyset &\emptyset & (01) &\emptyset \end{array} $$
接下来只要遍历所有的基矢，并在每个基矢上作用一遍哈密顿量的所有项，找到对应的基矢，把这两个基矢作为指标填上数字就行了。
def get_state_index(state: int, index: int) -> int:
"""获得一个整数某位置处的二进制值"""
mask = 1 << index
return (state & mask) >> index
def flip_state(state: int, index: int) -> int:
"""翻转一个整数某位置处的二进制值"""
mask = 1 << index
return state ^ mask
N = 4
H = np.zeros((2 ** N, 2 ** N))
# a is a basis
for a in range(2 ** N):
# i is position of this basis
for i in range(N):
# j is the nearest neighbor, mod N for periodic boundary conditions
j = (i + 1) % N
ai = get_state_index(a, i)
aj = get_state_index(a, j)
# Sz
b = a
if ai == aj:
H[a, b] += 1
else:
H[a, b] += -1
# SxSx + SySy
if ai != aj:
b = flip_state(a, i)
b = flip_state(b, j)
H[a, b] += 1 / 2
H112 = H
这样直接构造的结果与张量积乘出来的哈密顿量当然是一样的
eigvals, eigvecs = eigsh(H, k=1, which="SA")
eigvals, np.allclose(H112 - H111, 0)
利用对称性：组合型¶
从基矢构造哈密顿量的过程实际上分为两步，首先要确定基矢的总数和表象，之前我们取得表象就是最自然的全部排列，一共有 $2^N$ 个，然后把这些基矢按照哈密顿量的规则变换成另一组基矢，两组基矢作为哈密顿量矩阵的两个指标来填数字。
利用对称性主要就是改变前一步，不同的对称性子空间实际上限制的是基矢的数目，只要我们把某个子空间中的全部基矢找到，再用这组基矢构造哈密顿量，得到的就是这个对称性子空间中的哈密顿量了。只要把所有子空间都对角化一遍，自然就能得到全部的本征态。
这里先利用的对称性是最简单的一类，这些对称性只与基矢的总数有关，而与基矢内部的排列方式无关，例如总磁矩和总粒子数。总磁矩就是把每个格点上的磁矩全部加起来，总粒子数同样就是把每个格点上的粒子数加起来。总和一样的那些基矢就全部属于同一个子空间。反过来看，也可以是一个总的量子数在所有格点上的不同组合。
这里利用的是Heisenberg模型的总磁矩守恒 $$S^z_{\mathrm{tot}} = \sum_i s[i] $$ 例如 $|0011\rangle, |1010\rangle$ 这两个态的总磁矩就是一样的，都是 $1$。
所以找到总磁矩为 1 的所有基矢有两种方法，一是遍历一遍所有的态，二是按照组合数的方式找到N个格点中取2个的所有组合。
N = 4
# 这里的总自旋写2是因为总的因子 1/2 被省略了
Nz = 2
state_list = []
for a in range(2 ** N):
if bin(a).count("1") == Nz:
state_list.append(a)
state_list
from itertools import combinations
state_list = [sum([1 << i for i in c]) for c in combinations(range(N), Nz)]
state_list.sort()
state_list
接下来就是利用已经找到的这些基矢来构造哈密顿量，构造的过程完全一样，只是要注意，判断基矢规则时要用原来的基矢，之后要把这个基矢对应到新的子空间基矢里。
H = np.zeros((len(state_list),) * 2)
for new_a, a in enumerate(state_list):
for i in range(N):
j = (i + 1) % N
ai = get_state_index(a, i)
aj = get_state_index(a, j)
# Sz
new_b = new_a
if ai == aj:
H[new_a, new_b] += 1
else:
H[new_a, new_b] += -1
# SxSx + SySy
if ai != aj:
b = flip_state(a, i)
b = flip_state(b, j)
new_b = state_list.index(b)
H[new_a, new_b] += 1 / 2
可以发现对于4格点周期边界情况，基态是在总自旋为1的空间里的，因为这是个反铁磁的Heisenberg链。
eigvals, eigvecs = eigsh(H, k=1, which="SA")
eigvals
利用对称性：平移对称性¶
一维周期边界的链满足平移对称性，平移算符的本征态可以表示为 $$T | n \rangle = e^{ik}|n\rangle, k = m \frac{2\pi}{N}, m=0,1,\cdots,N-1$$ 平移算符的作用就是把态在格点上移动一个 $$T | S^z_1, S^z_2, \cdots, S^z_N \rangle = | S^z_N, S^z_1, \cdots, S^z_{N-1} \rangle $$
在平移算符的本征子空间中的基矢可以由投影算符构造 $$ P_k = \frac{1}{L} \sum_r^{L-1} e^{ i \frac{2\pi}{L} k r } T^r$$ 例如从 $|0011\rangle$ 生成的一组基矢为 $$ P_k |0011\rangle = \frac{1}{4} \left( |0011\rangle + e^{i \frac{2\pi}{4} k \cdot 1} |0110\rangle + e^{i \frac{2\pi}{4} k \cdot 2} |1100\rangle + e^{i \frac{2\pi}{4} k \cdot 3} |1001\rangle \right) $$ 这样就能够得到一组基矢，可以记作 $|r^a(k)\rangle=P_k|r^a\rangle$ ，对于其它的初始基矢还可以得到其它组基矢，记作 $|r^b(k)\rangle=P_k|r^b\rangle$ 等。
那么在这一组新的基矢 $\left\lbrace|r(k)\rangle\right\rbrace$ 下的哈密顿量一共有$L$ 个，分别对应 $k=0,1,\cdots,L-1$，它的矩阵元为 $$ H_{mn} = \langle r^m(k) | H | r^n(k) \rangle = \langle r^m| P_k H P_k | r^n \rangle = \langle r^m| P_k H | r^n \rangle$$ 其中用到 $[P_k, H]=0$ 因为 $[T,H]=0$ ，和 $P_k^2=P_k$ 因为 $P_k$ 是投影算符。注意：这里所有计算过程中保持归一化。
所以计算这个新的哈密顿量块只需要把每组新基矢在哈密顿量作用下进行变换，再求这个基矢与 $k=0$ 时的那个基矢的内积就行了。因为 $|r^a(k=0)\rangle = |r^a\rangle$ 。
def state_cycle_shift(state: int, shift: int, total: int) -> int:
"""循环移位操作"""
low = state >> (total - shift)
mask = 2 ** (total - shift) - 1
high = mask & state
return (high << shift) + low
首先要从所有态里面把属于同一组的找出来
state_list_temp = state_list.copy()
state_k_all = []
while len(state_list_temp) != 0:
state_k_list = []
s_first = state_list_temp[0]
s = s_first
while True:
state_k_list.append(s)
state_list_temp.remove(s)
s_next = state_cycle_shift(s, 1, N)
if s_next == s_first:
break
else:
s = s_next
state_k_list.sort()
state_k_all.append(state_k_list)
state_k_all
然后就可以构造哈密顿量了，构造哈密顿量的过程中需要计算哈密顿量与基矢的乘积，这里直接把基矢构造成矩阵，把最后一步整合成矩阵乘法了，否则的话要写嵌套4层循环才行。
k = 0
Hk = np.zeros((len(state_k_all),) * 2).astype(complex)
for new_b, b_list in enumerate(state_k_all):
# build state_b
state_new_b = np.zeros((1, len(state_list))).astype(complex)
for rb, b in enumerate(b_list):
ib = state_list.index(b)
state_new_b[0, ib] = 1
state_new_b = state_new_b / np.linalg.norm(state_new_b)
for new_a, a_list in enumerate(state_k_all):
state_new_a = np.zeros((len(state_list), 1)).astype(complex)
# build state_a
state_new_a_temp = np.zeros((len(state_list), 1)).astype(complex)
for ra, a in enumerate(a_list):
ia = state_list.index(a)
state_new_a_temp[ia, 0] = np.exp(1j * 2 * np.pi / N * k * ra)
state_new_a_temp = state_new_a_temp / np.linalg.norm(state_new_a_temp)
# H @ state_a
for ra, a in enumerate(a_list):
ia = state_list.index(a)
for i in range(N):
j = (i + 1) % N
ai = get_state_index(a, i)
aj = get_state_index(a, j)
# Sz
new_ia = ia
if ai == aj:
state_new_a[new_ia, 0] += state_new_a_temp[ia, 0]
else:
state_new_a[new_ia, 0] -= state_new_a_temp[ia, 0]
# SxSx + SySy
if ai != aj:
b = flip_state(a, i)
b = flip_state(b, j)
new_ib = state_list.index(b)
state_new_a[new_ib, 0] += 1 / 2 * state_new_a_temp[ia, 0]
# Hk
Hk[new_a, new_b] = state_new_b @ state_new_a
这里构造得到的 $k=0$ 子空间只有 $2\times 2$ 维，从这里计算得到的也是基态。一般来说，基态总是在 $k=0$ 的子空间中。
np.linalg.eigvalsh(Hk)
利用对称性：其它离散对称性¶
除了总磁矩和平移对称性之外还有一些其它的离散对称性，例如对于 $k=0,\pi$ 时有宇称对称性和自旋翻转对称性等等。但是这些对称性与平移对称性并不对易，所以不能从平移对称性的结果中继续细分，而是要联合平移对称性一起去找基矢。
比如联合使用平移对称性和宇称对称性，要同时考虑平移算符和宇称算符来构造基矢，例如 $$ \mathcal{P}_k = \frac{1}{L} \sum_r^{L-1} C^\sigma\left( \frac{2\pi}{L} k r \right) T^r (1+P) $$ $$ C^\sigma\left( \frac{2\pi}{L} k r \right) = \left\lbrace \begin{array}{ll} \cos( \frac{2\pi}{L} k r), & p = 1 \\ \sin( \frac{2\pi}{L} k r), & p = -1 \end{array} \right. $$
程序待写
Fermi-Hubbard 模型¶
$$ H = -t \sum_{i,\sigma} \left( c^\dagger_{i,\sigma} c_{i+1,\sigma} + c^\dagger_{i+1,\sigma} c_{i,\sigma} \right) + U \sum_{i} n_{i\uparrow}n_{i\downarrow}$$
首先定义哈密顿量作用到基矢上的结果，基本原则就是如果作用有效果就返回作用后的基矢，如果作用为0就返回 None。
from typing import Optional, Tuple
def check_U(state_up: int, state_down: int, i: int) -> Optional[Tuple[int, int]]:
state_up_i = get_state_index(state_up, i)
state_down_i = get_state_index(state_down, i)
if state_up_i == state_down_i == 1:
return state_up, state_down
else:
return None
def check_t(state: int, i: int, N: int) -> Optional[Tuple[int, int]]:
j = (i + 1) % N
state_i = get_state_index(state, i)
state_j = get_state_index(state, j)
if (state_i, state_j) == (0, 1) or (state_i, state_j) == (1, 0):
state_new = flip_state(state, i)
state_new = flip_state(state_new, j)
return state_new
else:
return None
先实现一个全对角化，以此作为检查标准。
N = 4
t = 1
U = 1
H = np.zeros((4 ** N, 4 ** N))
al = []
for a_up in range(2 ** N):
for a_down in range(2 ** N):
a = (a_up << N) + a_down
al.append(a)
for i in range(N):
# t
b_up = check_t(a_up, i, N)
if b_up is not None:
b = (b_up << N) + a_down
H[a, b] += -t
b_down = check_t(a_down, i, N)
if b_down is not None:
b = (a_up << N) + b_down
H[a, b] += -t
# U
b_result = check_U(a_up, a_down, i)
if b_result is not None:
b = (b_result[0] << N) + b_result[1]
H[a, b] += U
# eigvals, eigvecs = eigsh(H, k=1, which='SA')
# eigvals
eigvals_all = np.linalg.eigvalsh(H)
首先利用粒子数对称性，这里考虑的是粒子数的分别守恒，总数守恒就把两组整数一起组合就行了。
这里取半满且两种自旋数目相同。
N = 4
Nup = 2
Ndown = 2
t = 1
U = 1
from itertools import combinations
state_up_list = [sum([1 << i for i in c]) for c in combinations(range(N), Nup)]
state_up_list.sort()
state_down_list = [sum([1 << i for i in c]) for c in combinations(range(N), Ndown)]
state_down_list.sort()
state_total_list = [(i, j) for i in state_up_list for j in state_down_list]
state_up_list, state_down_list
接下来找平移对称性，因为有两组基矢，所以平移的时候要同时对它们都操作。
state_list_temp = state_total_list.copy()
state_k_all = []
while len(state_list_temp) != 0:
state_k_list = []
s_first = state_list_temp[0]
s = s_first
while True:
state_k_list.append(s)
state_list_temp.remove(s)
s0_next = state_cycle_shift(s[0], 1, N)
s1_next = state_cycle_shift(s[1], 1, N)
s_next = (s0_next, s1_next)
if s_next == s_first:
break
else:
s = s_next
state_k_list.sort()
state_k_all.append(state_k_list)
# state_k_all
接下来就构造哈密顿量。
这里取动量为$0$。
k = 0
Hk = np.zeros((len(state_k_all),) * 2).astype(complex)
for new_b, b_list in enumerate(state_k_all):
# build state_b
state_new_b = np.zeros((1, len(state_total_list))).astype(complex)
for rb, b in enumerate(b_list):
ib = state_total_list.index(b)
state_new_b[0, ib] = 1
state_new_b = state_new_b / np.linalg.norm(state_new_b)
for new_a, a_list in enumerate(state_k_all):
state_new_a = np.zeros((len(state_total_list), 1)).astype(complex)
# build state_a
state_new_a_temp = np.zeros((len(state_total_list), 1)).astype(complex)
for ra, a in enumerate(a_list):
ia = state_total_list.index(a)
state_new_a_temp[ia, 0] = np.exp(1j * 2 * np.pi / N * k * ra)
state_new_a_temp = state_new_a_temp / np.linalg.norm(state_new_a_temp)
# H @ state_a
for ra, a in enumerate(a_list):
a_up = a[0]
a_down = a[1]
ia = state_total_list.index(a)
for i in range(N):
# t
b_up = check_t(a_up, i, N)
if b_up is not None:
b = (b_up, a_down)
new_ib = state_total_list.index(b)
state_new_a[new_ib, 0] += -t * state_new_a_temp[ia, 0]
b_down = check_t(a_down, i, N)
if b_down is not None:
b = (a_up, b_down)
new_ib = state_total_list.index(b)
state_new_a[new_ib, 0] += -t * state_new_a_temp[ia, 0]
# U
b_result = check_U(a_up, a_down, i)
if b_result is not None:
b = b_result
new_ib = state_total_list.index(b)
state_new_a[new_ib, 0] += U * state_new_a_temp[ia, 0]
# Hk
Hk[new_a, new_b] = state_new_b @ state_new_a
# Hk
对角化一下新构造的哈密顿量，这个新哈密顿量的所有本征值都在全哈密顿量的谱里，这说明写对了。
从这个结果中我们知道，Fermi-Hubbard 模型的基态既是半满的，又是动量为$0$的。
eigvals = np.linalg.eigvalsh(Hk)
for x in eigvals:
index = np.abs(eigvals_all - x).argmin()
print(np.isclose(eigvals_all[index], x), eigvals_all[index], x)
评论
Comments powered by Disqus