脚本分享(四)——快速绘制能垒图

type

Post

status

Published

date

Feb 23, 2023

slug

summary

一个可以快速生成能垒图的脚本

简介：

催化计算相关的论文大多需要绘制能垒图（也称为台阶图）。我之前一直使用Origin中的两点图功能来绘制台阶图。这种方法有优点也有缺点。其中，优点是方便自定义绘图元素，例如可以随意修改图例位置和文字样式等。缺点则是速度慢，而且修改麻烦。我之前反复算了很多条路径，每个路径都要画台阶图，非常费时间。尤其是会遇到画完了，结果在组会上被老板否决的情况。因此，我写了这个脚本来解决这个问题。

这个脚本可以自动绘制出能垒图。其中，我定义了一个绘图函数，用于绘制不同状态点的能量和状态性质，例如稳态和过渡态。我还定义了一些其他参数，例如线宽、线条颜色和标签值等。

脚本的使用非常简单，只需要输入不同路径的能量和状态值即可。脚本还可以自动调整Y轴范围，使其适应不同路径的能垒差异。

最后，我还设置了一些绘图参数，例如字体、坐标轴范围和刻度大小等。绘制出来的能垒图清晰、美观，可以直接保存为图片。

这个脚本的作用不仅仅是解决了绘制能垒图的问题，同时也提高了工作效率，让我在组会和写论文时更加从容。此外，脚本的可重复性和可扩展性也让它变得非常实用。

脚本使用方法：

数据参数： 要绘制一条能量曲线，需要定义两个数组，energies数组和States数组。energies数组是不同状态点的能量，States数组表示不同状态点的状态。在States数组中，1代表稳态，0代表过渡态。如果需要绘制多条曲线，则可以设置多个参数数组。这些数据参数的正确设置是绘制能量曲线的前提，因为只有正确定义了能量和状态，绘制的曲线才能反映系统的真实情况。因此，在设置这些参数时，需要根据具体的实验或模拟数据进行合理的选择。如果数据的选择不当，绘制出来的曲线可能会偏离实际情况，导致分析和预测出现偏差。注意两个数组的长度应当一致



energies_1 = [0.000,-1,-0.1,-1.4,-2.8,-2.1,-2.7]
states_1 = [1, 1, 0, 1, 1, 0, 1]

绘图函数plot_states_energies ：该函数是用来绘制能垒图的重要函数。该函数的输入参数包括能量状态和能量值，线条宽度、线条颜色和标签值。该函数的主要作用是循环遍历所有能级，并对于稳态和过渡态分别进行处理。

对于稳态，该函数会判断当前状态是否是稳态以及前一个状态是否也是稳态。如果是，则使用虚线连接前一个稳态和当前稳态，并绘制一个实线段表示稳态；如果当前状态是稳态且能级值为0（即曲线起点），函数会绘制一个实线段表示稳态。这样，能量曲线中的稳态就会由水平线段表示出来。

对于过渡态，该函数会绘制一个平滑的曲线，通过 make_interp_spline 函数实现。该函数将过渡态的前一个稳态、当前过渡态和下一个稳态作为三个控制点，生成一条平滑的曲线。函数使用 prev_x 和 prev_y 记录前一个能级的位置，以便在画布上正确地放置每个能量点的线段。这样，能量曲线中的过渡态就会由平滑曲线表示出来。

最终，plot_states_energies 函数将生成一个能量曲线，其中稳态由水平线段表示，过渡态由平滑曲线表示。图例标签值将应用于稳态线段的标签。这个函数是能垒图绘制过程中不可或缺的一部分，通过调整其参数，可以让能垒图更符合绘图需求。。


