内聚能的定义是晶体的能量减去够成晶体的成分的能量,然后除以成分的“总个数”。因此对于分子晶体,是指分子晶体的能量减去各种分子的能量,然后除以各种分子的总个数,对于类似二氧化硅这样的晶体,Si和O就成为两个成分,和分子晶体的定义就完全类似了。公式上体现为:
$E_{coh}$ = $(E_{AmBn}$ - $mE_A$ - $nE_B)$/(m + n)
计算思路:用同样的参数,计算晶体、A、B,分别得到能量,然后带入上式。不过计算单体A、B,周期性可以设置为None,如此可以减小计算量的。
本文基于BAND模块,Quantum ESPRESSO的计算是完全类似的,参考:Quantum ESPRESSO内聚能的计算
本文使用软件版本:AMS2024.102。
例如萘晶体,单包包含2个萘分子,因此组分只有一种。晶体坐标(点击进入,复制内容,在AMSinput窗口Ctrl v粘贴即可完成建模):
看起来很碎,实际上就是正确的2个完整的萘分子,只不过分子处于Cell边界,因此根据周期边界条件,露出Cell边界的原子,出现在Cell另一侧。如果我们不习惯这样的显示方式,可以通过❄ → Crystal → Map Molecules Complete,改为这种方式显示:
这样就显示完整的2个萘分子了。实际上两种显示模式,可以看到右下角原子个数是一样的,实际上是同一个模型,因此不影响计算结果。
由于是分子间相互作用,因此色散矫正就非常重要了,因此需要使用-D类型的泛函,本例中我们使用PBE-D4(EEQ),首先我们计算晶体(即优化):
一般性参数设置(详细说明参考:【入门基础教程】单点计算与BAND的基本参数设置、计算效率):
由于是分子晶体,因此可以设置k点为gamma only并没有什么影响,普通晶体此处则需要改为Good,或点击后面的>按钮自行指定Number of points。 勾选晶格常数优化功能:
保存并运行作业(分子晶体优化晶格常数,收敛是非常慢的),另外晶体结构优化收敛标准中,梯度阈值默认为 0.001 Hartree/Å,如需修改,也在上图页面修改。
单个萘分子的参数设置,由于没有周期性,也就谈不上k点,因此参数就简单一些:
如果体系是原子,或者开壳层的分子,这里需要勾选Unrestricted。一般不需要设置具体Spin polarization,BAND一般会自行收敛到正确的Spin polarization值。
保存并运行作业。
计算过程中,应该跟踪SCM →logfile打开*.logfile文件,观察收敛情况。如果收敛则出现5个T:
<Jun16-2024> <21:04:27> current energy -4.35368845 Hartree <Jun16-2024> <21:04:27> energy change -0.00002594 0.00018000 T <Jun16-2024> <21:04:27> constrained gradient max 0.00062572 0.00100000 T <Jun16-2024> <21:04:27> constrained gradient rms 0.00027711 0.00066667 T <Jun16-2024> <21:04:27> gradient max 0.00062572 <Jun16-2024> <21:04:27> gradient rms 0.00027711 <Jun16-2024> <21:04:27> cart. step max 0.00289484 0.01000000 T <Jun16-2024> <21:04:27> cart. step rms 0.00121744 0.00666667 T
那么*.logfile尾部的能量就得到了(三种单位):
<Jun16-2024> <21:04:35> ENERGY OF FORMATION: -4.3537 A.U. <Jun16-2024> <21:04:35> -118.4699 E.V. <Jun16-2024> <21:04:35> -2731.9810 KCAL/MOL
这是萘分子的结果。类似可以得到晶体的能量:
<Jun17-2024> <13:45:00> ENERGY OF FORMATION: -8.8613 A.U. <Jun17-2024> <13:45:00> -241.1279 E.V. <Jun17-2024> <13:45:00> -5560.5427 KCAL/MOL
因此内聚能为: $E_{coh}$ = $(E_{AmBn}$ - $mE_A$ - $nE_B)$/(m + n) = [-5560.5427 - 2*(-2731.9810)]/2 = -48.29985 KCAL/MOL
这个值一般小于0,值的绝对值越大,表示稳定性越强,如果大于0则表示不可能结合为晶体。也有的地方定义内聚能的时候,把相减的顺序颠倒过来,因此正负的含义也相反。