作者：刘千锐，邮箱：terry_liu@pku.edu.cn
审核：陈默涵，邮箱：mohanchen@pku.edu.cn
最后更新时间：2024/10/15
对于物质的电导热导，一般可以分成电子和离子贡献，而对于电子贡献的部分可以使用 Kubo-Greenwood 的公式进行计算[Phys. Rev. B 83, 235120 (2011)]，该方法已被广泛应用于液态金属、温稠密(Warm Dense matter, WDM)等物质的计算。Kubo-Greenwood 通过直接计算含频的昂萨格系数$$L_{mn}(\omega)$$，来计算电子的电导与热导：
$$ L_{mn}(\omega)=(-1)^{m+n}\frac{2\pi e^2}{3\omega\Omega}\ \times\sum_{ij\alpha\mathbf{k}}W(\mathbf{k})\left(\frac{\epsilon_{i\mathbf{k}}+\epsilon_{j\mathbf{k}}}{2}-\mu\right)^{m+n-2} |\langle\Psi_{i\mathbf{k}}|\hat{v}\alpha|\Psi{j\mathbf{k}}\rangle|^2[f(\epsilon_{i\mathbf{k}})-f(\epsilon_{j\mathbf{k}})]\delta(\epsilon_{j\mathbf{k}}-\epsilon_{i\mathbf{k}}-\hbar\omega), $$
这里的 n,m 指标取值范围为 1, 2，
含频电导$$\sigma(\omega)=L_{11}(\omega)$$，直流电导$$\sigma_0=\lim_{\omega\to0}\sigma(\omega)$$
含频热导$$\kappa(\omega)=\frac{1}{e^2T}\left(L_{22}-\frac{L^2_{12}}{L_{11}}\right)$$, 热导$$\kappa_0=\lim_{\omega\to0}(\omega)$$
以上便是传统 KSDFT 计算的方法，通过 KSDFT 求解的波函数、本征能量，可以轻易的带入公式进行计算。
对于随机波函数密度泛函理论(sDFT)，其没有波函数、本征能量不能通过该方法进行计算。而应从原始的 Kubo 公式中出发：
电流、热流响应函数：$$C_{mn}(t)=-2\theta(t)\Im\left{Tr\left[\sqrt{\hat f}\hat{j}_m(1-\hat{f})e^{i\frac{\hat{H}}{\hbar}t}\hat{j}_ne^{-i\frac{\hat{H}}{\hbar}t}\sqrt{\hat f}\right]\right}$$, 这里$$\theta(x)$$为阶跃函数，$$\Im$$为取虚部，$$\hat{f}=\frac{1}{1+\exp(\frac{\hat{H}-\mu}{kT})}$$为费米狄拉克算符，$$\hat{j}_1=e\hat{v}$$为电流算符，$$\hat{j}_2=(\hat{v}\hat{H}+\hat{H}\hat{v})/2-\mu\hat{v}$$为热流算符，$$e^{-i\hat{H}t/\hbar}$$为演化算符。这里求迹便可通过随机波函数进行计算。
而昂萨格系数就可通过响应函数的傅里叶变换进行计算：$$L_{mn}=\frac{2e^{m+n-2}}{3\Omega\hbar\omega}\Im[\tilde{C}{mn}(\omega)]$$，这里计算$$\tilde{C}{mn}(\omega)$$可以使用不同的窗函数：
高斯型：$$\tilde{C}{mn}(\omega)=\int_0^\infty C{mn}(t)e^{-i\omega t}e^{-\frac{1}{2}s^2t^2}dt$$
洛伦兹型：$$\tilde{C}{mn}=\int_0^\infty C{mn}(t)e^{-i\omega t}e^{-\gamma t}dt$$
其分别与前面提到的高斯函数$$G(x)$$和洛伦兹函数$$L(x)$$对应。
下载链接：
https://github.com/MCresearch/abacus-user-guide/tree/master/examples/stochastic/cond_Si
cond_Si 文件夹：这是一个电子温度为 0.6 Ry（约 8.16 eV）的 单原子硅（Si）的电导热导计算算例，包含布里渊区 K 点的 KPT 文件和包含原子位置的 STRU 文件与传统的 KSDFT 计算并无区别，主要的不同在于输入文件 INPUT。为了计算电导热导，我们只需要把 cal_cond 参数打开即可，INPUT 文件如下：
INPUT_PARAMETERS
#Parameters (1.General)
suffix Si
calculation scf
esolver_type sdft
nbands 10
nbands_sto 40
nche_sto 120
seed_sto 20000
symmetry 1
kpar 1
bndpar 2
#Parameters (2.Iteration)
ecutwfc 40
scf_thr 1e-6
scf_nmax 100
#Parameters (3.Basis)
basis_type pw
#Parameters (4.Smearing)
smearing_method fd
smearing_sigma 0.6
#Parameters (5.Mixing)
mixing_type broyden
mixing_beta 0.4
#Parameters (6.Conductivity)
cal_cond 1
cond_smear 1
cond_fwhm 0.4
cond_wcut 20
cond_dw 0.02
cond_dt 0.1以上参数在 ABACUS 的线上文档中均有详细说明，这里再进行简单概述：
- cal_cond: 控制是否计算电子贡献的电导热导
- cond_smear: 设置展宽类型，1：高斯展宽（默认），2：洛伦兹展宽
-
cond_fwhm: 设置半高宽 FWHM，对于高斯展宽,
$$\mathrm{FWHM}=2\sqrt{2\ln2}s$$ ; 对于洛伦兹展宽$$\mathrm{FWHM}=2\gamma$$ . - cond_wcut: 计算频率的截断值，决定输出含频数据的频率范围, 单位: eV
- cond_dw: 频率的间隔，单位:eV
- cond_dt: 对响应函数积分时的积分间隔，原子单位，响应函数振荡过快需要较小积分间隔，一般取值在 0.01-0.2 a.u.，在屏幕输出中有推荐的取值：Recommended dt: 1.553e-01 a.u.
- cond_dtbatch: 其与 cond_dt 共同决定演化算符$$\exp(i\hat{H}* \mathrm{cond_dtbatch} * dt)$$的切比雪夫展开阶数，一般来说该值设的越大，内存消耗越多，但速度越快，当其达到一定值，速度就不会增加。当其设为 0 时，程序会自动决定 cond_dtbatch 以使切比雪夫展开阶数在 100 附近。
计算时间的说明：
- 使用 sDFT 计算电导热导还是比较慢的，计算时间正比于切比雪夫展开阶数和$$nt = t_{total}/(\mathrm{cond_dtbatch} * dt)$$
- 切比雪夫阶数由平面波截断能决定，阶数约正比于截断能
- nt 由$$t_{total}$$决定，而$$t_{total}$$由
cond_fwhm决定，其长度反比于cond_fwhm，因此该方法不能使用任意小的展宽。 - 增加$$\mathrm{cond_dtbatch} * dt$$会使切比雪夫展开阶数增加，会使 nt 减小，总体效应是其值越大，计算效率越高，但达到一定程度，效率就不再增高了。
注：结果仅供参考，以最新版本算出的结果为准
set cond_dtbatch to 69 自动设置cond_dtbatch参数
set N order of Chebyshev for KG as 113 演化算符的切比雪夫展开阶数
Calculating conductivity....
nw: 1000 ; dw: 2.000e-02 eV 总频率数，频率间隔
nt: 4862 ; dt: 1.000e-01 a.u.(ry^-1) 总时间数，时间间隔
Emin_KS(2): 1.373e+01 eV; Emax: 5.442e+02 eV; Recommended max dt: 1.553e-01 a.u. 推荐cond_dt
ik=0: (Time left 1.673e+02 s) 估计剩余计算时间
nt: 计算进度
69 138 207 276 345 414 483 552 621 690
759 828 897 966 1035 1104 1173 1242 1311 1380
1449 1518 1587 1656 1725 1794 1863 1932 2001 2070
2139 2208 2277 2346 2415 2484 2553 2622 2691 2760
2829 2898 2967 3036 3105 3174 3243 3312 3381 3450
3519 3588 3657 3726 3795 3864 3933 4002 4071 4140
4209 4278 4347 4416 4485 4554 4623 4692 4761 4830
DC electrical conductivity: 2.454778e+05 Sm^-1 直流电导值（直接线性外推的简单估计值，实际需要根据含频电导拟合外推）
Thermal conductivity: 5.968738e+02 W(mK)^-1 热导值
Lorenz number: 3.456403e+00 k_B^2/e^2 洛伦兹常数#t(a.u.) c11(t) c12(t) c22(t) decay
0 -0 -0 -0 1
0.1 -1.39981 6.61846 -38.6529 0.999999
0.2 -1.91869 6.94748 -25.4566 0.999997
0.3 -1.87453 5.46871 -15.4962 0.999993
0.4 -1.613 3.75367 -6.29163 0.999988
0.5 -1.36998 3.0228 -6.02995 0.999981
0.6 -1.1229 2.41076 -5.2121 0.999972
0.7 -0.84644 1.75119 -5.82033 0.999962
0.8 -0.529567 0.779473 -2.13103 0.99995储存响应函数的文件，各列分别为均为原子单位，时间 t, 响应函数$$C_{11}(t), C_{12}(t), C_{22}(t)$$, 窗函数
## w(eV) sigma(Sm^-1) kappa(W(mK)^-1) L12/e(Am^-1) L22/e^2(Wm^-1)
0.01 250820 604.318 -4.66249e+06 1.4392e+08
0.03 261503 619.205 -4.83345e+06 1.47997e+08
0.05 283291 649.006 -5.18082e+06 1.56228e+08
0.07 317019 693.79 -5.71545e+06 1.68767e+08
0.09 363915 753.684 -6.45332e+06 1.85835e+08
0.11 425574 828.87 -7.41519e+06 2.07723e+08
0.13 503909 919.555 -8.62601e+06 2.34774e+08
0.15 601083 1025.94 -1.01141e+07 2.67375e+08
0.17 719414 1148.17 -1.19101e+07 3.05944e+08储存昂萨格系数的文件，各列分别为频率，电导，热导，$$L_{12}(\omega), L_{22}(\omega)$$