基于应力-应变关系计算弹性常数—VASPKIT v1.2.0新功能

‍

先前王伟老师已经介绍了基于能量-应变关系计算弹性常数—VASPKIT v1.00新功能的教程，今天介绍VASPKIT 1.2版本最新功能—基于应力-应变关系计算弹性常数。

在材料的线性形变范围内(小应变的情况下)，整个体系的应力与应变满足胡克定律，

上面表达式中σi和εj分别表示应力和应变，Cij是弹性常数，其中1≤i≥6和1≤j≥6。通过给结构施加不同应变，分别计算出所对应的应力大小，然后利用上面公式拟合得到一次项系数，从而便可得到Cij。下图给出了立方晶系的具体形式：

下面我们以金刚石结构为例讲解如何采用应力-应变函数关系计算弹性常数，详见VASPKIT/examples/elastic/diamond_3D。由于金刚石具有立方晶体结构，一共有3个独立弹性常数C11、C12和C44 (不明白的请看原胞转化方法以及标准原胞在计算中的重要性)。

1. 准备优化好的POSCAR文件，注意通常采用具有标准基矢形式的原胞计算弹性常数(VASPKIT-603/604可以生成标准结构)，至于原因请看原胞转化方法以及标准原胞在计算中的重要性。

2. 运行VASPKIT-102生成KPOINTS (注意精度要稍高一些)

3. 运行VASPKIT-101-DC生成INCAR文件，并根据实际情况修改，以下仅供参考：

Global ParametersISTART = 0 (Read existing wavefunction; if there)LREAL = F (Projection operators: automatic)PREC = High (Precision level)LWAVE = F (Write WAVECAR or not)LCHARG = F (Write CHGCAR or not)ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic RelaxationISMEAR = 0 (Gaussian smearing; metals:1)SIGMA = 0.05 (Smearing value in eV; metals:0.2)NELM = 40 (Max electronic SCF steps)NELMIN = 4 (Min electronic SCF steps)EDIFF = 1E-08 (SCF energy convergence; in eV) # GGA = PS (PBEsol exchange-correlation)

Ionic RelaxationNELMIN = 6 (Min electronic SCF steps)NSW = 100 (Max electronic SCF steps)IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)EDIFFG = -1E-02      (Ionic convergence; eV/AA)

4. 准备VPKIT.in文件并设置第一行为1 (预处理)，运行VASPKIT并选择200, 将生成用于计算弹性常数的文件

1 ! 1 for prep-rocessing, 2 for post-processing3D ! 2D for slab, 3D for bulk7 ! number of strain-0.015 -0.010 -0.005 0.000 0.005 0.010 0.015                   ! magnitude of strain

运行VASPKIT会在屏幕输出以下信息：

-->> (01) Reading VPKIT.in File... +-------------------------- Warm Tips --------------------------+ See some examples in vaspkit/examples/elastic, Require the fully-relaxed and standard Conventional cell. The stress-strain method requires higher ENCUT and denser K-Mesh +---------------------------------------------------------------+ -->> (02) Reading Structural Parameters from POSCAR File... -> C11_C12_C44 folder created successfully! -> strain_-0.015 folder created successfully! -> strain_-0.010 folder created successfully! -> strain_-0.005 folder created successfully! -> strain_0.000 folder created successfully! -> strain_+0.005 folder created successfully! -> strain_+0.010 folder created successfully!-> strain_+0.015 folder created successfully!

5. 批量提交vasp作业，相应的脚本如下(根据实际情况修改)

#!/bin/bashroot_path=`pwd`for cij in `ls -F | grep /$`docd ${root_path}/$cijfor s in strain_*docd ${root_path}/$cij/$secho `pwd` cp ../../vasp.job . ./vasp.job# Add here your vasp_submit_job_scriptdonedone

6. 等VASP全部计算完成之后，再次修改VPKIT.in文件中第一行为2 (后处理)，然后再次运行VASPKIT并选择200，得到以下结果；

