材料的力学性质一直是DFT计算的重要方向，笔者在以往已有针对于静态结构的力学性质诸如弹性常数的相关计算，同时可通过VASPKIT借助相关方程导出力学性能。
bash+vasp+vaspkit能量应变计算弹性常数
vaspkit计算弹性常数的对称性指定
vaspkit计算弹性常数脚本补充
本文依据VASPKIT的EXAMPLE，简单介绍使用AIMD方法和VASPKIT软件计算含温力学性能。
材料选择为面心立方结构的Al为模型，空间群为 Fm-3m。
注意计算力学相关性质需要使用晶体学晶胞，而非原胞。
Generated by VASPKIT code1.0000004.0389299999999997 0.0000000000000000 0.00000000000000000.0000000000000000 4.0389299999999997 0.00000000000000000.0000000000000000 0.0000000000000000 4.0389299999999997Al4Direct0.0000000000000000 0.0000000000000000 0.0000000000000000 Al0010.0000000000000000 0.5000000000000000 0.5000000000000000 Al0020.5000000000000000 0.0000000000000000 0.5000000000000000 Al0030.5000000000000000 0.5000000000000000 0.0000000000000000 Al004
在计算AIMD时首先进行扩胞，这里使用2×2×2的超胞结构进行下一步计算。
手动生成VASPKIT计算所需要的设置文件INPUT.in，内容如下，参数详细说明附在后面，可根据需要自行更改。
1 ! 1 for prep-rocessing, 3 for finite-temperature MD3D ! 2D for slab, 3D for bulk4 ! number of strain-0.06 -0.03 0.03 0.06 ! magnitude of strain500 ! The fisrt 500 MD frames will be skipped
准备好VASP的输入文件后可运行VASPKIT来生成计算文件。
vaspkit -task 200
注意：通过AIMD计算弹性模量只支持应力应变方法，故这里使用vaspkit的200功能
vaspkit -task 200+---------------------------------------------------------------+| VASPKIT Standard Edition 1.5.1 (27 Jan. 2024) || Running VASPKIT Under Command-Line Mode |+---------------------------------------------------------------+-->> (01) Reading INPUT.in File.+---------------------------- Tip ------------------------------+| See some examples in vaspkit/examples/elastic. || Require the fully-relaxed and standard conventional cell. ||The stress-strain method requires higher ENCUT & denser K-Mesh.|+---------------------------------------------------------------+-> C11_C12_C44 folder created successfully!-> strain_-0.060 folder created successfully!-> strain_-0.030 folder created successfully!-> strain_+0.030 folder created successfully!-> strain_+0.060 folder created successfully!|---------------------------------------------------------------|| * CITATIONS * || When using VASPKIT in your research PLEASE cite the paper: || [1] 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, Computer Physics Communications || 267, 108033, (2021), DOI: 10.1016/j.cpc.2021.108033 |o---------------------------------------------------------------o
提交任务的推荐脚本为：
#!/bin/bashroot_path=`pwd`for cij in `ls -F | grep /$`docd ${root_path}/$cijfor s in strain_*docd ${root_path}/$cij/$secho `pwd`runvasp # Add here your vasp_submit_job_scriptdonedone
计算INCAR推荐
ENCUT = 400EDIFF = 1E-4ALGO = NormalIALGO = 48MAXMIX = 40IBRION = 0NSW = 2000NBLOCK = 1KBLOCK = 10POTIM = 2ISYM = 0# NVT ensembleISIF = 2SMASS = 2MDALGO = 2TEBEG = 400LREAL = FalseNELMIN = 4PREC = NormalISTART = 0ISMEAR = 2SIGMA = 0.2NPAR = 4NCORE = 1NSIM = 4NWRITE = 0LCHARG = .FALSE.LPLANE = .TRUE.LWAVE = .FALSE.IWAVPR = 11
计算完成后，更改INPUT.in文件：
再次运行VASPKIT200功能，获得结果数据
+---------------------------------------------------------------+| VASPKIT Standard Edition 1.5.1 (27 Jan. 2024) || Running VASPKIT Under Command-Line Mode |+---------------------------------------------------------------+-->> (01) Reading INPUT.in File.+---------------------------- Tip ------------------------------+| See some examples in vaspkit/examples/elastic. || Require the fully-relaxed and standard conventional cell. ||The stress-strain method requires higher ENCUT & denser K-Mesh.|+---------------------------------------------------------------+-->> (02) Calculating Fitting Coefficients of Stress vs. Strain.-->> (03) Reading Input Parameters From INCAR File.-->> Current directory: Fitting PrecisionC11_C12_C44 -> [XX]: 0.114E+03C11_C12_C44 -> [YY]: 0.103E+03C11_C12_C44 -> [ZZ]: 0.121E+03C11_C12_C44 -> [YZ]: 0.797E+01C11_C12_C44 -> [ZX]: 0.877E+00C11_C12_C44 -> [XY]: 0.000E+00+-------------------------- Summary ----------------------------+Crystal Class: m-3mSpace Group: Fm-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 C44Stiffness Tensor C_ij (in GPa):109.940 64.469 64.469 0.000 0.000 0.00064.469 109.940 64.469 0.000 0.000 0.00064.469 64.469 109.940 0.000 0.000 0.0000.000 0.000 0.000 31.301 0.000 0.0000.000 0.000 0.000 0.000 31.301 0.0000.000 0.000 0.000 0.000 0.000 31.301Compliance Tensor S_ij (in GPa^{-1}):0.016057 -0.005935 -0.005935 0.000000 0.000000 0.000000-0.005935 0.016057 -0.005935 0.000000 0.000000 0.000000-0.005935 -0.005935 0.016057 0.000000 0.000000 0.0000000.000000 0.000000 0.000000 0.031948 0.000000 0.0000000.000000 0.000000 0.000000 0.000000 0.031948 0.0000000.000000 0.000000 0.000000 0.000000 0.000000 0.031948Anisotropic mechanical properties of bulk singlecrystal:+---------------------------------------------------------------+| Mechanical Properties | Min | Max | Anisotropy||---------------------------------------------------------------|| Bulk Modulus B (GPa) | 79.626 | 79.626 | 1.000 || Young's Modulus E (GPa) | 62.279 | 82.993 | 1.333 || Shear Modulus G (GPa) | 22.738 | 31.301 | 1.377 || Poisson's Ratio v | 0.224 | 0.455 | 2.027 || Linear compressibility | 4.186 | 4.186 | 1.000 |+---------------------------------------------------------------+Linear compressibility beta is in unit of TPa^-1References:[1] Marmier A, Comput. Phys. Commun. 181, 2102–2115 (2010)[2] Gaillac R, J. Phys. Condens. Matter 28, 275201 (2016)Average mechanical properties of bulk polycrystal:+---------------------------------------------------------------+| Mechanical Properties | Voigt | Reuss | Hill ||---------------------------------------------------------------|| Bulk Modulus B (GPa) | 79.63 | 79.626 | 79.626 || Young's Modulus E (GPa) | 74.89 | 73.262 | 74.074 || Shear Modulus G (GPa) | 27.87 | 27.201 | 27.538 || Poisson's Ratio v | 0.34 | 0.347 | 0.345 || P-wave Modulus (GPa) | 116.79 | 115.895 | 116.344 || Pugh's Ratio (B/G) | 2.86 | 2.927 | 2.892 || Vickers Hardness (GPa)[6] | 1.09 | 0.922 | 1.007 || Vickers Hardness (GPa)[7] | 2.94 | 2.813 | 2.877 |+---------------------------------------------------------------+Pugh's Ratio (B/G): 2.89 --> Ductile region (> 1.75)Cauchy Pressure Pc (GPa): 33.2 --> Covalent-like bonding (< 0)Kleinman’s parameter: 0.98 --> Bond streching dominated (< 0.5)Universal Elastic Anisotropy: 0.12Chung-Buessem Anisotropy: 0.0Isotropic Poisson's Ratio: 0.34Longitudinal wave velocity (m/s): 6540.059Transverse wave velocity (m/s): 3181.823Average wave velocity (m/s): 3574.942Debye temperature (K): 418.3References:[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[6] Chen X-Q, Intermetallics 19, 1275 (2011)[7] Tian Y-J, Int. J. Refract. Hard Met. 33, 93–106 (2012)Eigenvalues of the stiffness matrix (in GPa):Eigenvalue lamda_1 = 31.301 > 0 meeted.Eigenvalue lamda_2 = 31.301 > 0 meeted.Eigenvalue lamda_3 = 31.301 > 0 meeted.Eigenvalue lamda_4 = 45.471 > 0 meeted.Eigenvalue lamda_5 = 45.471 > 0 meeted.Eigenvalue lamda_6 = 238.879 > 0 meeted.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.+---------------------------------------------------------------+-->> (04) Written ELASTIC_TENSOR File.|---------------------------------------------------------------|| * CITATIONS * || When using VASPKIT in your research PLEASE cite the paper: || [1] 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, Computer Physics Communications || 267, 108033, (2021), DOI: 10.1016/j.cpc.2021.108033 |o---------------------------------------------------------------o
计算过程中应力应变数据保存在的STRAIN_STRESS.dat中。
