FORCE_CONSTANTS中3阶力常数大小与原子间距的分析脚本

原创 2022-04-06 09:47:52 2.8k 阅读 AI写同款

odin_linux

93粉丝29篇原创

关注

从超胞构建到谱图绘制：VASP结合Phonopy的声子谱全流程解析本文详细解析了使用VASP结合Phonopy进行声子谱计算的全流程，从超胞构建到谱图绘制的关键步骤。重点介绍了VASP参数设置、力常数提取和声子谱计算技巧，帮助研究者避免常见错误并提升计算效率。适用于材料科学领域的声子谱分析需求。阅读详情

AI摘要

该博客介绍了一个Python脚本，用于读取和分析shengBTE所需的FORCE_CONSTANTS_3RD文件，展示力常数与原子间距的关系，并根据设定的新截断半径生成新的力常数文件。通过散点图描绘了力常数的最大值与原子间距的最大距离之间的关系，有助于理解热导率随截断半径的变化。

shengBTE所需的3阶力常数文件格式为FORCE_CONSTANTS_3RD, 想要得知所有力常数数值与原子间距的变化关系，或者输出更小截断半径下新的FORCE_CONSTANTS_3RD以验证热导率随截断半径的变化，需要分析FORCE_CONSTANTS_3RD内数据，这里提供一个PYTHON脚本，能够读取FORCE_CONSTANTS_3RD文件，并绘制力常数强度和原子间距的散点图，以及输出新的截断半径下的力常数文件，供大家参考，注意python的缩进。

