3.8.1 金刚石Si，300K下NVT系综的AIMD模拟#

POSCAR：常压下金刚石Si晶胞（8个原子）扩胞222，超胞64个原子，abc三个方向长度大于10Å

INCAR：

# Global Parameters
ISTART = 1
ISPIN = 1
LREAL = A
ENCUT = 400
NWRITE = 2

# Electronic Relaxation
ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-05

# Molecular Dynamics
ISYM = 0
IBRION = 0
NSW = 10000
POTIM = 1
SMASS = 0
TEBEG = 300
TEEND = 300
ISIF = 2
MDALGO = 2
#PSTRESS = 0.1
#LANGEVIN_GAMMA = 10 10 10
#LANGEVIN_GAMMA_L = 10
#PMASS = 1000

#ML_LMLFF = TRUE
#ML_MODE = train

前面两组参数是控制DFT计算的，就是普通的scf计算参数；LREAL=A时可以对大体系计算进行加速，参数具体含义见官网。注意，如果体系需要考虑SOC或者DFT+U，或者需要用其它泛函做AIMD模拟，则INCAR中需要加上相应的参数。

单独强调，对于大部分磁性体系，一般跑MD不需要开磁性参数。因为一般MD都是在300K及以上模拟，很多磁性材料在这个温度（磁性材料有个转变温度）之上是非磁的，因此不开自旋跑MD或许更合理。

第三组是控制AIMD计算，NVT参数设置要求见：
https://www.vasp.at/wiki/index.php/NVT_ensemble

IBRION=0即为开启MD计算，此时ISYM=0（关闭对称性）；

NSW=10000，模拟1万步，与POTIM配套。（模拟时长一般10ps以上）

POTIM=1，一步即1fs；如果体系有H原子（比如含H₂O），POTIM=0.5

ISIF=2 固定晶格和体积，移动离子，用于NVT模拟

MDALGO=2 选用Nose-Hoover热浴控温，适用于NVT模拟

SMASS参数对模拟过程影响较大，一般来说需要测试不同的SMASS。关于该参数的讨论，见本小节最后。

TEBEG和TEEND就是模拟的温度，单位K。前者是起始温度，后者是最后温度。二者一致时，就是恒温的模拟；如果不一致，就是升温或降温过程。升温/降温的模拟一般与NPT系综参数搭配，可以用来模拟体系的融化/结晶过程。

分子动力学中EDIFFG参数不起作用，即使你设置了，程序也会忽略掉。

INCAR中可以加入NCORE/NPAR/KPAR等并行参数，来提高计算效率。需要自行针对超胞做测试。如果不想测试，那就取NCORE=2或者4（任务调用的核数的因数，要取小一些，KPAR=2）下面是网上一些测试结果：

https://mp.weixin.qq.com/s/uiWuvILE1WeiB5SC3_Frvg

NWRITE=2（默认值）时，会在OUTCAR中输出每一步MD的原子受力，如果你后面需要提取AIMD的数据用于后续训练机器学习势，强烈建议设置为2，这样方便后面dpdata直接提取DFT数据。

有时候算MD时还会设置NBLOCK参数，默认为1，含义是把每一步的原子位置都输出到XDATCAR文件中，方便后续查看每一步动力学过程。如果你跑一个长时间的MD，比如几十万步，为避免输出的XDATCAR文件冗长，可以设置为10或者100，表示每隔10步把原子位置输出到XDATCAR。但建议取1.

AIMD计算的INCAR文件也可以用vaspkit生成，再自己添加/修改参数。再次强调，当你用vaspkit辅助生成INCAR时，务必清楚INCAR中每一个参数的含义以及作用！！

KPOINTS：原则上需要对超胞做能量收敛性测试得到合适的Kmesh；为了演示，我们这里只取Gamma点。

K-Spacing Value to Generate K-Mesh: 0.025
0
Monkhorst-Pack
1 1 1
0.0 0.0 0.0

POTCAR取默认的4价电子的赝势。

此时我们调用vasp的gam版本，可以对k111的计算进行加速。需要在任务提交脚本里面修改。之前我们算静态或者结构优化时，很少设置K111的情况，所以一般是调用vasp_std版本。这里调用gam版本。如果没有安装gam，找组里人或者工程师询问。

