约化密度矩阵理论简介


\(
\def \bra#1{{\langle #1 \rvert}}
\def \ket#1{{\lvert #1 \rangle}}
\def \braket#1#2{{\langle #1 | #2 \rangle}}
\def \expect#1{{\langle #1 \rangle}}
\)如何用约化密度矩阵方法计算量子体系能量等问题,最早在约六七十年前就已经提出。但由于N-representability问题的悬而未决,在过去的几十年中一直没有解决。自上世纪90年代来,随着控制理论中semidefinite programming优化方法的发展,虽然N-representability问题依然没有完全解决,但利用p-positivity条件已经可以由2-RDM计算得到较好的结果。另一种不同但与之相关的方法,即采用Contracted Schrödinger equation的SCF解来进行计算,在3-RDM和4-RDM近似展开在2-RDM上的条件下,也能得到不错的结果。

虽然这种方法在近年发展缓慢,甚至有可能是一种不太可行的方法(非常感谢刘文剑老师课上提醒,N-representable 条件可能是无穷多的!),但这种方法的formulation其实非常漂亮,其背后思想十分简洁。作为一种 alternative viewpoint,也许能对我们有一些启发。

这篇文章主要是结合一下这方面工作的代表性人物 Mazziotti 教授 (JFI, UChicago) 写的一些综述,以及他在07年出的一本书《Reduced Density Matrix Mechanics: with Application to Many-Electron Atoms and Molecules》(见Advances in Chemical Physics第134卷),简单介绍一些这方面的基础。

1. Introduction & History

密度算符之定义为

\[ \rho = \ket{\psi}\bra{\psi} = \sum_{i} |c_i|^2 \ket{\phi_i}\bra{\phi_i} \triangleq \sum_k |c_k|^2 \ket k \bra k \]

(von Neumann, 1927)

若考虑我们的系统由两个子系统\(A\), \(B\)构成,其对应的希尔伯特空间\(H_A\), \(H_B\)的规范正交基分别为\( \{ \ket{a_i} \} \), \( \{ \ket{b_j} \} \)

对密度矩阵求对 \( B \) 的偏迹(partial trace),得

\[ \begin{align}
\rho_A &= Tr_A[ \rho ] \\
&= Tr_A[ \ket{\psi}\bra{\psi} ] \\
&= Tr_A [\sum_k |c_k|^2 \ket k \bra k ] \\
&= \sum_j \bra{b_j}(\sum_k |c_k|^2 \ket k \bra k ) \ket{b_j} & , k>j \\
\end{align}\]

此即子系统A的约化密度算符。(Dirac, 1930)

若将A系统定义为我们的N电子系统中的任意两个电子,则可以得到一个二阶约化密度矩阵(second order reduced density matrix, 2-RDM),该矩阵就是所谓的二电子约化密度矩阵(two-electron reduced density matrix)。我们从以上定义可以看出,2-RDM的矩阵元实际上等价于对N电子全波函数进行部分的积分,但留下某2个电子不积完。

由于电子是不可分辨的,且电子之间的相互作用仅仅存在成对相互作用(pair-wise interactions),因此,系统的总能量(Energy)仅仅是2-RDM的线性泛函。

\[ E = \sum_{ij, kl} H^{ij}_{kl} \Gamma^{ij}_{kl} \]

\[ H^{ij}_{kl} = \frac{1}{N-1}(\delta_{ik} h_{jl} + \delta_{jl} h_{ik} – \delta_{il} h_{jk} – \delta_{jk} h_{il} ) + V_{ij, kl} \]

式中\(h\)是单电子算符,\(V\)是双电子算符,\(\Gamma\)是2-RDM。这意味着,或许我们可以仅仅采用2-RDM而非N电子体系的全波函数来计算各种分子和原子的能量,从而不必去找全波函数Schrödinger equation的解。

