Algorithm used in outwin in WIEN2k (CN)

背景

outwin.f在WIEN2k各程序中出现, 它包含例程outwin. 之前大略知道它是用来计算原子球内的径向波函数的, 但对于其算法一直很模糊, 一方面由于它涉及相对论方程, 另一方面它除了输入参数的德语注释外一句注释也没有. 最近在研究局域基组生成的问题, 而这个例程出现频率非常高, 因此准备多啃一下这块代码.

不同SRC文件下的 outwin.f 版本也不尽同. SRC_nmr, SRC_lapw7等仍然使用Adams-Moulton四阶算法, SRC_lapw2 对于第四个以外的格点允许用五阶算法. SRC_lapw7 中的注释更多一些, 但仍然用的是四阶算法. 这里尝试对 SRC_lapw7中的 outwin.f 源码进行解读.

原理

在 Rydberg 单位下, 具有量子数 $\kappa$ 的大分量波函数 $u(r)=G(r)/r$ 满足¹

$$\begin{array}{rl}\frac{d}{dr}G(r)& =\frac{G(r)}{r}+M(r)F(r)\\ \frac{d}{dr}F(r)& =-\frac{F(r)}{r}+(\frac{\kappa (\kappa +1)}{r^2}\frac{1}{M(r)}-(E-V(r)))G(r)\end{array}$$

其中

$$M(r)\equiv 1+{\left(\frac{\alpha }{2}\right)}^2(E-V(r))=1+\frac{E-V(r)}{c^2}$$

$\alpha$ 为精细结构常数, 在Rydberg单位下 $\alpha=2/c$, $c$ 为光速. 为在步长$h$的对数格点上进行数值计算, 作变量替换 $r=r_0 e^x$, $\mathrm{d}r=r\mathrm{d}x$, 得到关于$x$的方程组

$$\begin{array}{rl}{G}^{\prime }& =G+MrF\\ F^{\prime }& =-F+(\frac{\kappa (\kappa +1)}{rM}-r(E-V))G\end{array}$$

方便起见, 上式中略去了 G,F,M 的变量 $r\equiv r(x)$, 撇号表示关于x求导.

源码解读

参数

行66-81

行83-91

从81行开始是对第四个及以后的格点的循环. X为-h, DRDI为rh. 其他一些的中间量与式 $\text{(???)}$ 中量的关系是

$$\begin{array}{rl}\mathrm{PHI}& =rh\frac{E-V}{c}\\ \mathrm{U}& =rhc+\mathrm{PHI}=rhc[1+\frac{E-V}{c^2}]=rhcM\\ \mathrm{Y}& =-\kappa (\kappa +1)h^2/\mathrm{U}+\mathrm{PHI}=-\frac{h}{c}[\frac{\kappa (\kappa +1)}{rM}-r(E-V)]\end{array}$$

从而可以将式 $\text{(???)}$ 写成

$$\begin{array}{rl}{G}^{\prime }& =G+\frac{\mathrm{U}}{hc}F\\ F^{\prime }& =-F-\frac{c}{h}\mathrm{Y}G\end{array}$$

令$A=G, B=F/c$, $A’=G’, B’=F’/c$, 得到

$$\begin{array}{rl}{A}^{\prime }& =A+\frac{\mathrm{U}}{h}B\\ B^{\prime }& =-B-\frac{\mathrm{Y}}{h}A\end{array}$$

行92-96

由行列式解线性方程的知识可知, 这部分求解的是这样一个矩阵方程

$$\left[\begin{array}{cc}\frac{8}{3}+X& -U\\ Y& \frac{8}{3}-X\end{array}\right]\left[\begin{array}{c}{A}_c\\ B_c\end{array}\right]=\left[\begin{array}{c}B1\\ B2\end{array}\right]$$

这里下标c表示在代码(code)中的定义. B1 和 B2 在 93 和 94 行计算, 基于Adams-Moulton算法, 因为 8/3 来自于四阶算法

