1. Unrestricted Kohn-Sham 一阶梯度与中间矩阵
这一节我们相对系统地讨论一阶 GGA 梯度的相关性质，包括能量表达式、一阶梯度、U 矩阵的计算与关联；从而对后续的 UKS 的计算打下基础。
1.1. 准备工作
[1]:
%load_ext autoreload
%autoreload 2
%matplotlib notebook
from matplotlib import pyplot as plt
import numpy as np
from pyscf import gto, scf, dft, grad
from pyscf.scf import ucphf
from functools import partial
from pyxdh.DerivOnce import GradUSCF, GradSCF
from pyxdh.Utilities import GridHelper, KernelHelper
from pyxdh.Utilities import NucCoordDerivGenerator, NumericDiff
import warnings
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.allclose = partial(np.allclose, atol=1e-6, rtol=1e-4)
np.set_printoptions(5, linewidth=180, suppress=True)
warnings.filterwarnings("ignore")
为了简化计算量，我们大多数时候对格点积分不作非常精细的计算，因此使用非常小的格点 (50, 194)。使用的分子是非对称的 CH3 自由基。
[2]:
mol = gto.Mole()
mol.atom = """
C 0. 0. 0.
H 1. 0. 0.
H 0. 2. 0.
H 0. 0. 1.5
"""
mol.basis = "6-31G"
mol.spin = 1
mol.verbose = 0
mol.build()
# mol = gto.Mole()
# mol.atom = """
# N 0. 0. 0.
# H 1. 0. 0.
# H 0. 2. 0.
# H 0. 0. 1.5
# """
# mol.basis = "6-31G"
# mol.spin = 0
# mol.verbose = 0
# mol.build()
[2]:
<pyscf.gto.mole.Mole at 0x7fbfa1f4fcd0>
[3]:
grids = dft.Grids(mol)
grids.atom_grid = (50, 194)
grids.build()
[3]:
<pyscf.dft.gen_grid.Grids at 0x7fbfa1f4f7f0>
[4]:
scf_eng = dft.UKS(mol)
scf_eng.xc = "B3LYPg"
scf_eng.grids = grids
scf_eng.run()
scf_eng.e_tot
[4]:
-39.60377211830869
[5]:
gradh = GradUSCF({"scf_eng": scf_eng, "cphf_tol": 1e-10})
gradh.eng
[5]:
-39.60377211830869
我们顺便定义一下数值导数计算量 gradn。可以用它来进行若干矩阵的数值导数计算。
[6]:
def mol_to_grad_helper(mol):
g = dft.Grids(mol)
g.atom_grid = (50, 194)
g.build()
mf = dft.UKS(mol)
mf.xc = "B3LYPg"
mf.grids = g
return GradUSCF({"scf_eng": mf.run()})
gradn = NucCoordDerivGenerator(mol, mol_to_grad_helper)
[7]:
def plot_diff(anal_mat, num_mat):
fig, ax = plt.subplots(figsize=(2.4, 1.8)); ax.set_xscale("log")
ax.hist(abs(anal_mat.ravel() - num_mat.ravel()), bins=np.logspace(np.log10(1e-10), np.log10(1e-1), 50), alpha=0.5)
ax.hist(abs(num_mat.ravel()), bins=np.logspace(np.log10(1e-10), np.log10(1e-1), 50), alpha=0.5)
return fig.tight_layout()
1.2. 能量 UKS 计算与相关矩阵
1.2.1. 新的基础数据结构
由于 Unrestricted 计算会在同一分子中涉及到两套占据轨道信息，因此会产生 RKS 所不会出现的各种不便利，并且要求一种新的数据结构。我们会作简单的说明。
对于 GGA 而言，
cx\(c_\mathrm{x}\) 与xc泛函名称与 RKS 情形相同
[8]:
cx, xc = gradh.cx, gradh.xc
cx, xc
[8]:
(0.2, 'B3LYPg')
nmo\(n_\mathrm{MO}\)、nao\(n_\mathrm{AO}\)、natm\(n_\mathrm{Atom}\) 与 RKS 的情形相同
[9]:
nmo, nao, natm = gradh.nmo, gradh.nao, gradh.natm
nmo, nao, natm
[9]:
(15, 15, 4)
nocc\((n_\mathrm{occ}^\alpha, n_\mathrm{occ}^\beta)\),nvir\((n_\mathrm{vir}^\alpha, n_\mathrm{vir}^\beta)\) 则是Tuple[int]类型
[10]:
nocc, nvir = gradh.nocc, gradh.nvir
nocc, nvir
[10]:
((5, 4), (10, 11))
so,sv,sa作为占据、非占、全轨道的分割，其类型也变为了Tuple[slice]类型
[11]:
so, sv, sa = gradh.so, gradh.sv, gradh.sa
so, sv, sa
[11]:
((slice(0, 5, None), slice(0, 4, None)),
(slice(5, 15, None), slice(4, 15, None)),
(slice(0, 15, None), slice(0, 15, None)))
C\(C_{\mu p}^\sigma\),e\(e_p^\sigma\) 分别是轨道系数与轨道能，维度分别是 \((\sigma, \mu, p)\) 与 \((\sigma, p)\)一般来说，\(\sigma\) 在 np.einsum 程序中会用
x,y等表示\(\sigma\) 在以后的程序中，通常置于维度的第一位置，比被求导量 \(\mathbb{A}\) 优先
由于
C[0]\(C_{\mu p}^\alpha\) 与C[1]\(C_{\mu p}^\beta\) 维度相同，因此C使用np.ndarray储存；e同理
[12]:
C, e = gradh.C, gradh.e
C.shape, e.shape
[12]:
((2, 15, 15), (2, 15))
Co\(C_{\mu i}^\sigma\) 为占据轨道系数，维度是 \((\sigma, \mu, i)\)但是留意，由于 \(\alpha\) 与 \(\beta\) 自旋的占据轨道数不同，因此使用
Tuple[np.ndarray]储存；以后的矩阵通常也按照这两种方式分别处理
[13]:
Co = gradh.Co
Co[0].shape, Co[1].shape
[13]:
((15, 5), (15, 4))
eo\(e_i^\sigma\) 为占据轨道能，维度 \((\sigma, i)\)，类型同理Tuple[np.ndarray]
[14]:
eo = gradh.eo
eo[0].shape, eo[1].shape
[14]:
((5,), (4,))
上述的情况也可以用于定义非占轨道系数与非占轨道能：
[15]:
Cv, ev = gradh.Cv, gradh.ev
mo_occ\(\delta_{p \in \mathrm{occ}}^\sigma\) 表示轨道占据情况
[16]:
mo_occ = gradh.mo_occ
mo_occ
[16]:
array([[1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
1.2.2. 电子积分
H_0_ao\(h_{\mu \nu}\),S_0_ao\(S_{\mu \nu}\),eri0_ao\((\mu \nu | \kappa \lambda)\),H_1_ao\(h_{\mu \nu}^\mathbb{A}\),S_1_ao\(S_{\mu \nu}^\mathbb{A}\),eri1_ao\((\mu \nu | \kappa \lambda)^\mathbb{A}\) 与 RKS 没有区别
[17]:
H_0_ao, S_0_ao, eri0_ao, H_1_ao, S_1_ao, eri1_ao = gradh.H_0_ao, gradh.S_0_ao, gradh.eri0_ao, gradh.H_1_ao, gradh.S_1_ao, gradh.eri1_ao
H_0_mo\(h_{pq}^\sigma\),S_0_mo\(S_{pq}^\sigma\), dim: \((\sigma, p, q)\), type:np.ndarrayH_1_mo\(h_{pq}^{\mathbb{A}, \sigma}\),S_1_mo\(S_{\mu \nu}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, p, q)\), type:np.ndarray
[18]:
H_0_mo, S_0_mo = gradh.H_0_mo, gradh.S_0_mo
H_0_mo.shape, S_0_mo.shape
[18]:
((2, 15, 15), (2, 15, 15))
[19]:
H_1_mo, S_1_mo = gradh.H_1_mo, gradh.S_1_mo
H_1_mo.shape, S_1_mo.shape
[19]:
((2, 12, 15, 15), (2, 12, 15, 15))
以 H_1_mo \(h_{pq}^{\mathbb{A}, \sigma}\) 为例，
[20]:
np.allclose(np.einsum("Auv, xup, xvq -> xApq", H_1_ao, C, C), H_1_mo)
[20]:
True
eri0_mo\((pq|rs)^{\sigma \sigma'}\), dim: \((\sigma \sigma', p, q, r, s)\), type:np.ndarray上述的 \(\sigma \sigma'\) 所实际指代的是 \(\alpha \alpha, \alpha \beta, \beta \beta\)；\((pq|rs)^{\sigma \sigma'}\) 中，\(p, q\) 是 \(\sigma\) 自旋的，而 \(r, s\) 是 \(\sigma'\) 自旋的
[21]:
eri0_mo = gradh.eri0_mo
eri0_mo.shape
[21]:
(3, 15, 15, 15, 15)
[22]:
np.allclose(
np.array([
np.einsum("uvkl, up, vq, kr, ls -> pqrs", eri0_ao, C[0], C[0], C[0], C[0]),
np.einsum("uvkl, up, vq, kr, ls -> pqrs", eri0_ao, C[0], C[0], C[1], C[1]),
np.einsum("uvkl, up, vq, kr, ls -> pqrs", eri0_ao, C[1], C[1], C[1], C[1]),
]),
eri0_mo
)
[22]:
True
eri1_mo\((pq|rs)^{\mathbb{A}, \sigma \sigma'}\), dim: \((\sigma \sigma', \mathbb{A}, p, q, r, s)\), type:np.ndarray，与eri0_mo\((pq|rs)^{\sigma \sigma'}\) 同理
[23]:
eri1_mo = gradh.eri1_mo
eri1_mo.shape
[23]:
(3, 12, 15, 15, 15, 15)
1.2.3. 密度矩阵 \(D_{\mu \nu}^\sigma\)
D\(D_{\mu \nu}^\sigma\), dim: \((\sigma, \mu, \nu)\), type:np.ndarray
[24]:
D = gradh.D
D.shape
[24]:
(2, 15, 15)
[25]:
np.allclose(np.einsum("xup, xvp, xp -> xuv", C, C, mo_occ), D)
[25]:
True
另一种验证方式是
但需要留意到，两种自旋的占据轨道数量并不相同；因此，不能像上面一行代码即可验证。
[26]:
(
np.allclose(np.einsum("ui, vi -> uv", Co[0], Co[0]), D[0]),
np.allclose(np.einsum("ui, vi -> uv", Co[1], Co[1]), D[1])
)
[26]:
(True, True)
1.2.4. 库伦积分 \(J_{\mu \nu}\) 与交换积分 \(K_{\mu \nu}\)
该积分尽管常用，但以后不经常用该记号。通常以后会显式地写出如何从原子轨道与密度作张量缩并。
该积分的定义与 RKS 相同。如果我们定义影响库伦积分的密度矩阵为 \(X_{\kappa \lambda}\)，那么
\(J_{\mu \nu} = (\mu \nu | \kappa \lambda) X_{\kappa \lambda}\)
\(K_{\mu \nu} = (\mu \kappa | \nu \lambda) X_{\kappa \lambda}\)
[27]:
X = np.random.randn(nao, nao)
[28]:
np.allclose(np.einsum("uvkl, kl -> uv", eri0_ao, X), scf_eng.get_j(dm=X))
[28]:
True
[29]:
np.allclose(np.einsum("ukvl, kl -> uv", eri0_ao, X), scf_eng.get_k(dm=X, hermi=0))
[29]:
True
1.2.5. 轨道与密度格点、泛函核格点
轨道与密度格点仍然用
grdh记号表示，但在文档中，会使用Tuple[GridHelper]类型储存，长度为 2。两个GridHelper类型分别代表 \(\alpha, \beta\) 密度下的轨道与密度格点。在 pyxdh 程序中，使用zip[GridIterator]。
[30]:
grdh = (GridHelper(mol, grids, D[0]), GridHelper(mol, grids, D[1]))
[31]:
ngrid = grdh[0].ngrid
ngrid
[31]:
26836
泛函核格点仍然用
kerh记号表示
[32]:
kerh = KernelHelper(grdh, xc, deriv=3)
但由于使用了带自旋密度，因此许多量与 RKS 稍有不同。简单的情况中，在 RKS 中为向量的 \(f_\rho\) 在这里变成两条向量 \((f_{\rho^\alpha}, f_{\rho^\beta})\)。
[33]:
kerh.fr
[33]:
array([[-0. , -0. , -0. , ..., -0.01168, -0.01275, -0.01399],
[-0. , -0. , -0. , ..., -0.01172, -0.01275, -0.014 ]])
而 \(f_{\gamma}\) 则分为了 \((f_{\gamma^{\alpha \alpha}}, f_{\gamma^{\alpha \beta}}, f_{\gamma^{\beta \beta}})\)。
[34]:
kerh.fg
[34]:
array([[-0. , -0. , -0. , ..., -0.36551, -0.46675, -0.43452],
[ 0. , 0. , 0. , ..., 0.21322, 0.34198, 0.32137],
[-0. , -0. , -0. , ..., -0.36812, -0.43788, -0.4247 ]])
我们下面系统地整理一下 PySCF 中对泛函核导数的 说明:
exc\(f\) 与 RKS 较为接近
[35]:
kerh.exc.shape
[35]:
(26836,)
fr2:u, d\(f_{\rho^\alpha}, f_{\rho^\beta}\)
fg3:uu, ud, dd\(f_{\gamma^{\alpha \alpha}}, f_{\gamma^{\alpha \beta}}, f_{\gamma^{\beta \beta}}\)
[36]:
kerh.fr.shape[0], kerh.fg.shape[0]
[36]:
(2, 3)
frr3:u_u, u_d, d_d\(f_{\rho^\alpha \rho^\alpha}, f_{\rho^\alpha \rho^\beta}, f_{\rho^\beta \rho^\beta}\)
frg6:u_uu, u_ud, u_dd, d_uu, d_ud, d_dd\(f_{\rho^\alpha \gamma^{\alpha \alpha}}, f_{\rho^\alpha \gamma^{\alpha \beta}}, f_{\rho^\alpha \gamma^{\beta \beta}}, f_{\rho^\beta \gamma^{\alpha \alpha}}, f_{\rho^\beta \gamma^{\alpha \beta}}, f_{\rho^\beta \gamma^{\beta \beta}}\)
fgg6:uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd\(f_{\gamma^{\alpha \alpha} \gamma^{\alpha \alpha}}, f_{\gamma^{\alpha \alpha} \gamma^{\alpha \beta}}, f_{\gamma^{\alpha \alpha} \gamma^{\beta \beta}}, f_{\gamma^{\alpha \beta} \gamma^{\alpha \beta}}, f_{\gamma^{\alpha \beta} \gamma^{\beta \beta}}, f_{\gamma^{\beta \beta} \gamma^{\beta \beta}}\)
[37]:
kerh.frr.shape[0], kerh.frg.shape[0], kerh.fgg.shape[0]
[37]:
(3, 6, 6)
frrr4:u_u_u, u_u_d, u_d_d, d_d_d\(f_{\rho^\alpha \rho^\alpha \rho^\alpha}, f_{\rho^\alpha \rho^\alpha \rho^\beta}, f_{\rho^\alpha \rho^\beta \rho^\beta}, f_{\rho^\beta \rho^\beta \rho^\beta}\)
frrg9:u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd\(f_{\rho^\alpha \rho^\alpha \gamma^{\alpha \alpha}}, f_{\rho^\alpha \rho^\alpha \gamma^{\alpha \beta}}, f_{\rho^\alpha \rho^\alpha \gamma^{\beta \beta}}, f_{\rho^\alpha \rho^\beta \gamma^{\alpha \alpha}}, f_{\rho^\alpha \rho^\beta \gamma^{\alpha \beta}}, f_{\rho^\alpha \rho^\beta \gamma^{\beta \beta}}, f_{\rho^\beta \rho^\beta \gamma^{\alpha \alpha}}, f_{\rho^\beta \rho^\beta \gamma^{\alpha \beta}}, f_{\rho^\beta \rho^\beta \gamma^{\beta \beta}}\)
frgg12:u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd\(f_{\rho^\alpha \gamma^{\alpha \alpha} \gamma^{\alpha \alpha}}, f_{\rho^\alpha \gamma^{\alpha \alpha} \gamma^{\alpha \beta}}, f_{\rho^\alpha \gamma^{\alpha \alpha} \gamma^{\beta \beta}}, f_{\rho^\alpha \gamma^{\alpha \beta} \gamma^{\alpha \beta}}, f_{\rho^\alpha \gamma^{\alpha \beta} \gamma^{\beta \beta}}, f_{\rho^\alpha \gamma^{\beta \beta} \gamma^{\beta \beta}}, f_{\rho^\beta \gamma^{\alpha \alpha} \gamma^{\alpha \alpha}}, f_{\rho^\beta \gamma^{\alpha \alpha} \gamma^{\alpha \beta}}, f_{\rho^\beta \gamma^{\alpha \alpha} \gamma^{\beta \beta}}, f_{\rho^\beta \gamma^{\alpha \beta} \gamma^{\alpha \beta}}, f_{\rho^\beta \gamma^{\alpha \beta} \gamma^{\beta \beta}}, f_{\rho^\beta \gamma^{\beta \beta} \gamma^{\beta \beta}}\)
fggg10:uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd\(f_{\gamma^{\alpha \alpha} \gamma^{\alpha \alpha} \gamma^{\alpha \alpha}}, f_{\gamma^{\alpha \alpha} \gamma^{\alpha \alpha} \gamma^{\alpha \beta}}, f_{\gamma^{\alpha \alpha} \gamma^{\alpha \alpha} \gamma^{\beta \beta}}, f_{\gamma^{\alpha \alpha} \gamma^{\alpha \beta} \gamma^{\alpha \beta}}, f_{\gamma^{\alpha \alpha} \gamma^{\alpha \beta} \gamma^{\beta \beta}}, f_{\gamma^{\alpha \alpha} \gamma^{\beta \beta} \gamma^{\beta \beta}}, f_{\gamma^{\alpha \beta} \gamma^{\alpha \beta} \gamma^{\alpha \beta}}, f_{\gamma^{\alpha \beta} \gamma^{\alpha \beta} \gamma^{\beta \beta}}, f_{\gamma^{\alpha \beta} \gamma^{\beta \beta} \gamma^{\beta \beta}}, f_{\gamma^{\beta \beta} \gamma^{\beta \beta} \gamma^{\beta \beta}}\)
[38]:
kerh.frrr.shape[0], kerh.frrg.shape[0], kerh.frgg.shape[0], kerh.fggg.shape[0]
[38]:
(4, 9, 12, 10)
1.2.6. 能量计算
[39]:
(
+ np.einsum("uv, xuv -> ", H_0_ao, D)
+ 0.5 * np.einsum("uvkl, xuv, ykl -> ", eri0_ao, D, D)
- 0.5 * cx * np.einsum("ukvl, xuv, xkl -> ", eri0_ao, D, D)
+ np.einsum("g, g -> ", kerh.exc, grdh[0].rho_0)
+ np.einsum("g, g -> ", kerh.exc, grdh[1].rho_0)
)
[39]:
-47.22493668809642
[40]:
scf_eng.energy_elec()[0]
[40]:
-47.22493669052412
1.2.7. Fock 矩阵
F_0_ao\(F_{\mu \nu}^\sigma\), dim: \((\sigma, \mu, \nu)\), type:np.ndarray
需要注意，不像其他电子积分的原子轨道矩阵一样，Fock 矩阵是有自旋之分的。生成 Fock 矩阵的过程稍复杂，因此我们需要将 HF 贡献部分先计算，随后计算 GGA 贡献部分，最后加和。
[41]:
F_0_ao = gradh.F_0_ao
F_0_ao.shape
[41]:
(2, 15, 15)
[42]:
F_0_ao_ = (
+ H_0_ao
+ np.einsum("uvkl, ykl -> uv", eri0_ao, D)
- cx * np.einsum("ukvl, xkl -> xuv", eri0_ao, D)
)
[43]:
F_0_ao_GGA_ = np.zeros_like(F_0_ao_)
F_0_ao_GGA_[0] = (
+ 2 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg[0], grdh[0].rho_1, grdh[0].ao_1, grdh[0].ao_0)
+ np.einsum("g, rg, rgu, gv -> uv", kerh.fg[1], grdh[1].rho_1, grdh[0].ao_1, grdh[0].ao_0)
)
F_0_ao_GGA_[1] = (
+ 2 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg[2], grdh[1].rho_1, grdh[0].ao_1, grdh[0].ao_0)
+ np.einsum("g, rg, rgu, gv -> uv", kerh.fg[1], grdh[0].rho_1, grdh[0].ao_1, grdh[0].ao_0)
)
F_0_ao_GGA_ += F_0_ao_GGA_.swapaxes(-1, -2)
F_0_ao_GGA_ += np.einsum("xg, gu, gv -> xuv", kerh.fr, grdh[0].ao_0, grdh[0].ao_0)
[44]:
F_0_ao_ += F_0_ao_GGA_
[45]:
np.allclose(F_0_ao_, F_0_ao)
[45]:
True
F_0_ao\(F_{pq}^\sigma\), dim: \((\sigma, p, q)\), type:np.ndarray
[46]:
F_0_mo = gradh.F_0_mo
F_0_mo.shape
[46]:
(2, 15, 15)
[47]:
np.allclose(np.einsum("xuv, xup, xvq -> xpq", F_0_ao, C, C), F_0_mo)
[47]:
True
1.3. 一阶 Skeleton 导数
1.3.1. 一阶核坐标梯度
[48]:
E_1 = gradh.E_1
E_1
[48]:
array([[ 0.06565, -0.08058, -0.11324],
[-0.0919 , 0.01814, 0.02577],
[ 0.00907, 0.05126, 0.00839],
[ 0.01718, 0.01118, 0.07909]])
[49]:
E_1_ = (
+ np.einsum("Auv, xuv -> A", H_1_ao, D)
+ 0.5 * np.einsum("Auvkl, yuv, xkl -> A", eri1_ao, D, D)
- 0.5 * cx * np.einsum("Aukvl, xuv, xkl -> A", eri1_ao, D, D)
# - np.einsum("xApq, xpq, xp, xq -> A", S_1_mo, F_0_mo, mo_occ, mo_occ)
- np.einsum("Ai, i -> A", S_1_mo[0, :, so[0], so[0]].diagonal(axis1=-1, axis2=-2), eo[0])
- np.einsum("Ai, i -> A", S_1_mo[1, :, so[1], so[1]].diagonal(axis1=-1, axis2=-2), eo[1])
+ grad.rhf.grad_nuc(mol).reshape(-1)
).reshape((natm, 3))
[50]:
E_1_ += (
+ np.einsum("g, Atg -> At", kerh.fr[0], grdh[0].A_rho_1)
+ np.einsum("g, Atg -> At", kerh.fr[1], grdh[1].A_rho_1)
+ 2 * np.einsum("g, rg, Atrg -> At", kerh.fg[0], grdh[0].rho_1, grdh[0].A_rho_2)
+ 2 * np.einsum("g, rg, Atrg -> At", kerh.fg[2], grdh[1].rho_1, grdh[1].A_rho_2)
+ 1 * np.einsum("g, rg, Atrg -> At", kerh.fg[1], grdh[1].rho_1, grdh[0].A_rho_2)
+ 1 * np.einsum("g, rg, Atrg -> At", kerh.fg[1], grdh[0].rho_1, grdh[1].A_rho_2)
)
[51]:
E_1
[51]:
array([[ 0.06565, -0.08058, -0.11324],
[-0.0919 , 0.01814, 0.02577],
[ 0.00907, 0.05126, 0.00839],
[ 0.01718, 0.01118, 0.07909]])
1.3.2. Fock Skeleton 一阶导数：HF 部分
F_1_ao\(F_{\mu \nu}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, \mu, \nu)\), type:np.ndarray
[52]:
F_1_ao = gradh.F_1_ao
F_1_ao.shape
[52]:
(2, 12, 15, 15)
该式的推导方式比 RKS 的情形要复杂很多。我们首先需要定义变量
A_gamma_1\(\gamma^{\mathbb{A}, \sigma \sigma'}\), dim: \((\sigma \sigma', A, t, g)\), type:np.ndarray
[53]:
A_gamma_1 = np.zeros((3, natm, 3, grdh[0].ngrid))
A_gamma_1[0] = 2 * np.einsum("Atrg, rg -> Atg", grdh[0].A_rho_2, grdh[0].rho_1)
A_gamma_1[1] = (
+ np.einsum("Atrg, rg -> Atg", grdh[0].A_rho_2, grdh[1].rho_1)
+ np.einsum("Atrg, rg -> Atg", grdh[1].A_rho_2, grdh[0].rho_1)
)
A_gamma_1[2] = 2 * np.einsum("Atrg, rg -> Atg", grdh[1].A_rho_2, grdh[1].rho_1)
我们之后生成的 Fock 矩阵的 Skeleton 导数会放在变量 F_1_ao_ 中 (后面加上下划线以示区别)。HF 部分的贡献相当容易给出：
[54]:
F_1_ao_ = (
+ H_1_ao
+ np.einsum("Auvkl, ykl -> Auv", eri1_ao, D)
- cx * np.einsum("Aukvl, xkl -> xAuv", eri1_ao, D)
)
1.3.3. Fock Skeleton 一阶导数：GGA \(\alpha\) 自旋部分
GGA 部分的贡献会很复杂。我们先回顾一下 GGA 对 Fock 矩阵的贡献方式：
我们按照下述方式给出 Skeleton 导数。
0：\(f_{\rho^\alpha}\) 的所有 Skeleton 导数贡献
1：\(f_{\gamma^{\alpha \beta}}\) 关于 \(\rho^\sigma\) 的所有偏导部分的 Skeleton 导数贡献
2：\(f_{\gamma^{\alpha \beta}}\) 关于 \(\gamma^{\sigma \sigma'}\) 的所有偏导部分的 Skeleton 导数贡献
3：上式中出现的 \(\rho_r^\sigma\) 的 Skeleton 导数贡献
4：上式中出现的原子轨道格点的 Skeleton 导数贡献
由于两种自旋所取用泛函核不同，因此 \(\alpha, \beta\) 两种自旋对 Fock 矩阵 Skeleton 导数的贡献代码应当要分开写。我们将所有 \(\alpha\) 自旋下，GGA 的所有贡献写到 F_1_ao_GGA_alpha_ 中。
F_1_ao_GGA_alpha_, type:List[np.ndarray]; for each element in list, dim: \((A, t, \mu, \nu)\)之所以每个列表元素的维度不是 \((\mathbb{A}, \mu, \nu)\)，单纯地是因为在
GridHelper中的所有格点导数均采用 \((A, t)\) 即 \((n_\mathrm{Atom}, 3)\) 的维度大小的形式，而不是使用 pyxdh 项目中通常会用的 \((\mathbb{A}, )\) 即 \((n_\mathrm{Atom} \times 3, )\) 的形式。
[55]:
F_1_ao_GGA_alpha_ = [None for _ in range(5)]
[56]:
tmp = (
+ kerh.frr[0] * grdh[0].A_rho_1 + kerh.frr[1] * grdh[1].A_rho_1
+ kerh.frg[0] * A_gamma_1[0] + kerh.frg[1] * A_gamma_1[1] + kerh.frg[2] * A_gamma_1[2]
)
F_1_ao_GGA_alpha_[0] = np.einsum("Atg, gu, gv -> Atuv", tmp, grdh[0].ao_0, grdh[0].ao_0)
[57]:
tmp = (
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.frg[0] * grdh[0].rho_1 + kerh.frg[1] * grdh[1].rho_1, grdh[0].A_rho_1)
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.frg[3] * grdh[0].rho_1 + kerh.frg[4] * grdh[1].rho_1, grdh[1].A_rho_1)
)
F_1_ao_GGA_alpha_[1] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_alpha_[1] += F_1_ao_GGA_alpha_[1].swapaxes(-1, -2)
[58]:
tmp = (
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[0] * grdh[0].rho_1 + kerh.fgg[1] * grdh[1].rho_1, A_gamma_1[0])
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[1] * grdh[0].rho_1 + kerh.fgg[3] * grdh[1].rho_1, A_gamma_1[1])
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[2] * grdh[0].rho_1 + kerh.fgg[4] * grdh[1].rho_1, A_gamma_1[2])
)
F_1_ao_GGA_alpha_[2] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_alpha_[2] += F_1_ao_GGA_alpha_[2].swapaxes(-1, -2)
[59]:
tmp = 2 * kerh.fg[0] * grdh[0].A_rho_2 + kerh.fg[1] * grdh[1].A_rho_2
F_1_ao_GGA_alpha_[3] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_alpha_[3] += F_1_ao_GGA_alpha_[3].swapaxes(-1, -2)
[60]:
mol_slice = gradh.mol_slice
[61]:
tmp = 2 * kerh.fg[0] * grdh[0].rho_1 + kerh.fg[1] * grdh[1].rho_1
F_1_ao_GGA_alpha_[4] = np.zeros((natm, 3, nao, nao))
for A in range(natm):
sA = mol_slice(A)
F_1_ao_GGA_alpha_[4][A, :, sA, :] += - np.einsum("g, tgu, gv -> tuv", kerh.fr[0], grdh[0].ao_1[:, :, sA], grdh[0].ao_0)
F_1_ao_GGA_alpha_[4][A, :, sA, :] += - np.einsum("rg, trgu, gv -> tuv", tmp, grdh[0].ao_2[:, :, :, sA], grdh[0].ao_0)
F_1_ao_GGA_alpha_[4][A, :, sA, :] += - np.einsum("rg, tgu, rgv -> tuv", tmp, grdh[0].ao_1[:, :, sA], grdh[0].ao_1)
F_1_ao_GGA_alpha_[4] += F_1_ao_GGA_alpha_[4].swapaxes(-1, -2)
最后，我们将 \(\alpha\) 自旋的 GGA 贡献部分加到 HF 贡献部分，得到最终的 \(F_{\mu \nu}^{\mathbb{A}, \alpha}\)。
[62]:
F_1_ao_[0] += np.array(F_1_ao_GGA_alpha_).sum(axis=0).reshape((natm * 3, nao, nao))
[63]:
np.allclose(F_1_ao_[0], F_1_ao[0])
[63]:
True
1.3.4. Fock Skeleton 一阶导数：GGA \(\beta\) 自旋部分
但需要留意，我们尚没有实现 \(\beta\) 自旋的贡献部分。尽管思路是相同的，但代码应当要重新写一遍。
[64]:
F_1_ao_GGA_beta_ = [None for _ in range(5)]
[65]:
tmp = (
+ kerh.frr[1] * grdh[0].A_rho_1 + kerh.frr[2] * grdh[1].A_rho_1
+ kerh.frg[3] * A_gamma_1[0] + kerh.frg[4] * A_gamma_1[1] + kerh.frg[5] * A_gamma_1[2]
)
F_1_ao_GGA_beta_[0] = np.einsum("Atg, gu, gv -> Atuv", tmp, grdh[0].ao_0, grdh[0].ao_0)
[66]:
tmp = (
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.frg[2] * grdh[1].rho_1 + kerh.frg[1] * grdh[0].rho_1, grdh[0].A_rho_1)
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.frg[5] * grdh[1].rho_1 + kerh.frg[4] * grdh[0].rho_1, grdh[1].A_rho_1)
)
F_1_ao_GGA_beta_[1] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_beta_[1] += F_1_ao_GGA_beta_[1].swapaxes(-1, -2)
[67]:
tmp = (
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[2] * grdh[1].rho_1 + kerh.fgg[1] * grdh[0].rho_1, A_gamma_1[0])
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[4] * grdh[1].rho_1 + kerh.fgg[3] * grdh[0].rho_1, A_gamma_1[1])
+ np.einsum("rg, Atg -> Atrg", 2 * kerh.fgg[5] * grdh[1].rho_1 + kerh.fgg[4] * grdh[0].rho_1, A_gamma_1[2])
)
F_1_ao_GGA_beta_[2] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_beta_[2] += F_1_ao_GGA_beta_[2].swapaxes(-1, -2)
[68]:
tmp = 2 * kerh.fg[2] * grdh[1].A_rho_2 + kerh.fg[1] * grdh[0].A_rho_2
F_1_ao_GGA_beta_[3] = np.einsum("Atrg, rgu, gv -> Atuv", tmp, grdh[0].ao_1, grdh[0].ao_0)
F_1_ao_GGA_beta_[3] += F_1_ao_GGA_beta_[3].swapaxes(-1, -2)
[69]:
mol_slice = gradh.mol_slice
[70]:
tmp = 2 * kerh.fg[2] * grdh[1].rho_1 + kerh.fg[1] * grdh[0].rho_1
F_1_ao_GGA_beta_[4] = np.zeros((natm, 3, nao, nao))
for A in range(natm):
sA = mol_slice(A)
F_1_ao_GGA_beta_[4][A, :, sA, :] += - np.einsum("g, tgu, gv -> tuv", kerh.fr[1], grdh[0].ao_1[:, :, sA], grdh[0].ao_0)
F_1_ao_GGA_beta_[4][A, :, sA, :] += - np.einsum("rg, trgu, gv -> tuv", tmp, grdh[0].ao_2[:, :, :, sA], grdh[0].ao_0)
F_1_ao_GGA_beta_[4][A, :, sA, :] += - np.einsum("rg, tgu, rgv -> tuv", tmp, grdh[0].ao_1[:, :, sA], grdh[0].ao_1)
F_1_ao_GGA_beta_[4] += F_1_ao_GGA_beta_[4].swapaxes(-1, -2)
将 \(\beta\) 自旋的 GGA 贡献部分加到 HF 贡献部分，得到最终的 \(F_{\mu \nu}^{\mathbb{A}, \beta}\)。
[71]:
F_1_ao_[1] += np.array(F_1_ao_GGA_beta_).sum(axis=0).reshape((natm * 3, nao, nao))
最终，我们可以验证，我们计算所得到的 Fock 矩阵与 pyxdh 给出的结果一致：
[72]:
np.allclose(F_1_ao_, F_1_ao)
[72]:
True
需要指出的是，pyxdh 中计算 F_1_ao 的方式并不是直接执行类似于上述的代码，而是利用了 PySCF 的 hessian.uks.make_h1 方法。
1.4. A 张量
1.4.1. A 张量计算函数：前置工作
我们应当留意到 A 张量从定义上是从 Fock 矩阵 \(F_{\mu \nu}^\sigma\) 的 U 导数产生的。我们下面对当输入的矩阵 X \(X_{ai}^{\mathbb{A}, \sigma}\) 为任意矩阵时的情况作讨论。这里的 \(\mathbb{A}\) 可以是任意的，未必要指代原子核坐标分量。
[73]:
Ax0_Core = gradh.Ax0_Core
X = (np.random.randn(3, nvir[0], nocc[0]), np.random.randn(3, nvir[1], nocc[1]))
我们将 pyxdh 输出的 \(\mathtt{AX}_{ai}^{\mathbb{A}, \sigma}\) 储存在 Ax 中。由于 A 张量与 X 矩阵之间的缩并方式与 RKS 不太一样，因此我们不会采用类似于 \(A_{ai, bj} X_{bj}^\mathbb{A}\) 的写法。
X\(X_{ai}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, a, i)\), type:Tuple[np.ndarray]Ax\(\mathtt{AX}_{ai}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, a, i)\), type:Tuple[np.ndarray]
需要指出，Ax0_Core 函数不只可以实现非占-占据轨道的计算，也可以实现任意分子轨道下的分割计算。
[74]:
Ax = Ax0_Core(sv, so, sv, so)(X)
在实际进行 A 张量计算前，我们会对输入的 \(X_{ai}^{\mathbb{A}, \sigma}\) 作到原子轨道的转换：
dmX\(X_{\mu \nu}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, \mu, \nu)\), type:np.ndarray
[75]:
prop_dim = X[0].shape[0]
dmX = np.zeros((2, prop_dim, nao, nao))
dmX[0] = Cv[0] @ X[0] @ Co[0].T
dmX[1] = Cv[1] @ X[1] @ Co[1].T
dmX += dmX.swapaxes(-1, -2)
dmX.shape
[75]:
(2, 3, 15, 15)
1.4.2. A 张量计算函数：HF 贡献部分
ax_ao_HF_, dim: \((\sigma, \mathbb{A}, \mu, \nu)\), type:np.ndarray
[76]:
ax_ao_HF_ = (
+ np.einsum("uvkl, xAkl -> Auv", eri0_ao, dmX)
- cx * np.einsum("ukvl, xAkl -> xAuv", eri0_ao, dmX)
)
ax_ao_HF_.shape
[76]:
(2, 3, 15, 15)
Ax_HF_, dim: \((\sigma, \mathbb{A}, a, i)\)
[77]:
Ax_HF_ = (np.zeros((prop_dim, nvir[0], nocc[0])), np.zeros((prop_dim, nvir[1], nocc[1])))
Ax_HF_[0][:] = Cv[0].T @ ax_ao_HF_[0] @ Co[0]
Ax_HF_[1][:] = Cv[1].T @ ax_ao_HF_[1] @ Co[1]
1.4.3. A 张量计算函数：\(f_{\rho}\) 导数贡献部分
rho_X_0\(\varrho^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, g)\)rho_X_1\(\varrho_g^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, r, g)\)
[78]:
tmp_K = np.einsum("xAkl, gl -> xAgk", dmX, grdh[0].ao_0)
rho_X_0 = np.einsum("gk, xAgk -> xAg", grdh[0].ao_0, tmp_K)
rho_X_1 = 2 * np.einsum("rgk, xAgk -> xArg", grdh[0].ao_1, tmp_K)
gamma_X_0\(\gamma^{\mathbb{A}, \sigma \sigma'}\), dim: \((\sigma \sigma', \mathbb{A}, g)\)
需要留意到，这里单纯地因为符号不够用，因此与 Fock Skeleton 导数处使用了相同的 \(\gamma^{\mathbb{A}, \sigma \sigma'}\) 符号。但这两者的定义是不同的 (即使是相似的)。
[79]:
gamma_X_0 = np.zeros((3, prop_dim, grdh[0].ngrid))
gamma_X_0[0] = 2 * np.einsum("Arg, rg -> Ag", rho_X_1[0], grdh[0].rho_1)
gamma_X_0[1] = (
+ np.einsum("Arg, rg -> Ag", rho_X_1[0], grdh[1].rho_1)
+ np.einsum("Arg, rg -> Ag", rho_X_1[1], grdh[0].rho_1)
)
gamma_X_0[2] = 2 * np.einsum("Arg, rg -> Ag", rho_X_1[1], grdh[1].rho_1)
M_0_\(M^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, g)\)
[80]:
M_0_ = np.zeros((2, prop_dim, ngrid))
[81]:
M_0_[0] = (
+ kerh.frr[0] * rho_X_0[0]
+ kerh.frr[1] * rho_X_0[1]
+ kerh.frg[0] * gamma_X_0[0]
+ kerh.frg[1] * gamma_X_0[1]
+ kerh.frg[2] * gamma_X_0[2]
)
[82]:
M_0_[1] = (
+ kerh.frr[1] * rho_X_0[0]
+ kerh.frr[2] * rho_X_0[1]
+ kerh.frg[3] * gamma_X_0[0]
+ kerh.frg[4] * gamma_X_0[1]
+ kerh.frg[5] * gamma_X_0[2]
)
Ax_GGA_M_0_, dim: \((\sigma, \mathbb{A}, a, i)\), type:Tuple[np.ndarray]
[83]:
ax_ao_M_0_ = np.einsum("xAg, gu, gv -> xAuv", M_0_, grdh[0].ao_0, grdh[0].ao_0)
[84]:
Ax_GGA_M_0_ = (
np.einsum("Auv, ua, vi -> Aai", ax_ao_M_0_[0], Cv[0], Co[0]),
np.einsum("Auv, ua, vi -> Aai", ax_ao_M_0_[1], Cv[1], Co[1])
)
1.4.4. A 张量计算函数：\(f_{\gamma}\) 对 \(\rho\) 导数贡献部分
M_1_\(M_r^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, r, g)\)
[85]:
M_1_ = np.zeros((2, prop_dim, 3, ngrid))
[86]:
M_1_[0] = (
+ np.einsum("rg, Ag -> Arg", 2 * kerh.frg[0] * grdh[0].rho_1 + kerh.frg[1] * grdh[1].rho_1, rho_X_0[0])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.frg[3] * grdh[0].rho_1 + kerh.frg[4] * grdh[1].rho_1, rho_X_0[1])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[0] * grdh[0].rho_1 + kerh.fgg[1] * grdh[1].rho_1, gamma_X_0[0])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[1] * grdh[0].rho_1 + kerh.fgg[3] * grdh[1].rho_1, gamma_X_0[1])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[2] * grdh[0].rho_1 + kerh.fgg[4] * grdh[1].rho_1, gamma_X_0[2])
+ 2 * kerh.fg[0] * rho_X_1[0] + kerh.fg[1] * rho_X_1[1]
)
[87]:
M_1_[1] = (
+ np.einsum("rg, Ag -> Arg", 2 * kerh.frg[2] * grdh[1].rho_1 + kerh.frg[1] * grdh[0].rho_1, rho_X_0[0])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.frg[5] * grdh[1].rho_1 + kerh.frg[4] * grdh[0].rho_1, rho_X_0[1])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[2] * grdh[1].rho_1 + kerh.fgg[1] * grdh[0].rho_1, gamma_X_0[0])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[4] * grdh[1].rho_1 + kerh.fgg[3] * grdh[0].rho_1, gamma_X_0[1])
+ np.einsum("rg, Ag -> Arg", 2 * kerh.fgg[5] * grdh[1].rho_1 + kerh.fgg[4] * grdh[0].rho_1, gamma_X_0[2])
+ 2 * kerh.fg[2] * rho_X_1[1] + kerh.fg[1] * rho_X_1[0]
)
Ax_GGA_M_1_, dim: \((\sigma, \mathbb{A}, a, i)\), type:Tuple[np.ndarray]
[88]:
ax_ao_M_1_ = np.einsum("xArg, rgu, gv -> xAuv", M_1_, grdh[0].ao_1, grdh[0].ao_0)
ax_ao_M_1_ += ax_ao_M_1_.swapaxes(-1, -2)
[89]:
Ax_GGA_M_1_ = (
np.einsum("Auv, ua, vi -> Aai", ax_ao_M_1_[0], Cv[0], Co[0]),
np.einsum("Auv, ua, vi -> Aai", ax_ao_M_1_[1], Cv[1], Co[1])
)
最后我们可以验证一下计算过程的正确性：
[90]:
(
np.allclose(Ax_HF_[0] + Ax_GGA_M_0_[0] + Ax_GGA_M_1_[0], Ax[0]),
np.allclose(Ax_HF_[1] + Ax_GGA_M_0_[1] + Ax_GGA_M_1_[1], Ax[1])
)
[90]:
(True, True)
作者注意到了 A 张量的生成过程与 Fock Skeleton 导数生成的过程极其相似。在以后的程序版本中，可能会考虑引入对 M_0 \(M^{\mathbb{A}, \sigma}\) 与 M_1 \(M_r^{\mathbb{A}, \sigma}\) 的显式定义，避免程序重复并避免编写失误。
1.5. 一阶 U 矩阵
事实上，为一阶梯度所使用的所有变量都已经考虑完毕了 (因为一阶梯度不需要使用一阶 U 矩阵)。在此生成的一阶 U 矩阵有两个目的：其一是用来验证各种数值梯度；其二是熟悉 UKS 下 CP-HF 方程。
1.5.1. 数值一阶 U 矩阵
nd_C\(\partial_\mathbb{A} C_{\mu p}^\sigma\), dim: \((\sigma, \mathbb{A}, \mu, p)\)，通过数值导数获得C_inv\((\mathbf{C}^\sigma)^{-1}_{p \mu}\), dim: \((\sigma, p, \mu)\)
[91]:
nd_C = NumericDiff(gradn, lambda gradh: gradh.C).derivative.swapaxes(0, 1)
nd_C.shape
[91]:
(2, 12, 15, 15)
[92]:
C_inv = np.array([np.linalg.inv(C[0]), np.linalg.inv(C[1])])
C_inv.shape
[92]:
(2, 15, 15)
利用 \(\partial_\mathbb{A} C_{\mu p}^\sigma = C_{\mu m}^\sigma U_{mp}^{\mathbb{A}, \sigma}\)，我们能得到数值 U 矩阵
nd_U\(U_{mp}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, m, p)\)，通过数值导数获得
[93]:
nd_U = np.einsum("xmu, xAup -> xAmp", C_inv, nd_C)
我们能用数值 U 矩阵立即验证下述公式：
[94]:
plot_diff(nd_U + nd_U.swapaxes(-1, -2), - S_1_mo)
1.5.2. B 矩阵
B_1\(B_{pq}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, p, q)\)
[95]:
B_1 = gradh.B_1
B_1.shape
[95]:
(2, 12, 15, 15)
[96]:
F_1_mo = gradh.F_1_mo
B_1_ = (
+ F_1_mo
- np.einsum("xApq, xq -> xApq", S_1_mo, e)
- 0.5 * np.array(Ax0_Core(sa, sa, so, so)((S_1_mo[0, :, so[0], so[0]], S_1_mo[1, :, so[1], so[1]])))
)
[97]:
np.allclose(B_1_, B_1)
[97]:
True
1.5.3. CP-HF 方程验证
[98]:
plot_diff(
- (ev[0][:, None] - eo[0][None, :]) * nd_U[0][:, sv[0], so[0]]
- Ax0_Core(sv, so, sv, so)((nd_U[0][:, sv[0], so[0]], nd_U[1][:, sv[1], so[1]]))[0]
,
B_1[0][:, sv[0], so[0]]
)
上述等式还可以将 \(a, i\) 拓展到 \(p \neq q\) 的情形。
1.5.4. CP-HF 方程求解
在 PySCF 中，求解 Unrestricted CP-HF 方程需要使用的是 scf.ucphf 模块。为了讨论问题方便，即使是求取 U 矩阵，我们仍然不会使用到 solve_withs1 函数，而实际使用的是 solve_nos1 函数。当然，在 PySCF 中，两个函数都会被打包在 solve 函数中。
与 Restricted 情形完全不同的是，ucphf.solve 在计算过程中的传参方式与 cphf.solve 完全不同，尽管两个函数有非常类似的参数签名 (signature)。因此，我们传入的第一个参数不可以再像 Restricted CP-HF 一样使用 Ax0_Core(sv, so, sv, so)，而必须要额外定义一个函数 fx。传参过程中，U 矩阵会转化为维度 \((\mathbb{A}, n_\mathrm{vir}^\alpha n_\mathrm{occ}^\alpha + n_\mathrm{vir}^\beta n_\mathrm{occ}^\beta)\) 的中间矩阵；这些中间矩阵在下述函数中为 X 与 result 所指代。
[99]:
def fx(X):
prop_dim = X.shape[0]
X_alpha = X[:, :nocc[0] * nvir[0]].reshape((prop_dim, nvir[0], nocc[0]))
X_beta = X[:, nocc[0] * nvir[0]:].reshape((prop_dim, nvir[1], nocc[1]))
Ax = Ax0_Core(sv, so, sv, so, in_cphf=True)((X_alpha, X_beta))
result = np.concatenate([Ax[0].reshape(prop_dim, -1), Ax[1].reshape(prop_dim, -1)], axis=1)
return result
U_1_vo_\(U_{ai}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, a, i)\), type:Tuple[np.ndarray]
[100]:
U_1_vo_ = ucphf.solve(fx, e, mo_occ, (B_1[0, :, sv[0], so[0]], B_1[1, :, sv[1], so[1]]), max_cycle=100, tol=1e-10)[0]
[101]:
plot_diff(U_1_vo_[0], nd_U[0][:, sv[0], so[0]])
[102]:
plot_diff(U_1_vo_[1], nd_U[1][:, sv[1], so[1]])
1.5.5. 完整未“旋转”的 U 矩阵
U_1_\(U_{pq}^{\mathbb{A}, \sigma}\), dim: \((\sigma, \mathbb{A}, p, q)\)
[103]:
U_1_ = np.zeros((2, natm * 3, nmo, nmo))
U_1_[0, :, sv[0], so[0]] = U_1_vo_[0]
U_1_[1, :, sv[1], so[1]] = U_1_vo_[1]
[104]:
U_1_[0, :, so[0], sv[0]] = - S_1_mo[0, :, so[0], sv[0]] - U_1_vo_[0].swapaxes(-1, -2)
U_1_[1, :, so[1], sv[1]] = - S_1_mo[1, :, so[1], sv[1]] - U_1_vo_[1].swapaxes(-1, -2)
相同的公式可以用于非占-非占的情形。
[105]:
Ax_oo = Ax0_Core(so, so, sv, so)(U_1_vo_)
Ax_vv = Ax0_Core(sv, sv, sv, so)(U_1_vo_)
[106]:
U_1_[0, :, so[0], so[0]] = - (Ax_oo[0] + B_1[0, :, so[0], so[0]]) / (eo[0][:, None] - eo[0][None, :])
U_1_[1, :, so[1], so[1]] = - (Ax_oo[1] + B_1[1, :, so[1], so[1]]) / (eo[1][:, None] - eo[1][None, :])
U_1_[0, :, sv[0], sv[0]] = - (Ax_vv[0] + B_1[0, :, sv[0], sv[0]]) / (ev[0][:, None] - ev[0][None, :])
U_1_[1, :, sv[1], sv[1]] = - (Ax_vv[1] + B_1[1, :, sv[1], sv[1]]) / (ev[1][:, None] - ev[1][None, :])
[107]:
for p in range(nmo):
U_1_[:, :, p, p] = - S_1_mo[:, :, p, p] / 2
[108]:
plot_diff(U_1_, nd_U)
