进行QM计算时，软件在做什么？

Posted Oct 8, 2025 Updated Oct 8, 2025

By Bane Dysta

15 min read

进行QM计算时，软件在做什么？

理论背景

量子化学的核心目标是求解分子的薛定谔方程:

\[\hat{H}\psi = E\psi\]

对于一个包含 $N(N>1)$ 个电子的分子，哈密顿算符 $\hat{H}$ 包含电子动能、电子-核吸引能、电子-电子排斥能。但这个多体问题太复杂，无法精确求解。

Hartree-Fock方法

为避免直接求解多电子波函数，一个经典的思路是Hartree-Fock方法。HF的核心思想是平均场近似，即每个电子在其他所有电子产生的平均场中运动，如此我们可以用单电子波函数来近似描述多电子体系。在HF框架下，假设多电子波函数可以用单个Slater行列式表示：

\[\psi = \frac{1}{\sqrt{N!}}|\phi_1\phi_2\cdots\phi_N|\]

其中 $\phi_i$ 是分子轨道（MO），则待求解目标HF方程：

\[\hat{F}\phi_i = \varepsilon_i\phi_i\]

其中 $\hat{F}$ 是 Fock 算符，$\varepsilon_i$ 是轨道能量。

解决了多体问题，又产生了新的问题：Fock 算符本身依赖于我们要求解的分子轨道！在实际计算中，我们无法预知$\phi_i$的解析形式，这里需要用一种简单合理的数学形式来描述这个$\phi_i$。

基组

为了将抽象的$\phi_i$转化为可以看得到~~但看不懂~~的数字，通常会采用用一组基函数的线性组合来近似描述$\phi_i$：

\[\phi_i = \sum_{\mu=1} C_{\mu i}\chi_\mu\]

这里 $\chi_\mu$ 是基函数，$C_{\mu i}$ 是展开系数。基函数的集合称为基组（basis set）。

以最经典的STO-3G为例。STO代表Slater Type Orbital，其形式为:

\[\chi(r) = Nr^{n-1}e^{-\zeta r}Y(\theta,\phi)\]

但 STO 在计算双电子积分时非常困难。因此，STO-3G 用 3 个高斯型函数（GTO）的固定线性组合来拟合一个 STO:

\[\chi_{\text{STO}} \approx \sum_{i=1}^{3} d_i \chi_{\text{GTO},i}\]

对于 s 型轨道，每个高斯函数形式为：

\[\chi_{\text{GTO}} = Ne^{-\alpha r^2}\]

以氢原子 1s 轨道为例，STO-3G 使用以下三个高斯函数:

指数 $\alpha$	收缩系数 $d$
3.425250	0.154329
0.623914	0.535328
0.168855	0.444635

这意味着 H 原子的 1s 基函数实际上是:

\[\chi_{1s} = 0.154329 \times e^{-3.42525r^2} + 0.535328 \times e^{-0.623914r^2} + 0.444635 \times e^{-0.168855r^2}\]

又例如对于水分子 $\ce{H2O}$，使用 STO-3G 基组时:

O 原子：$1s, 2s, 2p_x, 2p_y, 2p_z$（5 个基函数）
H 原子：$1s$（每个 H 1 个基函数）

总共 7 个基函数。

SCF 迭代

将基函数展开式代入 HF 方程，左乘 $\langle\chi_\nu|$，得到： $\sum_\mu F_{\nu\mu}C_{\mu i} = \varepsilon_i\sum_\mu S_{\nu\mu}C_{\mu i}$

写成矩阵形式就是Roothaan 方程： $\mathbf{FC} = \mathbf{SC\varepsilon}$

其中:

①$\mathbf{F}$：Fock 矩阵（$K \times K$ 维，$K$ 是基函数数目）
\[F_{\mu\nu} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma}P_{\lambda\sigma}\left[(\mu\nu|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma)\right]\]
②$\mathbf{C}$：系数矩阵（每列是一个 MO 的展开系数）
- 每一列 $\mathbf{C}_i$ 是第 $i$ 个 MO 的展开系数
- 满足正交归一条件：$\mathbf{C}^\dagger\mathbf{SC} = \mathbf{I}$
③$\mathbf{S}$：重叠矩阵（$S_{ij} = \langle\chi_i \chi_j\rangle$）
\[S_{\mu\nu} = \langle\chi_\mu|\chi_\nu\rangle = \int\chi_\mu^*(\mathbf{r})\chi_\nu(\mathbf{r})d\mathbf{r}\]
- 如果基函数正交，$\mathbf{S} = \mathbf{I}$（单位矩阵）
- 实际基函数（如高斯函数）通常不正交，$S_{\mu\nu} \neq 0$
④$\boldsymbol{\varepsilon}$：轨道能量对角矩阵
对角矩阵，对角元是轨道能量 $\varepsilon_1, \varepsilon_2, \ldots$

