弹性常数
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 个参考例子：
- 硅晶体，采用 Stilinger-Weber-Weber 势 [SW85]
- <chem>SiO2</chem>（$\alpha$-石英），采用 “Pedone2006 Fe2” 势 [APGMMCM+06]
晶体硅
打开 VNL 的 Builder 
并点击 
从数据库中选择 “Silicon(alpha)” 添加到 Stash。然后发送结构到 Script Generator 
，做如下操作：
- 添加一个新的计算器
模块 - 选择 ATK-ForceField 和 “Stillinger_Weber_Si_1985” 势
注意
在 2017 年之前的 ATK 版本中，ATK-ForceField 计算器在 ATK-Classical 的下拉菜单中。
- 点击 Optimize
OptimizeGeometry 
选项，打开页面并设置：- 给最大力和最大应力设置足够的精确度，例如最大力为 0.001 eV/Å，最大应力为 0.001 eV/Å3；
- 不勾选 “Constrain cell” 以实现晶胞矢量的优化；
- 选择优化方法为 LBFGS；
- 返回到之前的页面，点击 Analysis ElasticConstants
模块，打开设置分析的参数。
第一组参数和选项可以参考应力/应变的计算：
参数：$\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 
进行设置。
为了获得应变构形间受到应力数值上的精确差异，我们通常不得不使用比默认值更严谨的计算器设置。
- 选择计算器 ATK-DFT，设置如下图所示的 Basic 参数
- 设置网格截断能量为 300 Hartree；
- 设置 K 点取样为 15×15×15。
- 在 Iteration control parameters 参数选项卡里
- 设置 Tolerance 参数为 1.0e-08；
- 设置 Damping factor（阻尼系数）为 0.5，历史步数为 10。
- 在 Basis set/exchange correlation 选项卡使用 PBEsol 泛函，它可以提供可靠的体模量及其相关量密度泛函理论假设。
- 交换关联函数类型选 GGA；
- 在 Predefined functionals 列表中选择 PBES。
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}
参考
- [APGMMCM+06] | (1, 2, 3) † Alfonso Pedone, † Gianluca Malavasi, † M. Cristina Menziani, ‡ Alastair N. Cormack, and † Ulderico Segre*. A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses. The Journal of Physical Chemistry B, 110(24):11780–11795, 2006. doi:10.1021/jp0611018.
- [Cow88] | E. Roger Cowley. Lattice dynamics of silicon with empirical many-body potentials. Phys. Rev. Lett., 60:2379–2381, Jun 1988. doi:10.1103/PhysRevLett.60.2379.
- [MBBT51] | H. J. McSkimin, W. L. Bond, E. Buehler, and G. K. Teal. Measurement of the elastic constants of silicon single crystals and their thermal coefficients. Phys. Rev., 83:1080–1080, Sep 1951. doi:10.1103/PhysRev.83.1080.
- [SW85] | F. H. Stillinger and T. A. Weber. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B, 31:5262–5271, Apr 1985. doi:10.1103/PhysRevB.31.5262.
- [YZY10] | R. Yu, J. Zhu, and H.Q. Ye. Calculations of single-crystal elastic constants made simple. Computer Physics Communications, 181(3):671 – 675, 2010. doi:10.1016/j.cpc.2009.11.017.