量子蒙特卡洛方法是什麼？

01-15

和蒙特卡洛方法有什麼區別？都有哪些應用？

嗯……本人物理本科，最近想開搞 NQS 遂開始補 QMC 的背景知識。然後開始調研文獻之後驚人地發現 QMC 這玩意居然沒什麼中文文獻。大多數計算物理教材（國外不清楚，國內比如馬文淦老師的計算物理學）中對 QMC 也是只講述了 Path Integral Quantum Monte Carlo 這種並不常用的形式。所以我就來把我做的一個個人學習總結放上來（不過我相信如果認真去翻我們老師們的教案，肯定能找到比我寫的更好的2333）。

首先明確一點：QMC 可以說是 MC 的應用，即也是利用（偽）隨機數去求解數值問題，至少可以說和 Quantum Computing 沒啥直接聯繫。QMC 的優勢主要在於它可以求解強關聯繫統波函數，這是 DFT 無法做到的。但是由於 QMC 方法的複雜度隨著體系的規模呈指數式上升，所以在大體系的求解上自然沒有 DFT 來的常用。

之後安利一本教材： Thijssen Computational Physics 3rd ed.(ISBN-13 978-0-521-83346-2). 作者介紹說這是一本「銜接計算物理的基礎課程和相關領域專著或 Review Papers 的計算物理教材」，書中對與幾種常用的計算物理方法均有比較系統（且易於理解）的介紹。我這裡介紹 QMC 主要也就是翻譯了其中關於 QMC 的幾個章節，權當貢獻一點中文翻譯了233。

問題的提出

求解我們感興趣狀態的電子波函數主要就是對於定態多體 Schr?dinger 方程的求解：

$E | psi angle = H | psi angle$ .

其中 $|psi angle$ 表示被考慮的（多體）量子系統的態，用坐標表象表示就是多體波函數 $psi(mathbf x_1, mathbf x_2, cdots, mathbf x_n)$ .

我們可以使用通常用以求解 PDE 的方法如 FDM 和 FEM 去求解上述方程，但是這類方法在處理上面的方程時往往會遇到稀疏程度非常高的迭代矩陣（往往這類矩陣是如此的稀疏以至於使用常用的 sparceBLAS 類方法也會造成大量的資源浪費），不易處理。因此，QMC 之類的基於隨機數的方法在很多情形下是求解 Schr?dinger 方程的更有效的手段。

首先使用變分的思想考慮我們的問題：量子態（或者說, 波函數）滿足定態 Schr?dinger 方程等價於系統能量在這個態下的的數學期望取得穩定值：

$E [psi] = frac{langle psi | H | psi angle}{langle psi | psi angle} = frac{int psi^*(mathbf x) H psi(mathbf x) mathrm d^3x} {int psi^*(mathbf x) psi(mathbf x) mathrm d^3x}$ .

i.e. 若這個態改變了一微小的量 $delta psi$ ，能量期望的改變數不存在一階項：

$delta E[psi] approx 0$ .

====下面開始可以跳過====

為了證明上面這一論斷, 我們定義 $P = langle psi | H | psi angle, Q = langle psi | psi angle$ ，並將 $Delta E$ 寫到 $Delta psi$ 的一階項：

$egin{aligned}delta E = frac{langle psi + delta psi | H | psi + delta psi angle}{langle psi + delta psi | psi + delta psi angle} - frac{langle psi | H | psi angle} {langle psi | psi angle} \ onumber approx frac{langle delta psi | H | psi angle - (P / Q) langle delta psi | psi angle}{Q} + \ onumber frac{langle psi | H | delta psi angle - (P / Q) langle psi | delta psi angle}{Q}end{aligned}$ .

結合上面幾個式子，注意到 $E = P / Q$ ，我們得到：

$H | psi angle = E | psi angle$

也就是定態 Schr?dinger 方程。

當然，我們可以對於定態 Schr?dinger 方程直接取一組規範正交基，然後在改規範正交基下把定態 Schr?dinger 方程表示成為矩陣形式 $mathbf {HC} = E mathbf C$ 來求解，但是正如開頭所說的，這種做法會導致太多零矩陣元出現，不合適具體進行計算。