Fock 矩阵的构建

如上所述，Fock 矩阵元素为:

\[F_{\mu\nu} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma} P_{\lambda\sigma}\left[(\mu\nu|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma)\right]\]

这里包含三部分:

1. 核心哈密顿矩阵 $H_{\mu\nu}^{\text{core}}$:

动能积分：$\langle\chi_\mu -\frac{1}{2}\nabla^2 \chi_\nu\rangle$
电子-核吸引：$\langle\chi_\mu -\sum_A\frac{Z_A}{r_A} \chi_\nu\rangle$

2. 密度矩阵 $\mathbf{P}$:

\[P_{\lambda\sigma} = 2\sum_{i}^{\text{occ}} c_{\lambda i}c_{\sigma i}\]

（注：上式适用于RHF）

3. 双电子积分 $(\mu\nu

\lambda\sigma)$:

\[(\mu\nu|\lambda\sigma) = \int\int\chi_\mu(\mathbf{r}_1)\chi_\nu(\mathbf{r}_1) \frac{1}{r_{12}} \chi_\lambda(\mathbf{r}_2)\chi_\sigma(\mathbf{r}_2) d\mathbf{r}_1d\mathbf{r}_2\]

这是最耗时的部分！对于 $K$ 个基函数，将计算 $K^4/8$ 个这样的积分（1/8来自双电子积分的对称性）。对于一个给定的分子和一套选定的基组，每个双电子积分 $(\mu\nu

\lambda\sigma)$ 的值都是唯一确定的，因此一些QC程序会选择将双电子积分存储下来，供后续反复读取使用。然而，若所用缓存目录挂载在HDD上，很可能读的速度还不如重新计算快，这也是QC计算对SSD需求强的原因。再加上双电子积分数量庞大，占用空间惊人，不少QC程序干脆不再存储双电子积分(direct SCF)，每次要用都重新计算。(PS:如果内存够大，其实把双电子积分全都存在内存中是最好的选择，但对于普通机子来说内存肯定存不下完整的双电子积分，顶多存一些比较重要的)

不难看出，为构建Fock矩阵，最关键的是得到密度矩阵$\mathbf{P}$。然而，这里有个”鸡生蛋”的问题：

构建 Fock 矩阵需要密度矩阵 $\mathbf{P}$
密度矩阵 $\mathbf{P}$ 需要从系数矩阵 $\mathbf{C}$ 计算
系数矩阵 $\mathbf{C}$ 需要通过求解 Fock 矩阵得到

初猜构建

为解决上述问题，我们需要一个$\mathbf{P}$的初猜，即初始密度矩阵 $\mathbf{P}^{(0)}$来启动整个SCF迭代。在 SCF 迭代开始前，我们还不知道电子分布，若姑且假设没有电子-电子相互作用，此时 Fock 矩阵简化为：

\[\mathbf{F}^{(0)} = \mathbf{H}^{\text{core}}\]

其中核心哈密顿矩阵：

\[H_{\mu\nu}^{\text{core}} = \langle\chi_\mu|-\frac{1}{2}\nabla^2|\chi_\nu\rangle + \langle\chi_\mu|-\sum_A\frac{Z_A}{r_A}|\chi_\nu\rangle\]

这种初猜方法是比较经典的Hcore初猜，相当于只考虑核的吸引，靠后续SCF迭代逐渐引入电子排斥效应。我们以Hcore初猜为例，介绍初猜构建流程(Hcore方法比较粗糙，并非很好的选择，现代QM软件一般都有更高级的初猜算法，如SAD等)。

考虑键长为0.74 Å的$\ce{H2}$ 分子，$\chi_1$为左侧 H 的 1s 基函数（位于 $z = -0.37$ Å），$\chi_2$为右侧 H 的 1s 基函数（位于 $z = +0.37$ Å）。以下数据使用PYSCF计算：

