使用pyWannier90计算局域化Wannier函数

用户7592569

发布于 2023-09-03 14:18:10

5720

发布于 2023-09-03 14:18:10

文章被收录于专栏：量子化学量子化学

Wannier函数简介

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

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

其中

是单胞数目（也等于

\vec{k}

的数目）。

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

|i\rangle = \sum_\mu |\mu\rangle C_{\mu i} \quad(2)

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

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

其中

|\mu,\vec{L}\rangle

是单胞

\vec{L}

中的

|\mu\rangle

。然后对每个

\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^{i\vec{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[\left\langle i,\vec{0}\left|r^2\right|i,\vec{0}\right\rangle-\left|\left\langle i,\vec{0}\left|\vec{r}\right|i,\vec{0}\right\rangle\right|^2\right] \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，即：

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

以lib模式编译wannier90:

make && make lib

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

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

编译pyWannier90:

make

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

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

python 47-pywannier90.py

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

计算局域化Wannier函数

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

#导入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_{\mu\nu}(\vec{L})=2\sum_{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_k\mathbf{P}(\vec{k})e^{i\vec{k}\vec{L}} \quad(11)

比较

(9)

和

(10)

是否相等：

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_{\mu\nu}(\vec{L}) \quad(12)

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