早在1955年开始,Mayer等人就进行了相关的计算尝试。但他们很快发现,对2-RDM的变分优化,事实上能把能量优化到比真实能量更低!这意味着,Rayleigh-Ritz 变分原理 (Rayleigh-Ritz variational principle) 对于2-RDM并不适用。为何会出现这种情况呢?很快,Tredgold, Coleman, Coulson 等人就已经认识到,有一些trail 2-RDM事实上并不可能由N电子全波函数积分得到。因此,对2-RDM的变分优化必须在某种限制条件\(\sigma\)下进行,即所有的trail 2-RDM,在该组限制条件下必须满足可由N电子全波函数积分得到。这一充分必要条件\(\sigma\),被称为N-representability conditions (N-可表条件)。1964年,Garrod和Percus等给出 了该条件的一个形式:

\[\forall H^{(2)},  tr[H^{(2)}\Gamma^{(2)}] \geq E_0(H^{(2)}) \]

\( H^{(2)} \)是任何有下确界的双电子算符,\(\Gamma^{(2)}\)是2-RDM,\(E_0(H^{(2)})\) 是 \( H^{(2)} \) 的最小本征值。这一形式显然无法用于实际的计算,因此我们还需要找到更具体的形式。

从以上的结果看来,我们似乎只要简单地解决N-representability问题,找到对于2-RDM合适的约束条件,以上的问题便迎刃而解。我们甚至可以马上找到4个十分强的条件:对于全同费米子的密度矩阵,至少需要满足

  1. normalized (归一性,保证粒子数守恒)
  2. Hermitian (厄米性)
  3. anti-symmetric under particle exchange (反对称性,费米子的要求)
  4. positive semi-definite (半正定的,保证概率的半正定性)

然而,计算表明,这4个条件还并不充分,但看上去已经十分接近。事实上,Coleman在他1951年提出N-representability问题的时候,也认为这一问题很快就会被某个数学家——也许Coleman自己完全地解决。然而,至今已经接近70年,N-representability依然是一个悬而未决的问题。

上世纪90年代初,随着一种新的方法被提出,即使用contracted Schrödinger equation  (CSE)的自洽解来直接计算基态2-RDM,人们重新开始对2-RDM和N-representability发生兴趣。在2000年左右,又基于 D、Q、G 条件发现了一组重要的约束条件,即 positivity conditions。该组条件很快被拓展成了一系列分级的 N-representability constraints,即p-positivity conditions。 Mazziotti等证明了每一级的条件均对应于要求一类算符的广义不确定性关系。在p-positivity conditions约束下,多电子原子和分子的基态能量能通过变分法准确计算。

Verstichel 将 N-可表条件推广到了非整数电子数目的情形,该方法比较复杂,可求解来自包含不同电子数的波函数的系综平均的情况。

以下,我们介绍两种具体的计算方法:变分能量优化,和收缩薛定谔方程(Contracted Schrödinger Equation, CSE)的非变分解。与传统的量子化学方法相比,变分法求解2-RDM的计算量现在还仍然很大。另一种方法,即求解CSE,面临的主要困难是 CSE 不仅依赖于 2-RDM,它还依赖于 3-RDM 和 4-RDM;即CSE方法的核心问题是如何用 2-RDM 近似计算 3-RDM 和 4-RDM。此外,CSE 方法同样面临求解方程时出现的收敛困难问题。

与许多方法不同,在给定有限的basis set下,2-RDM的限制条件下变分计算给出的是基态能量的下界。2-RDM的计算在过渡态和严重变形的结构等波函数非常难以参数化的结构中也能给出有用的准确度,因此,比较适合用于这些情况。

2. Mathematical Background

一个矩阵\(\mathbf M\)是positive semidefinite的,当且仅当\(\mathbf M\)的所有本征值是非负实数。记作

\[\mathbf M \geq 0 \]

显然,\(\mathbf M\)也是Hermitian matrix。

在positivity constraints下求2-RDM变分解需要用到所谓的 semidefinite programming (半定规划)的技巧。以下我们简单地介绍这一方法。

[TODO] semidefinite programming, see https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-251j-introduction-to-mathematical-programming-fall-2009/readings/MIT6_251JF09_SDP.pdf [/TODO]

3. p-positivity conditions & Variational 2-RDM GS methods (v-2-RDM)

考虑在基态波函数\( \ket{\psi} \)上作用一系列算符 \( \{ \hat{C}_i \} \)来生成一个基组\( \{ \ket{\phi_i} \} \),