# -*- coding: utf-8 -*-
"""
A concise PySCF script to demonstrate a Restricted Hartree-Fock (RHF)
calculation for the H2 molecule using a manually constructed Hcore initial guess.
"""
import numpy as np
from pyscf import gto, scf
from scipy import linalg

# --- 1. Define the molecule ---
mol = gto.Mole()
mol.atom = [
    ['H', (0, 0, -0.37)],
    ['H', (0, 0,  0.37)],
]
mol.basis = 'sto-3g'
mol.build()

# --- 2. Manually construct the Hcore initial guess ---
# Get overlap (S) and core Hamiltonian (H_core) matrices
S = mol.get_ovlp()
H_core = mol.get_hcore()

# Solve the generalized eigenvalue problem: H_core * C = S * C * E; This gives the initial guess for MO energies (E_0) and coefficients (C_0)
E_0, C_0 = linalg.eigh(H_core, S)

# Build the initial density matrix (P_0) from the occupied orbitals; For H2, there are 2 electrons, so we occupy the lowest-energy orbital.
occupied_mo_coeffs = C_0[:, :mol.nelectron // 2]
P_0 = 2 * np.dot(occupied_mo_coeffs, occupied_mo_coeffs.T)


# --- 3. Run the SCF calculation ---
# Initialize the RHF method
mf = scf.RHF(mol)

# Run the SCF cycle, providing our manually constructed density matrix as the start; The mf.kernel() function handles the entire iteration process.
mf.kernel(dm0=P_0)


# --- 4. Print the final, converged results ---
print(f"Calculation converged: {mf.converged}")
print(f"Final RHF Total Energy: {mf.e_tot:.10f} Hartree")
print("\nFinal Orbital Energies (Hartree):")
print(mf.mo_energy)
print("\nFinal MO Coefficients:")
print(mf.mo_coeff)
print("\nNote: The sign/phase of MO coefficient columns is arbitrary.")

首先，计算重叠矩阵 $\mathbf{S}$：

\[S_{\mu\nu} = \langle\chi_\mu|\chi_\nu\rangle\]

对于 H₂ 分子，计算得到的值为：

\[\mathbf{S} = \begin{pmatrix} 1.000 & 0.65987 \\ 0.65987 & 1.000 \end{pmatrix}\]

由于每个基函数自己与自己的重叠为100%，对对角元进行归一化，此时非对角元素体现了两个H原子基函数的重叠程度。以此为基础计算哈密顿量 $\mathbf{H}^{\text{core}}$：

\[\mathbf{H}^{\text{core}} = \begin{pmatrix} -1.12096 & -0.95938 \\ -0.95938 & -1.12096 \end{pmatrix}\]

对角元素$H_{11}^{\text{core}}$为基函数$\chi_1$局域在原子核1上的动能项加上两个核对该基函数的吸引项。最终值为负，说明体现吸引作用，即电子在原子1上是”能量有利的”。同理，电子在原子2上也是”能量有利的”。
非对角元素描述电子从原子1的轨道跳到原子2的轨道的耦合强度，包含动能耦合项和电子在原子1与原子2之间时感受到两个核的吸引积分，其值为负体现了成键作用。

现在，我们可以用 $\mathbf{H}^{\text{core}}$ 替代 Fock 矩阵，求解广义特征值问题：

\[\mathbf{H}^{\text{core}} \mathbf{C}^{(0)} = \mathbf{S} \mathbf{C}^{(0)} \boldsymbol{\varepsilon}^{(0)}\]

得到初猜系数矩阵（经过符号调整以方便观察）：

\[\mathbf{C}^{(0)} = \begin{pmatrix} 0.54884 & -1.21245 \\ 0.54884 & 1.21245 \end{pmatrix}\]

和初猜轨道能量：

\[\boldsymbol{\varepsilon}^{(0)} = \begin{pmatrix} -1.25331 & \\ & -0.47507 \end{pmatrix}\]

这给出了两个初猜分子轨道：

MO₁（成键，$\varepsilon = -1.25331$ Ha）：$\phi_1^{(0)} = 0.54884\chi_1 + 0.54884\chi_2$
MO₂（反键，$\varepsilon = -0.47507$ Ha）：$\phi_2^{(0)} = -1.21245\chi_1 + 1.21245\chi_2$

构建初猜密度矩阵 $\mathbf{P}^{(0)}$：

\[P_{\mu\nu}^{(0)} = 2\sum_{i}^{\text{occ}} c_{\mu i}^{(0)} c_{\nu i}^{(0)}\]

H₂ 有 2 个电子，占据能量最低的 MO₁：

\[\mathbf{P}^{(0)} = 2 \times \begin{pmatrix} 0.54884 \\ 0.54884 \end{pmatrix} \begin{pmatrix} 0.54884 & 0.54884 \end{pmatrix} = \begin{pmatrix} 0.60246 & 0.60246 \\ 0.60246 & 0.60246 \end{pmatrix}\]

有了 $\mathbf{P}^{(0)}$，就可以构建第一次真正的 Fock 矩阵 $\mathbf{F}^{(1)}$：

\[F_{\mu\nu}^{(1)} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma} P_{\lambda\sigma}^{(0)}\left[(\mu\nu|\lambda\sigma) - \frac{1}{2}(\mu\lambda|\nu\sigma)\right]\]

