Skip to content

Latest commit

 

History

History
138 lines (87 loc) · 8.98 KB

abacus-elastic.md

File metadata and controls

138 lines (87 loc) · 8.98 KB

ABACUS+pymatgen 计算弹性常数

作者:陈涛,邮箱:[email protected]

审核:陈默涵,邮箱:[email protected]

最后更新时间:2024/07/11

一、介绍

本教程旨在描述如何采用 ABACUS 来计算材料体系的弹性常数(elastic constants)。弹性常数是表征材料弹性的量,也是材料的重要力学性质,可以通过原子尺度模拟计算获得,例如通过密度泛函理论或者原子间的势函数方法。在晶体结构预测辅助的新材料计算设计中,弹性常数常用于检查结构的稳定性。

在晶体的线弹性范围内,应力$$\boldsymbol{\sigma}$$和相应应变$$\boldsymbol{\varepsilon}$$之间的关系符合胡克定律,$$\boldsymbol{\sigma=C\varepsilon},\sigma_{ij}=C_{ijkl}\varepsilon_{kl}$$,其中$$i,j,k,l$$是笛卡尔指标($$x,y,z$$),$$C_{ijkl}$$即为弹性常数。

由于$$\boldsymbol{\sigma}$$和$$\boldsymbol{\varepsilon}$$均是对称张量,基于 Voigt 表示法:$$xx\mapsto1,yy\mapsto2,zz\mapsto3,yz\mapsto4,xz\mapsto5,xy\mapsto6$$,则:$$\sigma_1=\sigma_{xx},\sigma_2=\sigma_{yy},\sigma_3=\sigma_{zz},\sigma_4=\sigma_{yz},\sigma_5=\sigma_{xz},\sigma_6=\sigma_{xy}$$

$$ \epsilon_1=\epsilon_{xx},\epsilon_2=\epsilon_{yy},\epsilon_3=\epsilon_{zz},\epsilon_4=\epsilon_{yz},\epsilon_5=\epsilon_{xz},\epsilon_6=\epsilon_{xy} $$

此时,$$C_{xxxx}\mapsto C_{11},C_{xxyy}\mapsto C_{12},.....$$

$$ \begin{bmatrix}\sigma_{1}\\sigma_{2}\\sigma_{3}\\sigma_{4}\\sigma_{5}\\sigma_{6}\end{bmatrix}=\begin{bmatrix}C_{11}&C_{12}&C_{13}&C_{14}&C_{15}&C_{16}\C_{12}&C_{22}&C_{23}&C_{24}&C_{25}&C_{26}\C_{13}&C_{23}&C_{33}&C_{34}&C_{35}&C_{36}\C_{14}&C_{24}&C_{34}&C_{44}&C_{45}&C_{46}\C_{15}&C_{25}&C_{35}&C_{45}&C_{55}&C_{56}\C_{16}&C_{26}&C_{36}&C_{46}&C_{56}&C_{66}\end{bmatrix}\begin{bmatrix}\epsilon_{1}\\epsilon_{2}\\epsilon_{3}\2\epsilon_{4}\2\epsilon_{5}\2\epsilon_{6}\end{bmatrix} $$

如果我们通过施加应变$$\varepsilon_i$$使晶体变形并计算相应的应力,则可以从方程中获得弹性常数。晶体晶胞上的变形矩阵为$$\boldsymbol{D=I+\varepsilon}$$,其中$$\boldsymbol{I}$$是 3*3 的单位矩阵,$$\boldsymbol{\varepsilon}$$是 Voigt 表示法中的应变矩阵。

在 3 维晶体中,应变矩阵为 $$\boldsymbol{\varepsilon}=\left[\begin{array}{lrr} \varepsilon_{1} & \varepsilon_{6} & \varepsilon_{5} \ \varepsilon_{6} & \varepsilon_{2} & \varepsilon_{4} \ \varepsilon_{5} & \varepsilon_{4} & \varepsilon_{3} \end{array}\right]$$,在2维平面晶体中,应变矩阵为$$\boldsymbol{\varepsilon}=\left[\begin{array}{lrr} \varepsilon_{1} & \varepsilon_{6} & 0\ \varepsilon_{6} & \varepsilon_{2} & 0 \ 0 & 0 & 0 \end{array}\right]$$。

