跳过正文

行列式驱动量子化学方法的编程实践(一)

·3830 字·8 分钟
Che
作者
Che
仅仅是个普通人
Determinant-based Quantum Chemistry - 这篇文章属于一个选集。
§ 1: 本文

概述
#

  在现代量子化学中,行列式驱动的方法,如 Hartree-Fock (HF) 及其后继的 post-HF 方法,因其在处理多电子系统时展现出计算效率与精度的卓越平衡而占据核心地位。这些方法的理论基础深深植根于 二次量子化 的表述。二次量子化作为标准量子力学在多体问题中的一种优雅而强大的语言,通过引入 产生算符 \(a_p^\dagger\) 和 湮灭算符 \(a_p\),能够灵活地构造 Fock 空间 中的任意算符,例如哈密顿量和波函数。这种框架不仅精简了复杂量子系统的数学描述,更关键的是,它通过算符之间的 反对易关系,自动确保了电子作为费米子所要求的波函数反对称性。这一特性在计算化学中至关重要,是行列式驱动方法能够高效实现的关键。

Slater 行列式的构造与反对称性
#

  考虑一个由 \(N\) 个电子组成的分子系统,其电子结构在一组 \(M\) 个正交 自旋轨道 \({\phi_p(\mathbf{x})}_{p=1}^M\) 下描述,其中 \(\mathbf{x} = (\mathbf{r}, \sigma)\) 表示电子的空间坐标 \(\mathbf{r}\) 和自旋 \(\sigma\)。这些自旋轨道被视为单粒子基函数,其具体形式在此不作深入讨论。根据 泡利不相容原理,多电子波函数必须在交换任意两个电子坐标时表现出反对称性。在 一次量子化 中,这一性质通过 反对称化算符 \(\hat{A}\) 实现:

$$ \hat{A} = \frac{1}{\sqrt{N!}} \sum_{P \in S_N} \epsilon(P) P, $$

其中 \(S_N\) 是 \(N\) 阶置换群,\(P\) 为置换操作,\(\epsilon(P)\) 表示置换的奇偶性(+1 为偶置换,-1 为奇置换)。将 \(\hat{A}\) 作用于 \(N\) 个自旋轨道的简单乘积态(即 Hartree 积):

$$ \Phi_{\text{Hartree}} = \phi_{p_1}(\mathbf{x}_1) \phi_{p_2}(\mathbf{x}_2) \cdots \phi_{p_N}(\mathbf{x}_N), $$

即可得到归一化的 Slater 行列式