def plot_states_energies(states, energies, linewidth, color, label):
    prev_x = 0
    prev_y = 0

    for i, (energy, state) in enumerate(zip(energies, states)):
        if state == 1:
            # 稳态
            if i > 0 and states[i-1] == 1:
                # 前一个状态也是稳态，使用虚线连接
                ax.plot([prev_x, prev_x+1], [prev_y, energy], linewidth=linewidth, linestyle='--', color=color)
                # 绘制稳态线段
                ax.plot([prev_x+1, prev_x+2], [energy, energy],linewidth = linewidth, color=color)
                prev_x += 2
                prev_y = energy
            elif energy == 0:
                ax.plot([prev_x, prev_x+1], [energy, energy],linewidth = linewidth, color=color,label=label)
                prev_x += 1
                prev_y = energy
            else:
                ax.plot([prev_x, prev_x+1], [energy, energy],linewidth = linewidth, color=color)
                prev_x += 1
                prev_y = energy
        else:
            # 过渡态
            # 过三个点画平滑曲线
            x = np.array([prev_x, prev_x+1.5, prev_x+3])
            y = np.array([prev_y, energy, energies[i+1]])
            xnew = np.linspace(x.min(),x.max(),300)
            smooth_curve = make_interp_spline(x, y, k=2)
            ynew = smooth_curve(xnew)
            plt.plot(xnew, ynew,linewidth = linewidth, color=color)
            prev_x += 3
            prev_y = energies[i+1]

其他绘图参数： 在这个部分，提供了一些常用的、需要调整的绘图参数，具体内容详见脚本注释。值得一提的是，这个脚本会根据提供的数据范围自动调整适合的坐标轴范围值。这一特性非常方便，它可以确保图形中的所有元素都能够完整呈现，并使得图像更加易于理解和解读。同时，我们也可以根据需要对这些绘图参数进行调整，以满足不同的绘图需求。这些参数的调整可以通过修改脚本中相应的数值来实现。在使用该脚本时，我们建议仔细查看这些参数，并根据需要进行相应的调整，以获得最佳的绘图效果。


font1 = {'family' : 'Arial','weight' : 'bold', 'size' : 30 }#设置字体
plt.rc('font', **font1)
ax.set_xlabel("Reaction coordinate",font = font1, labelpad = 25 )#设置X轴标签，字体，与坐标轴的距离
ax.set_ylabel("Energy (eV)", font = font1, labelpad = 25) #设置Y轴标签，字体，与坐标轴的距离
ax.set_xlim(0, 2*max_length+2) #设置X轴范围
ax.set_ylim(min_value-0.1*Energy_difference,max_value+0.1*Energy_difference)#设置y轴范围
ax.set_xticks([])#设置X轴刻度
ax.set_yticks([1.0, 0.0, -1.0, -2.0, -3.0, -4.0])#设置Y轴刻度
yticklabels = ['{:,.1f}'.format(y) for y in ax.get_yticks()]#y轴刻度保留一位小数
ax.set_yticklabels(yticklabels)
plt.tick_params(axis='both', labelsize = 28)#设置刻度大小
ax.tick_params(axis='y', which='both', length=10, width=3)#设置Y轴刻度线大小
plt.legend(labelcolor='linecolor', frameon=False, loc=1, fontsize = 27, ncol=2)#设置图例文字颜色与线条颜色一致。去除图框，右上角，字体大小，列数

# 设置图框
ax.spines['bottom'].set_linewidth(4)
ax.spines['left'].set_linewidth(4)
ax.spines['right'].set_linewidth(4)
ax.spines['top'].set_linewidth(4)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.spines['right'].set_color('black')
ax.spines['top'].set_color('black')
ax.tick_params(top=False, right=False)

完整脚本：


# By Zhang Yichen XMU
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline
import numpy as np
import sys

