4.1 力学性质

发表于 更新于 分类于

弹性常数

基本概念

弹性常数描述了晶体对外加应变的响应的刚度。在材料的线性变形范围内（应变较小的情况下），体系的应力与应变满足胡克定律。也就是说，对于足够的小的变形，应力与应变成正比，即应力分量(S)是应变分量(E)的线性函数，三维材料的弹性刚度常数矩阵是6×6的。公式中Cij就是我们通常所说的弹性常数。因为刚度矩阵是对称矩阵，因此，弹性常数的独立张量元数目至多只有21个。对不同的晶系的晶体，因为对称性的关系，其独立的弹性常数是确定的。因此，晶系的对称性越高，独立的张量元数目越少。
立方晶系 ——只有3个独立矩阵元（C11，C12，C44）
六角晶系 ——有5个独立矩阵元（C11，C12，C13，C33，C44）
三角晶系
a) 32，3m，-32/m——有6个独立矩阵元（C11，C12，C13，C14，C33，C44）
b) 3，-3,——有8个独立矩阵元（C11，C12，C13，C14，C15，C33，C44，C45）
四方晶系
a) 422，4mm，-42m，4/mmm——有6个独立矩阵元（C11，C12，C13，C33，C44，C66）
b) 4，-4，4/m——有7个独立矩阵元（C11，C12，C13，C16，C33，C44，C66）
正交晶系 ——有9个独立矩阵元（C11，C12，C13，C22，C23，C33，C44，C55，C66）
单斜晶系 ——有13个独立矩阵元
三斜晶系 ——有21个独立矩阵元

计算过程

采用第一性原理计算弹性常数有两种方法，第一种方法是应力-应变方法， 即通过给结构施加不同应变，分别计算出所对应的应力大小，然后利用公式（1）拟合得到一次项系数，从而得到。第二种方法是能量-应变方法， 即通过给结构施加不同应变后计算出体系总能相对于基态能量变化（应变能），再利用公式（2）进行二次多项式拟合，其中二次多项式系数是晶体的某个弹性常数或者弹性常数组合。第一种方法优点是每进行一次计算可一次得到六个独立分量，缺点是为了得到准确的应力大小，必须选择更高的截断能和更密集的K格点。在同样的精度下，能量-应变法相比应力-应变方法要求的截断能和K点数目相对较少，缺点是计算应变数目有所增加。VASPKIT中弹性常数计算基于能量-应变法，目前不支持三斜晶系。
一般计算弹性常数我们可以采用应变应力或者应变能量关系，进行拟合。或者可以直接通过VASP直接计算。

方法一：

首先介绍直接用VASP在输入文件中添加参数进行计算的方法。

  1. 当然首先要进行结构优化，

  2. 在INCAR中添加关键参数，进行弹性常数的计算

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    NSW=1
    EDIFFG=-1e-3
    ISIF=3
    IBRION=6 #计算弹性常数
    POTIM=0.015
    NFREE = 4

  3. 计算完成后可以查看OUTCAR得到弹性刚度矩阵。或者可以利用VASPKIT直接得到弹性常数。
    请注意OUTCAR中的刚度矩阵是未经过处理的刚度矩阵，二维材料需要乘以C轴长度的0.01，还需要了解Voigt简标，11>1 22>2 33>3 23>4 13>5 12>6。然后我们通过Origin等工具绘图即可得到材料的杨氏模量，泊松比等力学性质。

方法二：

