GSFE层错能计算（DFT）

1、在初始层错能结构中掺杂元素，由于层错能结构是“动态”的，所以可以选定某六个原子做八面体间隙掺杂，例如这里我选13-16个原子，底下就应该写为12-16

import os
from ase.io import read, write
from ase import Atom
import numpy as np# 输出文件夹
out_dir = "B_doped"
os.makedirs(out_dir, exist_ok=True)# 遍历当前目录所有 .vasp 文件
for file in os.listdir("."):if file.endswith(".vasp"):print(f"Processing {file} ...")atoms = read(file)# 找出 Ti 原子的索引和 z 坐标ti_indices = [i for i, a in enumerate(atoms) if a.symbol == "Ti"]ti_z = [(i, atoms[i].position[2]) for i in ti_indices]ti_z_sorted = sorted(ti_z, key=lambda x: x[1])# 选第13到第16高的 Ti 原子selected = ti_z_sorted[-18:-14]selected_indices = [i for i, _ in selected]# 输出它们的 z 坐标print("  Selected Ti atoms (index, z):")for i, z in selected:print(f"    index={i}, z={z:.6f}")# 计算这四个原子的几何中心positions = np.array([atoms[i].position for i in selected_indices])center = positions.mean(axis=0)# 添加B原子atoms.append(Atom("B", position=center))# 保存到新文件夹中out_path = os.path.join(out_dir, file)write(out_path, atoms, format="vasp", direct=True, vasp5=True)print(f"  Saved to {out_path}\n")print("✅ All .vasp files processed and saved in B_doped/")

2、固定c轴不动，间隙原子除外。分别运行下面两个代码relax.py和change.py