====上面為止可以跳過====

比較合理的計算方法是使用 Quantum Monte Carlo 方法（QMC），QMC 方法一共有三種（6/27/2017 修正：不是只有 3 種而是主要有 3 種）：變分 QMC （VMC）、擴散 QMC （DMC 或者叫 PMC：投影 QMC）和路徑積分 QMC。最後一個在好多基礎書上都有，我就不講了，主要寫一下前兩個：

Variational Monte Carlo

我們放棄考慮矩陣形式的方程，轉而用優化的思想討論如何找到使得能量 $E = frac{langle psi | H | psi angle}{langle psi | psi angle}$ 取得最小值的那個態。在正式開始 Quantum Monte Carlo 方法的介紹之前，我們先來類比一般的優化問題來簡要構造一下基於變分方法求解 Schr?dinger 方程的大致思路：

基於參數集 $mathcal S = (alpha_1, alpha_2, cdots, alpha_S)$ 構造試探多體波函數 $psi_{mathcal S} (mathbf x_1, mathbf x_2, cdots, mathbf x_n)$ ；
求能量的期望 $E [psi_{mathcal S}] = {langle psi_{mathcal S} | H | psi_{mathcal S} angle}/{langle psi_{mathcal S} | psi_{mathcal S} angle}$ ；
利用所得結果對 $mathcal S$ 進行評估, 循環。

評估能量期望和更新參數集的方法這裡暫時不予以討論，我們先來考慮上面步驟中求能量期望的這一步。在這裡，我們需要如下的積分式：

$E[psi_{mathcal S}] = frac{langle psi_{mathcal S} | H | psi_{mathcal S} angle}{langle psi_{mathcal S} | psi_{mathcal S} angle} = frac{int psi_{mathcal S}^*(mathbf x) H psi_{mathcal S}(mathbf x) mathrm d^3x} {int psi_{mathcal S}^*(mathbf x) psi_{mathcal S}(mathbf x) mathrm d^3x}$ .

這是一個高維積分，採用傳統的插值積分或者高維梯形積分的方法不僅構造繁瑣, 在電子數量 $n$ 上升時程序還不具有通用性。計算這個積分的最好方法自然就是使用 Monte Carlo 積分。

由於真實波函數在 Hilbert 空間的大部分區域的值都非常小，我們在求解積分式時又會不可避免地遇到稀疏問題。此時，如果採用傳統積分方法或者直接使用均勻隨機數進行 Monte Carlo 積分的話，依然不可避免地會造成大量的計算資源的浪費。在這種情形下，採用 Metropolis Monte Carlo （MMC）方法進行積分是比較可取的：使用大量粒子在波函數空間上進行隨機遊走，MMC 所對應的遊走規則會將它們「推向」波函數值較大的區域，達成較為高效的採樣積分的目的。

我們定義如下局部能量：

$E_{L mathcal S}(mathbf x) = frac{H psi_{mathcal S}(mathbf x)}{psi_{mathcal S}(mathbf x)}$ .

注意如果 $psi_{mathcal S}$ 與系統真實的基態 $psi(mathbf x)$ 相等, 這個局部能量是一個與 $mathbf x$ 無關的常數。這樣，系統能量的期望就能被寫成如下的形式：

$E[psi_{mathcal S}] = frac{int psi_{mathcal S}^2(mathbf x) E_{L mathcal S}(mathbf x) mathrm d^3 x}{int psi_{mathcal S}^2(mathbf x) mathrm d^3 x}.$

有了這樣形式的積分，我們就可以基於分布 $ho(mathbf x) = frac{psi_{mathcal S}^2 (mathbf x)}{int psi_{mathcal S}^2 (mathbf x ^prime) mathrm d^3 x ^prime}$ 進行 MMC 積分計算。

Diffuse/Projector Monte Carlo

