跳到主要内容
版本:2024.05-cpu

Appendix I: features Wiki

本节简要介绍 PWMLFF 中使用的特征。同时还列出了相关文献,供读者参考。

什么是特征?

特征(或描述符)是描述原子局部环境的量。它们需要保持平移、旋转和置换对称性。特征通常用作各种回归器(线性模型、神经网络等)的输入,这些回归器输出原子能量和力。

特征是空间坐标的可微函数,因此可以计算力:

Fi=dEtotdRi=j,αEjGj,αGj,αRi \mathbf{F_i} = - \frac{d E_{tot}}{d \mathbf{R_i}} = - \sum_{j,\alpha} \frac{\partial E_j}{\partial G_{j,\alpha}} \frac{\partial G_{j,\alpha}}{ \partial \mathbf{R_i}}

其中,jj 是在截断半径内的近邻原子的索引,α\alpha 是特征的索引。

注意:特征需要满足旋转、平移和置换不变性。

2-b and 3-b features with piecewise cosine functions (feature 1 & 2)

给定一个中心原子,利用分断余弦函数来描述其局部环境。通过下面的图表,可以大致了解它们的原理。

features

我们首先定义分段余弦函数,分别用于两体和三体特征。给定内部和外部截断 RinnerR_{inner}RouterR_{outer},基函数的阶数 MM,分段函数的宽度 hh,以及中心原子 ii 和近邻原子 jj 之间的原子间距 RijR_{ij},我们定义基函数为

