最近计算能量体积曲线的时候遇到一些问题，花了不少时间去解决，走了很多弯路。下面以Graphite的能量体积曲线的计算为例，来说明如何使用VASP计算能量体积曲线，并且使用orgin进行非线性曲线拟合BM状态方程，拟合得到体系平衡态时的各个参数。特别感谢轩师兄，他在我困惑的时候用他敏锐的思维点醒了我，向他学习。如果表达有不合适的地方还希望大家批评指正。

1.  能量体积曲线曲线的计算：简单的说就是通过连续改变体系的体积，来获取不同体积下对应的能量值；

以空间群为（P63/mmc）的石墨为例，在使用VASP计算能量体积散点值的时候，我们需要注意一下几个方面：

1）能量体积散点值计算的时候，我们需要连续改变晶格参数以获得不同的模型体积。通常的做法我们是直接在POSCAR中使用分数坐标，然后连续变化POSCAR中第2行的缩放因子来改变模型的体积。得到不同体积的模型后，我们需要做一次体积保持不变的离子弛豫的计算（relax），这个时候ISIF参数的选择比较关键，你可以选择ISIF=2或者ISIF=4。二者的区别在于：ISIF=4会改变原胞的形状，而使用2则不会改变原胞的形状，这里说的改变原胞的形状就是调整晶格中c/a或者b/a的比值。对于石墨这种靠层间范德华相互作用结合在一起的晶体，就不适合使用4去优化结构（金刚石使用2和4基本无差别），因为在使用4的时候会导致石墨的层间距发生大的变化，石墨又是由弱的层间范德华作用会结合在一起的，层间距的变化导致石墨结构转化为两层单独的石墨烯，而无法拟合得到能量体积曲线的抛物线形状，你使用4的时候会得到如下如所示的形状：

 2）直接从Material Studio里面导出来的Graphite模型的晶格参数对应的是实验值，并不适合我们直接导出POSCAR进行能量体积散点值的计算，如果你直接用Material Studio导出的模型进行计算，你就会发现根据能量体积曲线拟合得到的平衡体积总是与文献中的不一致。这个时候需要做的就是：先得到石墨的对应能量最低的晶格参数（你可以通过连续改变a和c的值，来拟合得到最低能量时候的晶格参数；也可以直接查文献获取平衡晶格参数 ），得到晶格参数后我们再按照上面的方法通过改变缩放因子来改变晶体的体积，最后计算出不同体积对应的能量，使用状态方程拟合得到平衡态时的体积、能量、体积模量等，这个时候得到结果就与文献非常一致，如下：

2.   使用origin进行B-M状态方程拟合体积模量、平衡态时的能量体积（非线性曲线拟合），以origin2018为例，其他版本方法一致：

B-M状态方程表达式：

1）确定拟合函数，步骤如下：

（1）origin工具栏下选择拟合函数生成器

（2）选择:创建新的函数,点击下一步：

（3）选择或新建类别，函数名称，文件名称，描述这四项都可以自己定义，函数模型中选择显函数，函数类型选择公式，选择完毕后点击下一步：

（4）自变量中填x，因变量中填y，x和y分别相当于BM状态方程中的V和E(也就是我们计算得到的能量和体积的散点数据)，参数一栏中我们填入B0,B1,E0,V0,分别代表平衡时体系的体积模量，体积模量的一阶导数，平衡时的能量和平衡时的体积，点击下一步：

（5）函数主体中把BM方程写入，点击完成：

y=E0+(9.0/16)*V0*B0*(((V0/x)^(2.0/3)-1)^(3.0)*B1)+(9.0/16)*V0*B0*((V0/x)^(2.0/3)-1)^(2.0)*(6-4*(V0/x)^(2.0/3))

2）上面5步我们创建了方程，下面我们使用我们刚刚编辑的BM方程进行能量体积曲线的拟合：

（1）选中计算得到的能量体积散点数据绘制散点图：

（2）绘制完散点图后选择分析下的拟合>非线性曲线拟合>打开对话框

（3）函数选择取中选择我们刚才创建的函数

（4）特别重要的是：由于BM方程是个多参数方程，所以我们在利用能量体积散点数据数据进行拟合的时候，需要预先给定要拟合值的初始参数（初始参数的给定可以根据散点图估计给出），下图中，我们给定了E0和V0的初始参数，给定后我们选择拟合就可以了：

（5）拟合结果如下(方框中的值即是我们拟合得到的结果)：

Origin能量体积曲线拟合.pdf