4. Psi4 的基础使用¶
在这一节中，我们将会给出一个实际的算例，了解 Psi4 的简单使用；并且通过 Restricted B3LYP 计算，了解最基础的量子化学相关矩阵的输出与使用方式．
4.1. 准备工作¶
4.1.2. 输出文件¶
Psi4 作为量子化学软件，它与普通的 Gaussian、NWChem 一样会有输出文件与暂存 (Scratch) 文件．下面的代码块则是决定了输出文件的名称．暂存文件的路径指定方式也列于下述代码块，但被注释．
In [2]:
# Output file path
# arg 0: Path of output file
# arg 1: Whether overwrite or not
psi4.set_output_file("output.dat", True)
# Scratch directory path
# psi4.core.IOManager.shared_object().set_default_path("/scratch")
提示
默认的输出文件是标准输出流 (Standard Output，有时简称 stdout)．一般来说，使用 Psi4 的 Python 接口是不会用到输出文件；不过有时我们仍然需要这些信息．如果不指定输出文件，通常会找不到输出信息，或者把 Jupyter Notebook 的命令行窗口搞得一团糟．一般来说，最好指定输出文件路径．
暂存文件的文件夹默认路径是
/tmp．如果要在服务器上执行大量大型任务，这个选项需要加上，否则会对主节点硬盘产生压力
(除非每个节点分了一块自己的硬盘给 /tmp)．
4.1.3. 设置内存¶
对于 Psi4，其内存大小通过 Byte 数来确定．如果我们希望分配约 0.5 GB 的空间，我们可以使用下述语句：
In [3]:
psi4.set_memory(int(5e8))
Out[3]:
500000000
4.1.4. 设置分子¶
Psi4 中，分子可以使用直角坐标来确定，也可以用内坐标确定．为了方便，我们始终使用下述的分子进行计算．为了避免使用对称性降低计算量所产生的对轨道与矩阵结构的复杂化，以后统一使用 \(C_1\) 对称性．
In [4]:
mol = psi4.geometry("""
O 0.000000000000 -0.000000000000 -0.079135765807
H 0.000000000000 0.707106781187 0.627971015380
H 0.000000000000 -0.707106781187 0.627971015380
symmetry c1
""")
4.1.5. 设置计算选项¶
Psi4 的选项一般是由全局字典 (Dictionary) 来控制．一般地，在这次教程中会使用下述选项．
对于下述选项，我们可以从 Psi4 的 官方文档关键词 了解详情．不过在最基础的简单应用中，我们不需要仔细检查每个选项，一般下面的选项足够了．其中一些选项的意义是：
scf_type：SCF 形式．pk选项代表使用精确电子积分 (ERI, Electron Repulsion Integral)，因此会使用四脚标双电子积分．另一个常用选项是dk，指使用 Density-fitting 解决复杂电子积分．若不使用复杂的、必须使用四脚标积分的 Post-HF 方法，则dk会是默认选项．详细信息参见 ERI Algorithms．在将来可能会在教程中引入 DF-MP2，届时需要将该选项改为dk．mp2_type: MP2 算法．默认为DF，即使用 Density-fitting．在最近的几个教程中不使用 Density-fitting，因此使用conv即传统算法来计算．dft_spherical_points与dft_radial_points：DFT 格点的球向与径向格点数．尽管可能会有所差别，但 (99, 590) 的格点选取应当与 Gaussian 的Int(UltraFine)非常类似．
In [5]:
psi4.set_options({
'basis': '6-31g',
'scf_type': 'pk',
'mp2_type': 'conv',
'e_convergence': 1e-8,
'd_convergence': 1e-8,
'dft_spherical_points': 590,
'dft_radial_points': 99,
})
4.2. 计算与输出¶
4.2.1. 分子计算与波函数导出¶
在基于轨道的 Post-HF 处理中，最关键的几部分是单电子积分、双电子积分、以及分子轨道系数和本征值．
Psi4 允许在计算后导出波函数对象
(psi4.core.Wavefunction)．波函数对象储存了分子的单电子矩阵、分子轨道、基组与占据轨道、非占轨道数等信息，都是在
Post-HF 计算中不可或缺的部分．它可以由下述的方式导出：
In [6]:
scf_e, scf_wfn = psi4.energy("B3LYP", molecule=mol, return_wfn=True)
通过下述方法，可以导出波函数对象中所储存的矩阵；其中小写 a 与 b 代表的是自旋信息，在 RHF 或 RKS 环境下将会是相同的．(表格可以参见 Psi4NumPy 文档，下表只是列举少数常用的方法)
| 数值量 | 方法 |
|---|---|
| 轨道系数 \(C\) | wfn.Ca(), wfn.Cb() |
| AO 基组密度 \(D\) | wfn.Da(), wfn.Db() |
| Fock 矩阵 \(F\) | wfn.Fa() |
| 基组 | wfn.basisset() |
| \(\alpha\) (\(\beta\)) 电子数 | wfn.nalpha(), wfn.nbeta() |
| 占据轨道数 | wfn.doccpi() |
| 体系能量 | wfn.energy() |
| 轨道能级 \(\boldsymbol{\varepsilon}\) | wfn.epsilon_a(), wfn.epsilon_b() |
不过，实际上上述代码的结果是波函数对象的子类：
In [7]:
type(scf_wfn)
Out[7]:
psi4.core.RHF
因此它还具有更多的方法．在 DFT 计算中，另一个重要的方法是
wfn.V_potential，它是 DFT 格点积分的引擎，在将来会对此做作解释．
4.2.2. Psi4 矩阵与 NumPy 矩阵转换¶
从波函数对象给出的轨道系数与积分等信息确实可以导出，但它并非是立即可用的 NumPy 对象．因此，我们需要了解这者之间的相互转换．
首先是从 Psi4 矩阵转换为 NumPy 矩阵．它可以直接由下述代码得到：
In [8]:
Da_psi4 = scf_wfn.Da()
Da_np = Da_psi4.np
Da_np = np.asarray(Da_psi4) # 与上一行代码等效
print("Da_Psi4 type: ", type(Da_psi4))
print("Da_np type: ", type(Da_np))
Da_Psi4 type: <class 'psi4.core.Matrix'>
Da_np type: <class 'numpy.ndarray'>
提示
从内存利用考虑，一般不使用 numpy.array 来从 Psi4 矩阵构建 NumPy
矩阵，因为这涉及深层复制 (Deep Copy)．
警告
在后面的教程中，可能会出现下述基本正确但偶尔会出错的语句来对 Da_np
进行赋值：
Da_np = scf_wfn.Da().np
Da_np = np.asarray(scf_wfn.Da())
之所以上述语句可能会出错是因为，我们获得的 Da_np 实际上是
scf_wfn.Da() 矩阵的向量视图
(View)，它本身其实没有数据，所有的数据引用自 scf_wfn
中所储存的保护变量 scf_wfn.Da_；我们必须通过方法 scf_wfn.Da()
来调取该矩阵．Da_np 表面上是 scf_wfn.Da_ 的向量视图，但
scf_wfn.Da_ 并不是我们所写的 Python 程序的任何一个变量，因此
Da_np 实际上是一个临时变量的向量视图，这个临时变量也是
scf_wfn.Da_ 的向量视图．而 Python 的垃圾回收机制 (Garbage
Collection)
可能会在内存空间不足时回收这个临时变量，导致 Da_np
被重定向到一个新的临时变量，从而产生未定义的值；也可能原来临时变量中储存了其它信息，导致
Da_np 所指向的值发生变化．
因此，比较稳妥的办法是用一个新的变量，譬如这里的 Da_psi4，储存
scf_wfn.Da()
变量，随后基于这个变量产生向量视图．同时，在程序执行的过程中始终不要重新定义变量
Da_psi4，除非你确定你不打算再使用该值了．
警告
正因为 Da_np 变量是向量视图，所以 Da_np 可以看成是对 AO
基组密度矩阵的引用；或者说，对 Da_np 的值的改变会直接影响
scf_wfn.Da_！特别是 Python
作为脚本语言，代码的安全性通常不能依靠声明变量类型解决；不像 C++
在这里可以声明 Da_np 为底层常量引用
const Matrix * Da_np = &scf_wfn.Da()
从而保护 scf_wfn.Da_
的数据不会被误改变．因此需要注意一般情况下，不应当对从波函数对象中提出的
Psi4
矩阵进行数值上的改变．下面两个代码块只是说明这个情况，但请不要在实际程序中写这类危险代码！
In [9]:
# Da_np 未改变
print(scf_wfn.Da().np[0, 0])
# 改变 Da_np
Da_np[0, 0] = 100.
# 从而底层数据被更改
print(scf_wfn.Da().np[0, 0])
# 如果发生这种更改底层数据，几乎唯一的解决方案是重新算一遍分子
scf_e, scf_wfn = psi4.energy("B3LYP", molecule=mol, return_wfn=True)
print(scf_wfn.Da().np[0, 0])
1.0431804467265933
100.0
1.0431804467265957
有时，我们也会碰到需要从 NumPy 矩阵转换为 Psi4 矩阵的情况．这种转换可以通过 Psi4 矩阵类的初始化直接完成：
In [10]:
mat_np = np.eye(3)
mat_psi4 = psi4.core.Matrix.from_array(mat_np)
print(type(mat_psi4))
print(mat_psi4.np)
<class 'psi4.core.Matrix'>
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
4.2.3. 积分引擎¶
Psi4 允许在分子计算前获得单电子与双电子积分．这些积分时通过 Psi4
的积分引擎 psi4.core.MintsHelper
给出的．通过下述语句，传入基组的信息，我们可以初始化一个积分引擎：
In [11]:
mints = psi4.core.MintsHelper(scf_wfn.basisset())
通过下述的方法，可以获得计算中所必要的电子积分．这些积分不需要经过实际的分子体系计算就能给出，是原子轨道基组下的积分．返回值是
psi4.core.Matrix 对象，因此需要转为 NumPy 对象再进行分析．
| 数值量 | 方法 |
|---|---|
| 原子轨道重叠积分 | mints.ao_overlap() |
| 动能积分 | mints.ao_kinetic() |
| 势能积分 | mints.ao_potential() |
| 电子排斥积分 | mints.ao_eri() |
| 偶极积分 | mints.ao_dipole() |
波函数对象也会给出单电子积分，我们可以简单验证一下这两者的积分是否相同：
In [12]:
# 动能积分 + 势能积分 = 单电子积分
np.allclose(mints.ao_kinetic().np + mints.ao_potential().np, scf_wfn.H().np)
Out[12]:
True
In [13]:
# 原子轨道重叠积分
np.allclose(mints.ao_overlap().np, scf_wfn.S().np)
Out[13]:
True
电子排斥积分 (ERI, Electron Repulsion Integral，有时称双电子积分) 按照化学惯例给予角标 \((\mu \nu \vert \kappa \tau)\)，且使用的是实空间轨道．因此，其具有八重的对称性．下面的代码会简单验证一下电子排斥积分的一些性质，同时可以简单了解一下张量的角标转换的方式．
In [14]:
eri_psi4 = mints.ao_eri()
eri = np.asarray(eri_psi4)
# 比较 ERI 的矩阵形状与基组大小
print("Shape of ERI tensor: ", eri.shape)
print("Shape of basis number: ", scf_wfn.nmo())
# ERI 积分的对称性
print("(uv|kl) = (vu|kl): ", np.allclose(eri, eri.transpose(1,0,2,3)))
print("(uv|kl) = (kl|uv): ", np.allclose(eri, eri.transpose(2,3,0,1)))
print("(uv|kl) = (uk|vl): ", np.allclose(eri, eri.transpose(0,2,1,3)))
# (uv|kl) -> (vu|kl) 简单的转置用 swapaxes 方便
print("transpose vs. swapaxes (uv|kl) -> (vu|kl): ", np.allclose(
eri.transpose(1,0,2,3), eri.swapaxes(0,1)
))
# (uv|kl) -> (kl|uv) 复杂的转置用 transpose 方便
print("transpose vs. swapaxes (uv|kl) -> (kl|uv): ", np.allclose(
eri.transpose(2,3,0,1), eri.swapaxes(0,2).swapaxes(1,3)
))
Shape of ERI tensor: (13, 13, 13, 13)
Shape of basis number: 13
(uv|kl) = (vu|kl): True
(uv|kl) = (kl|uv): True
(uv|kl) = (uk|vl): False
transpose vs. swapaxes (uv|kl) -> (vu|kl): True
transpose vs. swapaxes (uv|kl) -> (kl|uv): True
警告
电子排斥积分将会占用较大的资源．如果基组较大，就需要手动验证内存大小．其计算方法是，在知晓基组大小时，可以通过下面的代码块进行双电子积分内存大小预估．
如果计算过程中出现内存不足，一般 NumPy 会先报错；但也可能产生未定义的错误甚至系统崩溃！
In [15]:
nbf = scf_wfn.nmo()
I_size = (nbf**4) * 8.e-9
print("Size of the ERI tensor will be {:10.6f} GB.".format(I_size))
Size of the ERI tensor will be 0.000228 GB.
提示
一般地，基组数与分子轨道数没有必然的相等关系，且基组数不小于分子轨道数．然而在基组较小，或没有明显基组的线性依赖关系时，一般的量子化学软件都会使用与基组数相同的数来确定分子轨道数．这也是为何这里可以使用分子轨道数来表示原子积分的维度的原因．之后的教程中将始终认为基组数与分子轨道数的大小相同．