2.4 3D能带

发表于 2023-05-09 更新于 2025-12-18 分类于 DFT学习

vaspkit处理

开始之前可以做好正常的scf和band计算，确定需要第几条带。在二维能带里，我们只看沿着高对称路径的点上的能带。在三维能带图中，我们需要计算特定带的布里渊区的所有点。

新建文件夹用于3D计算，可以复制scf。vaspkit的231功能用于产生布里渊区撒点的KPOINTS。所生成的KP点包含1加权KP点和0加权KP点(注意，为了获得平滑的3D能带，用于生成KP点文件的分辨率值应该设置<0.01)。
提交计算，INCAR同scf。

后处理

vaspkit

vaspkit 232/233都可用于后处理，232可选择特定的带，233仅适用于非磁性半导体体系。完成之后会出现包含价带顶的BAND.HOMO.grd和导带底的BAND.LUMO.grd文件的能带。grep ^band PROCAR | head -n30
执行python how_to_visual.py，绘制3D能带图(python2.7环境)；这个文件在VASPKIT的安装包中，解压之后路径为：vaspkit/vaspkit.1.3.5/examples/3D_band；
特别的，如果能带不够美观，可通过能带插值进行改善。VASPKIT升级到1.3.0或更新版本，在~/.vaspkit文件中设置以下三个参数
how_to_visual.m对应于MATLAB代码， how_to_visual.py 对应于python代码，均可进行对BAND_HOMO.grd BAND_LUMO.grd KX.grd KY.grd 数据文件的读取和绘图。使用how_to_visual.m代码在MATLAB中进行绘制可实现角度拖动和放大缩放和自定义颜色等操作。

详细参考：https://mp.weixin.qq.com/s/CHnphDzUgeZA3hD0cENx5Q

https://mp.weixin.qq.com/s/4KWdZh6H3utVkF77DhuHFg

https://mp.weixin.qq.com/s/DxLl-mEZsdxq5gbft2cpzg

qvasp后处理

参考：https://mp.weixin.qq.com/s/A_Kb5IaSwY2U0UgpcsghiQ

qvasp -3dband处理提取VBM和CBM能带，
下载数据，放于origin中画图。步骤是全选数据，然后Plot->3D Surface-> Color Map Surface,点击ok
origin中精修（例如去除网格线之类的），并把价带顶和导带底merge到一起。并调整色彩美感。

脚本处理

绘图contour：就是类似等高线，如想看两个能带是否有nodal-line，可以将两个能带做差投影到一个平面。

from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
import sys
num=201
X = np.arange(0, num, 1)
Y = np.arange(0, num, 1)
X, Y = np.meshgrid(X, Y)

file=open("14u.dat",'r+')
Z=file.readlines()
Z=[Z[i].strip().split() for i in range(len(Z))]
Z=[float(Z[i][0]) for i in range(len(Z))]
Z=np.array(Z)
Z.resize(num,num)

file2=open("15u.dat",'r+')
Z2=file2.readlines()
Z2=[Z2[i].strip().split() for i in range(len(Z2))]
Z2=[float(Z2[i][0]) for i in range(len(Z2))]
Z2=np.array(Z2)
Z2.resize(num,num)
Z3 = Z- Z2print(np.min(np.abs(Z3)))

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1)
c=plt.contourf(Y,X,Z3,150,alpha=1,cmap=plt.cm.coolwarm)
plt.contour(Y,X,Z3,3,colors='black',linewidths=0.3)
plt.xlim(0,num)
plt.ylim(0,num)
plt.colorbar(c)
#plt.clabel(C, inline=True,fontsize=10)
plt.savefig("contour-3d",dpi=600)

脚本2

利用vasplit-232处理得到的数据绘制，输入对应vaspkit输出的能带index，生成2D和3D的指定能带，能带差值以及差值的绝对值，能带劈裂能的统计以及txt格式的数据格式。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
绘制VASPKIT 232功能输出的自旋劈裂图
修改版：调整3D图z轴标签位置，避免与标题重叠
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import os
import sys
import glob

# =============================================
# 1. 读取grd文件 - 修改版
# =============================================

