目录

弹性常数

QuantumATK 提供非常简单的方法计算 任意块体的弹性常数。该方法具有普遍适用性,极易与密度泛函理论(ATK-DFT)、半经验方法(ATK-SE)、经典势(ATK-ForceField)结合使用。在本教程中,您将学习如何借助 QuantumATK 计算弹性常数。

注意

块体的弹性常数与应力张量 $\boldsymbol{\sigma}$ 对体系外加应变 $\boldsymbol{\varepsilon}$ 的线性关系有关。因此,弹性常数描述了在特定变形下材料的定向刚度。其他的一般模量,如体模量、剪切模量、杨氏模量也都可以计算。

所有这些定量都表现的是固体的力学性质。除此之外,弹性常数和模量也是常被用来判断经典势的合适参数。

方法原理

应力和应变张量总是组成对称的 3×3 矩阵,因此它们可以通过 Voigt notation 表示法用 6 个矢量简洁地表达。

$$ \boldsymbol{\sigma} = ( \sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy} ) , $$

$$ \boldsymbol{\varepsilon} = ( \varepsilon_{xx},\varepsilon_{yy},\varepsilon_{zz},2\varepsilon_{yz},2\varepsilon_{xz},2\varepsilon_{xy} ). $$

应力张量和给定应变矢量的线性关系可以表达为:

$$ \boldsymbol{\sigma} = \boldsymbol{C} \cdot \boldsymbol{\varepsilon} , $$

这里的 $\boldsymbol{C}$ 是包含弹性常数的 6×6 对称矩阵。

随着晶体对称性的增高,矩阵 $\boldsymbol{C}$ 中独立项数目越少。例如,在立方晶体中只有 $C_{11}$、$C_{12}$ 和 $C_{44}$ 这 3 个独立项:

$$ \begin{split}\begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_4 \\ \sigma_5 \\ \sigma_6 \\ \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \\ \end{pmatrix} \cdot \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5 \\ \varepsilon_6 \\ \end{pmatrix}\end{split} $$

为了获得弹性常数,需要在模拟的晶胞上沿选定的应变矢量上施加一个小的形变 $\eta$,计算产生的应力张量。通过拟合每个 Vogit 应力和对应的应变矢量组成的 $\sigma_i(\eta)$ 曲线就可得到线性应力分布情况。然后将弹性常数作为计算线性方程组的最小二乘解,同时考虑晶体的对称性。

注意

为了计算弹性常数,ATK 采用拉格朗日应变和应力张量来代替相应的物理张量。

为使计算应力的数目尽量最少,ATK 选用通用线性无关耦合应变(ULICS)矢量(Ref.[YZY10])。对于每个应变矢量,通常有 3 个以参考构形($\eta=0$)为中心的形变($-\eta$, 0, $+\eta$)。

重要

对于晶胞中原子数多于 1 个的情况,每个被施加应变的晶胞的原子位置都要在计算应力前进行优化;ATK 会自动处理这个问题。但是在开始计算弹性常数之前,一定要确保参考构形的晶胞矢量(和原子位置)已完全优化。

利用经典势计算弹性常数

接下来您将基于经典势,采用 ATK-ForceField 模块完成弹性常数的计算。经典势的优点除了计算速度快以外,大多数的性质如能量、压力、应力等都是坐标的平滑函数。这使得为计算应变设置形变从而获得弹性常数的方法热别稳健。

以下为 2 个参考例子:

晶体硅

打开 VNL 的 Builder 并点击 从数据库中选择 “Silicon(alpha)” 添加到 Stash。然后发送结构到 Script Generator ,做如下操作:

注意

在 2017 年之前的 ATK 版本中,ATK-ForceField 计算器在 ATK-Classical 的下拉菜单中。

第一组参数和选项可以参考应力/应变的计算:

参数:$\eta_{max}$

设置施加在晶胞上的最大形变量。该参数要足够大以实现应力发生显著变化,也要小到使其保持在基本的线性范围内。多数情况下默认值 0.002 是一个不错的选择。用经典势计算弹性常数通常对这个参数不太敏感。

参数:$n_\eta$

设置每个应变矢量从 $-\eta_{max}$ 到 $+\eta_{max}$ 之间形变的数目。较高的该参数值配合较高的 “Fit order” 值(下一个要设置的参数),不仅可以有助于滤除可能的非线性部分影响,还会增加需要执行计算的数目。默认值 3 一般就可以保证运行良好。

参数:Fit order

设置应力函数对 $\eta$ 拟合的多项式最高阶数,其中 $\eta$ 应该要比 $n_\eta$ 小。为了评估弹性常数,仅使用线性贡献的部分。

选项:Enable symmetry

允许 ATK 检测晶体的对称性,只计算这种点阵对称的独立常数。否则,上三角矩阵 $\boldsymbol{C}$ 的所有 21 个常数都会被当做是独立的。

其余的设置与每个应变系统应力计算前调用的力的优化(对原子坐标)有关。

选项:Optimizer

设置优化方法。如果选择 “None”,应力计算前不优化力,这是速度最快但精确度最低的选项。若单胞中只有一个原子,就没有必要做内坐标的优化,优化程序即自动失效

最小化力和应力的参数

