All Hail Scarlet!

0%

A Minimal Introduction to Quantum Chemistry

Posted on Edited on
这篇文章介绍了量子化学和计算化学的一些基本知识，并说明了如何计算氢原子的基态能量。

在这篇文章中将介绍如何计算 \mathrm{H} 的基态能量。

哈密顿量 Hamiltonian

量子力学将所有的测量、变换和演化都用算符进行表示。整个系统的状态由一个右矢向量 |\alpha\rangle 表示；测量将右矢转换为测量操作的本征右矢（用矩阵乘法表示的话本征右矢也即特征向量）；变换即是对右矢的变换；演化即是右矢对时间的变化。

在推导并定义时间演化的算符的时候，借用了哈密顿量在经典力学中是时间演化生成元的概念，将哈密顿量的概念引入了时间演化。

哈密顿量究竟如何影响一个系统状态的时间演化在这里不作推导，有兴趣的读者可以自己参考量子力学相关书籍。我们直接给出如下结果：

  • 哈密顿量是与体系能量关联的算符。
  • 若哈密顿量不含时间，即随不时间变化而直接变化，那么哈密顿量代表了体系的总能量，哈密顿量的特征值的可能取值就是体系能量的可能取值。
  • 对粒子（分子/原子）的哈密顿量 \hat{H} 和粒子的波函数 |\psi\rangle ，其对应于薛定谔方程 \hat{H}|\psi\rangle = E|\psi\rangle ，其中 E 也即哈密顿量算符的本征值，波函数为一个本征函数。

因而，为了在量子力学的系统中计算 \mathrm{H} 的基态能量，我们需要求出其哈密顿量，并求解其特征值。

变分法 Variational Method

我们在方程 \hat{H}|\psi\rangle = E|\psi\rangle 两端同时乘上 |\psi\rangle 的共轭转置 \langle\psi| 得到
\langle\psi|\hat{H}|\psi\rangle = E\langle\psi|\psi\rangle
波函数代表粒子出现在某处的概率分布，因而满足总概率为 1 的归一化条件 \langle\psi|\psi\rangle = 1 。从而我们得到
\langle\psi|\hat{H}|\psi\rangle = E
变分原理(The variational principle)，也即最小能量原理，描述了如下事实：
如果我们用非能量最低的空间概率分布函数 |\psi\rangle 带入上式计算，得到的能量值 E 将会比实际基态情况要高。

因而，我们要求能量最低状态的 |\psi\rangle 和 E ，只需求 |\psi\rangle 使得上式中 E 最小。

氢原子的基态能量

上面的内容可能看起来比较抽象。为了进一步了解具体的计算方式，我们将具体介绍如何计算氢原子的基态能量（电子的结合能）。

首先我们需要一个具体表示能量的单位。为了避免在推导过程中引入一系列例如 \hslash 的常数，我们采用原子单位制(Atomic Units)，并略去其推导过程：

  • 长度 - 玻尔半径 1a_0 = 5.2918 \times 10^{-11} m
  • 质量 - 电子质量 1m_e = 9.1095 \times 10^{-31} kg
  • 电量 - 元电荷 1e = 1.6022 \times 10^{-19} C
  • 能量 - 哈特里 1\epsilon_a = 2625.5\ kJ/mol = 27.2eV

（注解：实际上这些单位是通过求解对应问题的时候将相关常数提取出来而得到，而长度和能量单位都和氢原子脱不开干系。可能从严谨的角度上来说存在一些问题。）

我们需要写出氢原子的哈密顿量。氢原子由原子核和电子组成。对于结合能而言，我们无需考虑原子核和电子单独的性质。

哈密顿量是关于能量的算符。氢原子中，电子围绕原子核转动具有动能；电子处在原子核的电场中具有势能。

我们有
H = E_k + E_p

E_k = \frac{1}{2} mv^2 = \frac{1}{2} m_e \times -\nabla_r^2 = -\frac{1}{2} m_e \times (\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2})

E_p = -\frac{kqQ}{r}

在上述单位制中书写得到
H = -\frac{1}{2} (\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}) - \frac{1}{r}
其中
r = \sqrt{x^2 + y^2 + z^2}

对上述 H ，我们要求 |\psi\rangle = |\psi(x, y, z)\rangle 使得 \langle\psi|H|\psi\rangle 最小。

