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

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