使用Pymatgen计算PhaseDiagram

1
lines, stable_entries_plot, unstable_entries_plot = plotter.pd_plot_data
2

3
energy_above_hull = {entry.entry_id: phase_diagram.get_e_above_hull(entry) for entry in entries}
4

5
entries_plot = {}
6
for (x, y), entry in stable_entries_plot.items():
7
    entries_plot[entry.entry_id] = {
8
        'reduced_formula': entry.reduced_formula,
9
        'energy': entry.energy,
10
        'energy_per_atom': entry.energy_per_atom,
11
        'x': x,
12
        'y': y,
13
        'energy_above_hull': energy_above_hull[entry.entry_id],
14
        'stable': True,
15
    }
16

17
for entry, (x, y) in unstable_entries_plot.items():
18
    entries_plot[entry.entry_id] = {
19
        'reduced_formula': entry.reduced_formula,
20
        'energy': entry.energy,
21
        'energy_per_atom': entry.energy_per_atom,
22
        'x': x,
23
        'y': y,
24
        'energy_above_hull': energy_above_hull[entry.entry_id],
25
        'stable': False,
26
    }
27

28
df_entries_plot = pd.DataFrame.from_dict(entries_plot, orient='index').reset_index(names='material_id')
29
df_lines = pd.DataFrame(
30
    [
31
        {'x1': line[0][0], 'y1': line[1][0], 'x2': line[0][1], 'y2': line[1][1]}
32
        for line in lines
33
    ]
34
)
35

36
df_entries_plot.sort_values(['stable', 'energy_above_hull'], ascending=[False, True])

1
def to_latex(formula_pretty: str) -> str:
2
    '''
3
    将化学式转换为LaTeX格式，适用于matplotlib的文本显示。
4
    例如：ScAl2 -> $\mathrm{ScAl_{2}}$
5
    '''
6
    formula = re.sub(r'\)(\d+)', r')_{\1}', formula_pretty)
7
    formula = re.sub(r'([A-Z][a-z]?)(\d+)(?![}])', r'\1_{\2}', formula)
8
    return f'$\\mathrm{{{formula}}}$'
9

10

11
def reconstruct_triangles(points: np.ndarray, lines_df: pd.DataFrame, tol: float = 1e-6) -> list[list[int]]:
12
    '''
13
    根据相图中的线段数据重建三角形连接关系，以便在matplotlib中使用Triangulation进行等高线绘制。
14
    '''
15
    if len(points) < 3 or lines_df.empty:
16
        return []
17

18
    tree = cKDTree(points)
19
    adjacency = [set() for _ in range(len(points))]
20

21
    p1 = lines_df[['x1', 'y1']].to_numpy(float)
22
    p2 = lines_df[['x2', 'y2']].to_numpy(float)
23
    d1, idx1 = tree.query(p1, k=1)
24
    d2, idx2 = tree.query(p2, k=1)
25

26
    for point_a, point_b, dist_a, dist_b in zip(idx1, idx2, d1, d2):
27
        if dist_a <= tol and dist_b <= tol and point_a != point_b:
28
            adjacency[point_a].add(point_b)
29
            adjacency[point_b].add(point_a)
30

31
    triangles = []
32
    for i in range(len(points)):
33
        neighbors = adjacency[i]
34
        for j in (value for value in neighbors if value > i):
35
            common_neighbors = neighbors.intersection(adjacency[j])
36
            for k in (value for value in common_neighbors if value > j):
37
                triangles.append([i, j, k])
38

39
    return triangles

1
df_stable = df_entries_plot[df_entries_plot['stable']].copy()
2
df_unstable = df_entries_plot[~df_entries_plot['stable']].copy()
3

4
fig, ax = plt.subplots(figsize=(15 * cm, 10 * cm))
5

6
if len(df_stable) >= 3 and not df_lines.empty:
7
    points = df_stable[['x', 'y']].to_numpy(float)
8
    triangles = reconstruct_triangles(points, df_lines)
9

10
    if triangles:
11
        triangulation = Triangulation(points[:, 0], points[:, 1], np.asarray(triangles, dtype=int))
12
    else:
13
        triangulation = Triangulation(points[:, 0], points[:, 1])
14

15
    contour = ax.tricontourf(
16
        triangulation,
17
        df_stable['energy_per_atom'].to_numpy(float),
18
        levels=50,
19
        cmap='viridis',
20
        alpha=0.7,
21
        zorder=1,
22
        extend='both',
23
    )
24
    cbar = plt.colorbar(contour, ax=ax, shrink=0.8)
25
    cbar.set_label('Formation Energy (eV/atom)', fontsize=10)
26
    cbar.ax.tick_params(labelsize=8)
27
    half_plot_set(cbar.ax)
28

29
for _, row in df_lines.iterrows():
30
    ax.plot([row['x1'], row['x2']], [row['y1'], row['y2']], c='#555555', linewidth=0.5, zorder=9, alpha=0.9)
31

32
if not df_stable.empty:
33
    ax.scatter(
34
        df_stable['x'],
35
        df_stable['y'],
36
        marker='o',
37
        color='green',
38
        s=30,
39
        label='Stable',
40
        zorder=30,
41
    )
42

43
if not df_unstable.empty:
44
    ax.scatter(
45
        df_unstable['x'],
46
        df_unstable['y'],
47
        marker='d',
48
        c='lightgray',
49
        s=10,
50
        label='Unstable',
51
        zorder=35,
52
    )
53

54
texts = []
55
for row in df_stable.itertuples(index=False):
56
    texts.append(
57
        ax.text(
58
            row.x,
59
            row.y,
60
            to_latex(row.reduced_formula),
61
            fontsize=10,
62
            ha='center',
63
            va='center',
64
            zorder=100,
65
            weight='regular',
66
            color='black',
67
        )
68
    )
69

70
# if adjust_text is not None and texts:
71
adjust_text(texts, ax=ax)
72

73
ax.set_xticks([])
74
ax.set_yticks([])
75
ax.spines['top'].set_visible(False)
76
ax.spines['right'].set_visible(False)
77
ax.spines['bottom'].set_visible(False)
78
ax.spines['left'].set_visible(False)
79
ax.set_xlim(-0.05, 1.05)
80
ax.set_ylim(-0.05, 0.95)
81
plt.tight_layout()
82

83
fig.savefig('Al-Cu-Sc_phase_diagram.webp', dpi=1200)
84
plt.show()