e/vasp.6.4.0/bin/vasp_std
vasp.6.4.0/bin/vasp_gam
e/vasp.6.4.0/bin/vasp_ncl

分子动力学的计算一般很耗时，因此在有些对精度要求不太高的情况下，通常会降低计算精度。比如ECUT和K点。ECUT一般取到不同元素的ENMAX中最大值多一些，除了Li之外，常见元素的ENAMX＜400。K点只取gamma点。

计算时，每一步的能量和温度信息会输出在OSZICAR文件

DAV:  13    -0.340518831657E+03    0.45111E-05   -0.38335E-06  212   0.423E-03
      1 T=   300  E= -.33807580E+03 F=-.34051883E+03E0=-.34051883E+03    EK= 0.24430E+01 SP= 0.00E+00 SK= 0.00E+00
       N       E                     dE            d eps       ncg     rms          rms(c)
DAV:   1    -0.340505398329E+03    0.13438E-01   -0.40360E-01  332   0.121E+00    0.329E-01
DAV:   2    -0.340507812474E+03   -0.24141E-02   -0.27415E-02  380   0.338E-01    0.119E-01
DAV:   3    -0.340507612268E+03    0.20021E-03   -0.88628E-04   380  0.692E-02    0.709E-02
DAV:   4    -0.340507578783E+03    0.33485E-04   -0.38057E-04  380   0.371E-02    0.225E-02
DAV:   5    -0.340507583906E+03   -0.51233E-05   -0.56890E-05  308   0.156E-02
      2 T=   299.E=-.33807558E+03F=-.34050758E+03 E0=-.34050758E+03      EK= 0.24320E+01 SP= 0.15E-06 SK= 0.74E-12
       N       E                     dE            deps        ncg     rms          rms(c)
DAV:   1     -0.340474888509E+03   0.32690E-01   -0.30478E-06  219   0.408E-01    0.294E-02
DAV:   2    -0.340474912497E+03   -0.23988E-04   -0.24059E-04  364   0.308E-02    0.126E-02
DAV:   3    -0.340474911772E+03    0.72542E-06   -0.46291E-06  219   0.523E-03
      3 T=   295.E=-.33807534E+03F= -.34047491E+03 E0=-.34047491E+03     EK= 0.23998E+01 SP=-.27E-03 SK= 0.24E-05
       N       E                     dE            d eps       ncg     rms          rms(c)
DAV:   1    -0.340422612114E+03    0.52300E-01   -0.15243E-04  419   0.412E-01    0.366E-02
DAV:   2    -0.340422650853E+03   -0.38739E-04   -0.50161E-04  356   0.419E-02    0.301E-02

每一步的结构信息在XDATCAR文件（可以借助OVITO软件可视化）

每一步的结构和原子速度以及Predictor-corrector coordinates信息在CONTCAR中 https://www.vasp.at/wiki/index.php/CONTCAR

我们一般关注OSZICAR和XDATCAR文件，前者反映温度和能量波动，后者反映结构的动力学演化情况。当然还有包含各种信息的OUTCAR

（1）OSZICAR中的能量

E即为体系总能量，是电子-离子体系的自由能F与离子动能EK的和，再加上恒温器的势能SP与动能SK。其中，离子动能Ek遵循经典统计力学的能量均分定理，等于3/2(N-1)kT（N为原子数），要扣除三个体系平动的自由度，因此要N-1。如果是NPT模拟，那F还包含了体积项pV。

https://www.vasp.at/wiki/index.php/OSZICAR

一般在评估AIMD模拟过程体系的总能量时，以F或者E0为准。在计算结束时，可以采用grep命令从OSZICAR文件中提取温度/能量随时间变化的数据：

grep E0 OSZICAR |awk '{print $1 , $3, $9}'> ./Energy.dat

当然，能量信息也可以从OUTCAR中抓取每一次离子步之后的能量。可自行查找。

（2）模拟过程中的温度波动

