引言

分子動(dòng)力學(xué)模擬是一套數(shù)值模擬方法。該方法通過(guò)數(shù)值求解牛頓運(yùn)動(dòng)方程的方式，來(lái)獲得系統(tǒng)中每個(gè)原子的位置及速度隨時(shí)間的變化，從而可以獲取我們所關(guān)心的物理量，如比熱容、擴(kuò)散系數(shù)、熱導(dǎo)率等。簡(jiǎn)單的說(shuō)，我們可以把分子動(dòng)力學(xué)模擬當(dāng)作在計(jì)算機(jī)上面完成的物理實(shí)驗(yàn)。但相比于真正的物理實(shí)驗(yàn)，分子動(dòng)力學(xué)模擬可以獲取更小空間尺度（nm）和時(shí)間尺度（ps）上的信息。

在分子動(dòng)力學(xué)模擬中，我們把每個(gè)原子當(dāng)作剛性粒子，并通過(guò)經(jīng)典作用勢(shì)函數(shù)來(lái)描述原子間的相互作用。通過(guò)選取一個(gè)足夠小的時(shí)間步長(zhǎng)Δt來(lái)對(duì)連續(xù)運(yùn)動(dòng)方程進(jìn)行離散化，從而可以跟蹤系統(tǒng)中每個(gè)原子的位置和速度隨時(shí)間變化的軌跡。在實(shí)際模擬中，一個(gè)系統(tǒng)可能包含幾千到幾十萬(wàn)個(gè)原子，同時(shí)可能需要模擬幾十萬(wàn)步或者更多，所以我們需要用專門(mén)的可進(jìn)行并行計(jì)算的仿真軟件對(duì)系統(tǒng)進(jìn)行模擬。

分子動(dòng)力學(xué)模擬

最常用的分子動(dòng)力學(xué)仿真軟件為L(zhǎng)AMMPS，即為L(zhǎng)arge-scale Atomic/Molecular Massively Parallel Simulator的縮寫(xiě)。LAMMPS為美國(guó)Sandia國(guó)家實(shí)驗(yàn)室所開(kāi)發(fā)的開(kāi)源仿真軟件。經(jīng)過(guò)多年的開(kāi)發(fā)及完善，LAMMPS如今已具備強(qiáng)大的仿真功能。但對(duì)于一些比較特殊的問(wèn)題，我們?nèi)匀恍枰约壕帉?xiě)代碼或?qū)ζ湓创a進(jìn)行修改才可以完成。一個(gè)典型的分析動(dòng)力學(xué)仿真流程如下圖所示：

需要給定的輸入信息為仿真系統(tǒng)的結(jié)構(gòu)文件（包含每個(gè)原子的初始位置及速度）以及原子之間的相互作用勢(shì)函數(shù)。有了這些輸入信息之后，就可以在模擬過(guò)程中跟蹤每個(gè)原子的位置及速度隨時(shí)間變化過(guò)程。需要注意的是在具體模擬過(guò)程中，我們可能需要對(duì)系統(tǒng)的狀態(tài)進(jìn)行控制，比如控制溫度、壓強(qiáng)、原子的受力等，以達(dá)到特定的仿真目的。有了系統(tǒng)中原子位置和速度隨時(shí)間的演變信息，我們就可以在后處理過(guò)程中運(yùn)用特定的理論或定律來(lái)求解我們所關(guān)心的物理量，如系統(tǒng)的比熱容，擴(kuò)散系數(shù)，熱導(dǎo)率等等。

在上述流程圖中，最關(guān)鍵的一步是求解牛頓運(yùn)動(dòng)方程，得到下一時(shí)刻原子的位置和速度。在分子動(dòng)力學(xué)模擬中最常用的積分算法為Velocity Verlet積分算法。因其具有精度高、穩(wěn)定性好的特點(diǎn)，而被廣泛的使用。其算法積分步驟如下圖所示：