那么， |\psi\rangle 到底是什么？在这里，它表示电子的概率分布。如果将它看作向量以与“右矢”的表现形式符合，那么它所处的空间维数是无穷大，每一个基向量代表着空间中的一个微小区域。
如果函数 f(x, y, z) 是电子的概率分布函数，那么 |\psi\rangle = \iiint_{\mathbb{R}^3} f(\mathrm{x}) |\mathrm{x}\rangle d^3\mathrm{x} ，其中 |\mathrm{x}\rangle 代表点 \mathrm{x} 附近的一个微小体积所对应的基向量。

为了方便计算，我们使用一系列高斯分布函数来拟合 |\psi\rangle 。氢原子中电子的波函数可以求解析解，但是其方法不能扩展到多电子情形因而不在此作介绍。

每个高斯分布函数 \phi_i 具有如下形式（归一化条件：平方在整个空间上的积分值是1）：
\phi_i = \left(\frac{2\alpha_i}{\pi}\right)^\frac{3}{4} e^{-\alpha_i r^2}

在z轴方向求导（x、y轴方向类似），我们有
\frac{\partial^2}{\partial z^2} \phi_i = \left(\frac{2\alpha_i}{\pi}\right)^\frac{3}{4} (e^{-\alpha_i (x^2 + y^2)}) (2\alpha_i) (e^{-\alpha_i z^2}) (2\alpha_i z^2 - 1)

拟合的过程如下
\psi = \sum k_i\phi_i

我们先使用一个高斯函数进行拟合求解。我们使用Nelder-Mead方法进行求解。

另外，在对高斯函数进行积分时进行了 5\sigma 截断。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import math
import numpy
import scipy.optimize as opt
import scipy.integrate as iint

def dst(x, y, z):
    return math.sqrt(x * x + y * y + z * z)

def gauss(a, r):
    return math.exp(-r * r * a)

def gauss_partial(a, var, fix1, fix2):
    var2 = var * var
    vexp = math.exp(-a * var2)
    c = 2 * a
    return math.exp(-a * (fix1 ** 2 + fix2 ** 2)) * c * vexp * (c * var2 - 1)

def gauss_d2(a, x, y, z):
    return (gauss_partial(a, x, y, z)
            + gauss_partial(a, y, z, x)
            + gauss_partial(a, z, x, y))

def hamilt_g(a, x, y, z, r):
    return -0.5 * gauss_d2(a, x, y, z) - gauss(a, r) / r

def variation(vec_x):
    opts = {'epsabs': 3e-4, 'epsrel': 3e-4}
    for kx in vec_x:
        if kx <= 1e-3:
            return 1e9
    a = vec_x[0]
    coe = (2 * a / math.pi) ** 1.5

    def _func_int(x, y, z):
        r = dst(x, y, z)
        left_sum = gauss(a, r)
        right_sum = hamilt_g(a, x, y, z, r)
        return right_sum * left_sum * coe

    L = 3e-4
    H = 5 / math.sqrt(a)
    return iint.nquad(_func_int, [(L, H)] * 3, opts=opts)[0] * 8

result = opt.minimize(variation,
                      numpy.array([0.3]),
                      tol=1e-3,
                      method='Nelder-Mead',
                      options={"adaptive": True})
print("Alpha:", result.x)
print("Energy (Hartrees):", variation(result.x))

运行结果：

1
2
3
Alpha: [0.2821875]
Energy (Hartrees): -0.42383029192249166
Wall time: 2.76 s

由上面的单位换算和基本实验数据可以计算出，氢原子基态能量是-0.5哈特里。这一结果的差异是由于一个高斯函数不能够很好地拟合氢原子外电子的概率分布。

下面我们用两个高斯函数进行拟合。作了一些trivial的调整：归一化由最后除以 \langle\psi|\psi\rangle 得到；积分截断到了 3\sigma 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import math
import numpy
import scipy.optimize as opt
import scipy.integrate as iint

def dst(x, y, z):
    return math.sqrt(x * x + y * y + z * z)

def gauss(a, r):
    return math.exp(-r * r * a)

def gauss_partial(a, var, fix1, fix2):
    var2 = var * var
    vexp = math.exp(-a * var2)
    c = 2 * a
    return math.exp(-a * (fix1 ** 2 + fix2 ** 2)) * c * vexp * (c * var2 - 1)