VASPKIT和VASP计算晶体的弹性常数具体计算步骤分为：

  1. 准备优化彻底的POSCAR文件，注意通常采用标准的惯用原胞计算弹性常数，如果不确信POSCAR文件中是否是标准的惯用原胞，可以用vaspkit-603/604生成标准结构；

  2. 运行vaspkit 选择102生成KPOINTS，由于计算弹性常数对K-mesh要求很高，因此对于半导体（金属体）体系，生成K点的精度应不小于 0.03(0.02)×2π {\AA}^{-1}

  3. INCAR参考设置如下；

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    Global Parameters
     ISTART =  0        
     LREAL  =  F      
     PREC   =  High  （截断能设置默认值1.5-2倍）
     LWAVE  = F       
     LCHARG = F     
     ADDGRID= .TRUE.    
    Electronic Relaxation
     ISMEAR =  0          
     SIGMA  =  0.05       
     NELM   =  40           
     NELMIN =  4           
     EDIFF  =  1E-08       
    Ionic Relaxation
     NELMIN =  6          
     NSW    =  100         
     IBRION =  2           
     ISIF   =  2    （切记选择2,如果选择3会把施加应变后原胞重新优化成平衡原胞）
     EDIFFG = -1E-02

  4. 准备VPKIT.in文件并设置第一行为1（预处理）；运行vaspkit并选择201, 将生成用于计算弹性常数文件；

    {.line-numbers}
    1
    2
    3
    4
    1                    ! 设置1将产生计算弹性常数的输入文件，2则计算弹性常数
    3D                   ! 2D为二维体系，3D为三维体系
    7                    ! 7个应变
    -0.015 -0.010 -0.005 0.000 0.005 0.010 0.015  ! 应变变化范围

  5. 批量提交vasp作业；
    总计有9个大文件夹，分别以该空间群所需要的独立弹性常数来命名；在每个大文件夹内则是以应变强度命名的实际计算执行的文件夹，文件夹的数量可根据预处理的VPKIT.in文件的最后一行的应变设置来更改。在每一个应变文件夹内则是施加了应变的POSCAR和事先准备好的INCAR、KPOINTS和POTCAR文件的链接，可直接进行vasp计算。
    这里通过一个简单的两重for循环命令依次序完成所需要的计算任务。
    脚本 :

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    #!/bin/bash
    path=`pwd`
    #记录当前文件夹位置并定义到变量 path中，称为初始目录
    pwd
    for cij in `ls -F | grep /$`
    #ls -F| grep /$ 命令为返回当前文件夹中所有子目录
    #将所有子目录添加到变量cij中，并准备进行for循环
    do
    cd ${path}/$cij
    #进入cij变量所代表的子目录中
    pwd
    for s in strain_*
    #将cij代表的子目录中所有strain文件夹添加到变量s中
    #此部添加变量s的方式可更改为添加变量cij的方式，但不建议这么做
    do
    cd ${path}/$cij/$s
    #进入初始目录下变量cij和变量s组成的目录，由于for循环的设置，将遍历所有可计算的strain文件夹
    pwd
     mpirun -np 16 vasp_std
     #执行vasp计算
    cd ${path}
    #返回初始目录
    pwd
    done
    done

  6. 再次修改VPKIT.in文件中第一行为2（后处理），然后再次运行vaspkit并选择201，得到以下结果；

  7. 如果在计算弹性常数时希望强制固定某个体系的空间群，可在SYMMETRY.in文件第二行设置空间群，这样晶系将通过指定空间群判断并计算相应的独立弹性常数，这对于掺杂或者合金体系比较有用。

根据弹性常数计算和三维可视化材料力学量

首先利用VASPKIT-204命令得到体弹性模量、杨氏模量、剪切模量及泊松比随角度的依赖关系，具体算法可参考文献【J. Phys. Condens. Matter, 28, 275201 (2016)和Comput. Phys. Commun. 181, 2102–2115 (2010)】。接下来我们以某单斜体系的弹性常数为例来演示如何进行材料力学量三维可视化。

  • 第一步新建ELASTIC_TENSOR.in并输入6x6弹性常数矩阵，注意第一行是注释行不可省略。如果是二维体系，则输入文件为ELASTIC_TENSOR_2D.in，注意二维体系弹性常数矩阵大小为3x3。
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    # comment line (in GPa)
       228.38   85.741   81.503        0   -0.737        0
       85.741   217.47   94.201        0  -20.213        0
       81.503   94.201   178.81        0   -9.472        0
            0        0        0   35.094        0  -17.851
       -0.737  -20.213   -9.472        0   37.778        0
            0        0        0  -17.851        0   42.708
  • 第二步运行VASPKIT 204命令得到MECHANICS_3D.dat
    打开MECHANICS_3D.dat可看到第一行给出了每一列数据所代表的物理量，第二行给出了球坐标系下分别划分仰角θ∈[0,180]和方位角ϕ∈[0,360]的格点大小。
  • 第三步从vaspkit/examples/angular_dependent_mechanics (VASPKIT ver. >= 1.3.2)文件夹中拷贝mechanics_3d_plot_matlab.m到当前目录，调用Matlab软件运行该脚本，得到以下信息.
    如果我们想可视化杨氏模量，则输入2回车即可得到。
  • 另外，以下两个软件也可以实现材料力学量三维可视化。
    Elate
    ElasticPOST