$$ \Phi_{p_1 p_2 \cdots p_N}(\mathbf{x}_1, \dots, \mathbf{x}_N) = |\phi_{p_1} \phi_{p_2} \cdots \phi_{p_N}| = \frac{1}{\sqrt{N!}} \begin{vmatrix} \phi_{p_1}(\mathbf{x}_1) & \phi_{p_2}(\mathbf{x}_1) & \cdots & \phi_{p_N}(\mathbf{x}_1) \\ \phi_{p_1}(\mathbf{x}_2) & \phi_{p_2}(\mathbf{x}_2) & \cdots & \phi_{p_N}(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_{p_1}(\mathbf{x}_N) & \phi_{p_2}(\mathbf{x}_N) & \cdots & \phi_{p_N}(\mathbf{x}_N) \\ \end{vmatrix}. $$

  Slater 行列式的代数结构天然满足反对称性:交换任意两个电子坐标 \(\mathbf{x}_i\) 和 \(\mathbf{x}_j\),波函数变号:

$$ \Phi(\cdots, \mathbf{x}_i, \cdots, \mathbf{x}_j, \cdots) = -\Phi(\cdots, \mathbf{x}_j, \cdots, \mathbf{x}_i, \cdots), $$

这与行列式交换两行引入负号的性质一致。

  从数学视角看,Slater 行列式是 \(N\) 个自旋轨道在 外代数(Exterior Algebra) 中的 楔积(Wedge Product) 在特定坐标表象下的实现:

$$ \phi_{p_1} \wedge \phi_{p_2} \wedge \cdots \wedge \phi_{p_N} = \frac{1}{\sqrt{N!}} \sum_{P \in S_N} \epsilon(P) \phi_{P(p_1)} \otimes \phi_{P(p_2)} \otimes \cdots \otimes \phi_{P(p_N)}. $$

楔积是构造反对称张量的普适方法,而 Slater 行列式则是其在电子坐标空间中的矩阵表达。在量子化学中,“行列式”、“反对称积”和“Slater 行列式”常互换使用,但需明确行列式是这一概念在特定基函数下的具体形式。

二次量子化与 Fock 空间
#

  二次量子化 将反对称性的处理提升至更高的抽象层次,并融入 Fock 空间 的代数结构。Fock 空间是由所有可能粒子数(包括真空态)的反对称态张成的希尔伯特空间。对于 \(N\) 电子系统,其状态空间是 Fock 空间的一个子空间。该空间的基矢(即 Slater 行列式)可通过 占据数(Occupation Number, ON)向量 \((n_1, n_2, \ldots, n_M)\) 唯一标识,其中 \(n_p = 1\) 表示自旋轨道 \(\phi_p\) 被占据,\(n_p = 0\) 表示空轨道,且 \(\sum_p n_p = N\)。例如,Slater 行列式 \(|\phi_{p_1} \cdots \phi_{p_N}|\) 对应 ON 态 \(|n_1, \ldots, n_M\rangle\),其中仅当 \(p \in {p_1, \ldots, p_N}\) 时 \(n_p = 1\),其余为 0。这种表示将多电子波函数从依赖 \(N\) 个显式坐标的函数简化为离散轨道占据信息的抽象态矢量,大幅简化了表示和运算。

  尽管一次量子化中的 Slater 行列式(坐标空间函数)与二次量子化中的 ON 态(Fock 空间基矢)形式不同,它们描述的是同一个 \(N\) 电子反对称量子态,具有等价的物理意义。在量子化学中,这两种表述(以及通过产生算符作用于真空态生成的 Fock 态)常交替使用,具体取决于上下文。

反对称性的算符代数实现
#

  在二次量子化中,费米子的反对称性由 产生算符 \(a_p^\dagger\) 和湮灭算符 \(a_p\) 的反对易关系 保证:

$$ { a_p, a_q^\dagger } = a_p a_q^\dagger + a_q^\dagger a_p = \delta_{pq}, $$

$$ { a_p, a_q } = a_p a_q + a_q a_p = 0, $$

$$ { a_p^\dagger, a_q^\dagger } = a_p^\dagger a_q^\dagger + a_q^\dagger a_p^\dagger = 0. $$

其中,\(a_p^\dagger\) 在真空态 \(|\text{vac}\rangle\) 上产生一个占据 \(\phi_p\) 的电子,\(a_p\) 则湮灭之。任意 \(N\) 电子 Slater 行列式可表示为:

$$ | \Phi_{p_1 p_2 \cdots p_N} \rangle = a_{p_1}^\dagger a_{p_2}^\dagger \cdots a_{p_N}^\dagger | \text{vac} \rangle. $$

由于 \({a_p^\dagger, a_q^\dagger} = 0\),即 \(a_p^\dagger a_q^\dagger = -a_q^\dagger a_p^\dagger\),交换两个产生算符引入负号:

$$ \cdots a_{p_j}^\dagger \cdots a_{p_i}^\dagger \cdots | \text{vac} \rangle = - \cdots a_{p_i}^\dagger \cdots a_{p_j}^\dagger \cdots | \text{vac} \rangle, $$

这与 Slater 行列式交换行列的性质一致。此外,\(a_p^\dagger a_p^\dagger = 0\) 直接反映了泡利原理——同一自旋轨道不可重复占据。因此,二次量子化通过算符代数隐式保证反对称性,无需显式构造行列式。

Slater行列式(ON向量)的计算机存储
#

  我们前面提到,Slater行列式和 ON 向量在Fock空间表述下的量子化学中是物理等价的两个概念。在接下来的内容中,我们将混合使用这两个名词。

  在Fock空间中,单个Slater行列式可以简洁地表示为一个长度为 \(N_{SO}\) 的向量,其中 \(N_{SO}\) 代表自旋轨道的总数。向量的每个元素取值为0(未占据)或1(占据)。例如,\( |1,1,0,0,0,0\rangle \) 表示一个包含6个自旋轨道的体系,其中前两个自旋轨道被电子占据。然而,在计算机中高效存储和操作这些行列式时,直接使用这种显式向量形式并不理想,通常会采用更紧凑且运算效率更高的**位串(bitstring)**编码方式。

位串编码(Bitstring Encoding)
#

  位串编码的基本思想是将每个自旋轨道映射为一个二进制位:自旋轨道被电子占据时对应位为1,未被占据时为0。这种表示方法不仅直观,还能显著节省存储空间。考虑到电子的自旋特性,为了便于计算(如处理自旋守恒的激发或应用自旋相关的算符),通常将自旋向上(\(\alpha\), \(\uparrow\))和自旋向下(\(\beta\), \(\downarrow\))的轨道占据情况分别存储在两个独立的位串中,即 alphabeta。因此,一个完整的Slater行列式在计算机中通常由这一对位串来表示。

  在量子化学中,首先需要定义一套包含 \(N_{mo}\) 个空间轨道的基组 \(\{\phi_k\}_{k=1}^{N_{mo}}\)。每个空间轨道 \(\phi_k\) 可与 \(\alpha\) 自旋函数结合形成 \(\alpha\) 自旋轨道 \(\psi_{k\alpha}\),或与 \(\beta\) 自旋函数结合形成 \(\beta\) 自旋轨道 \(\psi_{k\beta}\)。由此,\(\alpha\) 和 \(\beta\) 自旋轨道的数量均等于空间轨道数 \(N_{mo}\)。在限制性(Restricted)轨道表述中,第 \(k\) 个 \(\alpha\) 自旋轨道和第 \(k\) 个 \(\beta\) 自旋轨道共享相同的空间函数 \(\phi_k\)。因此,alpha 位串的第 \(j\) 位和 beta 位串的第 \(j\) 位分别表示第 \(j+1\) 个空间轨道的 \(\alpha\) 和 \(\beta\) 自旋的占据情况。

在非限制性(unrestricted)方法中,\(\alpha\) 和 \(\beta\) 电子可拥有不同的空间轨道,但位串编码的基本原理不变。而在广义(generalized)方法中,自旋轨道是 \(\alpha\) 和 \(\beta\) 自旋函数的线性组合,需采用不同的编码方式。

  现代计算机体系结构对固定长度整数类型(如64位无符号整数 uint64_t)的处理尤为高效。因此,当 \(N_{mo} \leq 64\) 时,alphabeta 各可用一个 uint64_t 变量存储。我们可以定义一个简单的结构体来表示行列式的占据情况:

struct Determinant {
    uint64_t alpha;  // Alpha自旋轨道占据情况
    uint64_t beta;   // Beta自旋轨道占据情况
};

若 \(N_{mo} > 64\),则需使用 uint64_t 数组(如 std::vector<uint64_t>)存储每个位串。例如,对于100个空间轨道,需 \(\lceil 100/64 \rceil = 2\) 个 uint64_t 整数来分别存储 alphabeta。不过,对于许多中小型体系,64位整数已足够。

相位因子(Phase Factor)与规范序(Canonical Ordering)
#

  仅靠轨道占据信息(即位串)无法完全定义一个Slater行列式。由于电子是费米子,交换任意两个电子(即交换产生式算符的顺序)会导致行列式符号反转:\(a_i^\dagger a_j^\dagger = -a_j^\dagger a_i^\dagger\)。因此,产生式算符的排列顺序对行列式的相位因子(通常为+1或-1)至关重要。

  为确保行列式表示的唯一性并正确处理相位,必须采用一个规范序。常见的约定包括:

  1. Alpha自旋块优先:所有 \(\alpha\) 电子的产生式算符 \(a_{i\uparrow}^\dagger\) 排列在所有 \(\beta\) 电子的产生式算符 \(a_{j\downarrow}^\dagger\) 之前。
  2. 块内轨道指数升序:在 \(\alpha\) 和 \(\beta\) 自旋块内部,产生式算符按轨道指数 \(i\) 或 \(j\) 的升序排列。

例如,对于行列式 \(|J\rangle = a_{k\uparrow}^\dagger a_{i\uparrow}^\dagger a_{j\downarrow}^\dagger |\text{vac}\rangle\)(其中 \(i < k\)),根据规范序应重排为 \(a_{i\uparrow}^\dagger a_{k\uparrow}^\dagger a_{j\downarrow}^\dagger |\text{vac}\rangle\)。若从原始形式到规范序需奇数次对换,则 \(|J\rangle = - |I_{\text{canonical}}\rangle\)。

  在实际计算中,通常以一个已知参考行列式(如Hartree-Fock行列式,相位定义为+1)为起点,通过激发操作生成新行列式。新行列式相对于参考态的相位需精确计算并存储,通常用0(+1)或1(-1)表示。可定义如下结构体:

struct DeterminantPhase {
    uint64_t alpha;
    uint64_t beta;
    uint8_t phase;  // 0 = +1相位, 1 = -1相位
};

phase 成员存储行列式 \((\text{alpha}, \text{beta})\) 相对于全局参考态(通常为HF态)的相位,具体计算方法常在辅助函数中实现。

存储与管理的优势
#

  使用 alphabeta 位串对并结合相位因子的存储方式具有以下优势:

  1. 存储紧凑:对于 \(N_{mo} \leq 64\) 的体系,每个行列式的轨道占据信息仅需两个 uint64_t 和一个字节(存储相位),比稀疏矩阵或字符串列表更高效。
  2. 快速比较与查找Determinant 结构体便于实现判等(operator==)和排序(operator<),适用于标准库的有序容器(如 std::set, std::map)或无序容器(如 std::unordered_set, std::unordered_map)。
  3. 高效位运算
    • 电子数统计:通过 __builtin_popcountllstd::popcount(C++20) 快速计算位串中1的个数,即 \(\alpha\) 或 \(\beta\) 电子数。
    • 激发操作:单激发、双激发等可通过位运算(如异或 ^、与 &、或 |)高效实现。例如,从轨道 \(i\) 到 \(p\) 的电子移动可通过 flip_bit(i)flip_bit(p) 操作完成。
    • 差异识别:两个行列式间的差异(激发类型)可通过位串异或操作快速确定。

  综上,将Slater行列式表示为 alphabeta 位串对,并显式存储其相对于规范序参考态的相位因子,是计算量子化学中兼顾存储与运算效率的成熟方案。

Determinant-based Quantum Chemistry - 这篇文章属于一个选集。
§ 1: 本文