def read_grd_file_simple(filename):
    """
    简单的.grd文件读取函数
    假设文件包含多个数值，每行可能有多个值
    """
    print(f"正在读取文件: {filename}")
    
    with open(filename, 'r') as f:
        lines = f.readlines()
    
    # 收集所有数值
    all_values = []
    for line in lines:
        line = line.strip()
        if line:
            # 分割数值，支持空格和制表符分隔
            parts = line.split()
            for part in parts:
                try:
                    all_values.append(float(part))
                except ValueError:
                    print(f"警告: 跳过非数值: {part}")
    
    # 转换为numpy数组
    data = np.array(all_values)
    
    # 尝试确定网格尺寸
    n_total = len(data)
    
    # 如果是KX或KY文件，我们可能需要从文件名推断
    if 'kx' in filename.lower() or 'ky' in filename.lower():
        # 对于KX/KY文件，尝试找到匹配的能带文件大小
        # 这里我们假设能带文件是正确格式的
        return data, 0, 0  # 暂时返回0尺寸，后面再处理
    else:
        # 对于能带文件，尝试推断尺寸
        nx = int(np.sqrt(n_total))
        ny = nx
        
        while nx * ny < n_total and nx <= ny:
            nx += 1
        
        while nx * ny > n_total and nx > 1:
            nx -= 1
        
        if nx * ny != n_total:
            ny = int(np.ceil(n_total / nx))
        
        print(f"简单读取: 数据点总数 = {n_total}, 推断网格 = {ny}×{nx}")
        
        # 重塑数据
        data_2d = data.reshape((ny, nx))
        
        return data_2d, nx, ny

# =============================================
# 2. 获取用户输入的能带编号
# =============================================

def get_band_index():
    """
    获取用户输入的能带编号
    """
    print("=" * 60)
    print("VASPKIT 232功能数据可视化")
    print("=" * 60)
    
    # 获取用户输入
    while True:
        try:
            band_input = input("请输入要分析的能带编号 (例如: 23): ").strip()
            band_index = int(band_input)
            if band_index > 0:
                return band_index
            else:
                print("错误: 能带编号必须是正整数!")
        except ValueError:
            print("错误: 请输入有效的数字!")
        except KeyboardInterrupt:
            print("\n程序已终止")
            sys.exit(0)

# =============================================
# 3. 查找文件
# =============================================

def find_files(band_index):
    """
    根据能带编号查找对应的文件
    """
    # 定义可能的文件名模式
    band_patterns = [
        f"BAND_B{band_index}_UP.grd",  # 默认格式
        f"BAND_B{band_index}_UP.GRD",  # 大写扩展名
        f"band_b{band_index}_up.grd",  # 全小写
        f"B{band_index}_UP.grd",       # 可能没有BAND前缀
        f"BAND_{band_index}_UP.grd",   # 可能没有B前缀
    ]
    
    dw_patterns = [
        f"BAND_B{band_index}_DW.grd",
        f"BAND_B{band_index}_DW.GRD",
        f"band_b{band_index}_dw.grd",
        f"B{band_index}_DW.grd",
        f"BAND_{band_index}_DW.grd",
    ]
    
    # 查找UP文件
    up_file = None
    for pattern in band_patterns:
        matches = glob.glob(pattern)
        if matches:
            up_file = matches[0]
            break
    
    # 查找DW文件
    dw_file = None
    for pattern in dw_patterns:
        matches = glob.glob(pattern)
        if matches:
            dw_file = matches[0]
            break
    
    # 查找KX和KY文件
    kx_patterns = ["KX.grd", "kx.grd", "KX.GRD", "KX.dat", "kx.dat"]
    ky_patterns = ["KY.grd", "ky.grd", "KY.GRD", "KY.dat", "ky.dat"]
    
    kx_file = None
    for pattern in kx_patterns:
        matches = glob.glob(pattern)
        if matches:
            kx_file = matches[0]
            break
    
    ky_file = None
    for pattern in ky_patterns:
        matches = glob.glob(pattern)
        if matches:
            ky_file = matches[0]
            break
    
    return kx_file, ky_file, up_file, dw_file

# =============================================
# 4. 主程序
# =============================================