AIMD模拟时，最开始20步以内，温度可能会从预设温度（300K）迅速下降到100K，后面又升高到500K，50步以后稳定在300K附近波动。对于2D材料，这个波动幅度可能有±100K；这是正常现象。

但有些体系，初始100步内，温度波动可能非常剧烈，比如预设是5000K，波动可能到9000K，这是因为初始结构原子位置太规整导致，这某些情况下可能会使得体系一开始温度太高而不稳定，比如变形或者融化，说明该模拟是不对的（我们需要的是在5000K附近的模拟，不应该到9000K），需要对初始的原子位置做一些微扰，可以考虑在低温下先NVT跑个100步得到带有原子位置扰动的结构，然后将CONTCAR拿去在高温下模拟；或者用vaspkit或是phonopy的原子位置扰动功能对原子位置做~0.01Å微小扰动。注意MD模拟过程中体系的对称性几乎都是P1，即完全没有晶格对称性了。

（3）AIMD模拟中的压强

AIMD模拟过程中的压强有两部分，一部分是电子-离子之间相互作用导致的内压，有的地方称为virial，另一部分是动能带来的压强；后者遵循PV=NRT，两部分压强可以在OUTCAR中找到：

FORCE on cell =-STRESS in cart. coord. units (eV):
Direction    XX          YY         ZZ          XY          YZ         ZX
Alpha Z  7689.59288  7689.59288 7689.59288
Ewald  -19475.06528-19247.11831-19859.91306   50.12228   149.64562   -27.55808
Hartree  2271.85963  2318.27494 1895.11861    -0.53742    42.27420   -17.66811
E(xc)   -3689.86335 -3690.99965 -3689.18958   -0.40733    -1.07063    -0.78619
Local    4162.37186  3930.80116 4982.23008   -28.14883  -166.17826   51.19869
n-local -2212.50882 -2210.92263 -2190.58966    3.88774     1.22559    -8.43636
augment  2500.97895  2490.13329 2479.45090    -6.76952    -6.64571    -2.62591
Kinetic  9836.31013  9827.90348 9822.70663   -11.10476     0.68046    9.33200
Fock        0.00000     0.00000    0.00000     0.00000     0.00000    0.00000
Total    1083.67599  1107.66516 1129.40679     7.04216    19.93127    3.45603
in kB    3024.54358  3091.49744 3152.17842    19.65468    55.62823    9.64579
external pressure :     3089.41 kB Pullay stress =        0.00 kB
kinetic pressure (ideal gas correction)     104.57 kB
total pressure      3193.98 kB
Total+kin.  3129.794    3197.161   3254.986      25.384      42.608     -0.244

因此，如果要读取MD过程中的压强，应该读取这里的total pressure或者Total+kin部分。另外，需要强调的是，in kB部分的内压强本身也和热波动有关，即，它并不等于你用一个完美格点的原子结构做一次静态得到的压强，而是会大一些。

（4）SMASS参数的影响

官网对SMASS参数的说明：https://www.vasp.at/wiki/index.php/SMASS

SMASS ≥ 0

For SMASS ≥ 0, a canonical ensemble is simulated using the algorithm of Nosé. The Nosé mass controls the frequency of the temperature oscillations during the simulation.^[1][2][3] For SMASS=0, a Nosé-mass corresponding to period of 40 time steps will be chosen. The Nosé-mass should be set such that the induced temperature fluctuation show approximately the same frequencies as the typical 'phonon'-frequencies for the specific system. For liquids something like 'phonon'-frequencies might be obtained from the spectrum of the velocity auto-correlation function. If the ionic frequencies differ by an order of magnitude from the frequencies of the induced temperature fluctuations, Nosé thermostat and ionic movement might decouple leading to a non-canonical ensemble. The frequency of the approximate temperature fluctuations induced by the Nosé-thermostat is written to the OUTCAR file.

实测同一体系，不同SMASS的影响，如下图；上面的图是温度变化，下面是能量。SMASS=2的模拟显然是不对的，因为温度没有控制好。因此在一个AIMD模拟结束后，需要确认温度的波动是否合理（右下）。

SMASS似乎会影响系统达到平衡所需的模拟时间。

https://journals.aps.org/prb/abstract/10.1103/PhysRevB.46.13756