$$y_{n+3}=y_{n+2}+h(\frac{9}{24}f(t_{n+3},y_{n+3})+\frac{19}{24}f(t_{n+2},y_{n+2})-\frac{5}{24}f(t_{n+1},y_{n+1})+\frac{1}{24}f(t_n,y_n))$$

其中 $f(t_{n}, y_{n})=y’_{n}$ 为第n格点上y的导数.

利用上式, $A’_{n+3}$ 可以写成

$$h{A}_{n+3}^{\prime }=\frac{8}{3}{A}_{n+3}-\frac{8}{3}{A}_{n+2}-\frac{19}{9}h{A}_{n+2}^{\prime }+\frac{5h}{9}{A}_{n+1}^{\prime }-\frac{h}{9}{A}_n^{\prime }$$

将式 $\text{(???)}$ 两边乘以$h$后, 在格点$n+3$处的表达式为

$$\begin{array}{rl}[\frac{8}{3}-h]A_{n+3}-U_{n+3}{B}_{n+3}& =\frac{8}{3}{A}_{n+2}+\frac{19}{9}h{A}_{n+2}^{\prime }-\frac{5}{9}h{A}_{n+1}^{\prime }+\frac{1}{9}h{A}_n^{\prime }\\ [\frac{8}{3}+h]B_{n+3}+Y_{n+3}{A}_{n+3}& =\frac{8}{3}{B}_{n+2}+\frac{19}{9}h{B}_{n+2}^{\prime }-\frac{5}{9}h{B}_{n+1}^{\prime }+\frac{1}{9}h{B}_n^{\prime }\end{array}$$

如果把上式看成未知数 $A_{n+3}$ 和 $B_{n+3}$ 的二元一次方程, 其系数矩阵和 $\text{(???)}$ 是相同的. 检查DGn和DFn的含义(见下面一节)以及系数Rn, 可以发现式子 $\text{(???)}$ 右侧与B1和B2也是一致的. 因此我们定义的$A$和$B$与outwin.f中的含义是一致的.

行98-103

更新最外的三个点的导数值, 以用于计算下一个格点上的B1和B2. 其中Dx1和Dx2分别用Dx2和Dx3替代, 即97, 98, 101和102行, 在格点迭代的语境下很好理解. DG3的更新表达式

DG3 = U*B(K) - X*A(K)

由式 $\text{(???)}$ 和X等于-h, 可得 $DG3=U_K{B}_K+h{A}_K=h(A_K+\frac{U_K}{h}{B}_K)=h{A}_K^{\prime }$ 因此DG3是$h$乘以A在格点K上的导数. 同理, DF3的更新表达式

DF3 = X*B(K) - Y*A(K)

意味着$DF3 = -h B_K - Y_K A_K = h B’_K$. 因此DF3是$h$乘以B在格点K上的导数.

行107-109

这一个循环对B进行了scaling, $B \to cB/2= F/2=\alpha F$. 报告¹指出F与小分量有关, 但具体关系暂未推导.

行111-113

由于$A=G=ru$, 因此VAL就是最后一个格点上的波函数值. 而SLO的数学表示

$$(h\frac{\mathrm{d}A}{hr\mathrm{d}x}-u)/r=\frac{\mathrm{d}(ru)}{r\mathrm{d}r}-\frac{u}{r}=\frac{r\mathrm{d}u}{r\mathrm{d}r}=\frac{\mathrm{d}u}{\mathrm{d}r}$$

为波函数在边界上的导数.

总结

作者阅读WIEN2k v16.1版本中例程outwin源码并试图理解其所用到的算法. outwin.f 中的outwin例程利用Adams-Moulton算法求解标量相对论方程 $\text{(???)}$, 在对数格点上得到量子数$\kappa$下的大分量波函数, 以函数乘矢径长的形式存储在$A$中.

---

1. 参考这一篇报告 ↩︎ ↩︎²