#!/usr/bin/env python3
"""
VASP文件坐标行修改脚本 - 修正版
解决添加Selective dynamics标记后的格式检查问题
"""import os
import globdef modify_vasp_file(filename):"""修改单个VASP文件，使用VASP能正确识别的约束标记格式"""try:# 读取原文件内容with open(filename, 'r', encoding='utf-8') as f:lines = f.readlines()# 检查文件格式if len(lines) < 9:print(f"警告: 文件 {filename} 行数不足，可能不是有效的VASP文件")return False# 重新分析文件结构，考虑可能已经添加了Selective dynamics行atom_counts_line = Nonecoord_type_line = Nonecoord_start_line = None# 查找关键行line_index = 0for i, line in enumerate(lines):stripped_line = line.strip()if i == 0:  # 第一行是元素continueelif i == 1:  # 第二行是缩放因子continueelif i >= 2 and i <= 4:  # 3-5行是晶格矢量continueelif atom_counts_line is None and len(stripped_line.split()) > 0:# 尝试解析原子数量行try:atoms = list(map(int, stripped_line.split()))if len(atoms) > 0 and all(a > 0 for a in atoms):atom_counts_line = icontinueexcept ValueError:passif atom_counts_line is not None and coord_type_line is None:# 下一行应该是坐标类型或Selective dynamicsif stripped_line in ['Cartesian', 'Direct', 'Selective dynamics']:coord_type_line = i# 如果这一行是Selective dynamics，那么下一行应该是坐标类型if stripped_line == 'Selective dynamics' and i + 1 < len(lines):next_line = lines[i + 1].strip()if next_line in ['Cartesian', 'Direct']:coord_start_line = i + 2else:coord_start_line = i + 1else:coord_start_line = i + 1breakif coord_start_line is None:print(f"警告: 无法确定文件 {filename} 的坐标起始行")return False# 计算总原子数try:atom_counts = list(map(int, lines[atom_counts_line].split()))total_atoms = sum(atom_counts)except ValueError:print(f"错误: 无法解析原子数量行: {lines[atom_counts_line]}")return False# 修改坐标行modified = Falseatoms_processed = 0for i in range(coord_start_line, len(lines)):if atoms_processed >= total_atoms:breakline = lines[i].strip()if line:  # 非空行# 移除可能存在的旧约束标记parts = line.split()if len(parts) >= 3:# 保留前三个数字（坐标）new_line = '   '.join(parts[:3])# 添加新的约束标记lines[i] = f'{new_line}   F   F   T\n'modified = Trueatoms_processed += 1if modified:# 备份原文件backup_name = filename + '.bak'if not os.path.exists(backup_name):with open(backup_name, 'w', encoding='utf-8') as f:f.writelines(lines)print(f"已创建备份文件: {backup_name}")# 写入修改后的内容with open(filename, 'w', encoding='utf-8') as f:f.writelines(lines)print(f"成功修改文件: {filename}，处理了 {atoms_processed} 个原子")return Trueelse:print(f"文件 {filename} 无需修改")return Trueexcept Exception as e:print(f"处理文件 {filename} 时出错: {str(e)}")return Falsedef check_and_fix_file_structure(filename):"""检查并修复文件结构，确保包含Selective dynamics标记返回修复后的文件内容"""try:with open(filename, 'r', encoding='utf-8') as f:lines = f.readlines()# 查找关键行索引element_line = 0scale_line = 1lattice_lines = [2, 3, 4]atom_species_line = 5  # 假设第6行是元素种类# 尝试找到原子数量行atom_counts_line = Nonefor i in range(5, min(10, len(lines))):  # 在第6-10行中查找try:atoms = list(map(int, lines[i].strip().split()))if len(atoms) > 0 and all(a > 0 for a in atoms):atom_counts_line = ibreakexcept ValueError:continueif atom_counts_line is None:print(f"错误: 无法找到原子数量行 in {filename}")return None# 检查坐标类型行coord_type_index = atom_counts_line + 1if coord_type_index >= len(lines):print(f"错误: 文件 {filename} 行数不足")return Nonecoord_line = lines[coord_type_index].strip()# 如果已经有Selective dynamics，直接返回if coord_line == 'Selective dynamics':return lines# 插入Selective dynamics行new_lines = lines[:coord_type_index]new_lines.append('Selective dynamics\n')# 如果当前行是坐标类型，保留它if coord_line in ['Cartesian', 'Direct']:new_lines.append(lines[coord_type_index])new_lines.extend(lines[coord_type_index + 1:])else:# 如果当前行不是坐标类型，可能需要调整new_lines.append('Cartesian\n')  # 默认使用Cartesiannew_lines.extend(lines[coord_type_index:])return new_linesexcept Exception as e:print(f"修复文件结构时出错: {str(e)}")return Nonedef process_vasp_files(file_list=None, directory='.'):"""处理VASP文件 - 增强版"""if file_list is None:vasp_files = glob.glob(os.path.join(directory, '*.vasp'))if not vasp_files:print(f"在目录 {directory} 中未找到.vasp文件")returnelse:vasp_files = file_listprint(f"找到 {len(vasp_files)} 个VASP文件")success_count = 0for vasp_file in vasp_files:if os.path.exists(vasp_file):print(f"\n处理文件: {vasp_file}")# 首先修复文件结构fixed_content = check_and_fix_file_structure(vasp_file)if fixed_content is not None:# 备份原文件backup_name = vasp_file + '.bak'if not os.path.exists(backup_name):with open(backup_name, 'w', encoding='utf-8') as f:with open(vasp_file, 'r', encoding='utf-8') as original:f.write(original.read())print(f"已创建备份文件: {backup_name}")# 写入修复后的结构with open(vasp_file, 'w', encoding='utf-8') as f:f.writelines(fixed_content)print("已修复文件结构")# 然后修改坐标行if modify_vasp_file(vasp_file):success_count += 1else:print(f"无法修复文件结构: {vasp_file}")else:print(f"文件不存在: {vasp_file}")print(f"\n处理完成！成功修改 {success_count}/{len(vasp_files)} 个文件")def main():"""主函数"""import argparseparser = argparse.ArgumentParser(description='VASP文件选择性动力学约束修改工具')parser.add_argument('files', nargs='*', help='要处理的.vasp文件列表')parser.add_argument('-d', '--directory', default='.', help='要搜索VASP文件的目录（默认当前目录）')args = parser.parse_args()if args.files:process_vasp_files(file_list=args.files)else:process_vasp_files(directory=args.directory)if __name__ == "__main__":main()

#!/usr/bin/env python3
import os
import globdef fix_last_line(filename):with open(filename, 'r', encoding='utf-8') as f:lines = f.readlines()if not lines:return False# 找到最后一行非空行for i in range(len(lines) - 1, -1, -1):if lines[i].strip():last_line = lines[i]breakelse:return False# 替换最后一行中 F F T 为 T T Tparts = last_line.strip().split()if len(parts) >= 6:# 只改后三个标志位parts[-3:] = ['T', 'T', 'T']lines[i] = '   '.join(parts) + '\n'with open(filename, 'w', encoding='utf-8') as f:f.writelines(lines)print(f"✅ 已修改: {filename}")return Trueelse:print(f"⚠️ 跳过: {filename} (最后一行格式不符合坐标+约束)")return Falsedef main():files = glob.glob("*.vasp")if not files:print("未找到 .vasp 文件")returnprint(f"找到 {len(files)} 个 .vasp 文件，开始处理...\n")count = 0for f in files:if fix_last_line(f):count += 1print(f"\n处理完成，共修改 {count}/{len(files)} 个文件。")if __name__ == "__main__":main()