def gauss_d2(a, x, y, z):
    return (gauss_partial(a, x, y, z)
            + gauss_partial(a, y, z, x)
            + gauss_partial(a, z, x, y))

def hamilt_g(a, x, y, z, r):
    return -0.5 * gauss_d2(a, x, y, z) - gauss(a, r) / r

def variation(vec_x):
    opts = {'epsabs': 3e-4, 'epsrel': 3e-4}
    for kx in vec_x:
        if kx <= 1e-3:
            return 1e9
    a1 = vec_x[0]
    a2 = vec_x[1]
    k1 = vec_x[2]
    k2 = vec_x[3]

    def _func_int(x, y, z):
        r = dst(x, y, z)
        left_sum = k1 * gauss(a1, r) + k2 * gauss(a2, r)
        right_sum = k1 * hamilt_g(a1, x, y, z, r) + k2 * hamilt_g(a2, x, y, z, r)
        return right_sum * left_sum

    def _func_norm(x, y, z):
        r = dst(x, y, z)
        left_sum = k1 * gauss(a1, r) + k2 * gauss(a2, r)
        return left_sum * left_sum

    L = 3e-4
    H = 3 / math.sqrt(min(a1, a2))
    norm = iint.nquad(_func_norm, [(L, H)] * 3)[0]
    return iint.nquad(_func_int, [(L, H)] * 3, opts=opts)[0] / norm

result = opt.minimize(variation,
                      numpy.array([0.3] * 4),
                      tol=1e-3,
                      method='Nelder-Mead',
                      options={"adaptive": True})
print("Alpha:", result.x)
print("Energy (Hartrees):", variation(result.x))

运行结果：

1
2
3
Alpha: [1.29184486 0.19795375 0.02100325 0.014888  ]
Energy (Hartrees): -0.48554126080356996
Wall time: 4min 52s

此时的拟合结果已经和精确值相去不远。事实上，计算时常常用三个高斯函数的线性组合来拟合氢原子中电子的波函数，已经可以达到相当的精度。用三个高斯函数近似得到的电子轨道在进一步的计算中常被使用，称为STO-3G Basis。

用两个函数去拟合比一个函数耗时大了许多。不过这只是最朴素的做法。首先，容易证明在这个归一化方法下参数k2是不必要的。如果愿意做一些小小的分析，k2可以通过解一元二次方程得到、而无需为归一化进行一次额外的积分。另外，容易看出，由于这个体积分的被积函数只和坐标到原点的距离有关，如果换元到球面坐标系将会更加容易计算。具体计算过程不在此展示。通过这样的方法可以在可以接受的时间范围内得到6个甚至更多高斯函数线性组合时的最优状态。

附：变分原理证明

考察对任意的向量 |\psi\rangle ， \frac{\langle\psi|H|\psi\rangle}{\langle\psi|\psi\rangle} 的取值范围。

用通常矩阵的语言，即考察对任意的向量 X ， \frac{X^\dagger H X}{X^\dagger X} 的取值范围。

根据量子力学的规定，哈密顿量 \hat{H} 是厄米的。写成矩阵形式即
H = H^\dagger

我们知道，这样的矩阵：一、特征值都是实数；二、可以正交相似于对角阵。

设
H = U^\dagger VU
其中
U^\dagger = U^{-1}

那么
\frac{X^\dagger H X}{X^\dagger X} = \frac{X^\dagger U^\dagger VUX}{X^\dagger (U^\dagger U) X} = \frac{(UX)^\dagger V (UX)}{(UX)^\dagger (UX)}

由于正交变换 X \to UX 是保长度且可逆的，不难证明取遍这里的 UX 时，原向量 X 也将被取遍，因而只需考虑 UX 取遍整个向量空间的情况。

设
UX = (y_1, y_2, …, y_n)

经整理可将上式化为
\frac{\lambda_1 |y_1|^2 + \lambda_2 |y_2|^2 + … + \lambda_n |y_n|^2}{|y_1|^2 + |y_2|^2 + … + |y_n|^2}

其中 \lambda_i 为 H 的各个特征值。而显然，上式
\geq \lambda_{min}

上述证明过程在极限意义下仍然成立，因而在无限维的情况依然适用。即证得了离散和连续情况下的变分原理。