cp2k能带处理脚本

对于不同版本的cp2k，能带数据格式不一样，这里提供两个不同版本的处理脚本：

（1）cp2k 7.1版本：

import matplotlib.pyplot as plt
import re

# 读取文件内容
file_path = './Graphene.bs'  # 请替换为您的文件路径
with open(file_path, 'r') as file:
    content = file.readlines()

# 初始化数据存储结构
data_refined = {}
current_key_refined = None

Fermi = -27.211386245988*0.07967043079508
# 处理文件内容，提取数据
for line in content:
    if "Nr." in line:
        # 使用正则表达式精确提取Nr.后的数字作为键
        current_key_refined = float(re.search(r'Nr\.\s+(\d+)', line).group(1))
        data_refined[current_key_refined] = []
    elif current_key_refined is not None:
        try:
            # 尝试将行内的每个分段转换为浮点数，并将其添加到当前键的列表中
            y_values = [float(y) - Fermi for y in line.split()]
            if y_values:  # 确保列表不为空
                data_refined[current_key_refined].extend(y_values)
        except ValueError:
            # 如果转换失败，则跳过此行
            continue

# 绘制图表
plt.figure(figsize=(10, 6))
for x, y_values in data_refined.items():
    plt.scatter([x] * len(y_values), y_values,color='b', label=f'Nr. {x}', s=10)

plt.xlabel('Band Index')
plt.ylabel('Energy (eV)')
plt.ylim(-5,5)
plt.title('Band Structure')
#plt.legend(markerscale=2)
plt.grid(True)
# 添加垂直虚线，统一使用蓝色
# plt.axvline(x=51, color='b', linestyle='--')  # 移除label以避免重复
# plt.axvline(x=101, color='b', linestyle='--')
plt.show()

print(data_refined.items())
# 输出数据到文件
with open('band.dat', 'w') as f:
    for x, y in data_refined.items():
        f.write(f'{x}'+' ')
        for y in y_values:
            f.write(f'{y}'+' ')
        f.write('\n')

（2）cp2k 2023.2版本：

# -*- coding: UTF-8 -*-
import numpy as np
import pylab as pl
import math
def readabc():
    nKpts=201  #这里自己修改！！！！！！！！！
    eigfile = 'WS2.bs'
    fIn = open(eigfile, 'r')
    lines = fIn.readlines()
    Nkpath = 7     ##这里自己修改！！！！！！！！！填写第一个“#  Point 1      Spin 1:”的行号
    # for i in range(3,Nkpath):
    #     tmp = lines[i].rstrip().split()
    #     nKpts = nKpts+int(tmp[0])
    # print(nKpts)
    # tmp = lines[0].rstrip().split()
    nBands = 1142   #这里自己修改！！！！！！！！！
    # print(nBands)
    tmp = lines[0].rstrip().split()
    E_Fermi=27.211386245988*(-0.16547759114603)
    print(E_Fermi)
    data = np.zeros((nKpts, nBands+1), dtype=np.float64)
    # print(lines[9])
    # exit()
    for i in range(0,nKpts):
        # tmp = lines[Nkpath+1+i*nBands].rstrip().split()    #自旋向上
        # tmp1 = lines[Nkpath+3+i*4].rstrip().split()    #自旋向下
        # tmp2 = lines[Nkpath+i*4].rstrip().split()      #存储k点路径
        data[i][0]=float(i)
        print('xxx')
        for j in range(0,nBands):
            tmp=lines[Nkpath+1+j+i*(nBands+2)].rstrip().split()
            data[i][j+1]=float(tmp[1])-E_Fermi
            
        # for k in range(1,4):
        #     data[i][k]=tmp2[k]      
        # for j in range(0,nBands):
        #     # data[i][j+4]=(float(tmp[j])-E_Fermi)*27.211
        #     # data[i][j+nBands+4]=(float(tmp1[j])-E_Fermi)*27.211
        #     # data[i][j+4]=(float(tmp[j])-E_Fermi)
        #     # data[i][j+nBands+4]=(float(tmp1[j])-E_Fermi)
        #     data[i][j+4]=float(tmp[j])
        #     data[i][j+nBands+4]=float(tmp1[j])
            


    file = open("./"+"cp2k"+".dat",'w')
    # file.write(str("iK"))
    # file.write(' ')
    # file.write(str("Kx"))
    # file.write(' ')
    # file.write(str("Ky"))
    # file.write(' ')
    # file.write(str("Kz"))
    # file.write(' ')
    # for i in range(1,nBands+1):
    #     file.write(str("Band")+str(i))
    #     file.write(' ')
    # for i in range(1,nBands+1):
    #     file.write(str("Band")+str(i))
    #     file.write(' ')
    # file.write("\n")
    for i in range(0,nKpts):
        for j in range(0,nBands+1):
            file.write(str(data[i,j]))
            file.write(' ')
        file.write("\n")
    file.close()
    # s=-0.1597417
    # for i in range(0,nKpts):
    #     for j in range(0,nBands*2+4):
    #         if math.fabs(data[i][j]-s)<0.000001:
    #             print(i,j-4+1)   #因为band是从1开始的，所以要加1
readabc()