3. 诸闭壳层量子化学方法的密度矩阵#

创建时间：2021-01-04；最后修改：2021-06-10

在这份简短笔记中，我们会回顾一些量子化学方法的密度矩阵，及其性质。大体的结论在下述表格中。

我们在这里只讨论闭壳层与实函数的情况。

方法	RHF 轨道基函数	能量关系	\(P_p^q\) 对称性	\(\Gamma_{pr}^{qs}\) 对称性	1-RDM 与电子数	\(\Gamma_{pr}^{qs}\) 与 \(P_p^q\) 的关系	\(\mathbf{P}\) 幂等性	\(P_i^a\) 为零	\(\Gamma_{ij}^{ab}\) 为零	\(F_p^q\) 对称性	1-RDM 偶极矩
RHF	√	√	√	√	√	√	√	√	√	√	√
Full-CI	√	√	√	√	√	√	×	×	×	√	√
MP2	√	√	√	√	√	×	×	√	×	×	×
CCSD	√	√	√	√	√	√	×	×	×	×	×
CISD	√	√	√	√	√	√	×	×	×	×	×
CASCI	×	√	√	√	√	√	×	N/A	N/A	×	×
CASSCF	×	√	√	√	√	√	×	N/A	N/A	√	√

之所以上面表格中 CASCI、CASSCF 方法不能说 \(P_i^a\) (密度矩阵的占据-非占) 与 \(\Gamma_{ij}^{ab}\) (2-RDM 的占据-非占)，是因为它们并是非基于 RHF 参考态的方法，不存在确切的占据与未占轨道。

3.1. 预定义#

from pyscf import gto, scf, mp, cc, ci, mcscf, fci
import numpy as np
from functools import partial

np.einsum = partial(np.einsum, optimize=True)
np.set_printoptions(precision=5, linewidth=150, suppress=True)

这里讨论的密度矩阵并非是 Full-CI 的情形矩阵，而是会因方法各异而不同的。

这里采用相对比较严格的 Einstein Summation Convention，即被求和角标必须是一个在上，一个在下。

这份文档中使用下述上下标：

\(p, q, r, s, m, n\) 分子轨道
\(i, j\) 分子占据轨道，\(a, b\) 分子未占轨道
\(\mu, \nu, \kappa, \lambda\) 原子轨道

分子轨道函数 \(\phi_p (\boldsymbol{r})\) 与原子轨道函数 \(\phi_\mu (\boldsymbol{r})\) 之间满足关系 (\(C_p^\mu\) 称为原子轨道系数)

\[ \phi_p (\boldsymbol{r}) = C_p^\mu \phi_\mu (\boldsymbol{r}) \]

我们假定了研究体系必然是实函数，但我们暂且定义函数的共轭记号如下：(不使用 Einstein Summation)

\[ \phi^p (\boldsymbol{r}) = \phi_p^* (\boldsymbol{r}) \]

分子轨道之间是正交归一的，但原子轨道需要用重叠积分：

\[ \int \phi_p (\boldsymbol{r}) \phi^q (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} = \delta_p^q, \; \int \phi_\mu (\boldsymbol{r}) \phi^\nu (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} = S_\mu^\nu \]

3.2. Full-CI 密度矩阵的定义与性质#

这里只作理论上的讨论。Full-CI 密度矩阵程序上的实现会在后面呈现。

密度矩阵与约化密度的定义有关。由于我们只讨论闭壳层情形，因此波函数可以安全地写成空间坐标的函数。

3.2.1. 1-RDM 与基转换关系#