def main():
    # 获取能带编号
    band_index = get_band_index()
    print(f"分析能带编号: B{band_index}")
    
    # 查找文件
    kx_file, ky_file, up_file, dw_file = find_files(band_index)
    
    # 检查文件是否存在
    files_found = []
    missing_files = []
    
    for fname, ftype in [(kx_file, "KX"), (ky_file, "KY"), (up_file, "UP"), (dw_file, "DW")]:
        if fname:
            files_found.append((fname, ftype))
            print(f"找到{ftype}文件: {fname}")
        else:
            missing_files.append(ftype)
    
    if missing_files:
        print(f"\n错误: 以下文件未找到: {', '.join(missing_files)}")
        print("请检查:")
        print(f"1. 当前目录下是否有KX.grd/KY.grd文件")
        print(f"2. 是否有BAND_B{band_index}_UP.grd和BAND_B{band_index}_DW.grd文件")
        print(f"3. 尝试不同的文件名格式 (如B{band_index}_UP.grd)")
        sys.exit(1)
    
    print(f"\n成功找到所有文件!")
    
    # 读取数据
    print("\n读取k点网格数据...")
    kx_data = read_grd_file_simple(kx_file)
    ky_data = read_grd_file_simple(ky_file)
    
    print("\n读取能带数据...")
    band_up_data = read_grd_file_simple(up_file)
    band_dw_data = read_grd_file_simple(dw_file)
    
    # 提取数据
    kx_2d, nx_kx, ny_kx = kx_data
    ky_2d, nx_ky, ny_ky = ky_data
    band_up_2d, nx_up, ny_up = band_up_data
    band_dw_2d, nx_dw, ny_dw = band_dw_data
    
    # 如果KX/KY是1D数组，尝试重塑
    if len(kx_2d.shape) == 1:
        # 推断网格尺寸
        n_total = len(kx_2d)
        nx = int(np.sqrt(n_total))
        ny = nx
        
        while nx * ny < n_total and nx <= ny:
            nx += 1
        
        while nx * ny > n_total and nx > 1:
            nx -= 1
        
        if nx * ny != n_total:
            ny = int(np.ceil(n_total / nx))
        
        print(f"重塑KX: {n_total}点 -> {ny}×{nx}网格")
        kx_2d = kx_2d.reshape((ny, nx))
        ky_2d = ky_2d.reshape((ny, nx))
    
    # 确保所有数组形状一致
    target_shape = kx_2d.shape
    print(f"\n目标网格形状: {target_shape}")
    
    # 调整能带数据形状
    if band_up_2d.shape != target_shape:
        print(f"调整UP能带形状: {band_up_2d.shape} -> {target_shape}")
        if band_up_2d.size == target_shape[0] * target_shape[1]:
            band_up_2d = band_up_2d.reshape(target_shape)
        else:
            print("警告: UP能带数据点数量不匹配!")
    
    if band_dw_2d.shape != target_shape:
        print(f"调整DW能带形状: {band_dw_2d.shape} -> {target_shape}")
        if band_dw_2d.size == target_shape[0] * target_shape[1]:
            band_dw_2d = band_dw_2d.reshape(target_shape)
        else:
            print("警告: DW能带数据点数量不匹配!")
    
    # 计算自旋劈裂
    spin_splitting = band_up_2d - band_dw_2d
    
    # 统计信息
    print("\n" + "=" * 60)
    print(f"能带 B{band_index} 自旋劈裂统计信息:")
    print(f"  最小值: {np.nanmin(spin_splitting):.6f} eV")
    print(f"  最大值: {np.nanmax(spin_splitting):.6f} eV")
    print(f"  平均值: {np.nanmean(spin_splitting):.6f} eV")
    print(f"  标准差: {np.nanstd(spin_splitting):.6f} eV")
    print(f"  绝对平均值: {np.nanmean(np.abs(spin_splitting)):.6f} eV")
    
    # =============================================
    # 5. 绘制2D映射图
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制2D映射图...")
    
    # 创建子图
    fig_2d, axes_2d = plt.subplots(2, 2, figsize=(14, 12))
    fig_2d.suptitle(f'2D Spin Splitting Maps (UP - DW) - Band B{band_index}', fontsize=16, fontweight='bold')
    
    # 定义颜色范围
    vmax = max(abs(np.nanmin(spin_splitting)), abs(np.nanmax(spin_splitting)))
    if vmax == 0:
        vmax = 0.01
    
    # 1. Spin-up 能带
    ax1 = axes_2d[0, 0]
    im1 = ax1.pcolormesh(kx_2d, ky_2d, band_up_2d, 
                        cmap='viridis', shading='auto')
    ax1.set_title(f'Spin-up Band Energy (B{band_index})', fontsize=12)
    ax1.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax1.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax1.set_aspect('equal')
    plt.colorbar(im1, ax=ax1, label='Energy (eV)')
    
    # 2. Spin-down 能带
    ax2 = axes_2d[0, 1]
    im2 = ax2.pcolormesh(kx_2d, ky_2d, band_dw_2d, 
                        cmap='viridis', shading='auto')
    ax2.set_title(f'Spin-down Band Energy (B{band_index})', fontsize=12)
    ax2.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax2.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax2.set_aspect('equal')
    plt.colorbar(im2, ax=ax2, label='Energy (eV)')
    
    # 3. 自旋劈裂 (UP - DW)
    ax3 = axes_2d[1, 0]
    im3 = ax3.pcolormesh(kx_2d, ky_2d, spin_splitting, 
                        cmap='RdBu_r', shading='auto',
                        vmin=-vmax, vmax=vmax)
    ax3.set_title(f'Spin Splitting (UP - DW) - B{band_index}', fontsize=12)
    ax3.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax3.set_aspect('equal')
    cbar3 = plt.colorbar(im3, ax=ax3, label='ΔE (eV)')
    
    # 4. 自旋劈裂绝对值
    ax4 = axes_2d[1, 1]
    im4 = ax4.pcolormesh(kx_2d, ky_2d, np.abs(spin_splitting), 
                        cmap='hot', shading='auto')
    ax4.set_title(f'Absolute Spin Splitting - B{band_index}', fontsize=12)
    ax4.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax4.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax4.set_aspect('equal')
    plt.colorbar(im4, ax=ax4, label='|ΔE| (eV)')
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_2d = f'spin_splitting_2d_maps_B{band_index}.png'
    plt.savefig(output_2d, dpi=600, bbox_inches='tight')
    print(f"2D映射图已保存为: {output_2d}")
    
    # =============================================
    # 6. 绘制3D曲面图 - 修改z轴标签位置
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制3D曲面图...")
    
    fig_3d = plt.figure(figsize=(14, 10))
    fig_3d.suptitle(f'3D Spin Splitting Surface (UP - DW) - Band B{band_index}', fontsize=16, fontweight='bold')
    
    # 创建3D子图
    # 1. 自旋劈裂3D曲面
    ax3d_1 = fig_3d.add_subplot(221, projection='3d')
    surf1 = ax3d_1.plot_surface(kx_2d, ky_2d, spin_splitting, 
                               cmap='RdBu_r', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_1.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_1.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离，避免与数字标签重叠
    ax3d_1.set_zlabel('ΔE (eV)', fontsize=10, labelpad=15)
    ax3d_1.set_title(f'Spin Splitting Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf1, ax=ax3d_1, shrink=0.5, aspect=10, label='ΔE (eV)', pad=0.12)
    
    # 2. 自旋劈裂3D线框图
    ax3d_2 = fig_3d.add_subplot(222, projection='3d')
    wire1 = ax3d_2.plot_wireframe(kx_2d, ky_2d, spin_splitting, 
                                 rstride=10, cstride=10, 
                                 color='blue', linewidth=0.5, alpha=0.7)
    ax3d_2.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_2.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_2.set_zlabel('ΔE (eV)', fontsize=10, labelpad=15)
    ax3d_2.set_title(f'Spin Splitting Wireframe - B{band_index}', fontsize=12)
    
    # 3. UP能带3D曲面
    ax3d_3 = fig_3d.add_subplot(223, projection='3d')
    surf3 = ax3d_3.plot_surface(kx_2d, ky_2d, band_up_2d, 
                               cmap='viridis', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_3.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_3.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_3.set_zlabel('E (eV)', fontsize=10, labelpad=15)
    ax3d_3.set_title(f'Spin-up Band Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf3, ax=ax3d_3, shrink=0.5, aspect=10, label='Energy (eV)', pad=0.12)
    
    # 4. DW能带3D曲面
    ax3d_4 = fig_3d.add_subplot(224, projection='3d')
    surf4 = ax3d_4.plot_surface(kx_2d, ky_2d, band_dw_2d, 
                               cmap='viridis', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_4.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_4.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_4.set_zlabel('E (eV)', fontsize=10, labelpad=15)
    ax3d_4.set_title(f'Spin-down Band Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf4, ax=ax3d_4, shrink=0.5, aspect=10, label='Energy (eV)', pad=0.12)
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_3d = f'spin_splitting_3d_surfaces_B{band_index}.png'
    plt.savefig(output_3d, dpi=600, bbox_inches='tight')
    print(f"3D曲面图已保存为: {output_3d}")
    
    # =============================================
    # 7. 绘制统计分布图
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制统计分布图...")
    
    fig_stats, (ax_hist, ax_scatter) = plt.subplots(1, 2, figsize=(14, 6))
    
    # 直方图
    spin_splitting_flat = spin_splitting.flatten()
    spin_splitting_flat = spin_splitting_flat[~np.isnan(spin_splitting_flat)]
    
    ax_hist.hist(spin_splitting_flat, bins=100, density=True, 
                alpha=0.7, color='steelblue', edgecolor='black')
    ax_hist.axvline(x=0, color='red', linestyle='--', linewidth=1.5, label='ΔE=0')
    ax_hist.set_xlabel('Spin Splitting ΔE (eV)', fontsize=12)
    ax_hist.set_ylabel('Probability Density', fontsize=12)
    ax_hist.set_title(f'Distribution of Spin Splitting Values - Band B{band_index}', fontsize=14)
    ax_hist.legend()
    ax_hist.grid(True, alpha=0.3)
    
    # 散点密度图
    scatter = ax_scatter.scatter(kx_2d.flatten(), ky_2d.flatten(), 
                               c=spin_splitting.flatten(), 
                               cmap='RdBu_r', s=1, alpha=0.7,
                               vmin=-vmax, vmax=vmax)
    ax_scatter.set_xlabel('$k_x$ (1/Å)', fontsize=12)
    ax_scatter.set_ylabel('$k_y$ (1/Å)', fontsize=12)
    ax_scatter.set_title(f'Scatter Density of Spin Splitting - Band B{band_index}', fontsize=14)
    ax_scatter.set_aspect('equal')
    plt.colorbar(scatter, ax=ax_scatter, label='ΔE (eV)')
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_stats = f'spin_splitting_statistics_B{band_index}.png'
    plt.savefig(output_stats, dpi=600, bbox_inches='tight')
    print(f"统计分布图已保存为: {output_stats}")
    
    # =============================================
    # 8. 保存数据
    # =============================================
    
    print("\n" + "=" * 60)
    print("保存数据文件...")
    
    # 保存自旋劈裂数据
    output_data = np.column_stack((
        kx_2d.flatten(),
        ky_2d.flatten(),
        band_up_2d.flatten(),
        band_dw_2d.flatten(),
        spin_splitting.flatten()
    ))
    
    # 生成带有能带编号的输出文件名
    data_filename = f'spin_splitting_data_B{band_index}.txt'
    np.savetxt(data_filename, output_data,
               header='kx(1/A) ky(1/A) E_up(eV) E_dw(eV) Delta_E(eV)',
               fmt='%.8f', comments='')
    
    print(f"数据已保存为: {data_filename}")
    
    # 保存统计摘要
    summary_filename = f'spin_splitting_summary_B{band_index}.txt'
    with open(summary_filename, 'w') as f:
        f.write("=" * 60 + "\n")
        f.write(f"Spin Splitting Analysis Summary - Band B{band_index}\n")
        f.write("=" * 60 + "\n\n")
        f.write(f"Grid dimensions: {kx_2d.shape}\n")
        f.write(f"Total k-points: {kx_2d.size}\n\n")
        f.write("Statistics:\n")
        f.write(f"  Minimum ΔE: {np.nanmin(spin_splitting):.6f} eV\n")
        f.write(f"  Maximum ΔE: {np.nanmax(spin_splitting):.6f} eV\n")
        f.write(f"  Mean ΔE: {np.nanmean(spin_splitting):.6f} eV\n")
        f.write(f"  Std Dev ΔE: {np.nanstd(spin_splitting):.6f} eV\n")
        f.write(f"  Mean |ΔE|: {np.nanmean(np.abs(spin_splitting)):.6f} eV\n")
        f.write(f"\nFiles generated:\n")
        f.write(f"  1. {output_2d}\n")
        f.write(f"  2. {output_3d}\n")
        f.write(f"  3. {output_stats}\n")
        f.write(f"  4. {data_filename}\n")
        f.write(f"  5. {summary_filename}\n")
    
    print(f"统计摘要已保存为: {summary_filename}")
    
    # =============================================
    # 9. 显示图形
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制完成! 显示所有图形...")
    print("=" * 60)
    
    plt.show()

if __name__ == "__main__":
    main()