VASP_AIMD+VASPKIT计算含温力学性质

材料的力学性质一直是DFT计算的重要方向，笔者在以往已有针对于静态结构的力学性质诸如弹性常数的相关计算，同时可通过VASPKIT借助相关方程导出力学性能。

bash+vasp+vaspkit能量应变计算弹性常数

vaspkit计算弹性常数的对称性指定

vaspkit计算弹性常数脚本补充

本文依据VASPKIT的EXAMPLE，简单介绍使用AIMD方法和VASPKIT软件计算含温力学性能。

材料选择为面心立方结构的Al为模型，空间群为 Fm-3m。

注意计算力学相关性质需要使用晶体学晶胞，而非原胞。

Generated by VASPKIT code 1.000000    4.0389299999999997    0.0000000000000000    0.0000000000000000    0.0000000000000000    4.0389299999999997    0.0000000000000000    0.0000000000000000    0.0000000000000000    4.0389299999999997   Al     4Direct    0.0000000000000000    0.0000000000000000    0.0000000000000000    Al001    0.0000000000000000    0.5000000000000000    0.5000000000000000    Al002    0.5000000000000000    0.0000000000000000    0.5000000000000000    Al003    0.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 /$`do  cd ${root_path}/$cij  for s in strain_*  do    cd ${root_path}/$cij/$s    echo `pwd`   runvasp # Add here your vasp_submit_job_script  donedone

    计算INCAR推荐

ENCUT   = 400EDIFF   = 1E-4ALGO    = NormalIALGO   = 48MAXMIX  = 40IBRION  = 0NSW     = 2000NBLOCK  = 1KBLOCK  = 10POTIM   = 2ISYM    = 0
# NVT ensembleISIF    = 2SMASS   = 2MDALGO  = 2TEBEG   = 400
LREAL   = FalseNELMIN  = 4PREC    = NormalISTART  = 0ISMEAR  = 2SIGMA   = 0.2
NPAR    = 4NCORE   = 1NSIM    = 4NWRITE  = 0
LCHARG  = .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 Precision     C11_C12_C44 -> [XX]:  0.114E+03    C11_C12_C44 -> [YY]:  0.103E+03    C11_C12_C44 -> [ZZ]:  0.121E+03    C11_C12_C44 -> [YZ]:  0.797E+01    C11_C12_C44 -> [ZX]:  0.877E+00    C11_C12_C44 -> [XY]:  0.000E+00 +-------------------------- Summary ----------------------------+ Crystal Class: m-3m Space Group: Fm-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):    109.940     64.469     64.469      0.000      0.000      0.000     64.469    109.940     64.469      0.000      0.000      0.000     64.469     64.469    109.940      0.000      0.000      0.000      0.000      0.000      0.000     31.301      0.000      0.000      0.000      0.000      0.000      0.000     31.301      0.000      0.000      0.000      0.000      0.000      0.000     31.301  Compliance 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.000000   0.000000   0.000000   0.000000   0.031948   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.031948   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.031948  Anisotropic 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^-1 References: [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.12 Chung-Buessem Anisotropy:     0.0 Isotropic Poisson's Ratio:   0.34 Longitudinal wave velocity (m/s):   6540.059 Transverse wave velocity (m/s):   3181.823 Average wave velocity (m/s):   3574.942 Debye temperature (K):   418.3 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 [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中。