\[ \ket{\phi_i} =   \hat{C}_i \ket{\psi} \]

则这些基组的度量矩阵(metric matrix)或者说重叠(overlap)矩阵\(\mathbf M\)的矩阵元可以表达为

\[ M_{ij} =  \braket{\phi_i}{\phi_j} = \bra{\psi} \hat{C}_i^\dagger \hat{C}_j \ket{\psi}\]

注意到度量矩阵必须满足半正定条件。如果某一个度量矩阵的矩阵元能够仅用p-RDM的信息表达,那么该度量矩阵的半正定条件就是一个p-positivity condition

例如,考虑1-positivity条件的问题。假设将每一个 \(\hat C_j\) 选为湮灭算符 \(\hat a_j\),那么得到的one-particle reduced density matrix (1-RDM) 矩阵元为

\[ ^1D_{ij} = \bra \psi \hat{a}_i^\dagger \hat a_j \ket \psi \]

同理,若选为产生算符,则得到one-hole reduced density matrix (1-HRDM,单空穴约化密度算符)

\[ ^1Q_{ji} = \bra \psi \hat{a}_j \hat a_i^\dagger \ket \psi \]

注意到费米子的二次量子化(数目表象)的反对易关系

\[ \hat{a}_i^\dagger \hat a_j + \hat{a}_j \hat a_i^\dagger = \delta_{ij} \]

i. e.

\[ {}^1D_{ij} + {}^1Q_{ji} = \delta_{ij} \]

显然其给出了\( {}^1\mathbf D , {}^1\mathbf Q\)的矩阵元间的一个线性映射。由于Kronecker \(\delta\)在酉变换下不变,知\( {}^1\mathbf D , {}^1\mathbf Q\)共享本征态。从而1-HRDM的第i个本征值可用1-RDM的第i个本征值\(n_i\)表示为\(1 – n_i\)。因此,我们得到一组1-positivity条件:

\[\begin{cases}
n_i \geq 0 & i=1, 2, … \\
1 – n_i \geq 0 & i=1, 2, … \\
\end{cases}\]

注意到\(n_i\)事实上就是第i个自旋轨道的占据数,因此,1-positivity条件实际上就等价于Pauli不相容原理,即某个自旋轨道上的占据数在0~1之间。

1963年,Coleman证明了以上的1-positivity conditions,加上之前所述的四条基本的约束,就是1-RDM N-representable的充分必要条件。

我们现在再来继续考虑2-positivity conditions。类似于以上的做法,按照产生/湮灭算符的数量,将算符\(\hat C\)分为3组:

\[ \hat C \in \{\hat a_i \hat a_j\}, \hat C \in \{\hat a_i^\dagger \hat a_j\}, \hat C \in \{\hat a_i^\dagger \hat a_j^\dagger\} \]

注意到以上分组之间的正交性(由于产生/湮灭算符的不平衡,不同组的算符产生的基矢必定正交,故度量矩阵必定为0),最终我们将得到3个度量矩阵,其矩阵元为

\[ {}^2D^{p, q}_{i, j} = \bra \psi \hat a_q^\dagger \hat a_p^\dagger \hat a_i \hat a_j \ket \psi \]

\[ {}^2Q^{p, q}_{i, j} = \bra \psi \hat a_q \hat a_p \hat a_i^\dagger \hat a_j^\dagger \ket \psi \]

\[ {}^2G^{p, q}_{i, j} = \bra \psi \hat a_q^\dagger \hat a_p \hat a_i \hat a_j^\dagger \ket \psi \]

得到的2-positivity conditions为

\[\begin{cases}
{}^2 \mathbf D \geq 0 \\
{}^2 \mathbf Q \geq 0 \\
{}^2 \mathbf G \geq 0 \\
\end{cases}\]

其物理意义是分别保证在两个轨道上找到两个粒子,找到两个空穴,找到一个粒子和一个空穴的概率均不是负数。反复交换升降算符,则依然由反对易关系可知这3个条件间存在线性相关关系,若注意到\( {}^1D^p_i = \frac{1}{N-1} \sum_j {}^2D^{p, j}_{i, j} \),则可以将该关系写为一组简单的张量方程

