通过materials project 中的pymatgen计算形成能formation energy和e above hull

if you ->会anaconda & 会 Python：then 拷贝代码，安装代码中的import包，运行代码->完事儿

前言

材料数据库Materials Project应该是大家用得最多的数据库，其中包含数十万的DFT计算数据。而其团队为数据库配套开发的python包Pymatgen也是功能丰富，本文将详细介绍如何利用pymatgen计算指定组分的energy above hull 以及形成能formation energy。
其实pymatgen的官网有给出利用materials project中的数据计算energy above hull的实例，但是有这么一种情况，我们做完DFT计算后，拥有若干化学式-静态总能（ev）这样的数据，要如何计算这些数据的e above hull 和形成能呢？毕竟官网给的实例，以及youtube上的实例是直接计算MP数据库中某种组分的相图和e above hull。其实通过看pymatgen的文档和e above hull 计算的源码可以发现，用自己计算的DFT数据也是可以计算其对应的e above hull 的。
其原理很简单，就是用总能-各元素基态能的加和，计算得到形成能，然后结合MP上的子成分形成能画相图，然后就可以计算e above hull了。举个例子，我手上有这么一组数据，CsPbI3-静态总能Etot，要计算CsPbI3的e above hull. 先从MP上下载CsPbI3包含的所有数据，如Cs，Pb,I,CsI,PbI2等，包含化学式和各种性质，每组数据被称为一个entry，如图1。可以看到每个entry有总能，计算细节等信息，其实还有形成能这些信息，没显示出来。有了这些子成分的信息（形成能，总能）后，就可以计算CsPbI3的形成能了，在子成分Cs，Pb，I中找能量最低的那个来计算就行了。有了形成能，也就可以画相图，计算e above hull了。如图2，先新建一个包含Cs4Pb4I12信息的entry，然后把这个entry和之前的子成分entries合并在一起，就可以画相图了。那个no ID的就是我们自己导入的数据，形成能是-1.226eV/atom，e above hull是0.039eV/atom。

图1 MP上Download数据

图2 CsPbI3相图计算

存在的问题

显然，这个CsPbI3的形成能和e above hull 是用自己计算的能量和MP数据库中子成分的能量计算出来的，如果两个能量的计算参数不同，或者DFT的泛函不同，那能量肯定就不在一个尺度上。这样计算出来的形成能和e above hull 就并不准确。有两种思路来解决这个问题，第一，自己DFT计算的时候严格保证和MP上的计算参数相同。第二，建立一个能量标准来划分稳定相和亚稳定相。比如我们计算的DFT数据中某个相是公认稳定的，我们就以这个相计算出来的e above hull 的值作为稳定相判据的上限。比如CsSnI3 是稳定的，它的e above hull是0.045eV/atom（瞎举的例子，我也不知道CsSnI3是否真的稳定！）。那么我们可以认为CsPbI3 也是稳定的，它的e above hull 小于0.045，就可以把0.045定为稳定相的上限值，同理，找个亚稳定相可以定一个亚稳相上限。第二种方法我也是偶然在某篇文献中看到的，不确定是否正确，文献也找不到了。

上代码

我把上述计算形成能和e above hull 的过程写成了函数，输入为excel表，输出为对应形成能和e above hull的excel表。大家可以利用这个函数进行批量计算。使用时代码需要修改两个地方，第一个为输入的excel路径，也就是第二个代码块中的tot_energy_excel_path=‘test.xlsx’，引号中为excel路径，像这样用相对路径时，需要把excel文件放到代码的同级目录，也可以用绝对路径。第二地方是MPAPI_KEY=‘xxxxx’，引号中的为MP提供的API KEY(2023年3月12日更正： 使用MP新版本中的API出现了返回的子成分能量与网站中的对应不上的问题，且有大量重复，导致ehull计算结果出现明显错误，建议使用老网站中的API_KEY，测试发现返回值正常，后续会向MP反映这个bug====2023年3月13日更正： MP开发人员反馈是mp-api版本原因，使用pip install mp-api --upgrade更新为最新版后问题解决，PS.只更新pymatgen不会同时更新mp-api)，是下载数据用的，需要在MP官网中拷贝。如图3，点击官网首页右上方的API。然后选中图中红色的API KEY 替换为代码中引号内的内容就行。
运行第一个代码块：运行函数
第二个代码块：调用函数开始计算
第三个代码块：保存结果为aboveHull_result.xlsx’，在代码文件的同级目录
第四个代码块：画第index个样本的相图
注意：导入的excel表数据格式如图4，共两列，第一列为化学式，第二列为总能，化学式一定是一个元素跟一个数字，1可以省略，不要有括号那些，总能是这个化学式下的总能量，千万不要用每原子能量。输出的excel表如图5，输出的文件名可以在第三个代码块的引号中更改。

from pymatgen.analysis.phase_diagram import PDEntry, PhaseDiagram
from pymatgen.analysis.phase_diagram import PDPlotter
from pymatgen.ext.matproj import MPRester
import re
def get_above_hull_and_formation_enrgy_by_total_energy(tot_energy_excel_path,MPAPI_KEY):
    import pandas as pd
    # 导入化学式——总能量数据表
    tot_energy_excel = pd.read_excel(tot_energy_excel_path)
    formu_energy_dic = dict(zip(tot_energy_excel.iloc[:,0],tot_energy_excel.iloc[:,1]))
    #批量处理分支
    #获取MP数据
    pd_list = []
    forma_energy_list = []
    e_above_hull_list = []
    with MPRester(MPAPI_KEY) as mpr:
        for key in formu_energy_dic:
            div = re.compile(r"([A-Z][a-z])|([A-Z])|(\(.*\))")
            ele_list = div.findall(key)
            formu_key=''
            for el in ele_list:
                for e in el:
                    if e != '':
                        formu_key+=e+'-'
            formu_key=formu_key.strip('-')
            print(key)
            entries = mpr.get_entries_in_chemsys(formu_key)#获取MP数据库所有相关相
            pdentry = PDEntry(composition=key,energy=formu_energy_dic[key])#创建自己导入的化学式-总能相
            entries.append(pdentry)
            # 画相图
            pd = PhaseDiagram(entries=entries)
           
            #计算导入相形成能
            forma_energy_list.append(pd.get_form_energy_per_atom(pdentry))
            #计算导入相above hull
            e_above_hull_list.append(pd.get_e_above_hull(pdentry))

#             plotter = PDPlotter(pd)
#             plotter.get_plot()
            pd_list.append(pd)
    tot_energy_excel['formation_energy_per_atom']=forma_energy_list
    tot_energy_excel['e_above_hull'] = e_above_hull_list
    res_excel = tot_energy_excel
    display(res_excel)
    return pd_list,formu_energy_dic,res_excel

pd_list,formu_energy_dic,res_excel = get_above_hull_and_formation_enrgy_by_total_energy(tot_energy_excel_path='test.xlsx',MPAPI_KEY='xxxxx')

# 将结果保存为Excel
res_excel.to_excel('aboveHull_result.xlsx')

#画相图，只能手动更改index值指定画第几个化学式的相图
index = 0 #  第index个样本的相图
display(res_excel.iloc[index,:])
plotter = PDPlotter(pd_list[index])
plotter.get_plot()

图3 MP上获取API KEY,红色那串就是

图4 导入的excel表数据格式

图5 输出结果

ok 以上就是全部内容，有问题可以留言。