除了基於評估和修改試探波函數的變分 QMC 方法，我們還可以將 Schr?dinger 的求解訴諸迭代法的思想，i.e. 通過一系列變換使得被研究系統的波函數 $(psi^{(i)}(mathbf x))$ 收斂到它的基態 $(psi^prime(mathbf x))$ 。

為了說明這一過程，我們先考慮單個粒子的自由隨機遊走：考慮一個位於一維離散時空間中的粒子：

$egin{array}{ll} x = n Delta l, n = cdots, -1, 0, 1, 2, cdots; \ t = m Delta t, m = 0, 1, 2, cdots, end{array}$

假設粒子向左或向右移動的概率都是 $(eta)$ ，而停留在原來位置的幾率為 $(1 - 2 eta)$ ，那麼在 $(t)$ 時刻找到粒子的概率服從 Markov 過程主方程：

$egin{aligned} p(x, t + Delta t) - p(x, t) = eta[p(x - Delta l, t) + p(x + Delta l, t) - 2 p(x, t)] onumber\ approx eta Delta x^2 frac{partial^2 p}{partial x^2}. onumber end{aligned}$

上式最左邊在 $(Delta t)$ 很小時也可以近似寫成 $(Delta t (partial p/partial t))$ ，再記 $(gamma = Delta x^2 eta / Delta t)$ ，我們可以得到上式的近似連續形式：

$frac{partial p(x, t)}{partial t} = gamma frac{partial^2 p(x, t)}{partial x^2}$ .

我們馬上發現這是一個含時擴散方程。於是我們可以就此寫出它的 Green 函數：

$G(x, x^prime, t) = frac{1}{sqrt{4 pi gamma t}} expleft(-frac{(x - x^prime)^2}{4 gamma t} ight),$

使得 $p(x, t) = int G(x, x^prime, t) p(x, t) mathrm dx^prime.$

在這裡我們注意到，可以通過如下方法構造另一個符合上述方程的 Markov 過程：取粒子 0 時刻的位置為 $(x_0)$ ，令它在 $(+ Delta t)$ 時刻遊走到 $(x_1)$ 的概率密度為 $(G(x_1, x_0, Delta t))$ ；完成這部遊走之後, 令它下一步 (時刻 $(x_2)$ ) 遊走到 $(+2 Delta t)$ 的概率密度為 $(G(x_2, x_1, Delta t))$ ，以此類推。用符號表示這一 Markov 過程就是：

$T_{Delta t}(x^prime ightarrow x) = G(x, x^prime, Delta t).$

該 Markov 過程和原先的隨機遊走的區別在於它利用了擴散方程的連續解的形式，使得它具有理論上更高的效率。這裡我們再使用兩個式子概括上面所描述的 Markov 過程：

$x(t + Delta t) = x(t) + eta sqrt{Delta t}$

式中 $(eta)$ 是一個滿足方差為 $(2 gamma)$ Gaussian 分布隨機變數：

$p(eta) = frac{1}{sqrt{4 pi gamma}} exp left(-frac{eta^2}{4 gamma} ight)$

擴散方程的一般形式是 $frac{partial ho(mathbf x, t)}{partial t} = mathcal L ho (mathbf x, t)$ ，其中 $mathcal L$ 是一個二階線性微分運算元。上式具有如下的形式解：

$ho(mathbf x, t) = exp (t mathcal L) ho(mathbf x, t).$

當然，由於上面形式解中含有微分算符的指數項，不能直接進行數值求解，但是從這個方程我們可以發現，上面關於 $mathcal L$ 方程的 Green 函數就是上面解中的微分算符 $(exp (t mathcal L))$ 在坐標表象下的矩陣元：

$G(mathbf x, mathbf x^prime, t) = langle mathbf x^prime | exp (t mathcal L) | mathbf x angle$

我們自然會想到使用 Green 函數的某種近似去構造 Markov 過程求解擴散方程。但是注意，Green 函數只有可以被歸一化， i.e. $(int G(mathbf x, mathbf x^prime, t) mathrm d^3 x = 1)$ 時才可以被用來構造 Markov 過程。我們會在下面對投影 QMC 的正式討論中遇到這個問題。