\[ \begin{cases}
{}^2Q^{p, q}_{i, j} / 2 = (\delta^p_i – 2 {}^1D^p_i) \wedge \delta^q_j + {}^2D^{p, q}_{i, j} / 2 \\
\\
{}^2G^{p, q}_{i, j} = \delta^i_p {}^1D^q_j – {}^2D^{i, q}_{p, j} \\
\end{cases}\]

其中\(\wedge\)是Grassmann wedge product。

由以上方法得到的2-positivity conditions是2-RDM N-representable的必要条件,但与1-RDM时的情况不同,对于2-RDM,2-positivity conditions并非N-representable的充分条件。不过,从计算实验结果来看,仅仅作2-positivity conditions的限制,对许多体系,无论是平衡还是非平衡构型,计算基态能量都已经比较有效了。

p-positivity conditions 需要通过 p-RDM 才能计算得到。Erdahl 在此基础上,提出了 \(T_1\) 和 \(T_2\) 条件,这两个条件是通过组合 3-positive 条件得到的:

\[ {}^3(T_1)^{p, q, r}_{i, j, k} = \bra \psi \{ \hat a_r^\dagger \hat a_q^\dagger \hat a_p^\dagger, \hat a_i \hat a_j \hat a_k\} \ket \psi \]

\[ {}^3(T_2)^{p, q, r}_{i, j, k} = \bra \psi \{ \hat a_r^\dagger \hat a_q^\dagger \hat a_p, \hat a_i^\dagger \hat a_j \hat a_k \} \ket \psi \]

Mazziotti 在此基础上给出了 (2, q) 条件,即类比以上的组合,对于 q-positive 条件也可以进行系统地组合来产生半正定条件。(2, q) 条件应当可以系统地改进计算结果,不过,随着 q 的增大,计算量会迅速增加。

如前文所述,目前找到的所有N-可表条件均是必要条件而非充分条件。因此,变分空间超出了N-可表空间。显然,加入更多条件,应当改善变分所得能量下界,即

\[ E_{DQ} \leq E_{DQG} \leq E_{DQGT_1} \leq  E_{DQGT_1T_2} \leq E_{DQGT_1T_2′} \leq E_{FCI} \]

4. Contracted Schrödinger equations & non-variational solutions

1976年,Nakasuji 证明了 CSE 与 Schrodinger 方程的等价性。数目表象下 CSE 为

\[\bra \psi \hat a_i^\dagger \hat a_j^\dagger \hat a_l \hat a_k H \ket \psi = E \Gamma^{ij}_{kl}\]

可以看出,CSE 不仅依赖于2-RDM,还依赖 3-RDM 和 4-RDM。用 2-RDM 近似计算 3-RDM 和 4-RDM 的方法被称为”重建”法 (reconstruction)。重建法的理论基础是 Rosina’s theorem,即对任何基态非简并的且仅存在pair-wise interaction的系统,2-RDM 与波函数存在一一映射。该定理可用反证法和 Rayleigh-Ritz 变分原理证明。

1993年,Valdemoro 等利用二次量子化产生和湮灭算符的对易和反对易关系构造了第一个近似重建方法。1996年,Nakasuji基于Green函数和微扰理论给出了一个重建方法。Mazziotti 在 1999 年用 cumulant theory 给出了另一个展开方法。3-RDM在此方法下的展开式为

\[ {}^3\mathbf{D}/6 = {}^1\mathbf{D} \wedge {}^1\mathbf{D} \wedge {}^1\mathbf{D} + 3 {}^2\mathbf{\Delta} \wedge {}^1\mathbf{D} + {}^3\mathbf{\Delta} \]

,其中,

\[ {}^2\mathbf{\Delta} = {}^2\mathbf{D}/2 – {}^1\mathbf{D} \wedge {}^1\mathbf{D} \]

求解 CSE 方程一般需要迭代方法。迭代过程中, CSE 方程会偏离 N-可表条件中的一些重要条件,因此类比 linear-scaling 1-RDM 计算,需要对 2-RDM 进行“纯化” (purification)。

