1 输入文件
输入脚本文件（命令文件？）要包含五个部分：
# 1) Initialization
# 2) System definition
# 3) Simulation settings
# 4) Visualization
# 5) Run
这五个部分不一定全都需要，也不一定完全按如下顺序。例如3和4可以翻过来，另外还可以多次重复。
1.1 初始化
# 1) Initialization
units lj
dimension 2
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
第一句采用lj归一化单位；第二句表示二维模型；第三句表示采用原子质点模型；第四句表示原子间作用力模型为LJ，截断值2.5；最后一行表示xyz三个方向采用周期性边界条件。
（二维模型三个周期性边界？It may seems strange to define the boundary conditions along the third dimension for a 2D simulation, but it is a LAMMPS requirement.）
cutoff是什么。
1.2 系统定义
region myreg block -30 30 -30 30 -0.5 0.5
create_box 2 myreg
create_atoms 1 random 1500 341341 myreg
create_atoms 2 random 100 127569 myreg
第一行构建模拟区域（几何模型），称为myreg。模拟区域为block，长方体，xyz范围如上。
第二行定义模拟盒子（？），其中含有两种原子。
第三行定义第一种原子随机布置，有1500个，341341是随机数种子。后续可以更改，以便进行不同的仿真。
第四行定义第二种原子，有100个。
1.3 仿真设置
# 3) Simulation settings
mass 1 1
mass 2 1
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 0.5 3.0
头两行定义原子质量，均为1。
后两行定义原子的作用势参数：势阱和势能为0的距离。注意，lammps用几何平均计算不同种类原子的相互作用参数，即：
$$\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}
\space
\space
\sigma_{ij} = \sqrt{\sigma_{ii} \sigma_{jj}}$$
这只是一种权益之计，不过它不影响仿真结果。
也可以采用pair_coeff 1 2明确指定相互作用参数。
1.4 可视化和运行
# 4) Visualization
thermo 10
# 5) Run
minimize 1.0e-4 1.0e-6 1000 10000
可视化时，thermo 表示输出热力学信息，如温度和能量等，10表示每10步输出一次。minimize表示进行能量最小化，约束条件是两次迭代之间的能量差值小于1e-4，两个原子之间的力的最大值小于1e-6，迭代次数大于1000次，力和能量计算次数超过10000次。
运行后输出结果如下
Step Temp E_pair E_mol TotEng Press
0 0 5.8997404e+14 0 5.8997404e+14 1.5732641e+15
10 0 56275037 0 56275037 1.5007118e+08
20 0 17731.329 0 17731.329 47458.738
30 0 350.68529 0 350.68529 972.14134
40 0 13.745948 0 13.745948 48.748312
50 0 0.5033657 0 0.5033657 8.3398718
60 0 -1.4427524 0 -1.4427524 1.1131474
70 0 -1.7136665 0 -1.7136665 -0.038162464
80 0 -1.7516642 0 -1.7516642 -0.15686171
81 0 -1.7518285 0 -1.7518285 -0.15730928
在步长零时，能量很高，得到1e+14量级。可能是因为粒子位置随机产生，有重叠，导致斥力大，能量高。随着计算进行，能量迅速减小。
lammps输出还有
Minimization stats: Stopping criterion = energy tolerance
说明满足了迭代中能量的条件。
lammps还会输出
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization
这是因为lammps认为构建粒子邻居列表比默认更为频繁，所以输出了warning。可忽略。
2 动力学模拟
# PART B - MOLECULAR DYNAMICS #
4) Visualization
thermo 1000
variable kinetic_energy equal ke
variable potential_energy equal pe
variable pressure equal press
fix myat1 all ave/time 10 1 10 v_kinetic_energy v_potential_energy v_pressure file energy.dat
# 5) Run
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917
fix myefn all enforce2d
timestep 0.005
run 10000
在前面写好的命令后面加上如上面命令。因为lammps执行命令顺序从上到下，所以这些命令是在能量最小化之后执行的。后面的thermo命令执行时，热力学信息输出改为1000步输出一次。
然后用variable命令定义了三个变量Kinetic_energy、potential energy, and the pressure，输出文件为energy.dat。
下面run部分，要求所有原子按照nve系综更新位置和速度。第二句采用郎之万热浴控温，温度为1，阻尼系数（？）为0.1，随机数种子1530917。然后enforce2D要求二维模拟。最后设置时间步长和模拟步数。
模拟输出
可以看到温度很快升到1。
此外，还输出了如下信息
可以看出，998个dangerous邻居列表。这说明两次计算邻居列表中原子移动的距离过大。在模拟设置部分需要增加如下命令：
neigh_modify every 1 delay 5 check yes
这样lammps就能够更频繁的构建邻居列表。
现在结果就正确了。
3 轨迹可视化
如果我们需要把轨迹可视化，需要定时输出dump原子的位置。在visualization一段加上如下命令：
dump mydmp all atom 1000 dump.lammpstrj
再次运行lammps。输出的文件为dump.lammpstrj，可用vmd或ovito可视化。
4 改进
比如，我们首先创建一个region，然后在其基础上创建一个box。再另外创建两个region，其中分别创建类型1和类型2的原子。
# 2) System definition
region mybox block -30 30 -30 30 -0.5 0.5
create_box 2 mybox
region mycylin cylinder z 0 0 15 INF INF side in
region mycylou cylinder z 0 0 15 INF INF side out
create_atoms 1 random 1000 341341 mycylou
create_atoms 2 random 150 127569 mycylin
side in和side out分别表示在圆柱内和圆柱外。
# 3) Simulation settings
mass 1 1
mass 2 1
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 0.5 3.0
neigh_modify every 1 delay 5 check yes
# 4) Visualization
thermo 10
dump mydmp all atom 10 dump.min.lammpstrj
# 5) Run
minimize 1.0e-4 1.0e-6 1000 10000
write_data minimized_coordinate.data
最后一行表示lammps将最终结果输出至minimized_coordinate.data文件。这个文件保存了能量最小化的状态，可以在以后重新模拟时使用。
更进一步，如果只想保留圆柱内类型2的原子和圆柱外类型1的原子，我们可以通过如下命令删除不符合条件的原子
region mycylin cylinder z 0 0 15 INF INF side in
region mycylou cylinder z 0 0 15 INF INF side out
group mytype1 type 1
group mytype2 type 2
group incyl region mycylin
group oucyl region mycylou
group type1in intersect mytype1 incyl
group type2ou intersect mytype2 oucyl
delete_atoms group type1in
delete_atoms group type2ou
运行时，可以发现lammps输出了相关的信息
为了输出相关信息，可以采用如下命令
variable Ntype1in equal count(mytype1,mycylin)
variable Ntype1ou equal count(mytype1,mycylou)
variable Ntype2in equal count(mytype2,mycylin)
variable Ntype2ou equal count(mytype2,mycylou)
fix myat1 all ave/time 1000 10 10000 v_Ntype1in v_Ntype1ou file population1vstime.dat
fix myat2 all ave/time 1000 10 10000 v_Ntype2in v_Ntype2ou file population2vstime.dat
compute coor12 type1 coord/atom cutoff 2.0 group type2
compute sumcoor12 all reduce ave c_coor12
fix myat3 all ave/time 1000 10 10000 c_sumcoor12 file coordinationnumber12.dat
这里定义了相关的变量，然后进行计算。fix ave/times求相关变量的平均，每10000步计算一次，10个步长的平均。后面是计算类型1和类型2原子间的配位数。
# 5) Run
velocity all create 1.0 4928459 mom yes rot yes dist gaussian
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917 zero yes
fix myefn all enforce2d
timestep 0.005
run 3000000
write_data data.mixed.lammps
这里，velocity create为每个原子施加初始速度。初始速度选取满足温度=1，随机产生，分布为高斯分布；mom yes 和rot yes表示该系统没有线动量和角动量（好奇怪的写法）。郎之万热浴中zero yes表示总的随机力为0。
