AlmaBTE使用和FORCE_CONSTANT原子排列转换

AlmaBTE可看作是材料热导率计算程序shengBTE的升级版本，不仅包含了shengBTE的全部公共能，而且能够实现超晶格、合金等复合结构热导率的计算。AlmaBTE的主页介绍链接如下：

Home - almaBTE - First-principles thermal transport simulation

AlmaBTE的计算速度相较于shengBTE要快很多，并且shengBTE安装如果没有intel编译包等，在安装比较困难，而且计算中格点大则速度很慢，而AlmaBTE安装使用cmake，对于无intel编译器或AMD的u非常方便，实测对于单晶非复合材料计算速度非常快。

对于单晶材料的计算，同样需要输入：

1. FORCE_CONSTANTS

即2阶力常数IFC文件，格式同phononpy

2. FORCE_CONSTANTS_3RD

即3阶力常数IFC文件，格式同shengBTE的thirdorder.py脚本制作格式。

3. _metadata文件

可参考安装后test_resources/中的示例，其实就是给出针对原胞的超胞扩胞倍数等信息。

4.BORN

对于极化材料，给出介电和BORN电荷，如果没有则认为无有极化，一般此项对热导率计算的影响很小。仅对色散中LO和TO的劈裂。

5. POSCAR

原胞晶格和原子坐标的描述，参考VASP介绍。只要有上述文件，AlmaBTE就能够产生“材料.h5“的结果输出，以此.h5文件就能够进一步计算声子谐振信息、散射空间、cumulative结果、热导率等。或者通过多种材料的.h文件，就能实现超晶格等热导计算。

上述POSCA文件中，一般将原胞内同种元素的原子放在一起，因此如果利用非VASP的第一性程序计算，则可能原胞内原子的排序与之不否，导致2阶和3阶力常数文件需要调整次序。因此，这里提供一个python脚本，用于实现原胞原子次序改变后，2结合3阶力常数文件的调整。

#############################################################