变形后,晶格矢量为$$\boldsymbol{A'=A \cdot D}$$,其中$$\boldsymbol{A}$$是初始的晶格矢量。

在实际计算中,6 种应变状态将逐一应用于初始弛豫(relaxation)后的结构,因此每次只考虑一种独立变形。对于 6 种应变状态中的每一种,应用 4 种不同的默认应变大小:$$\varepsilon_i\in{-0.01,-0.005,0.005,0.01}$$,将产生 24 种构型。接着对这 24 种构型分别使用 DFT 计算(固定晶格矢量,但允许原子弛豫,即 relax),对所获得的每种应变状态的 4 个应力和应变的集合使用线性拟合来计算相应的弹性常数。

更多细节可以参考:Elastic Constants | Materials Project Documentation

二、准备

1. 下载例子

首先可以下载一个 ABACUS 的计算实例,可以从 Gitee 上下载。具体来说,可以通过在网页右侧点击克隆/下载-> 下载 ZIP 得到算例,或者在 linux 终端执行如下命令得到算例:

git clone https://gitee.com/mcresearch/abacus-user-guide.git

下载后解压,之后进入 abacus-user-guide/examples/elastic 文件夹

里面有 relax 文件夹,为 8 原子金刚石的算例,周期性边界条件。

此外,gene_dfm.pycompute_dfm.py 为计算弹性常数使用的 python 脚本。

run_task.shsub.sh 为批量运行 abacus 计算的 shell 脚本,INPUTKPTC_ONCV_PBE-1.0.upfC_gga_7au_100Ry_2s2p1d.orb 为运行 abacus 计算所需的输入文件。

2. 下载并安装 pymatgen

Python Materials Genomics (pymatgen,https://pymatgen.org/)是一个 API 包,该软件可以与 materials project 结合进行高通量计算。该软件包是由加州大学圣地亚哥雅各布斯工程学院的纳米工程教授 Shyue Ping Ong 和他的材料虚拟实验室(Materials Virtual Lab)团队开发并维护的程序。

在本教程里会用到 pymatgen 来计算弹性常数,此外还会用到 dpdata,monty,numpy 这三个库,可以使用如下命令安装:

pip install monty numpy dpdata pymatgen

弹性常数所用方法相关文档:https://pymatgen.org/pymatgen.analysis.elasticity.html

三、金刚石的弹性常数计算

1. 结构弛豫

进入 relax 文件夹,运行 ABACUS 进行结构弛豫,完成后在 OUT.C8 文件夹下出现 STRU_ION_D 文件。之后的应变都将在这个构型文件基础上产生。此外,计算弹性常数所用的应力也会减去这个构型的应力,因此也需要 running_relax.log 输出应力

2. 产生应变构型

回到例子根目录,执行如下命令:

python gene_dfm.py abacus

默认应变的大小(代表晶格常数的倍数,例如 0.01,其形变量是 1%)如下,对应 gene_dfm.py 的 41-42 行:

norm_strains = [-0.010, -0.005, 0.005, 0.010]
shear_strains = [-0.010, -0.005, 0.005, 0.010]

将会产生 task.000task.023 共 24 个文件夹,分别对应第一节所说的 24 种构型。进入任意 task 文件夹,其中会有 INPUTKPTSTRUstrain.json 四个文件。

其中 INPUTKPT 拷贝自例子根目录下的 INPUTKPT,由于使用的赝势和轨道文件放在例子根目录下,所以 pseudo_dirorbital_dir 均设置为 ../。此外需要进行固定晶格矢量,但允许原子弛豫来计算应力,因此 calculation 设置为 relaxSTRUstrain.json 分别是生成的构型文件和相应的应变大小。

3. 计算应力

分别进入上述 task 文件夹,运行 abacus 计算相应构型的应力。也可以使用例子根目录下的 run_task.sh 脚本,但要注意根据实际修改其和 sub.sh 的内容。

4. 计算弹性常数

回到例子根目录,执行如下命令:

python compute_dfm.py abacus

屏幕输出如下:

# Elastic Constants in GPa
1043.31  107.39  107.39    0.00    0.00    0.00 
 107.39 1043.31  107.39    0.00    0.00    0.00 
 107.39  107.39 1043.31    0.00    0.00    0.00 
   0.00   -0.00   -0.00  557.05    0.00    0.00 
  -0.00   -0.00   -0.00    0.00  557.05    0.00 
   0.00    0.00    0.00    0.00    0.00  557.05 
# Bulk   Modulus BV = 419.37 GPa
# Shear  Modulus GV = 521.42 GPa
# Youngs Modulus EV = 1105.91 GPa
# Poission Ratio uV = 0.06

第 2 到 7 行即为计算的各种弹性常数(单位:GPa),分别是体弹性模量(Bulk Modulus)、剪切模量(Shear Modulus)、杨氏模量(Youngs Modulus)和泊松比(Poission Ratio)。精度更高的计算结果存储在 elastic.json 内。

以下是常用的几个弹性常数相关的名词解释:

弹性模量(Bulk Modulus):指当有力施加于物体或物质时,其弹性变形(非永久变形)趋势的数学描述。物体的弹性模量定义为弹性变形区的应力-应变曲线的斜率。

杨氏模量(Young's Modulus):是材料力学中的名词,一般将杨氏模量习惯称为弹性模量。弹性材料承受正向应力时会产生正向应变,在形变量没有超过对应材料的一定弹性限度时,定义正向应力与正向应变的比值为杨氏模量。杨氏模量的大小标志了材料的刚性,杨氏模量越大,越不容易发生形变。大部分金属在合金成分不同、热处理在加工过程中的应用,其杨氏模量值会有 5%或更大的波动。

剪切模量(Shear Modulus):又称刚度模量,描述材料在剪切力作用下抵抗形变的能力,即物体形变的大小与作用的剪切应力之间的关系。

泊松比(Poisson's Ratio):当材料在一个方向被压缩,它会在与该方向垂直的另外两个方向伸长,这就是泊松现象,泊松比是用来反映泊松现象的无量纲的物理量。泊松比一般是正值,表示在一方向拉伸后,在其他方向收缩。不过也存在泊松比为零(在一方向拉伸后,在其他方向的尺寸不变),其至为负的材料。