# 定义一个绘图函数，states数组，energies数组，线宽，线条颜色，标签值
def plot_states_energies(states, energies, linewidth, color, label):
    prev_x = 0
    prev_y = 0

    for i, (energy, state) in enumerate(zip(energies, states)):
        if state == 1:
            # 稳态
            if i > 0 and states[i-1] == 1:
                # 前一个状态也是稳态，使用虚线连接
                ax.plot([prev_x, prev_x+1], [prev_y, energy], linewidth=linewidth, linestyle='--', color=color)
                # 绘制稳态线段
                ax.plot([prev_x+1, prev_x+2], [energy, energy],linewidth = linewidth, color=color)
                prev_x += 2
                prev_y = energy
            elif energy == 0:
                ax.plot([prev_x, prev_x+1], [energy, energy],linewidth = linewidth, color=color,label=label)
                prev_x += 1
                prev_y = energy
            else:
                ax.plot([prev_x, prev_x+1], [energy, energy],linewidth = linewidth, color=color)
                prev_x += 1
                prev_y = energy
        else:
            # 过渡态
            # 过三个点画平滑曲线
            x = np.array([prev_x, prev_x+1.5, prev_x+3])
            y = np.array([prev_y, energy, energies[i+1]])
            xnew = np.linspace(x.min(),x.max(),300)
            smooth_curve = make_interp_spline(x, y, k=2)
            ynew = smooth_curve(xnew)
            plt.plot(xnew, ynew,linewidth = linewidth, color=color)
            prev_x += 3
            prev_y = energies[i+1]

# 定义状态点的能量与状态性质
# energies数组是不同状态点的能量
# States数组表明不同状态点的状态，1代表稳态，0代表过渡态
energies_1 = [0.000,-1,-0.1,-1.4,-2.8,-2.1,-2.7]
states_1 = [1, 1, 0, 1, 1, 0, 1]


max_length = 0 #获取最长的路径，
max_value = -sys.float_info.max#获取最大能量
min_value = sys.float_info.max#获取最小能量
global_vars = globals().copy()
for name, var in global_vars.items():
    if name.startswith('energies'):
        max_length = max(max_length, len(var))
        max_value = max(max_value, np.max(var))
        min_value = min(min_value, np.min(var))
        
# 获得能量极差，便于确定Y轴范围
Energy_difference = max_value - min_value

# 开始绘图
fig, ax = plt.subplots(figsize=(24,12), facecolor='white')#设置画布大小和底色
    
plot_states_energies(states_1, energies_1, 4, '#FA7F6F', r'path of Cu')

font1 = {'family' : 'Arial','weight' : 'bold', 'size' : 30 }#设置字体
plt.rc('font', **font1)
ax.set_xlabel("Reaction coordinate",font = font1, labelpad = 25 )#设置X轴标签，字体，与坐标轴的距离
ax.set_ylabel("Energy (eV)", font = font1, labelpad = 25) #设置Y轴标签，字体，与坐标轴的距离
ax.set_xlim(0, 2*max_length+2) #设置X轴范围
ax.set_ylim(min_value-0.11*Energy_difference,max_value+0.11*Energy_difference)#设置y轴范围
ax.set_xticks([])#设置X轴刻度
ax.set_yticks([1.0, 0.0, -1.0, -2.0, -3.0, -4.0])#设置Y轴刻度
yticklabels = ['{:,.1f}'.format(y) for y in ax.get_yticks()]#y轴刻度保留一位小数
ax.set_yticklabels(yticklabels)
plt.tick_params(axis='both', labelsize = 28)#设置刻度大小
ax.tick_params(axis='y', which='both', length=10, width=3)#设置Y轴刻度线大小
plt.legend(labelcolor='linecolor', frameon=False, loc=1, fontsize = 27, ncol=2)#设置图例文字颜色与线条颜色一致。去除图框，右上角，字体大小，列数

# 设置图框
ax.spines['bottom'].set_linewidth(4)
ax.spines['left'].set_linewidth(4)
ax.spines['right'].set_linewidth(4)
ax.spines['top'].set_linewidth(4)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.spines['right'].set_color('black')
ax.spines['top'].set_color('black')
ax.tick_params(top=False, right=False)

#保存图片
plt.savefig("energy.tiff",dpi=600)

绘制成品图，这里说明一下，这张图是因为示例数据太规整导致小刻度坐标恰好在边框处，正常的实验数据不会出现这种情况。

写在后面：

写这个脚本纯粹是为了更快速地组会出图，偷懒并节约时间。果然懒惰是第一生产力。。。。

新学期一周多和寒假，感觉自己啥也没干，痛苦后悔中。

振作起来，加油加油！！！

2022.2.23 于卢嘉锡253