#####################para_constants
bohr2ang=0.529177210
rhy2ev=13.60569193 #1 Ry = 1/2 Ha = 13.60569193 eV
#####################################
import os
import numpy as np
import matplotlib.pyplot as plt
######################
#输入1:力常数文件所在目录
filepath=R"填入总路径"
ifc3_std_filename=filepath+'/'+'FORCE_CONSTANTS_3RD'
#输入2：为了计算间距，需要输入uc内原子坐标索引(为XYZ下单位angstrom),同vesta用.xyz格式
uc_xyz=[[5.654830126,1.395610194,1.885231820], \
[7.273416028,4.199308702,1.885231820],\
[4.036048945,4.199421446,1.885231820],\
[5.654700060,5.133874941 , 3.114744400],\
[4.035983972,2.330176433, 3.114744400],\
[7.273416149 , 2.330176433 , 3.114744400],\
[1.885025911, 1.088320217 , 2.486000000],\
[9.424374089 , 5.441164917 , 2.514000000] \
]
#输入3：设置新的cutoff,输出相应的3阶力常数
ifc3_new_cut=6.0 #新的截断半径
ifc3_cut_filename=filepath+'/'+'FORCE_CONSTANTS_3RD_cut' #新的力常数文件
#############################
ifc3_std=[]
#3阶力常数的读取
ifc3_file=open(ifc3_std_filename,'r')
ifc3_line=ifc3_file.readlines()
ifc3_file.close()
#
#格式1行为总块数，空行，块(1序号2行3列超胞索引1行原胞原子索引27行4值)，每个块前有空行
#块数目
temp_line=ifc3_line[0].split()
ifc3_b_num=int(temp_line[0])
#读取块，转为同输入1的ifc_res列表
ifc3_sc=[] #字符串，1重为所有块，2重为2个相对晶格矢的行字符串，两个字符串
ifc3_name=[] #为原胞3原子索引，字符串
ifc3_val=[] #27数值
for m1 in range(len(ifc3_line)-1): #排除第一行的总数目
temp_line=ifc3_line[m1+1].split()
if temp_line and len(temp_line)==1:
#块的开始
temp_line1=ifc3_line[m1+1+3].strip()
ifc3_name.append(temp_line1)
temp_line2=ifc3_line[m1+1+1].strip()
temp_line3=ifc3_line[m1+1+2].strip()
ifc3_sc.append([temp_line2,temp_line3])
val_temp=[]
for m2 in range(27):
temp_line4=ifc3_line[m1+1+4+m2].split()
val_temp.append(float(temp_line4[3]))
ifc3_val.append(val_temp)
###
ifc3_std=[ifc3_b_num,ifc3_sc,ifc3_name,ifc3_val]
######################分析标准文件间距与
#ifc3_std=[ifc3_b_num,ifc3_sc,ifc3_name,ifc3_val]
ifc3_len=[] #1重为块索引，2重为2和3原子索引以及最大值，共3值
ifc3_max_val=[]
for k1 in range(ifc3_std[0]):
temp1=ifc3_std[1][k1][0].split()
temp2=ifc3_std[1][k1][1].split()
a2_sc=[float(temp1[0]),float(temp1[1]),float(temp1[2])]
a3_sc=[float(temp2[0]),float(temp2[1]),float(temp2[2])]
temp3=ifc3_std[2][k1].split()
a1_ind=int(temp3[0])-1
a2_ind=int(temp3[1])-1
a3_ind=int(temp3[2])-1
a21_xyz=[uc_xyz[a2_ind][0]-uc_xyz[a1_ind][0],uc_xyz[a2_ind][1]-uc_xyz[a1_ind][1],uc_xyz[a2_ind][2]-uc_xyz[a1_ind][2]]
a31_xyz=[uc_xyz[a3_ind][0]-uc_xyz[a1_ind][0],uc_xyz[a3_ind][1]-uc_xyz[a1_ind][1],uc_xyz[a3_ind][2]-uc_xyz[a1_ind][2]]
d21=[a2_sc[0]+a21_xyz[0],a2_sc[1]+a21_xyz[1],a2_sc[2]+a21_xyz[2]]
d31=[a3_sc[0]+a31_xyz[0],a3_sc[1]+a31_xyz[1],a3_sc[2]+a31_xyz[2]]
d21_len=(d21[0]**2+d21[1]**2+d21[2]**2)**0.5
d31_len=(d31[0]**2+d31[1]**2+d31[2]**2)**0.5
d_max=max([abs(d21_len),abs(d31_len)])
d_now=[d21_len,d31_len,d_max]
ifc3_len.append(d_now)
###找到每个块27值的最大值
val=[abs(x) for x in ifc3_std[3][k1]]
val_max=max(val)
ifc3_max_val.append(val_max)
###
######绘图
x_d=[x[2] for x in ifc3_len]
y_d=ifc3_max_val
cut_res=[x_d,y_d]
plt.scatter(x_d,y_d,s=5,c="r")
plt.xlabel('std_res:x_max(d12 vs d13) VS y_max(27-val)')
plt.show()
############筛选ifc3_std的三阶力常数，用ifc3_new_cut
#ifc3_std=[ifc3_b_num,ifc3_sc,ifc3_name,ifc3_val]
ifc3_cut_num=0
ifc3_sc_cut=[]
ifc3_name_cut=[]
ifc3_val_cut=[]
#ifc3_len--------ifc3_max_val
for k1 in range(ifc3_std[0]):
if ifc3_len[k1][2]<=ifc3_new_cut:
ifc3_sc_cut.append(ifc3_std[1][k1])
ifc3_name_cut.append(ifc3_std[2][k1])
ifc3_val_cut.append(ifc3_std[3][k1])
ifc3_cut_num=ifc3_cut_num+1
ifc3_cut=[ifc3_cut_num,ifc3_sc_cut,ifc3_name_cut,ifc3_val_cut]
#输出3阶力常数
ifc3_new_file=open(ifc3_cut_filename,'w')
#ifc3_tot=[ifc3_b_num,ifc3_sc,ifc3_name,ifc3_val]
print(" %d "%(ifc3_cut[0]),file=ifc3_new_file)
#print(" ",file=ifc3_new_file)
for kk1 in range(ifc3_cut[0]):
print(" ",file=ifc3_new_file)
print(" %d"%(kk1+1),file=ifc3_new_file)
print("%s"%(ifc3_cut[1][kk1][0]),file=ifc3_new_file)
print("%s"%(ifc3_cut[1][kk1][1]),file=ifc3_new_file)
print(" %s"%(ifc3_cut[2][kk1]),file=ifc3_new_file)
#27值
pp0=0
for pp1 in range(3):
for pp2 in range(3):
for pp3 in range(3):
be_str=' '+repr(pp1+1)+' '+repr(pp2+1)+' '+repr(pp3+1)+' '
be_val=ifc3_cut[3][kk1][pp0]
pp0=pp0+1
print("%s %.14f"%(be_str,be_val),file=ifc3_new_file)
ifc3_new_file.close()