ϕα(Rij)={12cos(RijRαhπ)+12,RijRα<h0,otherwise \phi_\alpha (R_{ij}) = \begin{cases} \frac{1}{2}\cos(\frac{R_{ij}-R_{\alpha}}{h}\pi) + \frac{1}{2} &, |R_{ij} - R_{\alpha}| < h \\ 0 &, \text{otherwise} \\ \end{cases}

其中

Rα=Rinner+(α1)h, α=1,2,...,M R_{\alpha} = R_{inner} + (\alpha - 1) h,\ \alpha = 1,2,...,M

中心原子 ii两体特征 表达式为

Gα,i=mϕα(Rij) G_{\alpha,i} = \sum_{m} \phi_{\alpha}(R_{ij})

三体特征 表达式为

Gαβγ,i=j,kϕα(Rij)ϕβ(Rik)ϕγ(Rjk) G_{\alpha\beta\gamma,i} = \sum_{j,k} \phi_{\alpha}(R_{ij}) \phi_{\beta}(R_{ik}) \phi_{\gamma}(R_{jk})

其中 m\sum_{m}m,n\sum_{m,n} 分别表示在中心原子 ii 的截断半径 RouterR_{outer} 内的所有原子的求和。

这两个特征通常是成对使用的。

参考文献

Huang, Y., Kang, J., Goddard, W. A. & Wang, L.-W. Density functional theory based neural network force fields from energy decompositions. Phys. Rev. B 99, 064103 (2019)

2-b and 3-b Gaussian feature (feature 3 & 4)

这两个特征是 Behler-Parrinello 神经网络中首次使用的特征。给定截断半径 RcR_c 及中心原子 ii 和近邻原子 jj 之间的原子间距 RijR_{ij},定义截断函数 fcf_c

fc(Rij)={12cos(πRijRc)+12,Rij<Rc0,otherwise f_c(R_{ij}) = \begin{cases} \frac{1}{2}\cos(\frac{\pi R_{ij}}{R_c}) + \frac{1}{2} &, R_{ij} < R_c \\ 0 &, \text{otherwise} \\ \end{cases}

中心原子 ii两体高斯 特征定义为

Gi=jie(η(RijRs)2)fc(Rij) G_i = \sum_{j \neq i} e^{(-\eta(R_{ij} - R_s)^2)} f_c (R_{ij})

其中 η\etaRsR_s 是用户定义的参数。

中心原子 ii三体高斯 特征定义为

Gi=21ζj,ki(1+λcosθijk)ζ eη(Rij2+Rik2+Rjk2)fc(Rij)fc(Rik)fc(Rjk) G_i = 2^{1-\zeta} \sum_{j,k \neq i} (1+\lambda \cos \theta_{ijk} )^\zeta\ e^{-\eta(R_{ij}^2 + R_{ik}^2 + R_{jk}^2)} f_c (R_{ij}) f_c (R_{ik}) f_c (R_{jk})

其中

cosθijk=RijRikRijRik \cos \theta_{ijk} = \frac{\mathbf{R_{ij}} \cdot \mathbf{R_{ik}}}{|\mathbf{R_{ij}}||\mathbf{R_{ik}}|}

η\etaζ\zetaλ=±1\lambda = \pm1 是用户定义的参数。

这两个特征通常是成对使用的。

参考文献

J. Behler and M. Parrinello, Generalized Neural-Network Representation of High Dimensional Potential-Energy Surfaces. Phys. Rev. Lett. 98, 146401 (2007)

Moment Tensor Potential (feature 5)

在 MTP 中,中心原子 ii 的局部环境由

ni=(zi,zj,rij) \mathbf{n_i} = (z_i, z_j, \mathbf{r_{ij}})

定义,其中 ziz_i 是中心原子的原子类型,zjz_j 是近邻原子的原子类型,rij\mathbf{r_{ij}} 是近邻原子的相对坐标。接下来,每个原子的能量贡献被展开为

Ei(ni)=αcαBα(ni) E_i(\mathbf{n_i}) = \sum_\alpha c_\alpha B_\alpha(\mathbf{n_i})

其中 BαB_\alpha 是用户选择的基函数,cαc_\alpha 是待拟合的参数。

为了构造基函数,我们引入矩张量 MμνM_{\mu\nu} 来定义基函数

Mμν(ni)=jfμ(rij,zi,zj)νrij M_{\mu\nu} (\mathbf{n_i}) = \sum_j f_\mu (|\mathbf{r_{ij}}|,z_i,z_j) \bigotimes_\nu \mathbf{r_{ij}}

这些矩张量包含径向和角度部分。径向部分可以展开为

fμ(rij,zi,zj)=βcμ,zi,zj(β)Q(β)(rij) f_\mu (|\mathbf{r_{ij}}|,z_i,z_j) = \sum_\beta c^{(\beta)}_{\mu,z_i,z_j} Q^{(\beta)}(|\mathbf{r_{ij}}|)

其中 Q(β)(rij)Q^{(\beta)}(|\mathbf{r_{ij}}|) 是径向基函数。具体地,

Q(β)(rij)={ϕ(β)(rij)(Rcut(rij))2,(rij)<Rcut0,otherwise Q^{(\beta)}(|\mathbf{r_{ij}}|) = \begin{cases} \phi^{(\beta)}(|\mathbf{r_{ij}}|) (R_{cut} - (|\mathbf{r_{ij}}|))^2 &, (|\mathbf{r_{ij}}|) < R_{cut} \\ 0 &,\text{otherwise} \end{cases}

其中 ϕ(β)\phi^{(\beta)} 是定义在区间 [Rmin,RcutR_{min},R_{cut}] 上的多项式(例如切比雪夫多项式)。

角度部分由 νrij\bigotimes_\nu \mathbf{r_{ij}} 给出,它表示对 rij\mathbf{r_{ij}} 进行 ν\nu 次张量积,包含了近邻 ni\mathbf{n_i} 的角度信息。ν\nu 决定了矩张量的秩。当 ν=0\nu=0 时,得到一个常数标量;当 ν=1\nu=1 时,得到一个向量(秩-1 张量);当 ν=2\nu=2 时,得到一个矩阵(秩-2 张量);以此类推。

最后,我们定义矩张量的级数为

lev(Mμν)=2+4μ+ν lev(M_{\mu \nu}) = 2 + 4\mu + \nu

这是一个经验公式。

参考文献

I.S. Novikov, etal, The MLIP package: moment tensor potential with MPI and active learning. Mach. Learn.: Sci. Technol, 2, 025002 (2021)

Spectral Neighbor Analysis Potential (feature 6)

在 SNAP 中,不使用高斯基函数。因此不计算两个原子局域环境图之间的距离和核函数。它首先定义一个原子局域环境,然后使用球谐函数(或 4D 球,带有旋转矩阵)来展开原子局域环境。然后使用双谱,使其具有旋转不变性。从某种意义上说,它类似于 MTP,但它使用一种特殊的方法来缩并方向指数,使其具有旋转不变性。它通常与线性回归一起使用。

首先,它定义位于 r\mathbf{r} 处的中心原子 ii 的邻近原子周围的原子局域环境为三维空间中的 δ\delta 函数之和:

ρ(r)=δ(r)+rki<RCfC(rki)ωkδ(rrki) \rho(\mathbf{r}) = \delta({\mathbf{r}}) + \sum_{\mathbf{r}_{ki}\lt R_C}f_C(\mathbf{r}_{ki})\omega_k\delta(\mathbf{r}-\mathbf{r}_{ki})

其中 rki\mathbf{r}_{ki} 是原子 ii 的第 kk 个近邻的位置,ωk\omega_k 是第 kk 个近邻的权重,径向函数 fC(rki)f_C(\mathbf{r}_{ki}) 确保每个近邻的贡献在截断半径 RCR_C 附近平滑地变为零:

fC(r)=0.5[cos(πrRC)+1]f_C(\mathbf{r}) = 0.5\left[\cos\left(\frac{\pi r}{R_C}\right)+1\right]

这个原子局域环境函数的角部分可以用球谐函数展开,球谐函数定义在 l=0,1,2,...l = 0, 1, 2, ...m=l,l+1,...,l1,lm = -l, -l+1, ..., l-1, l 的基础上。径向分布通常由一组径向基函数表示。然而,在这里,径向信息 r\mathbf{r} 被映射到 4D 超球面函数 Ummj(θ0,θ,ϕ)U^j_{mm^{'}}(\theta_0,\theta,\phi) 中,其中所有点(邻近原子)落入 3D 球面(在 4D 空间中),定向(角度)信息由三个角度给出:

r(xyz)ϕ=arctan(y/x)θ=arccos(z/r)θ0=34πr/rc\mathbf{r} \equiv \begin{pmatrix} x \\ y \\ z \end{pmatrix} \rightarrow \begin{matrix} \phi = \arctan(y/x) \\ \theta = \arccos(z/\mathbf{r}) \\ \theta_0 = \frac{3}{4} \pi \mathbf{r} / \mathbf{r}_{c} \end{matrix}

因此,上述原子局域环境函数可以用这些 4D 超球面函数 Ummj(θ0,θ,ϕ)U^j_{mm^{'}}(\theta_0,\theta,\phi) 展开,展开系数为 ummju^j_{mm^{'}}

ρ(r)=j=0,12,1,...m=j,j+1jm=j,j+1,...jummjUmmj(θ0,θ,ϕ) \rho(\mathbf{r}) = \sum_{j=0,\frac{1}{2},1,...}^\infin \sum_{m=-j, -j+1}^j \sum_{m^{'}=-j,-j+1,...}^j u^j_{mm^{'}} U^j_{mm^{'}}(\theta_0,\theta,\phi)

使用上述原子局域环境函数,可以计算 ummju^j_{mm^{'}}

ummj=Ummj(0,0,0)+rki<RCfC(rki)ωkUmmj(θ0(k),θ(k),ϕ(k)) u^j_{mm^{'}} = U^j_{mm^{'}}(0,0,0) + \sum_{\mathbf{r}_{ki}\lt R_C}f_C(\mathbf{r}_{ki})\omega_kU^j_{mm^{'}}(\theta_0(k),\theta(k),\phi(k))

其中,kk 是邻近原子的索引,θ0(k),θ(k),ϕ(k)\theta_0(k),\theta(k),\phi(k) 是原子 ii 的第 kk 个近邻的位置矢量的三个角度。注意,ummju^j_{mm^{'}} 是由于其指数 m,mm, m^{'} 而具有方向依赖性。它们可以与下面的缩并公式(使用三个 ummju^j_{mm^{'}})缩并:

F(j1,j2,j)=m1,m1=j1jm2,m2=j2jm,m=jj(ummj)um1m1j1um2m2j2×Cj1m1j2m2jmCj1m1j2m2jm F(j_1,j_2,j) = \sum^j_{m_1,m_1^{'}=-j_1} \sum^j_{m_2,m_2^{'}=-j_2} \sum^j_{m,m^{'}=-j} (u^{j}_{mm^{'}})^{*}u^{j_1}_{m_1m_1^{'}} u^{j_2}_{m_2m_2^{'}} \times C_{j_1 m_1 j_2 m_2}^{jm} C_{j_1 m_1^{'} j_2 m_2^{'}}^{jm}

这里,Cj1m1j2m2jmCj1m1j2m2jmC_{j_1 m_1 j_2 m_2}^{jm} C_{j_1 m_1^{'} j_2 m_2^{'}}^{jm} 是 Clebsch-Gordan 系数,最终的标量特征是 F(j1,j2,j)F(j_1,j_2,j)。通过设置不同的 j1,j2,jj_1,j_2,j,可以产生不同的特征。注意,在这些特征中,没有径向函数索引,而是有三个角动量索引。这是因为我们已经将径向距离信息转换为 3D 球面中的第三个角度信息。

DP-Chebyshev (feature 7)

这个特征类似于 DP 的嵌入网络。它使用切比雪夫多项式作为基础。

首先,我们将 S(rij)S(\mathbf{r}_{ij}) 定义为加权的距离的倒数函数:

S(r)=fC(r)r S(\mathbf{r}) = \frac{f_C(\mathbf{r})}{\mathbf{r}} fC(r)={1,r<RC212cos(πrRC2RcRC2)+12,RC2r<RC0,r>RC f_C(\mathbf{r}) = \Bigg\lbrace{\begin{matrix} 1, \qquad\qquad\qquad \mathbf{r} \lt R_{C_2}\\ \frac{1}{2} \cos(\pi \frac{\mathbf{r} - R_{C_2}}{R_c - R_{C_2}}) + \frac{1}{2}, \quad R_{C_2} \leq \mathbf{r} \lt R_C \\ 0, \qquad\qquad\qquad \mathbf{r} \gt R_C \end{matrix}}

这里,RC2R_{C_2} 是一个平滑的截断参数,它允许在由 RCR_C 定义的局部区域的边界上平滑地将 ri\mathbf{r_i} 的分量减小到零。这个平滑函数比之前使用的 RC2R_{C_2} 更复杂一些。S(rji)S(\mathbf{r}_{ji}) 减小了远离中心原子 ii 的原子的权重。然后,我们定义径向函数 gM(s)g_M(s)深度势能切比雪夫特征中的切比雪夫多项式 CMC_M

gM(s)=CM(2RminS1). g_M(s) = C_M(2R_{min} S - 1).

这里,RminR_{min} 是最小 r\mathbf{r} 的输入。

为了构造这样的特征,我们首先计算四个分量的向量:

TM(k)=rji<RCX^ji(k)S(rji)gM(S(rji)) T_M(k) = \sum_{\mathbf{r}_{ji} \lt R_C} \hat{X}_{ji}(k) S(\mathbf{r}_{ji}) g_M(S(\mathbf{r}_{ji}))

这里,k=0,1,2,3k = 0,1,2,3(四分量向量)。它们是由通常的 x,y,zx,y,z 分量构成的,再加上 SS 分量:

{xji,yji,zji}{S(rji),x^ji,y^ji,z^ji} \{ x_{ji}, y_{ji}, z_{ji}\} \rightarrow \{ S(\mathbf{r}_{ji}), \hat x_{ji}, \hat y_{ji}, \hat z_{ji} \}

其中 x^ji=Xjirji,y^ji=Yjirji,z^ji=Zjirji\hat x_{ji} = \frac{X_{ji}}{\mathbf{r}_{ji}}, \hat y_{ji} = \frac{Y_{ji}}{\mathbf{r}_{ji}}, \hat z_{ji} = \frac{Z_{ji}}{\mathbf{r}_{ji}}rji\mathbf{r}_{ji} 的单位向量。

从这些 4D 向量中,我们可以缩并分量索引以得到标量特征:

F(M1,M2)=k=03TM1(k)TM2(k) F(M_1,M_2) = \sum_{k=0}^3 T_{M_1}(k) T_{M_2}(k)

这里,M1M_1 也编码了除切比雪夫外的原子类型的数量。因此,如果最大切比雪夫阶数是 MM,特征的数量是 Mntype(Mntype+1)/2M \cdot n_{type} \cdot (M \cdot n_{type} +1 ) / 2。我们可以通过设置不同的 MM 来产生不同的特征。

DP-Gaussian (feature 8)

这个特征类似于 DP-Chebyshev,但我们使用高斯函数代替切比雪夫多项式,并且位置和宽度参数由用户指定。

类似于 DP-Chebyshev,4D 向量构造如下:

TM(k)=rji<RCX^ji(k)gM(rji) T_M(k) = \sum_{\mathbf{r}_{ji} \lt R_C} \hat {X}_{ji}(k) g_M(\mathbf{r}_{ji}) X^(0)=S(r),X^(1)=xr,X^(2)=yr,X^(3)=zr \hat X(0) = S(\mathbf{r}^{'}), \hat X(1) = \frac{x}{\mathbf{r}}, \hat X(2) = \frac{y}{\mathbf{r}}, \hat X(3) = \frac{z}{\mathbf{r}} gM(r)=fC(r) exp((rrM)ωM) g_M(\mathbf{r}) = f_C(\mathbf{r}) \ · \exp(-\frac{(\mathbf{r} - r_M)}{\omega M}) fC(r)=12cos(πrRC)+12 f_C(\mathbf{r}) = \frac{1}{2} \cos(\frac{\pi \mathbf{r}}{R_C}) + \frac{1}{2}

缩并过程如下:

F(M1,M2)=k=03TM1(k)TM2 F(M_1,M_2) = \sum_{k=0}^3 T_{M_1}(k) T_{M_2}

这里,M1M_1 也编码了除切比雪夫外的原子类型的数量。因此,如果最大切比雪夫阶数是 MM,特征的数量是 Mntype(Mntype+1)/2M \cdot n_{type} \cdot (M \cdot n_{type} +1 ) / 2。我们可以通过设置不同的 MM 来产生不同的特征。