分子动力学模拟与纳米尺度传热(一)：LAMMPS模拟软件入门

引言

分子动力学模拟是一套数值模拟方法。该方法通过数值求解牛顿运动方程的方式，来获得系统中每个原子的位置及速度随时间的变化，从而可以获取我们所关心的物理量，如比热容、扩散系数、热导率等。简单的说，我们可以把分子动力学模拟当作在计算机上面完成的物理实验。但相比于真正的物理实验，分子动力学模拟可以获取更小空间尺度（nm）和时间尺度（ps）上的信息。

在分子动力学模拟中，我们把每个原子当作刚性粒子，并通过经典作用势函数来描述原子间的相互作用。通过选取一个足够小的时间步长Δt来对连续运动方程进行离散化，从而可以跟踪系统中每个原子的位置和速度随时间变化的轨迹。在实际模拟中，一个系统可能包含几千到几十万个原子，同时可能需要模拟几十万步或者更多，所以我们需要用专门的可进行并行计算的仿真软件对系统进行模拟。

分子动力学模拟

最常用的分子动力学仿真软件为LAMMPS，即为Large-scale Atomic/Molecular Massively Parallel Simulator的缩写。LAMMPS为美国Sandia国家实验室所开发的开源仿真软件。经过多年的开发及完善，LAMMPS如今已具备强大的仿真功能。但对于一些比较特殊的问题，我们仍然需要自己编写代码或对其源代码进行修改才可以完成。一个典型的分析动力学仿真流程如下图所示：

需要给定的输入信息为仿真系统的结构文件（包含每个原子的初始位置及速度）以及原子之间的相互作用势函数。有了这些输入信息之后，就可以在模拟过程中跟踪每个原子的位置及速度随时间变化过程。需要注意的是在具体模拟过程中，我们可能需要对系统的状态进行控制，比如控制温度、压强、原子的受力等，以达到特定的仿真目的。有了系统中原子位置和速度随时间的演变信息，我们就可以在后处理过程中运用特定的理论或定律来求解我们所关心的物理量，如系统的比热容，扩散系数，热导率等等。

在上述流程图中，最关键的一步是求解牛顿运动方程，得到下一时刻原子的位置和速度。在分子动力学模拟中最常用的积分算法为Velocity Verlet积分算法。因其具有精度高、稳定性好的特点，而被广泛的使用。其算法积分步骤如下图所示：

1. 获取t时刻原子的位置、速度及受力。
2. 根据t时刻原子的速度及受力，计算t+Δt/2（半步长）时刻每个原子的速度。
3. 通过t时刻的原子的位置及t+Δt/2时刻原子的速度，计算在t+Δt（整步长）时刻原子的位置。
4. 根据原子在t+Δt时刻所在位置的受力情况，计算t+Δt（整步长）时刻的速度。
5. 更新原子在t+Δt时刻的所有信息，重复上述过程。

LAMMPS输入例程

在LAMMPS 模拟中，我们可以通过输入script 的方式来告诉LAMMPS我们想要具体执行的模拟过程。输入script中每一行表示一条指令，来告诉LAMMPS执行具体的操作。关于LAMMPS中具体指令的使用方法，可参考LAMMPS的官方网站。下面以一个具体的例子来简单介绍下LAMMPS 的使用。
（注：滑动屏幕可以显示一行中不完整的代码信息）

 1## 基本输入信息，一行中以# 开始的部分表示注释
 2units metal # 定义仿真过程中用到的单位系统  
 3boundary p p p # 定义边界条件，p p p 表示在三个方向上都使用周期性边界条件
 4atom_style atomic # 定义原子的类型及属性
 5# 读取结构文件“structure.d”，里面包含仿真结构模型的具体信息
 6read_data structure.d
 7## 定义原子间的作用势函数
 8pair_style tersoff #作用势类型，我们的仿真系统为Si, 所以用Tersoff三体作用势
 9pair_coeff * * SiC_Lindsay.tersoff Si Si # 所用势的具体参数
10## 定义不同的group，后面的仿真中可能会对不同group进行不同的操作
11group left_end id <> 1 288  # 我们将输入文件中ID 为1到288 号之间的原子定义成group left_end, 后面对该group会进行具体的操作
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 # 定义仿真步长， 单位为ps
17## 仿真过程的开始，一个"fix" 表示对系统进行某项控制，如控制能量、温度等
18# 给初始结果一定的初始速度分布，使得系统的初始温度为300 K
19velocity all create 300 34678 sum no mom yes rot yes dist gaussian
20# 用Langevin 控温器对系统的温度进行控制
21fix 1 all langevin 300 300 0.1 48279
22fix 2 all nve/limit 0.1
23# 对系统的位子信息进行输出，可以用后期可视化软件如VMD对其进行观看
24dump 1 all atom 100000 initial.lammpstrj
25# 对仿真过程中的热力学信息进行输出监测
26thermo_style custom step temp press pe ke etotal lx ly lz
27thermo 1000
28# 对系统进行仿真 1000000 步
29run 1000000
30unfix 1 # 每一个fix 过程完需要unfix
31unfix 2
32undump 1

在该系列的随后文章中，我们将介绍怎么用LAMMPS 来计算纳米尺度传热的一些具体问题，如怎么计算材料的热导率、界面热阻等，以及怎么通过分子动力学仿真的方式来获取更深层次的物理信息，敬请关注！