计算材料学案例--从总能到反应热

      在确定材料的初始结构后，通过结构优化得到体系对应的总能，我们可以进一步计算分子内部的键能、体材料的形成能，给出反应的热化学方程式等。我们给出相关的计算案例，讨论影响体系总能的因素，结合实验数据判断计算结果的合理性。（Benchmarks很重要！如果我们氢气、氧气都和实验没对上，那么含氢、氧元素的结构计算结果几乎是不可信的。）

       假如我们要计算H-H键能，E_{H-H}=2E_{H}-E_{H_2}，这里E_{H}和E_{H_2}分别是孤立的氢原子和氢分子的基态总能。由于体系是孤立分子，k点网格用1x1x1即可(VASP可以在KPOINTS文件中设置），而晶格常数不需要优化（在INCAR文件中设置ISIF=2），在SAGAR平台中设置kpts=1,1,1，ISIF=2即可。在计算孤立的氢原子时，体系存在未成对电子，需要使用自旋极化计算（在INCAR文件中设置ISPIN=2），在SAGAR平台中设置kpts=1,1,1，ISIF=2，ISPIN=2即可。

      从结果看，不设置ISPIN=2，孤立的H原子能量是-0.056eV，设置后，能量下降为-1.113eV，因此，E_{H}= -1.113 eV。而氢分子总能E_{H_2}= -6.737 eV，可知H-H键能为4.51eV，对应435.3kJ/mol，和实验值436kJ/mol符合。我们可以通过截断能测试的功能来考察计算精度对结果的影响，一般截断能越高，精度越高，但也更耗时。如下表所示，随着截断能的增加，H-H键能将收敛为4.54eV（438kJ/mol），和实验相符。

       对于单质体材料，内聚能计算公式：E_{coh}=E_{atom}-E_{bulk}，这里E_{atom}和E_{bulk}分别是孤立的原子基态总能和体结构总能对原子数目的平均。以Si的惯用晶胞为例，输入的晶格常数为实验值5.43，优化后为5.464，8个原子体系总能为-43.388 eV，而孤立的Si原子能量为 -0.885 eV，因此内聚能为 4.54 eV， 和实验值4.63 eV基本符合。类似的，可以计算金刚石结构，8个原子体系总能为 -72.750 eV，而孤立的C原子能量为 -1.321 eV，因此内聚能为 7.77 eV，和实验值7.37 有5%的偏差。优化后的晶格常数是3.573，和实验值3.567不超过0.2%的偏差。我们进一步可以考虑ENCUT和k点网格对结果的影响。

       要计算反应热，需要分别计算反应前后体系的总能，然后做差即可。

      比如氢气和氯气反应生成氯化氢分子，分别优化这三个分子，得到总能 -6.737 eV, -3.573 eV, -6.087 eV。计算可得，反应热绝对值是 1.86 eV， 右边总能更低，是放热反应， 对应179.5 kJ/mol, 和实验值184.6符合。

      再比如氢气和氧气反应生成水分子，水分子总能是 -14.221 eV，而氧气分子基态是自旋极化的，用ISPIN=2计算得到能量是 -9.861 eV， 不考虑自旋极化总能只有 -8.838 eV。计算方程两边总能差得到反应热是 2.55 eV（放热反应），对应246.6 kJ/mol，和实验值241.8 符合。

       注意：反应物和生成物的状态（固、气）和电子结构的基态对体系总能有影响，从而对反应热的数值有影响。 比如高温下木炭和水气反应得到一氧化碳和氢气，这里水、一氧化碳、氢气分子能量分别为 -14.221 eV，-14.777 eV ，-6.737 eV。如果C用金刚石的平均能量 -72.750/8 eV来计算，反应热数值为1.80 eV，对应174.1 kJ/mol。根据文献结果（Phys. Rev. B 67, 035423, 2003），无定形碳能量比金刚石高 0.3 eV左右，因此方程左边总能相对升高，反应热数值将变为 1.50 eV，对应145.1 kJ/mol，结果更合理。