计算弹性常数

背景

适用于大型周期性体系研究的密度泛函理论（DFT）方法取得的最新进展已使其成为解决材料设计和加工问题的重要手段。DFT方法可用于指导和引领新材料的设计，使研究人员能够了解反应和变化过程的基本化学和物理原理。

介绍

在本教程中，将学习使用CASTEP模块计算弹性常数和其他力学性能的方法。首先要优化立方BN晶体的结构，然后计算其弹性常数。本教程包括如下部分：

开始

优化立方BN晶体的结构

计算BN的弹性常数

分析弹性常数计算结果文件

注意：为了和本教程中的参数保持一致，可以使用Settings Organizer对话框将项目中所有参数都设置为BIOVIA的默认值。

1、开始

首先启动Materials Studio并创建一个新项目。

打开New Project对话框，输入BN_elastic作为项目名，单击OK按钮。

新项目将以BN_elastic为项目名显示于Project Explorer中。下一步是导入BN晶体结构。

在菜单栏中选择File | Import...，打开Import Document对话框。导航至Structures/semiconductors/文件夹，选择BN.xsd文件。单击Open按钮。

显示BN的晶体结构。

立方BN的结构

这是BN结构的惯用胞表示。为缩短计算时间，应该将其转换为初级胞表示。

在菜单栏中选择Build | Symmetry | Primitive Cell。

2、优化立方BN晶体的结构

在计算弹性常数之前并不一定要进行几何优化，可以利用实验观测到的结构计算出Cij数据。尽管如此，如果对晶胞进行完整的几何优化过程（包括晶胞参数的优化），可以获得相容性更好的结果，进而计算与理论基态对应的结构的弹性常数。

弹性常数，尤其是剪切常数的精确度，主要取决于SCF计算的精度，特别是布里渊区取样和波函数收敛程度的精度。所以应设置SCF、k点取样和FFT格子的精度为Fine。

现在设置几何优化参数。

在Modules工具条中单击CASTEP按钮

后面的下拉箭头，然后从下拉列表中选择Calculation，或从菜单栏中选择Modules / CASTEP / Calculation。

打开CASTEP Calculation对话框。

在Setup选项卡中，把计算任务Task设置为几何优化Geometry Optimization，将精度Quality设置为Fine，并且将泛函Functional设置为GGA和PBESOL。

单击More...按钮，打开CASTEP Geometry Optimization对话框。从晶胞优化Optimize cell的下拉列表中选择Full，关闭对话框。

在CASTEP Calculation对话框中，选择Job Control选项卡，选择想要运行CASTEP计算任务的网关Gateway。

单击Run按钮。

几何优化完成后，结构的晶胞参数应为a=b=c=2.553 Å，即对应于惯用胞中晶格常数为3.610 Å（实验值为3.615 Å）。

在3D Viewer中单击鼠标右键，然后从快捷菜单中选择Lattice Parameters。

将显示晶胞参数。现在可以继续计算优化后结构的弹性常数。

3、计算BN的弹性常数

在CASTEP Calculation对话框中选择Setup选项卡。从计算任务Task下拉列表中选择弹性常数Elastic Constants，单击More...按钮。

打开CASTEP Elastic Constants对话框。

将“每次应变加载步数Number of steps for each strain”由4更改为6，关闭对话框。确保BN CASTEP GeomOpt/BN.xsd为当前文档，单击CASTEP Calculation对话框中的Run按钮。

注意：如果在打开CASTEP Elastic Constants对话框之前，已经将BN CASTEP GeomOpt/BN.xsd设置为当前文档，则该对话框上的应变模式Strain pattern列表将包含数值。

CASTEP的弹性常数计算任务的结果以一组.castep输出文件的形式给出。其中每一个文件代表对于给定的应变模式和应变幅度，使用固定单元进行几何优化的结果。

4、分析弹性常数计算结果文件

注意：算得的结果可能和所展示的结果有细微不同，是由于输入构型的微小差异导致的。

对于这种点阵类型，仅需要一种应变模式。从相应的.castep文件中提取的计算得到的应力汇总如下：

所有应变模式的结果总结如下：

仅当使用两个以上的应变振幅值时，才可以计算误差，因为如果只以两个点拟合直线，则不存在统计不确定性。

然后，弹性常数以常规的6×6张量形式表示，随后是柔度系数的类似6×6张量形式：

上述算得的数值大致在DFT计算与实验测量值（B=396 GPa，C11=820 GPa，C12=190 GPa，C44=457 GPa）典型误差（10%）的范围内。

晶体弹性性质的计算相比于晶格参数和原子坐标，对电子结构计算的准确性明显更加敏感。因此有必要检查计算得到的性质相对于以下参数的收敛性：

k点密度（最重要）

截断能

增强密度比例因子，当使用超软赝势的情况下，动态生成或列于表中

弹性性质计算中需要考虑的另一个因素是交换关联泛函的选择。PBESOL和Wu-Cohen等现代泛函，比传统LDA或PBE泛函能更精确地再现固体性质；建议使用这些泛函计算固体的弹性常数。