使用pyWannier90计算局域化Wannier函数

Wannier函数简介

Wannier函数是周期性体系里和分子轨道对应的概念。很多固体物理教材都详细介绍了Wannier函数，如南京大学教材《固体理论》[1]的第八章。Wannier函数定义为Bloch函数的一个傅立叶变换：

|i,vec{L}rangle = frac{1}{sqrt{N}}sum_{vec{k}}|i,vec{k}rangle e^{-ivec{k}vec{L}} quad(1)

其中

是单胞数目（也等于

vec{k}

的数目）。

本文以量子化学LCAO的角度简要介绍Wannier函数，应当会更便于研究分子电子结构的人理解。分子轨道写成LCAO形式：

|irangle = sum_mu |murangle C_{mu i} quad(2)

周期性体系的HF是类似的，但正则轨道需满足Bloch定理。为此，首先对AO做如下变换，让AO满足Bloch定理：

|mu,vec{k}rangle = frac{1}{sqrt{N}}sum_{vec{L}}e^{ivec{k}vec{L}}|mu,vec{L}rangle quad(3)

其中

|mu,vec{L}rangle

是单胞

vec{L}

中的

|murangle

。然后对每个

vec{k}

，把正则轨道写成LCAO：

|i,vec{k}rangle = sum_{mu}|mu,vec{k}rangle C_{mu i}{(vec{k})} quad(4)

C_{mu i}(vec{k})

是每个

vec{k}

的LCAO系数，通过PBC-HF获得。

接下来，我们把Wannier函数也写成LCAO形式。由于各单胞的Wannier函数之间的平移关系，我们只需知道一个单胞的Wannier函数。在

(1)

里令

L=0

：

|i,vec{0}rangle = frac{1}{sqrt{N}}sum_{vec{k}}|i,vec{k}rangle quad(5)

其他单胞的Wannier函数通过

(5)

平移即可。把

(3)(4)

代入

(5)

，即可推出Wannier函数的LCAO：

|i,vec{0}rangle = sum_{mu,vec{L}}|mu,vec{L}rangle C_{mu i}(vec{L}) quad(6)

C_{mu i}(vec{L}) = frac{1}{N}sum_{vec{k}}C_{mu i}(vec{k})e^{ivec{k}vec{L}} quad(7)

类似分子计算里的轨道局域化，Wannier函数也可以做局域化，得到最大局域化Wannier函数（Maximally localized Wannier functions, MLWF）。[2][3][4]

程序简介和安装

Wannier90

Wannier90是计算MLWF的程序。[5]详情见https://wannier.org/。与多个电子结构程序有接口：QE，VASP，Wien2K，PySCF等。Wannier90既可以独立运行，也可以作为库使用（注意lib模式使用时只能串行）。

简要介绍Wannier90计算MLWF的原理。对每一个

vec{k}

，对Bloch函数做一个酉变换：

|i^{prime}, vec{k}rangle=sum_{i} U_{i^{prime}i}^{(vec{k})}|i,vec{k}rangle quad(8)

然后将新得到的

|i^{prime}, vec{k}rangle

按照

(5)

变为Wannier函数。酉变换

U_{i^{prime}i}^{(vec{k})}

通过最小化Wannier函数位置的方差来确定（即Boys方法）：

Omega = sum_{i} left[leftlangle i,vec{0}left|r^2right|i,vec{0}rightrangle-left|leftlangle i,vec{0}left|vec{r}right|i,vec{0}rightrangleright|^2right] quad(9)

具体的实现细节比分子复杂很多，详情见文献。[2][3][5]

pyWannier90

pyWannier90是PySCF和Wannier90的接口，调用Wannier90（lib模式），由PySCF计算的轨道生成MLWF。详情见https://github.com/hungpham2017/pyWannier90。可参考作者Hung Q. Pham周期性DMET的工作[6]。

编译和安装

本文使用wannier90-3.1.0和intel编译器。pyWannier90需要先安装好pybind11和PySCF。

解压wannier90-3.1.0和pyWannie90，用pyWannier90/src/wannier90_lib.F90替换掉wannier90-3.1.0/src/wannier90_lib.F90

把wannier90-3.1.0/config/make.inc.ifort拷贝为wannier90-3.1.0/src/make.inc，把COMMS和MPIF90两行去掉（因为lib模式只支持串行），给FCOPTS加上-fPIC，即：

代码语言：javascript复制

F90 = ifort
FCOPTS= -O2 -fPIC
LDOPTS= -O2

以lib模式编译wannier90:

代码语言：javascript复制

make && make lib

进入pyWannier90/src，修改Makefile。修改W90DIR为wannier90-3.1.0的路径，修改LIBDIR为intel的路径，例如：

代码语言：javascript复制

W90DIR =[path-to-]/wannier90-3.1.0
LIBDIR =[path-to-]/intel

编译pyWannier90:

代码语言：javascript复制

make

得到一个.so文件，例如：libwannier90.cpython-310-x86_64-linux-gnu.so

把.so文件放到[path-to]/pyscf/examples/pbc里，运行47号例子（我的PySCF版本是2.2.1）：

代码语言：javascript复制

python 47-pywannier90.py

如果成功的话会产生8个.xsf文件MLWF-[0-7].xsf，用VESTA打开可以看到Si晶体的8个轨道，例如

计算局域化Wannier函数

如下是一个例子展示pyWannier90的用法，计算金刚石价带（每个碳原子2个轨道，共4个）的MLWF：