#####此脚本注意python的缩进，

    #shengBTE的2阶和3阶force_constant文件中，改变原胞内原子次序后的相应文件转换
    #220318
    #输入1：所在总路径
    filepath=R"./alamode_cp2k_InP3sc441_odin3_220311"
    #输入2：原力常数文件名
    ifc2_filename=filepath+'/'+'FORCE_CONSTANTS_2ND'
    ifc3_filename=filepath+'/'+'FORCE_CONSTANTS_3RD'
    #输入3：更改后文件名
    ifc2_new_filename=filepath+'/'+'FORCE_CONSTANTS_2ND'+'_vaspnew'
    ifc3_new_filename=filepath+'/'+'FORCE_CONSTANTS_3RD'+'_vaspnew'
    #输入4：超胞原子次序更改:索引为输入2中原胞内的原子次序
    #######原胞原子排列次序 
    uc_ind0=[0,1,2,3,4,5,6,7]
    uc_ind1=[0,1,2,5,6,7,3,4]
    #输入5:2阶力常数的扩胞倍数
    sc_mul=[4,4,1]
    ########
    #扩胞次序：uc_atom*x*y*z
    uc_atom=len(uc_ind0)
    sc_atom=uc_atom*sc_mul[0]*sc_mul[1]*sc_mul[2]
    #1读取
    ifc2_file=open(ifc2_filename,'r')
    ifc2_line=ifc2_file.readlines()
    ifc2_file.close()
    ifc3_file=open(ifc3_filename,'r')
    ifc3_line=ifc3_file.readlines()
    ifc3_file.close()
    #2超胞对应list: uc_atom*x*y*z 按照标签匹配
    ifc2_ind0=[]
    ifc2_ind1=[]
    for k1 in range(uc_atom):
        for k2 in range(sc_mul[0]):
            for k3 in range(sc_mul[1]):
                for k4 in range(sc_mul[2]):
                    ifc2_ind0.append([uc_ind0[k1],k2,k3,k4])
                    ifc2_ind1.append([uc_ind1[k1],k2,k3,k4])
    #索引对应list
    sc_tran=[]
    for k1 in range(sc_atom):
        s_v=-1  #默认值
        for k2 in range(sc_atom):
            if ifc2_ind0[k1][0]==ifc2_ind1[k2][0] and ifc2_ind0[k1][1]==ifc2_ind1[k2][1] \
                and ifc2_ind0[k1][2]==ifc2_ind1[k2][2] and ifc2_ind0[k1][3]==ifc2_ind1[k2][3]:
                #
                s_v=k2
        sc_tran.append(s_v)
    #3更改2阶力常数文件
    #格式：第一行为注释128  128,接着为每个原子对 1 2 以及3行3列值    
    #初始化最终结果：2重列表，1重为原子对，2重为3*3=9个取值
    ifc2_new_val=[]
    for k1 in range(sc_atom):
        for k2 in range(sc_atom):
            ifc2_new_val.append([0 for x in range(9)])
    #读取ifc2原有文件中每条数据，放入新的ifc2_new_val值中：
    r_num=0 #当前读取块数
    r_beg=1 #当前行索引，0行为注释行
    for k1 in range(sc_atom):
        for k2 in range(sc_atom):
            r_now=r_beg+r_num*4
            templine=ifc2_line[r_now].split()
            atom1=int(templine[0])
            atom2=int(templine[1])
            templine1=ifc2_line[r_now+1]
            templine2=ifc2_line[r_now+2]
            templine3=ifc2_line[r_now+3].strip()
            r_val=templine1+templine2+templine3
            #放入新ifc2_new_val
            m1=sc_tran[k1]
            m2=sc_tran[k2]
            m_new=m1*sc_atom+m2
            ifc2_new_val[m_new]=r_val
            r_num=r_num+1
    #4输出
    ifc2_new_file=open(ifc2_new_filename,'w')
    print("%d "%(sc_atom),file=ifc2_new_file)
    k0=0
    for k1 in range(sc_atom):
        for k2 in range(sc_atom):
            print("%d %d"%(k1+1,k2+1),file=ifc2_new_file)
            print("%s"%(ifc2_new_val[k0]),file=ifc2_new_file)
            k0=k0+1
    ifc2_new_file.close()
    #5ifc3的修改
    #格式1行为总块数，空行，块(1序号2行3列超胞索引1行原胞原子索引27行4值)
    #块数目
    temp_line=ifc3_line[0].split()
    force1_b_num=int(temp_line[0])
    #读取块，转为同输入1的ifc_res列表
    ifc3_res_name1=[] #为原胞3原子索引
    ifc3_res_sc_ang1=[]  #2重列表，每个sc索引为1*3,共有2重
    ifc3_res_val1=[]  #改为string
    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].split()
            ifc3_res_name1.append([int(temp_line1[0]),int(temp_line1[1]),int(temp_line1[2])])
            temp_line2=ifc3_line[m1+1+1].split()
            temp_line3=ifc3_line[m1+1+2].split()
            ifc3_res_sc_ang1.append([[float(temp_line2[0]),float(temp_line2[1]),float(temp_line2[2])],[float(temp_line3[0]),float(temp_line3[1]),float(temp_line3[2])]])
            val_temp=''
            for m2 in range(27):
                temp_line4=ifc3_line[m1+1+4+m2]
                val_temp=val_temp+temp_line4
            ifc3_res_val1.append(val_temp)
    ###转换ifc3_res_name1内原胞索引，使用uc_ind1，默认原始原胞为uc_ind1的索引
    #输出ifc3
    ifc3_new_file=open(ifc3_new_filename,'w')
    #
    print("%d "%(force1_b_num),file=ifc3_new_file)
    print(" ",file=ifc3_new_file)
    for k1 in range(force1_b_num):
        print("%d"%(k1+1),file=ifc3_new_file)
        print("%.14f %.14f %.14f"%(ifc3_res_sc_ang1[k1][0][0],ifc3_res_sc_ang1[k1][0][1],ifc3_res_sc_ang1[k1][0][2]),file=ifc3_new_file)
        print("%.14f %.14f %.14f"%(ifc3_res_sc_ang1[k1][1][0],ifc3_res_sc_ang1[k1][1][1],ifc3_res_sc_ang1[k1][1][2]),file=ifc3_new_file)
        #uc索引输出
        uc_new=[uc_ind1[ifc3_res_name1[k1][0]-1]+1,uc_ind1[ifc3_res_name1[k1][1]-1]+1,uc_ind1[ifc3_res_name1[k1][2]-1]+1]
        print("%d  %d  %d"%(uc_new[0],uc_new[1],uc_new[2]),file=ifc3_new_file)
        #27值输出字符串
        print("%s"%(ifc3_res_val1[k1]),file=ifc3_new_file)
        #print("",file=ifc3_new_file)
    ifc3_new_file.close()