我们回顾一阶约化密度 \(\rho(\boldsymbol{r}; \boldsymbol{r}')\)：

\[ \rho(\boldsymbol{r}; \boldsymbol{r}') = \idotsint \Psi^* (\boldsymbol{r}, \boldsymbol{r}_2, \boldsymbol{r}_3, \cdots, \boldsymbol{r}_{n_\mathrm{elec}}) \Psi (\boldsymbol{r}', \boldsymbol{r}_2, \boldsymbol{r}_3, \cdots, \boldsymbol{r}_{n_\mathrm{elec}}) \, \mathrm{d} \boldsymbol{r}_2 \, \mathrm{d} \boldsymbol{r}_3 \cdots \, \mathrm{d} \boldsymbol{r}_{n_\mathrm{elec}} \]

但是现在只有有限的基函数展开一阶约化密度；如果这组基函数是 RHF 分子轨道 \(\{ \phi_{p} (\boldsymbol{r}) \}\)，那么定义下述分子轨道基一阶约化密度矩阵 \(P_p^q\) (One-Order Reduced Density Matrix, 1-RDM)

\[ \rho(\boldsymbol{r}; \boldsymbol{r}') = P_p^q \phi^p (\boldsymbol{r}) \phi_q (\boldsymbol{r}') \]

如果是原子轨道 \(\{ \phi_\mu (\boldsymbol{r}) \}\)，那么它称为原子轨道基 1-RDM \(P_\mu^\nu\)

\[ \rho(\boldsymbol{r}; \boldsymbol{r}') = P_\mu^\nu \phi^\mu (\boldsymbol{r}) \phi_\nu (\boldsymbol{r}') \]

依据分子轨道与原子轨道间的关系，有

\[ \rho(\boldsymbol{r}; \boldsymbol{r}') = C_\mu^p P_p^q C_q^\nu \phi^\mu (\boldsymbol{r}) \phi_\nu (\boldsymbol{r}') \]

因此，原子轨道基与分子轨道基的 1-RDM 间存在关系

\[ P_\mu^\nu = C_\mu^p P_p^q C_q^\nu \]

3.2.2. 1-RDM 迹#

当 \(\boldsymbol{r}, \boldsymbol{r}'\) 相同时，我们会将一阶约化密度简记为电子密度 \(\rho(\boldsymbol{r}) = \rho(\boldsymbol{r}; \boldsymbol{r})\)。

一阶约化密度具有积分为电子数的性质：

\[ \int \rho(\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} = n_\mathrm{elec} \]

在分子轨道基的表示下，上式可以写为

\[ P_p^q \int \phi^p (\boldsymbol{r}) \phi_q (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} = P_p^q \delta^p_q = \mathrm{tr} (\mathbf{P}) = n_\mathrm{nelec} \]

3.2.3. 1-RDM 对称性#

首先我们可以证明下述等式：

\[\begin{split} \begin{align} &\quad\ \iint \phi_p (\boldsymbol{r}) \rho(\boldsymbol{r}; \boldsymbol{r}') \phi^q (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' \\ &= \iint \phi_p (\boldsymbol{r}) P_r^s \phi^r (\boldsymbol{r}) \phi_s (\boldsymbol{r}') \phi^q (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' \\ &= P_r^s \delta_p^r \delta_s^q = P_p^q \end{align} \end{split}\]

对于上式，如果我们交换被积元变量 \(\boldsymbol{r}, \boldsymbol{r'}\)、并对表达式取共轭，得到 (不使用 Einstein Summation)

\[ \iint \phi_p (\boldsymbol{r}) \rho(\boldsymbol{r}; \boldsymbol{r}') \phi^q (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' = \iint \phi_q (\boldsymbol{r}) \rho^*(\boldsymbol{r}'; \boldsymbol{r}) \phi^p (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' \]

如果我们再利用实数情形下，根据一阶约化密度的定义，有 \(\rho(\boldsymbol{r}; \boldsymbol{r}') = \rho(\boldsymbol{r}'; \boldsymbol{r}) = \rho^*(\boldsymbol{r}'; \boldsymbol{r})\)，那么可以立即得到 (不使用 Einstein Summation)

\[ P_p^q = \iint \phi_p (\boldsymbol{r}) \rho(\boldsymbol{r}; \boldsymbol{r}') \phi^q (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' = \iint \phi_q (\boldsymbol{r}) \rho(\boldsymbol{r}; \boldsymbol{r}') \phi^p (\boldsymbol{r}') \, \mathrm{d} \boldsymbol{r} \, \mathrm{d} \boldsymbol{r}' = P_q^p \]

即 1-RDM 矩阵 \(P_p^q\) 是对称矩阵。

3.2.4. 2-RDM#

与 1-RDM 相同地，依据二阶约化密度 \(\gamma (\boldsymbol{r}_1, \boldsymbol{r}_2; \boldsymbol{r}'_1, \boldsymbol{r}'_2)\) 的定义：

\[ \gamma (\boldsymbol{r}_1, \boldsymbol{r}_2; \boldsymbol{r}'_1, \boldsymbol{r}'_2) = \idotsint \Psi^* (\boldsymbol{r}_1, \boldsymbol{r}_2, \boldsymbol{r}_3, \cdots, \boldsymbol{r}_{n_\mathrm{elec}}) \Psi (\boldsymbol{r}'_1, \boldsymbol{r}'_2, \boldsymbol{r}_3, \cdots, \boldsymbol{r}_{n_\mathrm{elec}}) \, \mathrm{d} \boldsymbol{r}_3 \cdots \, \mathrm{d} \boldsymbol{r}_{n_\mathrm{elec}} \]

用分子轨道基作展开，可以定义分子轨道基的二阶约化密度矩阵 \(\Gamma_{pr}^{qs}\) (Two-Order Reduced Density Matrix, 2-RDM)

\[ \gamma (\boldsymbol{r}_1, \boldsymbol{r}_2; \boldsymbol{r}'_1, \boldsymbol{r}'_2) = \Gamma_{pr}^{qs} \phi^p (\boldsymbol{r}_1) \phi_q (\boldsymbol{r}'_1) \phi^r (\boldsymbol{r}_2) \phi_s (\boldsymbol{r}'_2) \]

原子与分子轨道基转换也与 1-RDM 类似：

\[ \Gamma_{\mu \kappa}^{\nu \lambda} = C_\mu^p C^\nu_q \Gamma_{pr}^{qs} C_\kappa^r C^\lambda_s \]

3.2.5. 2-RDM 与 1-RDM 的关系#

出于全同粒子的性质，2-RDM 与 1-RDM 之间存在关系：

\[ \rho(\boldsymbol{r}_1; \boldsymbol{r}'_1) = \frac{1}{n_\mathrm{elec} - 1} \iint \gamma (\boldsymbol{r}_1, \boldsymbol{r}_2; \boldsymbol{r}'_1, \boldsymbol{r}_2) \, \mathrm{d} \boldsymbol{r}_2 \]

对上式展开并作一部分积分后，可以得到

\[ P_p^q \phi^p (\boldsymbol{r}_1) \phi_q (\boldsymbol{r}'_1) = \frac{1}{n_\mathrm{elec} - 1} \Gamma_{pr}^{qm} \phi^p (\boldsymbol{r}_1) \phi_q (\boldsymbol{r}'_1) \delta_m^r \]

由于上式要在任意的 \(\boldsymbol{r}_1, \boldsymbol{r}'_1\) 的取值下成立，因此可以认为

\[ P_p^q = \frac{1}{n_\mathrm{elec} - 1} \Gamma_{pr}^{qm} \delta_m^r \]

注意上式要对等式右边作关于 \(r, m\) 角标的求和。

3.2.6. 2-RDM 对称性#

分析 2-RDM 对称性相对比较麻烦。这里就略过讨论了。我们仅指出，实数闭壳层下的 2-RDM 应当具有二重对称性，而不具有更高的对称性：

\[ \Gamma_{pr}^{qs} = \Gamma_{rp}^{sq} \]

3.2.7. 密度矩阵与能量#

这是最关键的一个性质。密度矩阵可以用来表示电子态的能量。

现在记原子轨道基组下的单电子算符积分为 \(h_\mu^\nu\)、双电子算符积分为 \(g_{\mu \kappa}^{\nu \lambda}\)，其中单电子算符包含动能、原子核-电子库伦势能、电场势能等贡献，双电子算符包含电子-电子库伦势能贡献。那么，体系单点能为

\[ E_\mathrm{tot} = E_\mathrm{elec} + E_\mathrm{nuc} = h_\mu^\nu P_\nu^\mu + \frac{1}{2} g_{\mu \kappa}^{\nu \lambda} \Gamma_{\nu \lambda}^{\mu \kappa} + E_\mathrm{nuc} \]

3.2.8. 1-RDM 与偶极矩#

Full-CI 的 1-RDM 可以直接用以计算偶极矩；以 \(z\) 轴方向施加的电场为例：

\[ d_z = - z_\mu^\nu P_\nu^\mu \]

其中 \(z_\mu^\nu = \langle \mu | z | \nu \rangle\)。

3.2.9. 广义 Fock 矩阵#

广义 Fock 矩阵定义为

\[ F_p^q = h_p^r P_r^q + g_{pr}^{ms} \Gamma_{ms}^{qr} \]

特别地，在 RHF 下，它是对角矩阵。而在 Full-CI 下，它会是对称矩阵。

3.3. RHF 密度矩阵特有性质#

3.3.1. 幂等性#

在 Conanical-HF 下，幂等性是几乎显然的：\(P_p^q\) 一定是对角矩阵；如果 \(p\) 所代表的轨道是占据轨道，那么一定填了因为填了两个电子而值为 2，否则为零。因此，Conanical-HF 下一定满足

\[ P_p^m P_m^q = 2 P_p^q \]

一般程序都只会给出 Conanical-HF 的结果。但若讨论 Nonconanical-HF 时 \(P_p^q\) 未必是对角矩阵，但上述结论应仍然成立。

3.3.2. 非占-占据部分为零#

Hartree-Fock 方法严格地将轨道分为占据与非占据。因此，Canonical 或 Nonconanical HF 方法都会保证 1-RDM 是块状对角化的；即在占据-非占 \(P_i^a\)，非占-占据 \(P_a^i\)，非占-非占 \(P_a^b\) 均严格为零。对于 2-RDM 是类似的。

由于 Hartree-Fock 方法没有考虑非占轨道的贡献，因此任何 Post-HF 方法均一定程度上有激发态的贡献。一般来说，非占-非占的 \(P_a^b\) 贡献总是存在的；但占据-非占或非占-占据的 \(P_i^a\) 与 \(P_a^i\) 则未必存在。

3.4. 通用计算函数#

下面的文档仅仅是验证开头的表格的代码。

3.4.1. 分子定义：水分子#

mol 水分子实例；
nelec \(n_\mathrm{elec}\) 电子数；
nocc \(n_\mathrm{occ}\) 占据轨道数；
h 原子轨道基 \(h_\mu^\nu\)，维度 \((\mu, \nu)\)；
g 原子轨道基 \(g_{\mu \kappa}^{\nu \lambda}\)，维度 \((\mu, \nu, \kappa, \lambda)\)；
S 原子轨道基 \(S_\mu^\nu\)，维度 \((\mu, \nu)\)；
mf_rhf RHF 实例。

mol = gto.Mole()
mol.atom = """
O  0. 0. 0.
H  0. 0. 1.
H  0. 1. 0.
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()

<pyscf.gto.mole.Mole at 0x7feae8ff9580>

nelec = mol.nelectron
nocc = mol.nelec[0]
nelec, nocc

(10, 5)

h = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
g = mol.intor("int2e")
S = mol.intor("int1e_ovlp")

mf_rhf = scf.RHF(mol).run()

3.4.2. 验证能量表达式#

验证

\[ E_\mathrm{tot} = h_\mu^\nu P_\nu^\mu + \frac{1}{2} g_{\mu \kappa}^{\nu \lambda} \Gamma_{\nu \lambda}^{\mu \kappa} + E_\mathrm{nuc} = h_p^q P_q^p + \frac{1}{2} g_{pr}^{qs} \Gamma_{qs}^{pr} + E_\mathrm{nuc} \]

eng_nuc = mol.energy_nuc()

def verify_energy_relation(eng, eng_nuc, rdm1, rdm2, h_mo, g_mo):
    return np.allclose(np.einsum("pq, qp ->", h_mo, rdm1) + 0.5 * np.einsum("pqrs, qpsr ->", g_mo, rdm2) + eng_nuc, eng)

3.4.3. 验证 1-RDM 对称性#

验证 \(P_p^q = P_q^p\)。

def verify_rdm1_symm(rdm1):
    # Output: 1-RDM symmetric property
    return np.allclose(rdm1, rdm1.T)

3.4.4. 验证 2-RDM 对称性#

验证 \(\Gamma_{pr}^{qs} = \Gamma_{rp}^{sq}\)。

def verify_rdm2_symm(rdm2):
    return np.allclose(rdm2, np.einsum("pqrs -> rspq", rdm2))

3.4.5. 验证 1-RDM 的迹#

验证 \(P_p^r \delta_r^p = n_\mathrm{elec}\)。

def verify_rdm1_tr(rdm1):
    return np.allclose(rdm1.trace(), nelec)

3.4.6. 验证 1-RDM 与 2-RDM 的关系#

验证 \(P_p^q = (n_\mathrm{elec} - 1)^{-1} \Gamma_{pr}^{qm} \delta_m^r\)。

def verify_rdm12_relation(rdm1, rdm2):
    return np.allclose(rdm1, (nelec - 1)**-1 * rdm2.diagonal(axis1=-1, axis2=-2).sum(axis=-1))

3.4.7. 验证 1-RDM 幂等性#

验证 \(P_p^m P_m^q = 2 P_p^q\)。

def verify_rdm1_idomp(rdm1):
    return np.allclose(rdm1 @ rdm1, 2 * rdm1)

3.4.8. 验证 \(P_i^a\) 为零#

这里实际上同时验证 \(P_a^i\) 是否为零。

def verify_rdm1_ov(rdm1):
    mat1 = rdm1[nocc:, :nocc]
    mat2 = rdm1[:nocc, nocc:]
    return np.allclose(mat1, np.zeros_like(mat1)) and np.allclose(mat2, np.zeros_like(mat2))

3.4.9. 验证 \(\Gamma_{ij}^{ab}\) 为零#

def verify_rdm2_ovov(rdm2):
    mat = rdm2[:nocc, nocc:, :nocc, nocc:]
    return np.allclose(mat, np.zeros_like(mat))

3.4.10. 验证广义 Fock 矩阵对称性#