分子動(dòng)力學(xué)模擬分為兩種，一種是基于力場(chǎng)的經(jīng)典分子動(dòng)力學(xué)模擬。通過(guò)施加力場(chǎng)以使體系在勢(shì)能面上游走。其常用軟件為（GROMACS和 LAMMPS）。另一種是基于量子化學(xué)、第一性原理的從頭算分子動(dòng)力學(xué)模擬,通用軟件為AIMD（Ab-initio Molecular Dynamics）。本期我們將講述AIMD的基本操作。

在AIMD中，我們首先要確定原子坐標(biāo)。這一般通過(guò)結(jié)構(gòu)優(yōu)化得到。其次，初始的POSCAR必須是擴(kuò)得足夠大的超胞，含有足夠多的原子數(shù)。相對(duì)于在勢(shì)能面上尋找極小點(diǎn)的結(jié)構(gòu)優(yōu)化，AIMD根據(jù)牛頓定律使體系隨機(jī)地在勢(shì)能面上運(yùn)動(dòng)，進(jìn)行退火模擬。

因此，我們可為原子提供初始速度，它們?cè)赪ycoff位置后的一個(gè)獨(dú)立段落中。如果沒(méi)有提供初始速度，系統(tǒng)則會(huì)在計(jì)算開(kāi)始時(shí)假定隨機(jī)速度。但要注意的是，由于初始速度的隨機(jī)性，從不同的計(jì)算中得到的軌跡是很難比較的。

以下是在AIMD中，INCAR的主要參數(shù)設(shè)置：

IBRION=0：通過(guò)將IBRION設(shè)置為0，可以啟用分子動(dòng)力學(xué)計(jì)算。

POTIM：設(shè)置分子動(dòng)力學(xué)運(yùn)行的時(shí)間步長(zhǎng)（fs）。

NSW：設(shè)置分子動(dòng)力學(xué)運(yùn)行的步驟數(shù)，單位為fs。ps級(jí)別的退火模擬，一般設(shè)在NSW=5000~20000，最長(zhǎng)20ps就足夠了。

TEBEG：如果使用恒溫器，則定義分子動(dòng)力學(xué)計(jì)算運(yùn)行的初態(tài)溫度。

TEEND：分子動(dòng)力學(xué)計(jì)算運(yùn)行的末態(tài)溫度。

SMASS：由于AIMD計(jì)算的尺度較小，體系的溫度并不是一個(gè)確定值。SMASS確定了耦合強(qiáng)度，以此控制體系溫度的漲落范圍。強(qiáng)度越大，溫度波動(dòng)越小。通常使用SMASS=2，Nosé算法控制溫度變化。

ISIF：在分子動(dòng)力學(xué)計(jì)算中，這個(gè)標(biāo)簽用于選擇NVT系綜或NpT系綜（NVE系綜為特殊情況）。對(duì)于ISIF=2，體積保持不變，使用NVT系綜。使用這個(gè)標(biāo)簽可以計(jì)算出應(yīng)力張量，因此可以監(jiān)測(cè)體系壓力。對(duì)于ISIF=3，應(yīng)力張量（壓力）保持不變，使用NpT系綜。使用這個(gè)標(biāo)簽可以計(jì)算出體積，并對(duì)之進(jìn)行監(jiān)測(cè)。NVE系綜是特殊情況。它可以通過(guò)選擇Andersen恒溫器、與熱浴的非碰撞（ANDERSEN_PROB=0）來(lái)實(shí)現(xiàn)。

MDALGO：決定AIMD計(jì)算在哪個(gè)恒溫器下執(zhí)行。對(duì)于常規(guī)的AIMD計(jì)算，恒溫器可通過(guò)數(shù)字選取（例如，1：Andersen，2：Nosé-Hoover等）。對(duì)于Biased Molecular Dynamics， Metadynamics等，恒溫器的選擇與VASP 6或更高版本相同。在VASP 5.x中，它是由一個(gè)兩位數(shù)的數(shù)字選擇的，其中第一位數(shù)字對(duì)應(yīng)于恒溫器，類似于常規(guī)分子動(dòng)力學(xué)，第二位數(shù)字對(duì)應(yīng)于分子動(dòng)力學(xué)類型（例如，11：Andersen恒溫器下的Metadynamics，21：Nosé-Hoover恒溫器下的Metadynamics等）。

以下為AIMD系綜、恒溫器及恒壓器的可能選?。?/strong>

The following combinations of thermostats and barostats is possible:

Thermostat

Ensemble

Andersen

Nose-Hoover

Langevin

Multiple Andersen

NVE

MDALGO=1,?ANDERSEN_PROB=0.0

NVT

MDALGO=1

MDALGO=2

MDALGO=3

MDALGO=13

ISIF=2

ISIF=2

ISIF=2

ISIF=2

NpT

not available

not available

MDALGO=3

not available

ISIF=3

NpH

MDALGO=3,?ISIF=3,?LANGEVIN_GAMMA_L=0.0

AIMD計(jì)算耗時(shí)較長(zhǎng)，運(yùn)行個(gè)一兩周左右是正常的，所以需要耐心等待。計(jì)算完成后，在linux窗口輸入指令grep E= OSZICAR |awk '{print $1,$9}' > 1.dat，即可得到數(shù)據(jù)。第一列是退火模擬時(shí)長(zhǎng)，以fs為單位，通常轉(zhuǎn)換成ps為橫坐標(biāo)單位。第二列是退火過(guò)程中體系的能量漲落，除以體系原子數(shù)，轉(zhuǎn)換成縱坐標(biāo)單位eV/atom。

以下示圖為PdSeO3在500K下的退火模擬，能量曲線在合理的波動(dòng)范圍內(nèi)。

關(guān)于第一性原理的分子動(dòng)力模擬，更完整的理論請(qǐng)參考：

https://www.vasp.at/wiki/index.php/Category:Molecular_dynamics及以下文獻(xiàn)。

參考文獻(xiàn)：

[1] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, London, 1996).

[2] M. P. Allen and D. J. Tildesley, Computer simulation of liquids (Oxford university press: New York, 1991).

[3] M. Qiao, J. Liu, Y. Wang, et al. J. Am. Chem. Soc., 2018, 140, 12256-12262.

標(biāo)簽：