3、runvasp，然后把结果提取出来绘制层错能曲线。注意文件命名方式要与代码中一致。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
DFT GSFE Calculator
放在包含 GSF_* 文件/文件夹 的目录下运行即可。
生成：- GSFE_DFT_<tag>.csv- GSFE_DFT_<tag>.png
"""import os, re, csv
import numpy as np
import matplotlib.pyplot as plt
from glob import glob# ================= 可配置区 =================
TAG = 'gammaTiAl_111_112'                          # 只用于输出文件名/图标题
DFT_DIRS_GLOB = 'GSF_gammaTiAl_111_112_*'  # 读取这些目录下的 OUTCAR
CSV_OUT = f'GSFE_DFT_{TAG}.csv'
PNG_OUT = f'GSFE_DFT_{TAG}.png'
J_PER_M2 = 16.02177                               # eV/Å^2 -> J/m^2
# ===========================================def read_outcar_energy(outcar_path: str) -> float:"""读取 VASP OUTCAR 最终能量(eV)。"""E = Nonewith open(outcar_path, 'r', errors='ignore') as f:for line in f:if 'free  energy   TOTEN' in line:try:E = float(line.split()[-2])  # "... =  -XXXX eV"except Exception:passif E is None:with open(outcar_path, 'r', errors='ignore') as f:for line in f:if 'energy  without entropy' in line:try:E = float(line.split()[-1])except Exception:passif E is None:raise RuntimeError(f'未在 {outcar_path} 找到能量记录。')return Edef surface_area_from_poscar(poscar_path: str) -> float:"""从 POSCAR/CONTCAR 计算表面积，假设 z 为法向，表面积 = |a×b| （Å^2）。"""from ase.io import readat = read(poscar_path)a, b, _ = at.get_cell()return float(np.linalg.norm(np.cross(a, b)))def parse_disp_from_name(name: str):"""从文件/目录名解析 (disa1, disa2)：1) ..._0.123_0.000.vasp / 目录名同理2) ..._0167(_SP)?(.vasp)?  -> (0.167, 0.0)"""m = re.search(r'_(\d+\.\d{3})_(\d+\.\d{3})', name)if m:return float(m.group(1)), float(m.group(2))m = re.search(r'_(\d{4})(?:_SP)?(?:\.vasp)?$', name)if m:d = int(m.group(1)) / 1000.0return d, 0.0raise ValueError(f'无法解析位移: {name}')def collect_dft_gsfe_data(root_glob: str):"""收集 DFT 层错能数据：返回 {'data': {(d1,d2):(E, area)}, 'E0': Eref, 'area0': area_ref}"""dmap = {}dirs = sorted(glob(root_glob))if not dirs:raise RuntimeError(f'未匹配到 DFT 目录：{root_glob}')print(f"找到 {len(dirs)} 个DFT目录")# 首先找到参考结构 (0,0)ref_dir = Nonefor d in dirs:try:d1, d2 = parse_disp_from_name(d)if abs(d1) < 1e-6 and abs(d2) < 1e-6:ref_dir = dbreakexcept ValueError:continue# 如果没有(0,0)，找位移最小的作为参考if ref_dir is None:min_disp = float('inf')for d in dirs:try:d1, d2 = parse_disp_from_name(d)disp_norm = np.sqrt(d1**2 + d2**2)if disp_norm < min_disp:min_disp = disp_normref_dir = dexcept ValueError:continueprint(f"参考结构目录: {ref_dir}")# 读取参考能量和面积ref_outcar = os.path.join(ref_dir, 'OUTCAR')ref_E0 = read_outcar_energy(ref_outcar)ref_poscar = os.path.join(ref_dir, 'CONTCAR')if not os.path.isfile(ref_poscar):ref_poscar = os.path.join(ref_dir, 'POSCAR')ref_area = surface_area_from_poscar(ref_poscar)print(f"参考能量: {ref_E0:.6f} eV")print(f"参考面积: {ref_area:.6f} Ų")# 收集所有数据点for d in dirs:try:d1, d2 = parse_disp_from_name(d)except ValueError as e:print(f"[WARN] 跳过目录 {d}: {e}")continueoutcar = os.path.join(d, 'OUTCAR')if not os.path.isfile(outcar):print(f'[WARN] 缺少 OUTCAR：{outcar}')continueE = read_outcar_energy(outcar)# 优先使用CONTCAR，其次POSCAR计算面积poscar = os.path.join(d, 'CONTCAR')if not os.path.isfile(poscar):poscar = os.path.join(d, 'POSCAR')if os.path.isfile(poscar):area = surface_area_from_poscar(poscar)else:print(f'[WARN] 目录 {d} 缺少CONTCAR/POSCAR，使用参考面积')area = ref_areadmap[(d1, d2)] = (E, area)return {'data': dmap, 'E0': ref_E0, 'area0': ref_area}def calculate_gsfe(dft_data):"""计算层错能并生成结果"""keys = sorted(dft_data['data'].keys(), key=lambda k: (k[0], k[1]))if not keys:raise RuntimeError('没有有效的DFT数据点。')rows = []for k in keys:d1, d2 = kE_dft, area_dft = dft_data['data'][k]g_dft = (E_dft - dft_data['E0']) / area_dft * J_PER_M2rows.append((d1, d2, E_dft, area_dft, g_dft))return rowsdef save_gsfe_csv(rows, csv_out):"""保存CSV结果"""with open(csv_out, 'w', newline='') as f:w = csv.writer(f)w.writerow(['disa1', 'disa2', 'Energy_DFT(eV)', 'Area(Å^2)', 'GSFE_DFT(J/m^2)'])for r in rows:w.writerow([f'{r[0]:.3f}', f'{r[1]:.3f}', f'{r[2]:.6f}', f'{r[3]:.6f}', f'{r[4]:.6f}'])print(f'[OK] 写出：{csv_out}（{len(rows)} 点）')def plot_gsfe_curve(rows, png_out, tag):"""绘制层错能曲线"""# 判断主要位移方向use_a1 = sum(1 for r in rows if abs(r[1]) < 1e-9) >= len(rows) * 0.9use_a2 = sum(1 for r in rows if abs(r[0]) < 1e-9) >= len(rows) * 0.9if use_a1:xs = [r[0] for r in rows]xlabel = 'Displacement along a1'title_suffix = 'along a1'elif use_a2:xs = [r[1] for r in rows]xlabel = 'Displacement along a2'title_suffix = 'along a2'else:# 二维情况，用位移模量xs = [np.sqrt(r[0]**2 + r[1]**2) for r in rows]xlabel = 'Displacement magnitude'title_suffix = '2D'gsfe = [r[4] for r in rows]# 按x值排序sorted_data = sorted(zip(xs, gsfe))xs_sorted, gsfe_sorted = zip(*sorted_data)plt.figure(figsize=(10, 6))plt.plot(xs_sorted, gsfe_sorted, 'bo-', linewidth=2, markersize=6, label='DFT')plt.xlabel(xlabel)plt.ylabel('GSFE (J/m²)')plt.title(f'DFT Generalized Stacking Fault Energy: {tag} {title_suffix}')plt.grid(True, alpha=0.3)plt.legend()plt.tight_layout()plt.savefig(png_out, dpi=300, bbox_inches='tight')plt.show()print(f'[OK] 写出：{png_out}')def print_statistics(rows):"""打印统计信息"""gsfe_values = [r[4] for r in rows]print("\n=== DFT层错能统计 ===")print(f"数据点数: {len(rows)}")print(f"层错能范围: {min(gsfe_values):.3f} ~ {max(gsfe_values):.3f} J/m²")print(f"最大层错能: {max(gsfe_values):.3f} J/m²")print(f"最小层错能: {min(gsfe_values):.3f} J/m²")print(f"平均层错能: {np.mean(gsfe_values):.3f} J/m²")# 找到不稳定层错能和本征层错能zero_point = next((r for r in rows if abs(r[0]) < 1e-6 and abs(r[1]) < 1e-6), None)if zero_point:print(f"完美晶体层错能: {zero_point[4]:.3f} J/m² (应为0)")# 不稳定层错能通常是最大值max_idx = np.argmax(gsfe_values)unstable_point = rows[max_idx]print(f"不稳定层错能: {unstable_point[4]:.3f} J/m² at ({unstable_point[0]:.3f}, {unstable_point[1]:.3f})")def main():"""主函数"""print("开始计算DFT层错能...")# 收集DFT数据dft_data = collect_dft_gsfe_data(DFT_DIRS_GLOB)# 计算层错能gsfe_rows = calculate_gsfe(dft_data)# 保存结果save_gsfe_csv(gsfe_rows, CSV_OUT)# 绘制曲线plot_gsfe_curve(gsfe_rows, PNG_OUT, TAG)# 打印统计信息print_statistics(gsfe_rows)print("\n[完成] DFT层错能计算完毕！")if __name__ == '__main__':main()