VASP笔记之:计算德拜温度,杨氏模量,弹性矩阵
最近需要计算杨氏模量,但是上面三个量都是一起算出来的,so,一起记录一下笔记。
使用版本为VASP5.4.4,为了方便计算使用的是vaspkit1.2.1软件辅助自动生成的脚本进行的计算,微信公众号:学术之友,原文链接为:这里
下面我们以金刚石结构为例讲解如何采用应力-应变函数关系计算弹性常数,详见VASPKIT/examples/elastic/diamond_3D。由于金刚石具有立方晶体结构,一共有3个独立弹性常数C11、C12和C44 (不明白的请看原胞转化方法以及标准原胞在计算中的重要性)。
-
准备优化好的POSCAR文件,注意通常采用具有标准基矢形式的原胞计算弹性常数(VASPKIT-603/604可以生成标准结构),至于原因请看原胞转化方法以及标准原胞在计算中的重要性。
-
运行VASPKIT-102生成KPOINTS (注意精度要稍高一些)
3. 运行VASPKIT-101-DC生成INCAR文件,并根据实际情况修改,以下仅供参考:
Global Parameters
ISTART = 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 Relaxation
ISMEAR = 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 Relaxation
NELMIN = 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)
- 准备VPKIT.in文件并设置第一行为1 (预处理),运行VASPKIT并选择200, 将生成用于计算弹性常数的文件
1 ! 1 for prep-rocessing, 2 for post-processing
3D ! 2D for slab, 3D for bulk
7 ! 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!
- 批量提交vasp作业,相应的脚本如下(根据实际情况修改)
#!/bin/bash
root_path=`pwd`
for cij in `ls -F | grep /$`
do
cd ${root_path}/$cij
for s in strain_*
do
cd ${root_path}/$cij/$s
echo `pwd`
cp ../../vasp.job .
./vasp.job
# Add here your vasp_submit_job_script
done
done
- 等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-3m
Space Group: Fd-3m
Crystal System: Cubic system
Including Point group classes: 23, 2/m-3, 432, -43m, 4/m-32/m
There are 3 independent elastic constants
C11 C12 C12 0 0 0
C12 C11 C12 0 0 0
C12 C12 C11 0 0 0
0 0 0 C44 0 0
0 0 0 0 C44 0
0 0 0 0 0 C44
Stiffness Tensor C_ij (in GPa):
1050.316 126.488 126.488 0.000 0.000 0.000
126.488 1050.316 126.488 0.000 0.000 0.000
126.488 126.488 1050.316 0.000 0.000 0.000
0.000 0.000 0.000 559.816 0.000 0.000
0.000 0.000 0.000 0.000 559.816 0.000
0.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.000000
0.000000 0.000000 0.000000 0.001786 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.001786 0.000000
0.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.193
Cauchy Pressure (GPa): -433.328
Universal Elastic Anisotropy: 0.044
Chung-Buessem Anisotropy: 0.004
Isotropic Poisson's Ratio: 0.073
Longitudinal wave velocity (in m/s): 17942.831
Transverse wave velocity (in m/s): 12176.403
Average wave velocity (in m/s): 13279.977
Debye temperature (in K): 2212.733
References:
[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)。