|
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()
|