1. 獲取t時(shí)刻原子的位置、速度及受力。
2. 根據(jù)t時(shí)刻原子的速度及受力，計(jì)算t+Δt/2（半步長(zhǎng)）時(shí)刻每個(gè)原子的速度。
3. 通過(guò)t時(shí)刻的原子的位置及t+Δt/2時(shí)刻原子的速度，計(jì)算在t+Δt（整步長(zhǎng)）時(shí)刻原子的位置。
4. 根據(jù)原子在t+Δt時(shí)刻所在位置的受力情況，計(jì)算t+Δt（整步長(zhǎng)）時(shí)刻的速度。
5. 更新原子在t+Δt時(shí)刻的所有信息，重復(fù)上述過(guò)程。

LAMMPS輸入例程

在LAMMPS 模擬中，我們可以通過(guò)輸入script 的方式來(lái)告訴LAMMPS我們想要具體執(zhí)行的模擬過(guò)程。輸入script中每一行表示一條指令，來(lái)告訴LAMMPS執(zhí)行具體的操作。關(guān)于LAMMPS中具體指令的使用方法，可參考LAMMPS的官方網(wǎng)站。下面以一個(gè)具體的例子來(lái)簡(jiǎn)單介紹下LAMMPS 的使用。
（注：滑動(dòng)屏幕可以顯示一行中不完整的代碼信息）

 1## 基本輸入信息，一行中以# 開(kāi)始的部分表示注釋
 2units metal # 定義仿真過(guò)程中用到的單位系統(tǒng) ?
 3boundary p p p # 定義邊界條件，p p p 表示在三個(gè)方向上都使用周期性邊界條件
 4atom_style atomic # 定義原子的類(lèi)型及屬性
 5# 讀取結(jié)構(gòu)文件“structure.d”，里面包含仿真結(jié)構(gòu)模型的具體信息
 6read_data structure.d
 7## 定義原子間的作用勢(shì)函數(shù)
 8pair_style tersoff #作用勢(shì)類(lèi)型，我們的仿真系統(tǒng)為Si, 所以用Tersoff三體作用勢(shì)
 9pair_coeff * * SiC_Lindsay.tersoff Si Si # 所用勢(shì)的具體參數(shù)
10## 定義不同的group，后面的仿真中可能會(huì)對(duì)不同group進(jìn)行不同的操作
11group left_end id <> 1 288 ?# 我們將輸入文件中ID 為1到288 號(hào)之間的原子定義成group left_end, 后面對(duì)該group會(huì)進(jìn)行具體的操作
12group left_heat id <> 289 864
13group right_end id <> 29089 29376
14group right_heat id <> 28513 29088
15group mobile subtract all left_end right_end # 定義group mobile
16timestep 0.0005 # 定義仿真步長(zhǎng)， 單位為ps
17## 仿真過(guò)程的開(kāi)始，一個(gè)"fix" 表示對(duì)系統(tǒng)進(jìn)行某項(xiàng)控制，如控制能量、溫度等
18# 給初始結(jié)果一定的初始速度分布，使得系統(tǒng)的初始溫度為300 K
19velocity all create 300 34678 sum no mom yes rot yes dist gaussian
20# 用Langevin 控溫器對(duì)系統(tǒng)的溫度進(jìn)行控制
21fix 1 all langevin 300 300 0.1 48279
22fix 2 all nve/limit 0.1
23# 對(duì)系統(tǒng)的位子信息進(jìn)行輸出，可以用后期可視化軟件如VMD對(duì)其進(jìn)行觀看
24dump 1 all atom 100000 initial.lammpstrj
25# 對(duì)仿真過(guò)程中的熱力學(xué)信息進(jìn)行輸出監(jiān)測(cè)
26thermo_style custom step temp press pe ke etotal lx ly lz
27thermo 1000
28# 對(duì)系統(tǒng)進(jìn)行仿真 1000000 步
29run 1000000
30unfix 1 # 每一個(gè)fix 過(guò)程完需要unfix
31unfix 2
32undump 1

在該系列的隨后文章中，我們將介紹怎么用LAMMPS 來(lái)計(jì)算納米尺度傳熱的一些具體問(wèn)題，如怎么計(jì)算材料的熱導(dǎo)率、界面熱阻等，以及怎么通過(guò)分子動(dòng)力學(xué)仿真的方式來(lái)獲取更深層次的物理信息，敬請(qǐng)關(guān)注！