lammps模拟气体扩散

1.lammps简介

1.什么是lammps

• 分子动力学软件包
• 全称：Large-scale Atomic/Molecular Massively Parallel Simulation
• 官方网址：https://www.lammps.org/
• 开发单位：美国圣地亚国家实验室
• 主要开发人员：Steve Plimpton
• 性质：开源

2.lammps的特点

• 可以串行或并行计算，支持GPU加速
• 开源，C++语言编写，高移植性，高性能
• 可以方便的扩展，增加新功能
• 自定义变量和方程
• 一个输入脚本实现一个或多个模拟任务

3.lammps的功能

• 模拟对象：原子，金属，聚合物，生物分子，粒状和粗粒化体系
• 计算体系：小至几个粒子，大到上百万甚至上亿个粒子
• 运行平台：单处理器/多处理器计算机，Linux(推荐)/Windows
• 大多数使用者将其用作独立代码，也可将其用作库：
      1.直接耦合到另一个代码，如多物理场或多尺度模拟
      2.在其他代码的工作流中使用lammps，在机器学习应用中越来越普遍
      3.从Python脚本启动和使用lammps，使用mpi4py模块并行运行Python

2.LJ约化单位制

1.计算模拟的单位制

• 所有的计算模拟都需要首先定义单位制，它决定了所有输入脚本、数据文件和所有输出到屏幕、日志文件以及dump文件中物理量的单位。该命令通常用在输入脚本最开始的位置

  style = lj or real or metal or si or cgs or electron
  

• 物理量	si单位制 (国际单位制)	metal单位制
  Mass	公斤 kilograms	grams/mole
  Distance	米 meters	Angstroms
  Time	秒 seconds	picoseconds
  Energy	焦耳 Joules	eV
  Velocity	米/秒 meters/second	Angstroms/picosecond
  Force	牛顿 Newtons	eV/Angstrom
  Torque	牛顿·米 Newton-meters	eV
  Temperature	开尔文 Kelvin	Kelvin
  Pressure	帕斯卡 Pascals	bars

2.约化Lennard-Jones势

• 原始公式：$$V(r)=4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$$
• 约化公式：$$V^{\prime }(r)=V(r)/\epsilon r_i^{\prime }=r_i/\sigma$$ $$V^{\prime }(r)=4[\frac{1}{r^{\prime 12}}-\frac{1}{r^{\prime 6}}]$$

3.其他约化标度

• 质量单位 $$m_i^{\prime }=m_i/m$$
• 长度单位 $$r_i^{\prime }=r_i/\sigma$$
• 能量单位 $$V_i^{\prime }=V_i/\epsilon$$
• 时间单位 $$\frac{1}{2}m(\frac{\sigma }{\tau }{)}^2=\epsilon \tau =\sigma \sqrt{m/\epsilon }$$ $$t_i^{\prime }=t_i/\tau$$

4.lammps中应用

• # a = 1 埃米
  units   metal
  lattice   hex   1.0
  
  # 产生等效密度为1.0的晶格
  units    lj
  lattice   hex   1.0
  

• 实际原子密度(单位面积内原子数)与等效原子密度的换算：$$\rho =\frac{N}{S}=\frac{N}{L^2}{\rho }^{\prime }=\frac{N}{(L/\sigma {)}^2}$$ $${\rho }^{\prime }=\rho {\sigma }^2=(\frac{N}{S}){\sigma }^2$$ $$1.0=(\frac{1}{\frac{\sqrt{3}}{2}{a}^2}){\sigma }^2$$
• 原子间距为：$$a=\sqrt{\frac{2}{\sqrt{3}}}\approx 1.07\sigma$$

3.气体扩散模拟

1.模拟思路

• 将 Ar 原子在二维的模拟盒子里一分为二，在左半部分放入 Ar 原子，右半部分为空。
• 原子初始速度随机分布，但保持体系总动量为零。从时间 t=0 开始体系演化 10000 个时间步