現在我們正式介紹投影 QMC 方法：考慮如下的擴散方程

$frac{partial ho(mathbf x, au)}{partial au} = frac12 abla^2 ho(mathbf x, au) - V(mathbf x) ho(mathbf x, au)$ ,

注意到, 如果我們在上式中令 $( au = mathrm i t)$ ，這個方程就成為了單位質量粒子的含時 Schr?dinger 方程。利用上面的討論，我們可以把解寫成：

$ho(mathbf x, au) = exp[ au (-K - V)] ho(mathbf x, 0)$ .

其中 $(K = p^2/2 = -1/2 abla^2)$ 是動能算符。由於 $(K, V)$ 不對易，上式的指數項不能直接被展開。但是，如果我們令 $( au)$ 取一個很小的值，那麼指數項中的 Campbell-Baker-Hausdoff （CBH）對易子就可以被忽略：

$exp[- au (K + V)] = exp(- au K) exp(- au V) + mathcal O ( au^2). onumber$

現在, 我們的主要任務就成為了找到上式右邊各指數項的顯式表達式。Green 函數勢能項的指數展開是平凡的，問題仍然出在如何找到動能項所對應的矩陣元：

$G_mathrm{kin}(mathbf x, mathbf x^prime, au) = langle mathbf x^prime | exp(- au K) | mathbf x angle.$

為了做到這件事，我們在上式中插入完備正交歸一化關係 $int | mathbf p angle langle mathbf p | mathrm d^3 p = 1$ ，並利用基本投影關係 $langle mathbf x | mathbf p angle = frac{1}{sqrt{2 pi}} exp(mathrm i mathbf p cdot mathbf x)$ ，由於動能算符 $(K)$ 在動量表象下是對角化的，我們可以算得：

$G_mathrm{kin}(mathbf x, mathbf x^prime, au) = frac{1}{sqrt{2 pi au}} exp left(-frac{(mathbf x - mathbf x^prime)^2}{2 au} ight)$ .

我們發現，上式作為虛時間 Schr?dinger 方程在不含勢能項 $(V)$ 情況下的 Green 函數實際上就是不含時擴散方程的解，這也使我們確信了它作為正確解的地位。現在我們把虛時間 Schr?dinger 方程的 Green 函數再寫一下：

$G(mathbf x, mathbf x^prime, au) = G_mathrm{kin}(mathbf x, mathbf x^prime, au) exp(- au V) + mathcal O ( au^2). onumber$

但不幸的是，勢能項因子 $(exp(- au V))$ 破壞了上面 Green 函數的歸一性。對上式進行歸一化的方法是在其上乘一個收斂因子 $(exp( au E_T))$ (但是注意, 我們預先不知道這個 $(E_T)$ 的值)。乘上了收斂因子的新 Green 函數當然已經不再是原方程的解，而是存在著一個移位：

$frac{partial ho(mathbf x, au)}{partial au} = frac12 abla^2 ho(mathbf x, au) - (V - E_T) ho(mathbf x, au)$

如果我們選取 $(E_T)$ 使得 Green 函數歸一化，那麼虛時間 Schr?dinger 方程就會存在一個不隨 $( au)$ 演化的解：

$E_T ho(mathbf x) = -frac12 abla^2 ho(mathbf x) + V(mathbf x) ho(mathbf x)$ .

這時我們就會發現，上式正是定態 Schr?dinger 方程。那麼接下來我們就需要求解它。

接下來我們考慮 Green 函數的性質：

$egin{aligned} G = langle mathbf x^prime | exp(- au(H - E_T)) | mathbf x angle\ label{eqn:dmc18} approx frac{exp(- au(V - E_T)) exp left(-frac{(mathbf x - mathbf x^prime)^2}{2 pi au} ight)}{sqrt{2 pi au}} end{aligned}$ .

我們先把算符 $(exp(- au(H - E_T)))$ 按它的本徵態 $(| phi_n angle)$ 展開：

$exp(- au(H - E_T))= sum_n | phi_n angle exp(- au(E_n - E_T)) langle phi_n |. onumber$

