xDH 在 Psi4 下的简易教程¶
在这份文档中，我们将会介绍 xDH 型函数的典型：XYG3，其在 Psi4 下的实现过程。作者希望，读者可以借助于量子化学软件 Psi4 所提供的 Python API 所得到的电子积分与底层函数，以及 NumPy 对矩阵、张量计算的强大支持，可以较为轻松地在 Post-HF 或 DFT 方法理论推导，与上机实验的结果进行相互比对，并最终享受亲手实现新方法的乐趣。
动机¶
在量子化学领域中，计算方法的开发瓶颈除了最初的灵感与公式推导之外，表达式的程序实现与其正确性也非常重要。事实上，往往后者左右了方法的开发效率、完成度，以及方法的推广流传程度。通常，应用一个量子化学方法或者会在已有程序上作改动，或者自己重新构建一个程序。
对于程序改动，方法开发者通常需要面临理解复杂程序结构、阅读大量冗余代码、了解编译流程、缺少充足的文档等困难，在程序上所消耗的时间不亚于、甚至远高于方法开发的耗时。成功的程序改动决定于方法开发者的程序调试能力。
而对于从头构建程序，耗时将会因人而异，这取决于方法开发者的综合程序能力，包括编译、调试、框架搭建、底层工具设计、公式的程序表达等等。同时，由于大多数方法开发者只关心实现而不关心代码质量，因此这类程序的复用性较差，通常不会有充足的文档，也不适合进行稍复杂的改动。在量子化学发展的早期，正因为程序的开发难以跟上量化方法的发展，同时课题组的交流与合作受限于没有因特网，从而催生了数量众多的量化软件或程序。
这份文档更为偏向于从头构建程序。早期的计算机发展主要受限于计算效率与内存，而即使是很小的量子化学的任务，也对这两者有较高的消耗。然而，现在的微机已经可以快速处理中小体系的量子化学计算；从方法开发的角度上看，硬件不应当成为瓶颈——方法的开发只需要较小的模型体系通常即可。而对硬件要求的降低，催生了高级语言的流行。对于 Python，它对方法开发者的限制会显著减小：不存在编译问题，调试可以通过实时交互界面 (interactive interface，例如 Python consoles 或 Jupyter Notebook) 或 IDE 完成。对于底层工具，Python 社区提供众多的库函数可以胜任；在量子化学领域，底层工具绝大部分与矩阵计算有关，因此 NumPy 可以胜任。
对于量子化学中的 Post-HF 与 DFT 领域，剩下的问题则是获得量子化学积分、量子化学相关的底层工具、程序框架设计与复用性、以及公式与程序的对应。对于前两者，现有的软件，例如 Psi4 与 PySCF 均提供了不少对象接口 (API)，对它们的活用将会简化程序书写上的困难。对于第三个问题，我们暂且不考虑。
因此，也许看起来从头构建程序是很可怕的事情，但在现在众多程序便利的环境下，方法开发者所真正面临的问题已经可以回归到最为本质的问题，即如何将公式写成程序。在这份教程中，一个目标便是对于我们需要计算的公式中的大多数项 (譬如张量乘积)，尽可能在 1 行代码内解释，至少一般可以在 5 行以内的代码块内说明。并且，我们不希望涉及非必要的算法细节；整个教程中，绝大多数代码开始将不会多于 3 个 Tab，尽少使用循环与判断，让程序的书写的目的回归于公式表达本身。
这份教程的很多内容借鉴或直接从 Psi4 的一份官方简易实现 Psi4NumPy 获得。后者包含了许多高级 Post-HF 方法简单但有效的实现。
环境搭建¶
Anaconda 环境¶
这份教程的所有文档、计算都通过 Python 实现；因此需要安装一个 Python 环境。Anaconda 作为 Python 的一种非官方发行版，它集成了众多科学计算中所必须与经常使用的库；它至少包含可以实现矩阵计算的 NumPy、绘图实现 MatPlotLib、交互笔记本 Jupyter、仿 Matlab IDE 的 Spyder 等库与工具，以及文档工具 Pandoc、Python 管理工具 conda 和 pip 等管理工具。Anaconda 在普遍的科学计算领域足够完备，同时库的依赖关系可以通过 pip 与 conda 等方便地管理，使得程序员不必耗费许多精力在准备环境上。因此它也是普遍推荐的 Python 安装的解决方案。
由于 Psi4 不具备 Windows 版本，因此这里介绍 Linux 下的安装。
注意
WSL (Windows Subsystem of Linux) 尽管也是良好的 Linux 环境，但由于 Windows 系统不区分文件名中的大小写，因此在 WSL 下 Anaconda 时可能会遇到问题。因此，推荐在双系统或者服务器上安装 Anaconda。
前往 清华开源镜像 下，找到最新版本的 Anacodna 的 Linux 版本并下载；若在服务器上不太容易打开网页，则可以使用下述命令行下载 5.3 版本的 Anaconda：
$ wget https://mirror.tuna.tsinghua.edu.cn/anaconda/archive/Anaconda3-5.3.0-Linux-x86_64.sh
以
Anaconda3-5.3.0-Linux-x86_64.sh文件为例，在下载文件夹下直接执行：$ chmod +x Anaconda3-5.3.0-Linux-x86_64.sh $ ./Anaconda3-5.3.0-Linux-x86_64.sh
随后按照指示进行安装即可。
小技巧
安装后，需要前往
$HOME/.bashrc检查是否将该版本的 Anaconda 加入了环境变量；参考安装日志中出现下述输出的上下文：You may wish to edit your $HOME/.bashrc to setup Anaconda3:
对于服务器用户，可能需要改动的文件不是
$HOME/.bashrc而是$HOME/.bash_profile。若在中国，可以使用清华镜像的仓库加速 Anaconda 的库的下载速度与更新速度；其使用方式是 (引用自 Anaconda 镜像使用帮助)
$ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/ $ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/ $ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/ $ conda config --set show_channel_urls yes
若要使用 Anaconda 发行版的 Python，退出当前 Terminal 再重新进入即可。
Psi4 环境¶
在这份教程中，我们不需要改动 Psi4 的程序，也不需要进行编译，只需要获得其 Python 版本的二进制可执行文件即可。这里的安装过程主要参考 Psi4NumPy 上的说明。
在安装完 conda 或 Anaconda 后，执行
$ conda create -n p4env psi4 -c psi4/label/dev
注意
一方面，我们需要使用 DFT 模块，因此需要下载
psi4/label/dev而并非psi4的 Psi4 版本；另一方面，上述的命令是创建了一个虚拟环境，它是专门为 Psi4 创建的环境。这么做是因为避免与最新版本的 Anaconda 产生库的依赖冲突，保证默认的 Python 比较干净。因此，这里没有直接在默认的 Python 环境下安装 Psi4。这样做多少会对使用产生不便，但避免库依赖关系混乱可能导致的更严重的问题。
在每次需要使用 Psi4 或维护其库依赖关系时，需要在 Bash 下执行
$ source activate p4env当 Terminal 前有提示
(p4env)时，即意味着进入 Psi4 的虚拟环境了。以后我们假设所有的命令都在该虚拟环境下执行。Psi4 的 Python 二进制文件已经可以使用了；但 Jupyter 与 MatPlotLib 并不在其依赖关系中；而这些库是我们需要的。因此，我们需要在 Psi4 的虚拟环境下执行
(p4env) $ conda install jupyter matplotlib
至此我们已经完成了 Psi4 的安装。Psi4 可以作为一个量化软件，也可以作为 Python API 使用。对于前者，我们可以简单地使用一个输入文件作测试：
input.dat在 Bash 下使用下述命令进行测试：
(p4env) $ psi4 input.dat
如果能正常地看到
output.dat且有正常的输出信息，即表明安装正常。我们也可以尝试在 Python 下做一个小测试；如果看到与下述输出一样的信息，则表明 Python API 可以正常调用：
>>> import psi4 >>> mol = psi4.geometry(""" ... O 0.000 -0.000 -0.079 ... H 0.000 0.707 0.628 ... H 0.000 -0.707 0.628 ... symmetry c1 ... """) >>> psi4.set_options({'basis': '6-31g'}) >>> psi4.core.set_output_file('output.dat', False) >>> scf_e, scf_wfn = psi4.energy('B3LYP', return_wfn=True) >>> scf_e -76.3771897718305
Jupyter 服务器环境¶
我们的操作系统通常设为 Linux (Mac 亦可)。通常这会不太方便，因为主要操作系统一般是 Windows，因此我们会期望将 Psi4 部署在远程 Linux 服务器上。而我们又会大量使用 Jupyter Notebook，其默认使用的地址是服务器的本地地址 (127.0.0.1:8080)；而这对于本地电脑而言是不可访问的。因此，需要对 Jupyter Notebook 的地址进行更改，才能让本地电脑访问服务器所启动的 Jupyter Notebook。这里参考 Jupyter Notebook 的官方文档 Running a notebook server 讲述如何配置 Jupyter Notebook。
小技巧
如果 Psi4 可以部署在本地电脑，这一节的内容可以跳过。
首先执行下述语句：
(p4env) $ jupyter notebook --generate-config
这将产生 Jupyter Notebook 的配置文件
$HOME/.jupyter/jupyter_notebook_config.py在 Jupyter Notebook 配置文件中，你将看到下述语句：
#c.NotebookApp.ip = 'localhost'对上述语句取消注释，并将其中的
localhost更改为服务器的 IP 地址。Jupyter Notebook 的服务器环境就设立好了。我们可以试一下 Jupyter Notebook 了。在 Bash 下执行
(p4env) $ jupyter notebook --no-browser
将会弹出一些输出。我们关心下述输出
Copy/paste this URL into your browser when you connect for the first time, to login with a token:
后面一行的地址；将该地址复制到本地计算机的浏览器中，就可以使用服务器的 Jupyter Notebook 了。
Define molecule¶
In [1]:
# ==> Import statements & Global Options <==
import psi4
import numpy as np
# Set Psi4 memeory as about 2GB
psi4.set_memory(int(0.5e9))
# Set output file as output.dat
psi4.core.set_output_file('output.dat', False)
In [2]:
# ==> Molecule & Psi4 Options Definitions <==
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
""")
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})
In [3]:
# Get the SCF wavefunction & energies
scf_e, scf_wfn = psi4.energy('B3LYP', return_wfn=True)
In [4]:
# B3LYP Energy
scf_e
Out[4]:
-76.37718285084564
MP2 part of XYG3¶
In [5]:
# ==> ERIs <==
# Create instance of MintsHelper class
mints = psi4.core.MintsHelper(scf_wfn.basisset())
In [6]:
# ==> Get orbital information & energy eigenvalues <==
# Number of Occupied orbitals & MOs
ndocc = scf_wfn.nalpha()
nmo = scf_wfn.nmo()
# Get orbital energies, cast into NumPy array, and separate occupied & virtual
eps = np.asarray(scf_wfn.epsilon_a())
e_ij = eps[:ndocc]
e_ab = eps[ndocc:]
In [7]:
# ==> Construct MO integrals (ia|jb) = <ij|ab> <==
Co = scf_wfn.Ca_subset('AO','OCC')
Cv = scf_wfn.Ca_subset('AO','VIR')
MO = np.asarray(mints.mo_eri(Co, Cv, Co, Cv))
In [8]:
# ==> Compute MP2 Correlation & MP2 Energy <==
# Compute energy denominator array
e_denom = 1 / (e_ij.reshape(-1, 1, 1, 1) - e_ab.reshape(-1, 1, 1) + e_ij.reshape(-1, 1) - e_ab)
# Compute SS & OS MP2 Correlation with Einsum
mp2_os_corr = np.einsum('iajb,iajb,iajb->', MO, MO, e_denom)
mp2_ss_corr = np.einsum('iajb,iajb,iajb->', MO, MO - MO.swapaxes(1,3), e_denom)
# Total MP2 Energy
MP2_corr = mp2_os_corr + mp2_ss_corr
In [9]:
# ==> Compute Scaled MP2 Correlation <==
MP2_corr * 0.3211
Out[9]:
-0.06342111430728188
B3LYP Energy Decomposition¶
XC Part¶
XC part of B3LYP energy can be calculated when forming XC potential $ V_{
}^
^
In [10]:
# Number of basis functions
nbf = scf_wfn.nmo()
# Empty matrix; XC potential (contribution to Fock matrix) should be stored in this variable
# once `compute_V' is called
V = psi4.core.Matrix(nbf, nbf)
# DFT potential calculation engine
Vpot = scf_wfn.V_potential()
# SCF AO density
D = scf_wfn.Da()
# Set AO density to DFT potential engine
Vpot.set_D([D])
# Calculate XC potential, meanwhile energy contribution of XC part to total B3LYP energy
# is also obtained
Vpot.compute_V([V])
# XC energy can be dumped by the following method
Vpot.quadrature_values()["FUNCTIONAL"]
Out[10]:
-7.537715046264638
Comparing to the energy calculated by Psi4
In [11]:
scf_wfn.get_energies("XC")
Out[11]:
-7.537715054583692
One Electron Part¶
Energy of core Hamiltonian can be obtained by $ E_{
H^
} = H_{
}^
} $
In [12]:
# Obtain core Hamiltonian in AO basis
H = np.asarray(scf_wfn.H())
In [13]:
# Calculate
np.einsum("ij, ij ->", H, D) * 2
Out[13]:
-122.34603005652934
Two Electron Part¶
Coulomb and exchange contribution to energy can be obtained by
$ J = P_{
} (
|
) P_{
} $
$ K = P_{
} (
|
) P_{
} $
Since $ c_
= 0.2 $ when using B3LYP, so two electron contribution to energy is $ 2 J - c_
J Part
In [14]:
# Atomic integral in (pq|rs)
I = np.asarray(mints.ao_eri())
# two electron part: 2 J - 0.2 K
2 * np.einsum("pqrs, pq, rs ->", I, D, D, optimize=True) \
- np.einsum("pqrs, pr, qs ->", I, D, D, optimize=True) * 0.2
Out[14]:
44.66554208711975
Total Energy¶
All energy contributions have been calculated. We can compare the following result to the B3LYP energy.
In [15]:
scf_wfn.get_energies("Nuclear") \
+ np.einsum("ij, ij ->", H, D) * 2 \
+ 2 * np.einsum("pqrs, pq, rs ->", I, D, D, optimize=True) \
- np.einsum("pqrs, pr, qs ->", I, D, D, optimize=True) * 0.2 \
+ Vpot.quadrature_values()["FUNCTIONAL"]
Out[15]:
-76.37718285084576
XYG3 Non-Consistent Part¶
First, we need to define XYG3 non-consistent DFT energy form:
In [16]:
# ==> Define Functional <===
def build_xyg3_nc_superfunctional(name, npoints, deriv, restricted):
# Build a empty superfunctional
sup = psi4.core.SuperFunctional.blank()
# No spaces, keep it short and according to convention
sup.set_name('XYG3NC')
sup.set_description(' XYG3 Non-Consistent Functional without MP2 Part\n')
# -0.0140 * LDA Exchange
# -0.0140 = 1 - 0.8033 - 0.2107
lda_x = psi4.core.LibXCFunctional("XC_LDA_X", restricted)
lda_x.set_alpha(-0.0140)
sup.add_x_functional(lda_x)
# 0.2107 * B88 Exchange
gga_x = psi4.core.LibXCFunctional("XC_GGA_X_B88", restricted)
gga_x.set_alpha(0.2107)
sup.add_x_functional(gga_x)
# 0.6789 * LYP Correlation
# 0.6789 = 1 - 0.3211
lyp_c = psi4.core.LibXCFunctional("XC_GGA_C_LYP", restricted)
lyp_c.set_alpha(0.6789)
sup.add_c_functional(lyp_c)
# 0.8033 Exact Exchange
sup.set_x_alpha(0.8033)
return sup
Although we do not need to make calculation using XYG3 non-consistent functional. However, since wavefunction can be obtained after an SCF job, I just do one energy calculation for convenience.
This energy calculation is only to obtain wavefunction and thus DFT energy and potential calculation engine.
In [17]:
nscf_e, nscf_wfn = psi4.energy("SCF", dft_functional=build_xyg3_nc_superfunctional, return_wfn=True)
In [18]:
# Although we can get XC part of energy
# !!! THIS VALUE CAN NOT BE USED !!!
# since we are doing SCF calculation where XC defined as XYG3 non-consistent
nscf_wfn.get_energies("XC")
Out[18]:
-2.003929060426492
In [19]:
# Non-consistent potential
Vn = psi4.core.Matrix(nbf, nbf)
# DFT potential calculation engine of XYG3 XC form
Vnpot = nscf_wfn.V_potential()
# Set **SCF** AO density to XYG3 **non-SCF** potential engine
Vnpot.set_D([D])
Vnpot.compute_V([Vn])
# The following value is non-consistent XC contribution to XYG3 energy
Vnpot.quadrature_values()["FUNCTIONAL"]
Out[19]:
-2.0048836531553964
Then Total energy can be obtained by the following code, using SCF AO density:
In [20]:
xyg3_e = scf_wfn.get_energies("Nuclear")
xyg3_e += np.einsum("ij, ij ->", H, D) * 2
xyg3_e += 2 * np.einsum("pqrs, pq, rs ->", I, D, D, optimize=True)
xyg3_e -= np.einsum("pqrs, pr, qs ->", I, D, D, optimize=True) * 0.8033
xyg3_e += Vnpot.quadrature_values()["FUNCTIONAL"]
xyg3_e += MP2_corr * 0.3211
In [21]:
# Compare that to Gaussian within 6 decimals
psi4.compare_values(xyg3_e, -0.76282393305943E+02, 6, 'XYG3 Energy')
XYG3 Energy.......................................................PASSED
Out[21]:
True
In [22]: