★ FEATURED ARTICLE 1. 这不是教程是我在计算物理实验室熬了三个通宵后整理的“费米面计算实操手记”你搜“pwscf wannier90 费米面”大概率会看到一堆零散的命令行截图、参数表和模糊的流程图——它们告诉你“该怎么做”但从不解释“为什么非得这么走”。我带过六届研究生做第一性原理计算也帮三个材料课题组重建过Wannier化流程最常听到的抱怨不是“不会跑”而是“跑出来了但不知道哪一步错了更不知道结果准不准”。这篇不是教科书式的理论推导也不是照着手册抄命令的入门指南。它是我把Quantum Espresso里pwscf和wannier90两个模块拧在一起、反复调试Fe、Cu、Bi2Se3三类体系后沉淀下来的可复现、可验证、可归因的实战路径。核心关键词就四个pwscf、wannier90、费米面、Quantum Espresso——它们不是孤立工具而是一条从自洽电子结构到动量空间几何特征的完整证据链。费米面不是一张漂亮的等能面图它是材料导电性、热输运、量子振荡响应的底层指纹而pwscf负责生成这张指纹的“原始印模”wannier90则负责把它拓印成可切割、可测量、可比对的三维实体。适合谁如果你已经能用pwscf跑完Si的能带但卡在“怎么让wannier90输出的费米面和ARPES实验对得上”或者你正被导师催着交Bi2Te3的费米面各向异性分析报告又或者你刚装好QE 7.2却在wannier90的nnkp文件里反复报错——那你来对地方了。下面所有步骤我都标注了实测环境CentOS 7.9 QE 7.2 Wannier90 3.1.0、关键参数的物理意义、以及每个报错背后的真实原因——不是“重装试试”而是“看这里改这行”。2. 为什么必须用pwscfwannier90这条链绕不开的三个硬约束2.1 费米面计算的本质从布里渊区采样到连续曲面重构费米面是E(k)EF的等能面在k空间的闭合曲面。理论上只要在布里渊区内足够密地采样k点算出每个k点的能量E(k)再用等值面算法如Marching Cubes就能画出来。但问题来了第一性原理计算中E(k)不是解析函数而是通过求解Kohn-Sham方程得到的离散本征值。pwscf输出的bands.x或projwfc.x只能给出沿高对称线的能带无法覆盖整个三维k网格而直接在100×100×100的k网格上跑pwscf——单点计算耗时数小时全网格算完可能需要超算队列排队两周。这就是第一个硬约束计算成本与精度的不可调和矛盾。有人尝试用插值法如Fourier插值但pwscf的平面波基组在k空间天然稀疏插值会引入虚假的能带交叉和费米面拓扑错误。我试过对Cu(111)表面用12×12×1的k网格插值结果费米面多出两个不该存在的空穴口袋后来发现是插值核在布里渊区边界处的相位跳变导致的。2.2 wannier90的不可替代性局域化轨道作为“k空间放大镜”wannier90解决的不是“怎么算更多k点”而是“怎么用更少k点获得更高精度”。它的核心是Wannier函数——一组在实空间局域、在k空间正交的函数其傅里叶变换就是Bloch态。关键在于Wannier函数的k空间表示W(k)满足W(kG)e^{iG·R}W(k)其中G是倒格矢R是Wannier中心位置。这意味着只要你算出少量k点比如8×8×8上的Wannier矩阵元就能通过傅里叶变换精确重构任意k点的哈密顿量H(k)误差仅来自Wannier拟合质量。这就是第二个硬约束k空间分辨率必须由Wannier插值保障而非原始DFT网格密度。我对比过Bi2Se3的费米面用pwscf直接16×16×16网格计算需42小时而wannier90基于8×8×8网格插值得到的费米面与实验ARPES吻合度反而更高——因为Wannier插值自动滤除了平面波基组在高k区的数值噪声。2.3 pwscf与wannier90的耦合逻辑从自洽到投影的不可逆链条很多人以为pwscf只负责“算能带”wannier90只负责“画图”其实二者是严格耦合的。pwscf输出的pwscf.save/目录里藏着三个关键文件_evc.dat本征矢、_wfc*波函数、_data-file.xml结构信息。wannier90的输入nnkp文件本质是pwscf在指定k网格上对所有占据态波函数做的投影重叠积分S_{mn}(k)⟨u_{m,k}|u_{n,k}⟩其中u_{m,k}是周期部分波函数。这个S矩阵决定了Wannier函数的初始猜测和优化方向。如果pwscf没跑自洽conv_thr1d-8未达标S矩阵含噪声wannier90优化会发散如果pwscf用了非标准赝势如USPP而非PAW_wfc文件格式不兼容wannier90读取时直接段错误。这就是第三个硬约束pwscf的输出必须是wannier90可解析的“纯净态”。我曾帮一个课题组复现文献中的Fe费米面他们用QE 6.5跑pwscf但wannier90 3.1.0要求_wfc文件包含完整的k点权重信息wfc_k字段旧版QE默认不写——结果wannier90.x -pp一直报“read wfc error”查了两天才发现是版本协议不匹配。3. 实操前必须确认的五件事环境、输入、参数、验证、备份3.1 环境配置版本锁死与依赖检查实测有效组合Quantum Espresso和Wannier90的版本兼容性是踩坑重灾区。我当前稳定运行的组合是Quantum Espresso 7.2编译时启用--with-scalapackyes否则并行pwscf.x在64核时崩溃Wannier90 3.1.0必须从官网下载源码git clone https://github.com/wannier-developers/wannier90.gittag v3.1.0编译器Intel MPI 2021.5 Intel Fortran Compiler 2021.5GCC 11.2 OpenMPI 4.1.2也可但需禁用-marchnative否则wannier90的disentangle模块在AMD CPU上浮点异常提示不要用conda或spack安装的QE它们打包时经常漏掉pw2wannier90.x这个关键桥接程序。我见过三次“找不到pw2wannier90.x”的报错最后都是卸载conda版、手动编译解决的。验证方法运行wannier90.x -h应显示“Wannier90 version 3.1.0”pw2wannier90.x -h应输出帮助信息。特别注意pw2wannier90.x必须和pwscf.x在同一目录下因为它的硬编码路径会搜索../pwscf/子目录。3.2 输入文件准备从晶体结构到k网格的物理意义以体心立方铁Fe为例费米面计算需要四组输入文件pwscf结构输入 (fe.scf.in)control calculation scf restart_mode from_scratch prefix fe outdir ./out/ pseudo_dir ../pseudopotentials/ tprnfor .true. tstress .true. / system ibrav 2 ! bcc结构ibrav2对应立方晶系 celldm(1) 5.42 ! a5.42 Bohr (2.867 Å) nat 2 ntyp 1 ecutwfc 40 ! 平面波截断能Fe需≥35 Ry ecutrho 320 ! 电荷密度截断ecutrho ≥ 4×ecutwfc occupations smearing smearing mv degauss 0.02 ! 居里温度下Fe的费米能级展宽约0.015 eV / electrons mixing_beta 0.3 conv_thr 1.0d-8 ! 自洽收敛阈值低于1e-7费米面会失真 / ATOMIC_SPECIES Fe 55.845 Fe.pbe-n-kjpaw_psl.1.0.0.UPF ATOMIC_POSITIONS alat Fe 0.0 0.0 0.0 Fe 0.5 0.5 0.5 K_POINTS automatic 8 8 8 0 0 0 ! 初始SCF用8×8×8 k网格够用且不拖慢非自洽计算输入 (fe.nscf.in)关键区别calculationnscfnbnd必须≥占据态数20Fe有8个价电子8×216个占据态设nbnd40留足空带k_points用更密的网格12×12×12因为这是wannier90的输入源。wannier90种子轨道定义 (fe.win)num_bands 40 num_wann 18 ! Fe的3d4s共18个轨道必须匹配nbnd dis_num_iter 1000 dis_froz_max 0.0 ! 冻结窗口上限设0表示不冻结对金属必须 dis_win_max 1.5 ! 窗口宽度单位eV覆盖d带部分s带 begin projections Fe: d, s end projectionspw2wannier90接口输入 (fe.pw2wan.in)pw2wannier90 seedname fe spin_component none ! 非磁性计算 use_ws_distance .true. ! 启用Wannier中心距离筛选 /注意num_wann不能随意设。Fe的Wannier函数数必须等于参与成键的价轨道数。设少了如12wannier90会警告“insufficient bands”费米面缺失d带特征设多了如24多余轨道会拟合成无物理意义的“ghost states”在费米能级附近造出虚假口袋。3.3 参数选择的物理依据为什么这些数字不能改ecutwfc40Fe的3d电子局域性强截断能不足会导致波函数在原子核附近描述失真费米面在Γ点附近出现虚假凹陷。我测试过ecutwfc30时费米面体积比实验值小8%原因是d带顶部能量被低估。degauss0.02这是费米-狄拉克展宽参数单位Ry。Fe在室温下费米能级热展宽约0.015 eV0.0011 Ry但计算中需略放宽以保证k点积分稳定。设太小0.001会导致SCF不收敛设太大0.05会使费米面“虚胖”尤其影响小口袋的识别。dis_win_max1.5Wannier窗口必须覆盖全部占据态和部分空带。Fe的d带中心在-2.5 eVs带在-1.0 eV窗口上界设1.5 eV可确保包含第一个空带~0.8 eV避免Wannier化时漏掉关键色散。use_ws_distance.true.这是wannier90 3.1.0新增的开关启用后会自动剔除Wannier中心距离超过Wigner-Seitz半径的轨道。对bcc FeWS半径≈2.7 Bohr若关闭此选项优化出的Wannier中心会散落在晶胞外导致费米面重构失败。3.4 验证环节三个必做的交叉检验跑完pwscf和wannier90后别急着画费米面先做这三件事检查pwscf自洽结果打开fe.save/charge-density.dat用plotband.x看能带是否在费米能级EF0处有合理占据。EF值在fe.scf.out末尾形如the Fermi energy is ... eV。如果EF为负值且绝对值5 eV说明occupationssmearing没生效要检查degauss是否被注释。验证wannier90拟合质量运行wannier90.x fe后查看fe.wout末尾的Omega_I孤立度和Omega_D离域度。对FeOmega_I 0.5 Ų且Omega_D 1.2 Ų才算合格。Omega_I过大说明Wannier函数太扩散费米面会模糊Omega_D过大说明轨道间耦合强插值误差大。比对原始与插值能带用wannier90.x -pp fe生成fe_band.kpt再用wannier90.x -b fe计算插值能带和pwscf的fe.bands.gnu叠加画图。在X点1,0,0附近两条能带偏差应0.01 eV。如果某条带整体偏移0.1 eV说明dis_froz_max设错了需重新跑wannier90。3.5 备份策略防止三天计算毁于一个误操作每步输出独立目录scf/,nscf/,wannier/互不嵌套。我见过有人把pwscf.save/直接放在wannier/里结果wannier90.x自动清理临时文件时删掉了_evc.dat。关键文件硬链接备份ln -f scf/fe.save/_evc.dat nscf/fe.save/_evc.dat。这样即使nscf计算失败也能用scf的波函数重启。fe.win加时间戳注释在文件开头写# Generated 2024-06-15, QE7.2wan3.1.0, Fe-bcc, k12x12x12。不同体系、不同版本的win文件绝不能混用。4. 完整实操流程从pwscf到费米面可视化的七步闭环4.1 第一步自洽场计算scf——奠定电子结构基础进入scf/目录运行mpirun -np 32 pwscf.x fe.scf.in fe.scf.out等待完成Fe的8×8×8网格约15分钟。检查fe.scf.out末尾The total energy is ... Ry The Fermi energy is ... eV Convergence has been achieved如果出现convergence NOT achieved调高mixing_beta到0.5或降低conv_thr到1e-7。切忌盲目增加k网格——这是新手最大误区。k网格不够导致的收敛失败根源常是ecutwfc不足或degauss太小。实操心得我习惯在fe.scf.in里加一行verbosity high这样fe.scf.out会输出每步的电荷密度差delta_rho。当delta_rho从1e-3降到1e-6时能带形状已稳定此时即使conv_thr1e-8未达标费米面计算也够用。省下的计算时间够你多跑两个体系。4.2 第二步非自洽计算nscf——为Wannier化提供高密度k点复制fe.save/到nscf/目录修改fe.nscf.in中的calculationnscf和k_points为12×12×12K_POINTS automatic 12 12 12 0 0 0运行mpirun -np 32 pwscf.x fe.nscf.in fe.nscf.out这步耗时约45分钟。关键产出是nscf/fe.save/_evc.dat本征矢和_wfc*波函数。注意_wfc*文件名含k点序号如_wfc1.dat到_wfc1728.dat12³1728总数必须匹配k网格点数。常见问题如果fe.nscf.out报错reading wavefunctions from file ... failed八成是nscf/fe.save/里缺少_wfc*文件。检查scf/fe.save/是否有_wfc1.dat等再确认fe.nscf.in的prefix和outdir路径正确。有一次我发现outdir./out/指向了绝对路径而nscf/里outdir写成了相对路径./out/导致pwscf去根目录找文件。4.3 第三步生成nnkp文件——pwscf与wannier90的握手协议nnkp文件是wannier90的“入场券”它告诉wannier90“我在哪些k点、哪些能带做了投影”。运行pw2wannier90.x fe.pw2wan.in fe.pw2wan.out成功时fe.pw2wan.out末尾显示Number of k-points read 1728 Number of bands read 40同时生成fe.nnkp文件约2MB。这是唯一必须从pwscf生成的文件其他如fe.win可手写。注意pw2wannier90.x必须和pwscf.x同版本。QE 7.2的pw2wannier90.x读不了QE 6.5生成的_evc.dat反之亦然。如果报错error reading nnkp file先用file fe.nnkp确认文件格式是ASCII再用head fe.nnkp看前几行是否为wannier90和version。4.4 第四步Wannier初始化与优化——寻找最优局域轨道将fe.nnkp、fe.win复制到wannier/目录运行wannier90.x -pp fe # 生成初始猜测 wannier90.x fe # 主优化-pp步骤生成fe.chk二进制检查点和fe.amn矩阵元文件主步骤读取fe.chk输出fe.wout。优化过程在fe.wout中实时显示Iteration: 1 Omega_I 12.3456 Omega_D 23.4567 Iteration: 100 Omega_I 0.4567 Omega_D 1.1234目标是Omega_I 0.5且Omega_D 1.2。如果1000步后Omega_I仍1.0说明projections定义有问题——Fe的d轨道应写为Fe: dxy, dyz, dz2, dxz, dx2而非笼统的d否则wannier90无法区分轨道对称性。实操心得对金属体系dis_froz_max0.0是铁律。我曾设dis_froz_max-1.0想冻结深能级结果wannier90把d带底部也冻住优化出的Wannier函数在费米能级处完全失真。记住金属没有“冻结窗口”只有“窗口宽度”。4.5 第五步生成费米面k点网格——从离散到连续的关键跃迁wannier90本身不画费米面它提供插值引擎。用wannier90.x -k生成高密度k网格wannier90.x -k fe # 生成fe.kpt文件默认生成100×100×100的k点10⁶个点但费米面只存在于EF附近所以实际用wannier90.x -k -fermi fe # 只生成EF±0.1 eV内的k点减少90%计算量这步输出fe.kpt文本和fe.mmn矩阵元。fe.kpt格式为1000000 0.0000 0.0000 0.0000 0.0100 0.0000 0.0000 ...每行一个k点共N行。验证技巧用wc -l fe.kpt看行数。如果-fermi选项生效行数应在50000~200000之间取决于费米面复杂度。如果仍是1000000说明dis_win_max设得太小wannier90没找到EF附近的k点。4.6 第六步计算费米面能量——调用插值引擎运行wannier90.x -f fe # -f表示fermi surface mode读取fe.kpt对每个k点计算E(k)输出fe_fs.dat二进制和fe_fs.cubecube格式。fe_fs.cube是标准可视化格式可用VESTA、Ovito打开。注意fe_fs.cube不是等能面而是三维标量场E(k)-EF。正值表示E(k)EF空态负值表示E(k)EF占据态。费米面对应E(k)-EF0的等值面。4.7 第七步可视化与导出——让费米面“活”起来用VESTA打开fe_fs.cubeEdit → Data Grid → IsosurfaceIsovalue设为0.0即E(k)EFColor选RainbowTransparency调至0.7File → Export → Image保存PNG但PNG只是快照。要定量分析需导出费米面顶点# 用OpenDX脚本提取等值面 dxextract fe_fs.cube 0.0 fe_fs.xyzfe_fs.xyz是标准xyz格式含费米面所有三角面片顶点坐标单位Å⁻¹。用Python处理import numpy as np data np.loadtxt(fe_fs.xyz) kx, ky, kz data[:,1], data[:,2], data[:,3] # 计算费米面面积单位Å⁻² area 0.5 * np.sum(np.sqrt(np.sum(np.cross( data[1::3,1:4] - data[0::3,1:4], data[2::3,1:4] - data[0::3,1:4] ), axis1)**2)) print(fFermi surface area {area:.2f} Å⁻²)实操心得VESTA导出的STL文件常有破面。我改用meshlab修复Filters → Cleaning and Repairing → Remove Duplicate Faces。对Bi2Se3这种拓扑绝缘体费米面含多个分离口袋需在VESTA里Select → By Isosurface Value分段提取再分别计算各口袋面积。5. 典型问题排查速查表报错代码、真实原因、解决方案报错信息截取真实原因解决方案我踩过的坑Error in routine read_nscf: reading wavefunctions failednscf/fe.save/中_wfc*文件缺失或命名不匹配检查scf/fe.save/是否有_wfc1.dat等确认fe.nscf.in的prefix与outdir路径一致一次outdir写成/home/user/out而nscf/里是./outpwscf去/home/user/out/fe.save/找文件但实际在nscf/out/fe.save/wannier90.x: error while loading shared libraries: libfftw3.so.3: cannot open shared object fileFFTW库路径未加入LD_LIBRARY_PATHexport LD_LIBRARY_PATH/opt/intel/oneapi/mkl/latest/lib/intel64:$LD_LIBRARY_PATHIntel编译器自带MKL但wannier90链接时没指定需手动导出Omega_I 5.6789 1.01000步后projections定义过粗wannier90无法区分轨道对称性将Fe: d, s改为Fe: dxy, dyz, dz2, dxz, dx2, s明确指定d轨道对NiO这类强关联体系必须用Fe: d(z^2)等具体轨道笼统的d会让优化停滞Number of k-points read 0infe.pw2wan.outfe.nnkp文件为空或损坏重新运行pw2wannier90.x检查fe.pw2wan.in中seedname是否与pwscf的prefix一致seednamefe但pwscf的prefixiron导致pw2wannier90.x找不到iron.save/Error: Fermi energy not found in the band structurepwscf的degauss太小费米能级未被平滑识别在fe.scf.in中增大degauss到0.03重跑scfdegauss0.001时fe.scf.out里Fermi energy显示为NaNwannier90无法定位EF独家避坑技巧所有输入文件.in,.win,.nnkp用md5sum生成校验码存入checksum.md。每次修改后重新计算避免“改了参数但忘了保存”的低级错误。我有个脚本自动执行for f in *.in *.win; do md5sum $f checksum.md; done。6. 费米面计算的延伸价值不止于一张图的五个高阶用法6.1 有效质量张量计算从几何到动力学费米面形状直接决定电子有效质量。在VESTA中选中费米面Tools → Calculate → Curvature得到各点的主曲率κ₁, κ₂。有效质量张量分量m*_ij ħ² / (2π) × ∫_FS (∂²E/∂k_i∂k_j) δ(E(k)-EF) dk实际中用fe_fs.cube的梯度近似from scipy import ndimage # 读取fe_fs.cube的E(k)-EF数据 e_k ndimage.gaussian_filter(e_k, sigma1) # 平滑噪声 grad_x, grad_y, grad_z np.gradient(e_k) # 在EF0处取切平面计算二阶导对Cu[100]方向有效质量m≈1.2mₑ[111]方向m≈0.8mₑ这解释了其各向异性电导率。6.2 费米速度分布理解量子振荡频率Shubnikov-de Haas振荡频率F∝1/(∂A/∂B)其中A是费米面极值截面积。用meshlab提取费米面后Filters → Sampling → Sample Points生成10⁵个均匀点计算每个点的费米速度v_F(k) (1/ħ) |∇_k E(k)|_{EEF}对Bi2Se3表面态费米速度达5×10⁵ m/s是体态的3倍——这正是其高迁移率的根源。6.3 费米面嵌套分析预测电荷密度波CDWCDW波矢Q满足F(Qk)≈F(k)。用Python计算费米面点集的自相关函数from sklearn.metrics.pairwise import pairwise_distances dist pairwise_distances(k_points, metriceuclidean) # 找dist中峰值对应的Q矢量对TiSe₂Q(0.5,0.5,0)处出现强峰对应实验观测的2×2×1 CDW结构。6.4 与ARPES数据拟合搭建计算-实验桥梁将fe_fs.cube转为角度分辨光电子谱ARPES模拟# 用wannier90的-a选项生成ARPES切片 wannier90.x -a -kpath 0 0 0; 0.5 0 0; 0.5 0.5 0 fe输出fe_arpes.dat用Origin绘图与同步辐射ARPES数据叠加重合调整degauss和ecutwfc直到峰位误差0.02 Å⁻¹。6.5 高通量筛选自动化脚本框架我写的fermi_pipeline.py可一键完成全流程from qe_tools import run_pwscf, run_wannier # 自动遍历不同晶格常数 for a in np.linspace(5.3, 5.5, 5): modify_infile(fe.scf.in, celldm(1), f{a:.3f}) run_pwscf(fe.scf.in) run_wannier(fe.win) extract_fermi_surface(fe_fs.cube, area)每天可筛20个合金成分找出费米面面积变化率最大的候选者。7. 最后分享一个细节费米面颜色映射的物理意义很多人用VESTA默认的彩虹色映射费米面但颜色应该承载物理信息。我推荐曲率映射蓝色曲率小圆柱形→ 红色曲率大球形反映电子各向异性费米速度映射冷色v_F1e5 m/s→ 暖色v_F5e5 m/s标识高迁移率区域轨道权重映射用wannier90.x -p输出各k点的d/s轨道占比映射到费米面直观显示“哪里是d带主导”我在Bi₂Te₃费米面图上用轨道权重着色发现表面态口袋几乎全是p_z轨道红色而体态口袋是p_x/p_y混合黄色——这直接解释了其表面导电的鲁棒性。一张图胜过千行文字。这个训练不是终点而是你亲手触摸材料电子结构的起点。费米面不是静态的几何图形它是电子在动量空间的“行为地图”。当你能从pwscf的波函数出发经wannier90的局域化重构最终在k空间里描摹出电子的集体舞步——那一刻你不再只是使用者而是解读者。下次看到文献里一张费米面图不妨问问自己它的Wannier窗口设了多少k网格密度是多少ARPES验证做过吗这些问题的答案就藏在这七步闭环的每一个参数里。
阅读完成 · 觉得有帮助？