2.模拟脚本

# 初始化
units   lj      # 设置lj单位制
dimension   2           # 设置模拟维度
boundary   p p p          # 设置边界条件
atom_style   atomic         # 定义原子的类型，它会决定原子包括哪些属性
# 原子模型创建
lattice   hex 1.0          # 定义晶格类型
region   mybox block 0 20 0 10 -0.1 0.1       # 定义空间几何区域
create_box   1 mybox          # 创建模拟盒子
region   2 block 0 10 0 10 -0.1 0.1
create_atoms   1 region 2           # 在定义的晶格点阵上创建原子
mass   1 1.0          # 设置原子质量
velocity   all create 0.5 78598           # 设置原子速度
# 参数设置
pair_style   lj/cut 2.5           # 设置原子间相互作用势
pair_coeff   1 1 1.0 1.0 2.5        # 设置对势的参数
neighbor   0.3 bin             # 定义近邻列表
neigh_modify   every 20 delay 0 check no           # 设置近邻列表算法的参数
fix   1 all nvt temp 0.5 0.5 0.01          # 定义分子动力学时间步中的操作
fix   2 all enforce2d
# 动力学模拟
dump   1 all custom 100 toEquil.lammpstrj id type x y z vx vy vz      # 输出分子动力学模拟的信息
thermo 500           # 输出热力学信息的间隔
run 10000          # 运行MD

3.模拟涉及命令

• 初始化：在原子被定义前，设置一些基本的环境参数
      units，dimension，boundary，atom_style
• 原子定义：创建或读入模型
      lattice，region，create_box，create_atoms，read_data或read_restart
• 参数设置：指定一系列的计算参数：力场、模拟参数、输出控制选项等
      力场：pair_style，pair_coeff
      模拟：neighbor，neigh_modify，group，timestep，min_style
      边界条件、积分：fix
      计算：compute、compute_modify、variable
      输出控制：thermo，dump，restart
• 启动模拟：minimize，run

4.多组分气体模拟

# 初始化
units   lj          # 使用LJ单位制
dimension   2          # 使用二维空间
boundary   p p p          # 使用周期性边界条件
atom_style   atomic          # 使用普通原子
variable   t equal 0.5          #定义系统温度
# 创建原子模型
lattice   sq 1.0          # 使用简单正方晶格，晶格常数为1.0
region   box block 0 100 0 100 -0.5 0.5          # 定义盒子的范围
create_box   2 box          # 在盒子中创建两种类型的原子
create_atoms   1 random 2500 23546 box          # 随机创建2500个类型为1的原子
create_atoms   2 random 2500 56145 box          # 随机创建2500个类型为2的原子
mass   1 1.0          # 设置原子1的质量
mass   2 1.0          # 设置原子2的质量
# 参数设置
pair_style   hybrid lj/cut 2.5 soft 5.0          # 使用hybrid混合两种势
pair_coeff   1 1 lj/cut 1.0 1.0 2.5          # 类型为1的原子之间为LJ势
pair_coeff   2 2 lj/cut 1.0 1.0 2.5          # 类型为2的原子之间为LJ势
pair_coeff   1 2 soft 5.0          # 类型为1和2的原子之间为soft势
# 动力学模拟  
compute   eng all pe/atom          # 计算每个原子的势能
compute   eatoms all reduce sum c_eng          # 计算所有原子的总势能
thermo_style   custom step temp epair etotal press c_eatoms          # 定义输出的热力学量
thermo   1000          # 梅1000步输出一次热力学量
dump   id all atom 100 dump.lammpstrj          # 每100步将所有原子的状态输出到dump.lammpstrj文件中
minimize   1e-4 1e-6 1000 10000          # 防止随机生成的原子重叠能量过大
velocity   all create $t 23561          # 赋予所有原子初始速度
fix   nvt all nvt temp $t $t 0.01          # 使用nvt系综，温度为t，时间阻尼为0.01
run   50000          # 运行50000步