晶界能计算
直接上公式，我计算了两种晶界有关的能量。
晶界形成能是展示稳定性的，我觉得和偏聚功差别并不大。
一、TiAl-TiB2界面
| 晶界能TiAl-TiB2 | |||||||
| type | 堆叠方式 | EGB+B | A | bulk-TiAl | bulk-TiB2 | γGB | 稳定性 |
| B端 | 1ABC | -216.910 | 34.2717 | -58.51969191 | -108.543 | -0.72723 | 稳定 |
| 2ACB | -214.415 | -59.8739562 | -108.543 | -0.67108 | |||
| 3BAC | -217.876 | -58.17774813 | -108.543 | -0.74631 | |||
| 5CAB | -216.853 | -59.87370293 | -108.543 | -0.70665 | |||
| 4BCA | -216.679 | -58.17779715 | -108.543 | -0.72885 | |||
| 6CBA | -214.112 | -58.52002092 | -108.543 | -0.68641 | |||
| Ti端 | ABC | -191.996 | -58.52040265 | -99.1439 | -0.50088 | 稳定 | |
| ACB | -191.987 | -58.17828086 | -99.1443 | -0.50573 | |||
| BAC | -201.359 | -59.87465277 | -99.1439 | -0.61771 | |||
| CAB | -200.617 | -59.87458076 | -99.1438 | -0.60689 | |||
| BCA | -201.635 | -58.17821694 | -99.144 | -0.6465 | |||
| CBA | -198.692 | -58.52062103 | -99.1437 | -0.59856 | |||
根据论文，这种界面分为Ti端和B端终止，这里还对他们的堆叠方式做出了区分。（B端的123456代表文件夹里的文件序号qwq）。EGB+B是拼接后界面的能量，A用下面代码即可计算：
import numpy as np
import osdef calculate_gb_area(poscar_file, gb_plane='001'):with open(poscar_file, 'r') as f:lines = f.readlines()a = np.array([float(x) for x in lines[2].split()])b = np.array([float(x) for x in lines[3].split()])c = np.array([float(x) for x in lines[4].split()])if gb_plane == '001':a1, a2 = a, belif gb_plane == '110':a1, a2 = a, (b + c)/2else:raise ValueError("Unsupported gb_plane: use '001' or '110'")area = np.linalg.norm(np.cross(a1, a2))return area# 获取当前目录下所有.vasp文件
vasp_files = [f for f in os.listdir('.') if f.endswith('.vasp')]# 打开输出文件
with open('gb_areas.txt', 'w') as out_file:out_file.write("Filename\t\tInterface Area (Ų)\n")out_file.write("--------------------------------------\n")for vasp_file in vasp_files:try:area = calculate_gb_area(vasp_file, gb_plane='001')out_file.write(f"{vasp_file}\t{area:.4f}\n")except Exception as e:out_file.write(f"{vasp_file}\tError: {str(e)}\n")
事实证明，不同的堆叠方式并不会影响晶面面积。
接着，因为我当初用materials studio直接扩胞建模，所以我并没有晶界两端bulk的模型，所以这里直接将生成的晶界模型劈成两端，只要具有周期性那就没有问题（我也不知道为啥）。分成两部分的代码在这里：
import numpy as np
import osdef read_poscar(poscar_file):with open(poscar_file, 'r') as f:lines = f.readlines()# 读取标题和缩放因子title = lines[0].strip()scale = float(lines[1].strip())# 读取晶格向量lattice = []for i in range(2, 5):lattice.append([float(x) for x in lines[i].split()])lattice = np.array(lattice)# 读取元素种类和数量elements = lines[5].split()counts = [int(x) for x in lines[6].split()]# 确定坐标类型coord_type = lines[7].strip().lower()# 读取所有原子坐标coords = []for line in lines[8:8+sum(counts)]:coords.append([float(x) for x in line.split()[:3]])coords = np.array(coords)return title, scale, lattice, elements, counts, coord_type, coordsdef write_poscar(filename, title, scale, lattice, elements, counts, coord_type, coords):with open(filename, 'w') as f:f.write(title + '\n')f.write(f"{scale:20.16f}\n")for vec in lattice:f.write(f"{vec[0]:20.16f} {vec[1]:20.16f} {vec[2]:20.16f}\n")f.write(' '.join(elements) + '\n')f.write(' '.join(map(str, counts)) + '\n')f.write(coord_type + '\n')for coord in coords:f.write(f"{coord[0]:20.16f} {coord[1]:20.16f} {coord[2]:20.16f}\n")def get_composition_name(elements):"""根据元素组成返回适当的文件名部分"""elem_set = set(elements)if elem_set == {'Ti', 'B'}:return 'TiB2'elif elem_set == {'Ti', 'Al'}:return 'TiAl'elif elem_set == {'Ti', 'Al', 'B'}:return 'TiAlB'elif 'Ti' in elem_set and 'Al' in elem_set:return 'TiAl'elif 'Ti' in elem_set and 'B' in elem_set:return 'TiB2'else:return '_'.join(sorted(elements))def split_poscar(poscar_file, output_prefix):try:title, scale, lattice, elements, counts, coord_type, coords = read_poscar(poscar_file)# 创建原子类型列表atom_types = []for elem, count in zip(elements, counts):atom_types += [elem] * count# 分离z坐标>0.5的原子mask = coords[:, 2] > 0.5coords_z_gt_05 = coords[mask]types_z_gt_05 = [atom_types[i] for i in range(len(atom_types)) if mask[i]]coords_z_le_05 = coords[~mask]types_z_le_05 = [atom_types[i] for i in range(len(atom_types)) if not mask[i]]# 统计元素数量并保持原始顺序def count_elements(types, elements_order):counts = []for elem in elements_order:counts.append(types.count(elem))return counts# 保持原始元素顺序elements_order = elements# 统计上部分(z>0.5)counts_z_gt_05 = count_elements(types_z_gt_05, elements_order)filtered_elements_z_gt_05 = [elem for elem, cnt in zip(elements_order, counts_z_gt_05) if cnt > 0]filtered_counts_z_gt_05 = [cnt for cnt in counts_z_gt_05 if cnt > 0]# 统计下部分(z<=0.5)counts_z_le_05 = count_elements(types_z_le_05, elements_order)filtered_elements_z_le_05 = [elem for elem, cnt in zip(elements_order, counts_z_le_05) if cnt > 0]filtered_counts_z_le_05 = [cnt for cnt in counts_z_le_05 if cnt > 0]# 重新排序坐标以匹配元素顺序def reorder_coords(types, coords, elements_order):reordered_coords = []for elem in elements_order:for i, atom_type in enumerate(types):if atom_type == elem:reordered_coords.append(coords[i])return np.array(reordered_coords)coords_z_gt_05 = reorder_coords(types_z_gt_05, coords_z_gt_05, filtered_elements_z_gt_05)coords_z_le_05 = reorder_coords(types_z_le_05, coords_z_le_05, filtered_elements_z_le_05)# 生成智能文件名base_name = os.path.splitext(os.path.basename(poscar_file))[0]comp_name_upper = get_composition_name(filtered_elements_z_gt_05)comp_name_lower = get_composition_name(filtered_elements_z_le_05)output_upper = f"{output_prefix}/{base_name}_{comp_name_upper}.vasp"output_lower = f"{output_prefix}/{base_name}_{comp_name_lower}.vasp"# 创建输出目录(如果不存在)os.makedirs(output_prefix, exist_ok=True)# 写入文件write_poscar(output_upper, title, scale, lattice, filtered_elements_z_gt_05, filtered_counts_z_gt_05, coord_type, coords_z_gt_05)write_poscar(output_lower, title, scale, lattice, filtered_elements_z_le_05, filtered_counts_z_le_05, coord_type, coords_z_le_05)return True, base_name, comp_name_upper, comp_name_lowerexcept Exception as e:print(f"处理文件 {poscar_file} 时出错: {str(e)}")return False, os.path.basename(poscar_file), "", ""def process_all_vasp_files(input_dir=".", output_dir="split_results"):# 获取当前目录下所有.vasp文件和POSCAR文件vasp_files = [f for f in os.listdir(input_dir) if f.endswith('.vasp') or f == 'POSCAR']if not vasp_files:print("当前目录下未找到.vasp文件")returnprint(f"找到 {len(vasp_files)} 个VASP文件需要处理")success_count = 0for vasp_file in vasp_files:file_path = os.path.join(input_dir, vasp_file)success, filename, comp_upper, comp_lower = split_poscar(file_path, output_dir)if success:print(f"已处理: {filename} -> {comp_upper}.vasp 和 {comp_lower}.vasp")success_count += 1else:print(f"处理失败: {filename}")print(f"\n处理完成! 成功处理 {success_count}/{len(vasp_files)} 个文件")print(f"结果保存在 {output_dir} 目录中")if __name__ == "__main__":process_all_vasp_files()然后，将分成的vasp文件分别跑一下单点，得到表中两列bulk的能能量，根据公式计算即可。表中的结果看来，TiB2-TiAl的界面比单胞更稳定，而且B端终止稳定性优于Ti端终止。
二、TiAl-Ti3Al界面
这个界面模型是师兄建立的，我在界面的基础上掺杂了B。
掺杂位置主要选择了S2、S4、S6三种，计算掺杂后的总能量为E-TiAlB。break成上图的三段后，分别计算top、middle、low的能量，面积A计算方法也同上。即可得到晶界能。
| top | middle | low | A | E-TiAl | E-TiAlB | γGB晶界能 | γseg偏析能 | ||
| PST-300_gama-Ti3Al-relaxed | S2 | -106.304 | -105.44297 | -92.983 | 28.3977 | -317.774 | -324.3302 | -0.3451 | -6.55643 |
| S4 | -105.05091 | -324.0399 | -0.34689 | -6.26608 | |||||
| S6 | -104.16193 | -324.2835 | -0.36683 | -6.50971 | |||||
| PST-240_gama-Ti3Al-relaxed | S2 | -94.982667 | 28.3313 | -318.527 | -325.3086 | -0.54778 | -6.78167 | ||
| S4 | -105.46357 | -324.8508 | -0.35473 | -6.32382 | |||||
| S6 | -105.40963 | -325.0233 | -0.35873 | -6.49635 | |||||
| PST-180_gama-Ti3Al-relaxed | S2 | -104.54207 | 28.3174 | -318.457 | -323.9392 | -0.35508 | -5.48205 | ||
| S4 | -105.35087 | -324.7065 | -0.35435 | -6.2493 | |||||
| S6 | -104.31086 | -324.4751 | -0.36863 | -6.01792 | |||||
| PST-120_gama-Ti3Al-relaxed | S2 | -98.866308 | 28.323 | -318.523 | -325.2518 | -0.47838 | -6.72918 | ||
| S4 | -105.33213 | -324.7574 | -0.35551 | -6.23483 | |||||
| S6 | -104.44461 | -324.5445 | 0.032671 | -6.02188 | |||||
| PST-60_gama-Ti3Al-relaxed | S2 | -102.53886 | 28.3228 | -318.456 | -325.0617 | -0.4102 | -6.60592 | ||
| S4 | -102.82955 | -324.6415 | -0.39765 | -6.18573 | |||||
| S6 | -102.17671 | -324.8985 | -0.41371 | -6.44274 | |||||
| PST-0_gama-Ti3Al-relaxed | S2 | -104.94677 | 28.3995 | -317.681 | -323.9896 | -0.34782 | -6.3082 | ||
| S4 | -104.1744 | -324.3511 | -0.36778 | -6.6697 | |||||
| S6 | -94.60448 | -323.656 | -0.52403 | -5.97457 | |||||
| gammaTiAl-111-120-relaxed | S2 | -98.894862 | 28.0138 | -296.501 | -303.183 | -0.08926 | -6.68245 | ||
| S4 | -87.290235 | -302.0732 | -0.27658 | -5.5727 | |||||
| S6 | -86.968389 | -302.8875 | -0.29685 | -6.38696 | |||||
| gamma-gama-111-0-relaxed | S2 | -99.329188 | 27.9281 | -297.36 | -303.9991 | -0.09637 | -6.63927 | ||
| S4 | -87.756831 | -303.9992 | -0.30355 | -6.63932 | |||||
| S6 | -86.846781 | -303.6142 | -0.31295 | -6.25432 | |||||
| gamma-gama-60-relaxed | S2 | -98.886153 | 28.0247 | -296.48 | -303.1616 | -0.089 | -6.68125 | ||
| S4 | -88.365886 | -303.2357 | -0.27802 | -6.75533 | |||||
| S6 | -86.345828 | -302.9114 | -0.30827 | -6.43109 |
偏析能就更加简单了，用E-TiAlB-E-TiAl即可。由结果可见，B的掺杂使界面更加稳定。