代码语言：javascript复制

#导入pywannier90和PBC-HF计算需要的模块
import numpy as np
from pyscf import lib 
from pyscf.pbc import gto, scf 
from pyscf.pbc.tools import pywannier90

#---------------- 建立一个金刚石单胞，并定义基组和赝势--------------
cell = gto.Cell()
cell.atom='''
C     4.462500000000E-01  4.462500000000E-01  4.462500000000E-01
C    -4.462500000000E-01 -4.462500000000E-01 -4.462500000000E-01
'''
cell.basis = 'gth-szv'; cell.pseudo = 'gth-pade'
cell.a = ''' 
0.000000000000E 00   0.178500000000E 01   0.178500000000E 01
0.178500000000E 01   0.000000000000E 00   0.178500000000E 01
0.178500000000E 01   0.178500000000E 01   0.000000000000E 00
'''
cell.build()

#-------------------------- PBC-HF计算 --------------------------
kmesh = [2,2,2]
kpts = cell.make_kpts(kmesh);nkpts = len(kpts)
mf = scf.KRHF(cell, kpts)
ehf = mf.kernel()

#-------------- ---------用pyWannier90计算MLWF --------------------
nocc = int(sum(mf.mo_occ[0]))//2；ncore = 0
keywords = 
"""
exclude_bands : 5,6,7,8 
"""
# 用exclude_bands关键字把虚轨道去除（1-4为占据轨道，5-8为虚轨道）
# 如果是全电子计算的话核轨道也可以去除

w90 = pywannier90.W90(mf, cell, kmesh, nocc-ncore, other_keywords=keywords)
# 初始化W90对象，五个参数依次为：1) PBC-HF; 2) 单胞; 3)k点; 4) Wannier函数的数目; 5) 其他参数
w90.use_bloch_phases = True # 定义初猜：此选项直接使用正则轨道作为初猜
w90.kernel() # 计算MLWF，得到公式(7)中的U矩阵

mo_coeff_kpts = [w90.mo_coeff_kpts[ik][:,w90.band_included_list] for ik in range(nkpts)] #每个k点的正则占据轨道
rotated_mo_coeff_kpts = lib.einsum('kui,kji->kuj', mo_coeff_kpts, w90.U_matrix) # 酉变换，公式(7)

Ts = lib.cartesian_prod((range(kmesh[0]), range(kmesh[1]), range(kmesh[2])))
phase = 1/nkpts*np.exp(1j*2*np.pi*np.dot(Ts, w90.kpt_latt_loc.T))
C_Lui = lib.einsum("kuv,Rk->Ruv", rotated_mo_coeff_kpts, phase) # 转换为Wannier函数的LCAO系数，公式(6)

w90.plot_wf(supercell=kmesh, grid=[20,20,20]) # 绘制Wannier函数

展示其中一个Wannier函数的图像：

注意这是个比较粗糙的用法，有一些问题没考虑，如不同

vec{k}

的轨道数目可能不相等，entangled bands等[3]。更精细的用法可参阅pyscf/pbc/tools/pywannier90.py源代码及Wannier90手册。

用Wannier函数的LCAO系数可以算出密度矩阵：

P_{munu}(vec{L})=2sum_{R}sum_i C_{mu i}(vec{R})C_{nu i}(vec{R} vec{L}) quad(10)

然后可以用密度矩阵来检验Wannier函数是否正确，如：

检验与HF密度矩阵是否相等（前提是计算MLWF时没冻核）

把HF的密度矩阵变换到实空间：

mathbf{P}(vec{L})=frac{1}{N_k}sum_kmathbf{P}(vec{k})e^{ivec{k}vec{L}} quad(11)

比较

(9)

和

(10)

是否相等：

代码语言：javascript复制

iLs = lib.cartesian_prod([range(kmesh[0]), range(kmesh[1]), range(kmesh[2])])
Ls = np.dot(iLs, cell.lattice_vectors())
expkLd = np.array(np.exp(1j*np.dot(kpts, Ls.T)), order='C')
dm_L_HF = np.einsum("kuv,kL->Luv", dm_k_hf, expkLd)/nkpts
assert np.allclose(dm_L_HF, dm_L)

检验电子数

N = sum_{L}sum_{mu nu}P_{mu nu}(vec{L})S_{munu}(vec{L}) quad(12) 代码语言：javascript复制

nelec_by_dm = lib.einsum("Luv,Luv->", dm_L, ovlp_L)
assert abs(nelec_by_dm - cell.nelectron) < 1e-6

另外，注意Wannier90一些可能的风险：

迭代圈数是固定的（通过num_iter设定），而不是迭代到收敛为止。所以算完之后最好打开.wout文件检查。
以正则轨道做初猜可能很难收敛，最好给一个合适的投影做初猜。详见Wannier90手册。

参考文献

[1]:李正中，固体理论，第二版，ISBN: 978-704011576-5 [2]:PRB, 1997, 56, 12847 [3]: PRB, 2001, 65, 035109 [4]: RMP, 2012, 84, 4, 1419 [5]: J. Phys.: Condens. Matter 2020, 32, 165902 [6]: JCTC, 2020, 16, 130

使用pyWannier90计算局域化Wannier函数

Wannier函数简介

程序简介和安装

Wannier90

pyWannier90

编译和安装

计算局域化Wannier函数

参考文献

相关阅读