其他的参数是控制优化程序的设置。这些选项实际上跟 OptimizeGeometry 模块的一样。而且这些参数中的大多数只在一些特殊情况中有相关性,“Max forces” 这个参数的选择就要足够精确。默认值为 0.05 eV/Å(比通常 ATK 的默认值低 10倍)代表了准确性和效率之间达到了很好的平衡。由于案例中的使用了快速的 ATK-ForceField 计算,您可以轻松地将此值进一步降低至 0.001 V/Å以实现更高的精度。

通过把脚本提交到 Job Manager 执行计算,然后再点击 按钮。也许会出现被要求再次保存脚本的情况,点击 按钮启动计算。在本机上完成计算大概需要的 5 秒的时间。

结果分析

您可以在日志类文件中找到计算和结果的信息。如您所见,ATK 可以正确地识别 FCC 晶体硅的立方对称性,而且只需要 3 个独立的弹性常数。

  

+------------------------------------------------------------------------------+
| Calculating elastic constants for the given bulk configuration               |
+------------------------------------------------------------------------------+
| Detected space group number:    227                                          |
| Detected lattice symmetry:      Cubic                                        |
| This lattice symmetry has 3 independent elastic constants:                   |
| C11  C12  C44                                                                |
+------------------------------------------------------------------------------+

日志文件底部显示有弹性常数矩阵,且与使用了 Stillinger-Weber 势的文献中的结果高度吻合 [Cow88]。 $$ \begin{split}\begin{align} C_{11} &= 151.4\, \mathrm{GPa}, \\ C_{12} &= 76.4\, \mathrm{GPa}, \\ C_{44} &= 56.4\, \mathrm{GPa}. \end{align}\end{split} $$

提示

日志文件中还记录了跟弹性常数矩阵有关的更多详细分析。它包含了弹性常数矩阵的逆矩阵—弹性柔度矩阵。还有很多常见模量,如体模量、剪切模量和杨氏模量,都是由弹性常数计算得到的。需要注意的是,根据 Voigt 和 Reuss 计算公式的细微差别,关于体模量和剪切模量有 3 种不同的描述。为了方便比较,这些结果都罗列在输出日志文件中。

SiO2石英

您可以使用同样的方法去计算更复杂的晶体的弹性性质,SiO2 ɑ 相石英。按照上个例子中相同的步骤设置计算,选择参考文献 Ref.[APGMMCM+06] 中的 “Pedone2006_Fe2” 势。

精确的长程静电

在设置了 Script Generator 中的脚本后,Script Detail 处在 Global IO 场中将 “Minimal” 更改为 “Show defaults”,然后将脚本发送到 Editor 做一个小改动。在计算器模块中,您将会发现所有的对势都被单独罗列出并添加到势的总集(这是设置了 “Show defaults” 的结果)。

Pedone 势包含长程静电相互作用。您可以手动地增加这些交互的截断距离(r_cut)。为了实现应力的高精度计算,您可以把这个参数从默认值 9.0 Å 增大到 15 Å,以考虑库仑势的长程性质。

potentialSet.setCoulombSolver(CoulombDSF(r_cut=15.0*Angstrom, alpha = 0.2))
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setInternalOrdering("default")

这样做会影响计算速度,但弹性常数的计算普遍模拟时间都不是很长,所以您几乎不会察觉到影响。

结果

计算完成后可以得到如下的结果:

  
+------------------------------------------------------------------------------+
| Elastic Constants in GPa                                                     |
+------------------------------------------------------------------------------+
|    86.55       8.71      11.16     -18.36       0.00       0.00              |
|               86.55      11.16      18.36       0.00       0.00              |
|                         106.67       0.00       0.00       0.00              |
|                                     49.41       0.00       0.00              |
|                                                49.41     -18.36              |
|                                                           38.92              |
+------------------------------------------------------------------------------+

由于举例的结构是斜方六面体,因此有 6 个独立常数 $C_{11}, C_{12}, C_{13}, C_{14}, C_{33}$ 和 $C_{44}$。所有的值都跟原始出版物 Ref.[APGMMCM+06] 中的结果一致。

运用DFT计算弹性常数

最后,您将运用密度泛函理论计算硅的弹性常数。可以按照上面列出的计算框架进行,只需要替换计算器,即使用 ATK-DFT 和 GGA 交换关联

再次用 Builder 新建一个硅的结构并添加到 Script Generator。然后依次点击 New Calculator OptimizeGeometry ElasticConstants 进行设置。

为了获得应变构形间受到应力数值上的精确差异,我们通常不得不使用比默认值更严谨的计算器设置。

OptimizeGeometryElasticConstants 模块在经典势的选择设置上基本相同,除了 ElasticConstants 的 $\eta_{max}$ 应该设为 0.003(意味着施加更大的应变,导致应力发生更显著的变化)。

这次计算只用 10 分钟就可以完成,得到的弹性常数为:

  
+------------------------------------------------------------------------------+
| Elastic Constants in GPa                                                     |
+------------------------------------------------------------------------------+
|   161.02      65.22      65.22       0.00       0.00       0.00              |
|              161.02      65.22       0.00       0.00       0.00              |
|                         161.02       0.00       0.00       0.00              |
|                                     75.68       0.00       0.00              |
|                                                75.68       0.00              |
|                                                           75.68              |
+------------------------------------------------------------------------------+

计算结果与如下所列的文献中实验数据 Ref. [MBBT51] 相当一致:

\begin{split}\begin{align} C_{11} &= 167.4\, \mathrm{GPa} \\ C_{12} &= 65.2\, \mathrm{GPa} \\ C_{44} &= 79.6\, \mathrm{GPa} \end{align}\end{split}

参考