-->> (01) Reading VPKIT.in File...+-------------------------- Warm Tips --------------------------+See some examples in vaspkit/examples/elastic,Require the fully-relaxed and standard Conventional cell.The stress-strain method requires higher ENCUT and denser K-Mesh+---------------------------------------------------------------+-->> (02) Reading Structural Parameters from POSCAR File...-->> (03) Calculating fitting coefficients of stress vs strain.+-------------------------- Summary ----------------------------+Based on the Strain versus Energy method.Crystal Class: m-3mSpace Group: Fd-3mCrystal System: Cubic systemIncluding Point group classes: 23, 2/m-3, 432, -43m, 4/m-32/mThere are 3 independent elastic constantsC11 C12 C12 0 0 0C12 C11 C12 0 0 0C12 C12 C11 0 0 00 0 0 C44 0 00 0 0 0 C44 00 0 0 0 0 C44

Stiffness Tensor C_ij (in GPa):1050.316 126.488 126.488 0.000 0.000 0.000126.488 1050.316 126.488 0.000 0.000 0.000126.488 126.488 1050.316 0.000 0.000 0.0000.000 0.000 0.000 559.816 0.000 0.0000.000 0.000 0.000 0.000 559.816 0.0000.000 0.000 0.000 0.000 0.000 559.816

Compliance Tensor S_ij (in GPa^{-1}):0.000977 -0.000105 -0.000105 0.000000 0.000000 0.000000-0.000105 0.000977 -0.000105 0.000000 0.000000 0.000000-0.000105 -0.000105 0.000977 0.000000 0.000000 0.0000000.000000 0.000000 0.000000 0.001786 0.000000 0.0000000.000000 0.000000 0.000000 0.000000 0.001786 0.0000000.000000 0.000000 0.000000 0.000000 0.000000 0.001786

Elastic stability criteria as seen in PRB 90, 224104 (2014).Criteria (i) C11 - C12 > 0 meeted.Criteria (ii) C11 + 2C12 > 0 meeted.Criteria (iii) C44 > 0 meeted.This Structure is Mechanically Stable.

Average mechanical properties for polycrystalline:+---------------------------------------------------------------+| Scheme | Voigt | Reuss | Hill |+---------------------------------------------------------------+| Bulk modulus K (GPa) | 434.431 | 434.431 | 434.431 || Shear modulus G (GPa) | 520.655 | 516.065 | 518.360 || Young's modulus E (GPa) | 1116.095 | 1109.045 | 1112.574 || P-wave modulus (GPa) | 1128.638 | 1122.517 | 1125.577 || Poisson's ratio v | 0.072 | 0.075 | 0.073 || Bulk/Shear ratio | 0.834 | 0.842 | 0.838 |+---------------------------------------------------------------+Pugh Ratio: 1.193Cauchy Pressure (GPa): -433.328Universal Elastic Anisotropy: 0.044Chung-Buessem Anisotropy: 0.004Isotropic Poisson's Ratio: 0.073Longitudinal wave velocity (in m/s): 17942.831Transverse wave velocity (in m/s): 12176.403Average wave velocity (in m/s): 13279.977Debye temperature (in K): 2212.733References:[1] Voigt W, Lehrbuch der Kristallphysik (1928)[2] Reuss A, Z. Angew. Math. Mech. 9 49-58 (1929)[3] Hill R, Proc. Phys. Soc. A 65 349-54 (1952)[4] Debye temperature J. Phys. Chem. Solids 24, 909-917 (1963)[5] Elastic wave velocities calculated using Navier's equation+---------------------------------------------------------------+

以下是利用VASPKIT结合VASP计算得到的结果与实验结果的对比，通过比较发现采用VASPKIT结合VASP得到的理论计算弹性常数与实验值符合较好。

‍‍

最后如果大家在研究中使用了VASPKIT请帮忙引用，格式如下：

V. WANG, N. XU, J.C. LIU, G. TANG, W.T. Geng, VASPKIT: A User-Friendly Interface Facilitating High-Throughput Computing and Analysis Using VASP Code, arXiv:1908.08269 (2019)。