目录

弹性常数

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 进行设置。

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

OptimizeGeometry 和 ElasticConstants 模块在经典势的选择设置上基本相同，除了 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}

参考