在 $( au)$ 很大的情況下基態 $(| phi_0 angle, E_0)$ 在上式以 $(exp[- au(E_1 - E_G)])$ 的優勢佔據主導地位。所以上式起到了一個將任意初態投影到基態的投影算符。而由於我們在 $( au)$ 很小的情況下的 $(G)$ 的近似表達式，我們必須多次將 Green 函數作用在波函數 $(psi^{(i)}(mathbf x))$ 上. 這就是投影 Monte Carlo。

在實際的計算模擬中，一批粒子被模擬在位形空間中進行隨機遊走，每個遊走過程分為擴散和分支兩步。在擴散步中，每個粒子按 Green 函數的擴散項進行隨機取點（由於 Gauss 函數的抽樣技術非常簡單，這裡應該直接使用解析的抽樣公式）；而在分支步中，每個粒子按勢能項（權函數）之比的概率 $expleft[frac{-(V(mathbf x^{(i + 1)}) - E_T)}{-(V(mathbf x^{(i)}))} ight]$

來決定是否真的遊走到所取的點 $(mathbf x^{(i + 1)})$ 。這個過程的效率依然不是很高，因為粒子很容易隨機遊走到勢能權重非常低的地方去，停留在這類地方會造成計算資源的浪費。為了緩解這個問題，我們可以引入出生--死亡法：每當粒子遊走到一個點 $(mathbf x^{(i + 1)})$ ，我們計算 $(mathbf x^{(i + 1)})$ 的權函數：

$q = exp[- au(V(mathbf x^{(i + 1)}) - E_T)]$ .

如果 $(q < 1)$ ,這個粒子就有 $(1 - q)$ 的幾率死亡；而如果 $(q > 1)$ ，這個粒子就會使 $(mathrm{floor}(q - 1))$ 或 $(mathrm{floor}(q))$ 個粒子在 $(mathbf x^{(i + 1)})$ 點、出生"，出生 $(mathrm{floor}(q))$ 個粒子的幾率由 $(q - mathrm{floor}(q))$ 給出。

我們注意到，出生--死亡法會使得進行隨機遊走的全部粒子數發生變化，而 $(E_T = E_0)$ 時隨機遊走的粒子數目在統計上是恆定的。所以我們可以在進行了一段時間的模擬之後對 $(E_T)$ 作如下調整：

$E_T = E_0 + alpha log ( ilde M/M),$

式中 $( ilde M)$ 是我們期望的粒子數而 $(M)$ 是真實觀察到的粒子數。在經過一段時間的調整之後，最終 $(E_T)$ 的平均值會收斂到系統的基態能量 $(E_0)$ ，而隨機遊走的粒子則會收斂到系統的基態波函數。

打完了。。原文在此： http://home.ustc.edu.cn/~rbh1998/QMC_intro.pdf

最後吐槽下編輯器有夠難受。。

既然已經知道了經典的蒙卡，那經典蒙卡是求解經典問題的，那麼量子蒙卡當然就是求解量子問題的了，經典問題可以直接寫出概率積分公式（就是求一個量在某個概率分布下的期望值，統計物理中，的經典模型，它的概率大都是玻爾茲曼分布），但是量子問題不能直接寫出概率分布（因為連能級都不知道），所以也不能直接寫出積分公式，所以從原來的求均值（就是含有狄拉克符號的那個求期望值的公式）的公式出發，變變變，想方設法變，到最後，仍然可以寫成經典那種形式，就可以直接算了，這就叫量子蒙卡，和經典蒙卡的區別就是多了中間的「變變變」那一步。

所謂「變變變」，說起來簡單，做起來可難了，針對不同的問題提出了不同變變變的方法，如果求基態，有投影蒙卡，求有限溫度的性質，有隨機序列展開方法（SSE，部分問題），還有行列式蒙卡，變分蒙卡，這些都是將量子模型轉化為經典模型的方法，但是蒙卡的核心沒變，就是得到概率分布，寫出積分公式。不同的量子蒙卡，所不同的，就是如何「變變變」。