Mazziotti 在 2006 年进一步发展了反厄米CSE方法 (anti-Hermitian CSE)。ACSE 形式为

\[\bra \psi [\hat a_i^\dagger \hat a_j^\dagger \hat a_l \hat a_k, H ] \ket \psi = E \Gamma^{ij}_{kl}\]

显然,若CSE成立,ACSE也成立。反之不然。由于用 3-RDM 近似 2-RDM 的误差一般要大于用 ACSE 近似 CSE 的误差,故此问题可能影响不大。

5. Applications

2-RDM在一些体系的计算中取得了成功。继 Garrod 等在1975年的计算之后,Nakata 等在2001年基于20世纪90年代发展起来的 SPD 算法采用 D, Q, G 条件计算了一系列分子的能量。与 FCI 相比,相关能误差在120%以内。其后 Mazziotti 用类似的方法计算了几个分子的能量,发现在远离平衡构型的情况下与 FCI 结果吻合得仍然很好,相关能误差在 110% 左右。Zhao 等进一步采用 D, Q, G, \(T_1\), \(T_2\) 条件求解 2-RDM,对一系平衡位置附近的分子得到的结果精度与 CCSD(T) 相当。Mazziotti 等采用 D, Q, G, \(T_1\), \(T_2’\) 以及 3-positive 条件计算了一些分子的能量,结果显示 3-positive 得到的结果精度比\(T_1\), \(T_2’\) 有显著改善。

Figure 1. The shapes of the potential energy curves of the OH radical from the 2-RDM methods with DQG and DQGT2 conditions

Nakata 在2001年以及 Mazziotti 在2002 年采用的是 primal-dual interior point SDP 算法(PDSDP),采用 DQG 条件的计算量为 \(N^{12}\), N为基函数数目。这个方法虽然计算量很大,但它是二阶方法,收敛性非常好,而且能够提供计算结果的误差。 Mazziotti 其后发展了一阶的 RRSDP 算法,计算量降为\(N^{6}\),但可能会影响收敛性。Verstichel 等通过采用迭代求解线性方程组的方法把 PDSDP算法的计算量降为\(N^{6}\),但是在 DQG 条件下,这个方法的计算量仍为 CCSD 方法的 100~1000 倍。

最近 Mazziotti 又进一步发展了零阶的 boundary point SDP 算法(BPSDP),他报道 BPSDP 算法的计算速度比 RRSDP 快 10 ~ 20 倍。但 BPSDP 算法在接近收敛标准时,收敛速度很慢。虽然 BPSDP 速度明显比 PDSDP 更快,但是与其他 SDP 算法如 MBSDP (modified barrier SDP)方法相比计算效率相差不大。

另一方面,2-RDM也能与传统量子化学方法结合。例如,Mazziotti 提出可以在 CASSCF 计算中,将在活性空间中的 FCI 计算用 2-RDM 计算代替,从而可能可用于较大活性空间的计算。不过,这个方法对激发态能量计算存在问题。

6. Conclusions

与传统方法相比,约化密度矩阵方法提供了一个新的思路,并且在最近的20年内取得了很大的进步。2-RDM理论在近年的主要发展包括给出了一系列可系统改进的 N-可表条件,利用SDP方法减小了变分方法计算量,CSE方面发展了更高精度的重建方法,找到了一些克服收敛困难的算法,包括演化方程法等。

显然,虽然RDM方法计算精度不错并且有系统改进的方法,在目前来看,距离其完全实用尚有许多重要问题需要解决,包括寻找更有效的 N-可表条件;高效的 SDP 算法;更精确的重建和纯化方法;含时问题等。

References

[1] David A. Mazziotti. Acc. Chem. Res. 2006, 39, 207-215

[2] David A. Mazziotti. Reduced-Density-Matrix Mechanics – With Application to Many-Electron Atoms and Molecules (Advances in Chemical Physics, Volume 134), 2007, Wiley-Interscience

[3] 黎乐民, 中国学科发展战略 · 理论与计算化学, 2016, 科学出版社

2 thoughts on “约化密度矩阵理论简介”

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.