SCF 迭代循环

进入迭代循环:

flowchart TD
    Start([开始scf循环(n=0)]) --> Step1[用 P<sup>n</sup> 计算 Fock 矩阵 F<sup>n</sup>]
    
    Step1 --> Step2[计算双电子积分贡献<br/>最耗时!]
    
    Step2 --> Step3[求解 Roothaan 方程<br/>得到新的系数矩阵 C<sup>n+1</sup>]
    
    Step3 --> Step4[用 C<sup>n+1</sup> 的占据轨道构建新密度矩阵 P<sup>n+1</sup>]

    Step4 --> Step5{"能量变化:<br/>ΔE < 10<sup>-8</sup><br/>且<br/>密度变化:<br/>ΔP < 10<sup>-10</sup>"}
    
    Step5 -->|是| End([计算结束])
    Step5 -->|否, n=n+1| Start
    
    style Step2 fill:#ffcccc
    style End fill:#ccffcc
    style Start fill:#cce5ff

计算结果

SCF 收敛后，软件得到的是什么？本质上是系数矩阵 $\mathbf{C}$(当然，还有能量，不过要在前述基础上加上核排斥项)。每个分子轨道都是基函数的线性组合，可以告诉我们每个分子轨道由哪些原子轨道以多大比例混合而成。

$\ce{H2}$ 分子示例

让我们以 $\ce{H2}$ 分子（键长 0.74 Å，STO-3G）为具体案例:

基函数设置:

$\chi_1$：$\ce{H}$（左）的 1s 基函数
$\chi_2$：$\ce{H}$（右）的 1s 基函数

收敛后的系数矩阵与轨道能量:

	$\text{MO}_1(\sigma)$	$\text{MO}_2(\sigma^*)$
$\chi_1$	0.54884	-1.21245
$\chi_2$	0.54884	1.21245
$\varepsilon$	-0.57855 Hartree	0.67114 Hartree

成键轨道（$\text{MO}_1$，占据）:

\[\phi_\sigma = 0.54884\chi_1 + 0.54884\chi_2\]

展开具体形式:

\[\begin{aligned} \phi_\sigma(\mathbf{r}) &= 0.54884 \times \left[0.154e^{-3.425r_1^2} + 0.535e^{-0.624r_1^2} + 0.445e^{-0.169r_1^2}\right] \\ &\quad + 0.54884 \times \left[0.154e^{-3.425r_2^2} + 0.535e^{-0.624r_2^2} + 0.445e^{-0.169r_2^2}\right] \end{aligned}\]

其中 $r_1$ 是到左侧 H 的距离，$r_2$ 是到右侧 H 的距离。

反键轨道（$\text{MO}_2$，未占据）:

\[\phi_{\sigma^*} = -1.21245\chi_1 + 1.21245\chi_2\]

可视化时，软件可以读取这些系数，在空间格点上计算:

\[\phi_i(x,y,z) = \sum_\mu c_{\mu i} \chi_\mu(x,y,z)\]

然后绘制等值面，形成我们看到的分子轨道图。

总结

量子化学软件计算时，本质上是在求解一个大型的非线性特征值问题。通过迭代过程，软件最终得到一组展开系数 $C_{i \mu}$，可以组合计算时所用基组来实现对分子波函数的描述。

Quantum Chemistry

This post is licensed under CC BY 4.0 by the author.