LAMMPS案例

1.4. LAMMPS案例#

LAMMPS入门案例——水中的PEG

https://2301032.xyz/peginwater-6.webp

教程来自法国格勒诺布尔物理研究所 (LIPhy), 由Simon Gravelle 开发和维护： https://lammpstutorials.github.io

教程用到的文件打包在此: 百度网盘

环境准备: lammps软件的安装参见 安装GROMACS 节

1.4.1. 准备水盒子#

  • 首先创建并平衡一个矩形水盒子

  • 使用 SPC/Fw 水模型：这是刚性 SPC（简单点电荷）模型的变体。

1.4.1.1. 输入文件准备#

创建一个名为 water.lmp 的文件，并将以下行复制到其中：

units real
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
pair_style lj/cut/coul/long 10
kspace_style ewald 1e-5
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 angle yes

region box block -30 30 -15 15 -15 15
create_box 8 box bond/types 7 angle/types 8 dihedral/types 4 extra/bond/per/atom 3 extra/angle/per/atom 6 extra/dihedral/per/atom 10 extra/special/per/atom 14

include parameters.inc

molecule h2omol water.mol
create_atoms 0 random 700 87910 NULL mol h2omol 454756 overlap 1.0 maxtry 50

group H2O type OW HW
minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

fix mynpt all npt temp 300 300 100 iso 1 1 1000

dump viz all image 250 myimage-*.ppm type type &
shiny 0.1 box no 0.01 view 0 90 zoom 3 size 1000 600
dump_modify viz backcolor white &
acolor OW red acolor HW white &
adiam OW 3 adiam HW 1.5

variable myvol equal vol
variable myoxy equal count(H2O)/3
variable NA equal 6.022e23
variable Atom equal 1e-10
variable M equal 0.018
variable rho equal ${myoxy}*${M}/(v_myvol*${NA}*${Atom}^3)
thermo 500
thermo_style custom step temp etotal v_myvol v_rho

timestep 1.0
run 15000

write_restart water.restart

  • unit命令设定模拟使用的单位制

    • real单位制：能量为 kcal/mol, 距离为 Å（埃）, 质量 amu, 时间 fs, 电荷 e

  • atom_style full命令定义了原子数据的存储格式。full格式支持的原子属性最为全面(当模拟的体系包含复杂分子结构时，通常会选用这种格式)，涵盖：

    • 原子 ID 和类型

    • 分子 ID（用于识别原子所属分子）

    • 电荷

    • 键、角、二面角等拓扑结构信息

  • bond_style 、 angle_style 和 dihedral_style 命令定义模拟中使用的键、角和二面角的势函数形式

  • pair_style命令设置的是原子间的非键相互作用势函数，具体包含两部分：

    • Lennard-Jones 势：用于描述短程(lj/cut)的排斥和吸引作用，其截断距离为 10Å。

    • 库仑势：采用长程作用处理方式（coul/long），同样设置截断距离为 10Å。长程库仑作用会在后续通过kspace_style进一步处理。

  • kspace_style命令设置的是长程静电相互作用的计算方法。

    • ewald代表使用 Ewald 求和法，这是一种高精度计算周期性体系中静电相互作用的方法。

    • 1e-5则表示计算的精度要求（相对误差）。

  • special_bonds命令用于处理分子内特殊原子对之间的非键相互作用，具体规则如下：

    • Lennard-Jones 势：

      • 1-2 键相连原子对：相互作用被完全屏蔽（系数为 0.0）。

      • 1-3 键角相连原子对：相互作用也被完全屏蔽（系数为 0.0）。

      • 1-4 二面角相连原子对：相互作用减弱为正常值的 0.5 倍。

    • 库仑势：

      • 1-2 和 1-3 原子对：相互作用被屏蔽（系数为 0.0）。

      • 1-4 原子对：相互作用保持完整（系数为 1.0）。

    • angle yes：表明在计算中要考虑键角对原子位置的约束作用。

  • region和create_box创建一个尺寸为 \(6×3×3 nm^3\) 的 3D 模拟盒子, 并设置拓扑结构参数（如键、角、二面角类型的数量）及原子属性上限。

    • LAMMPS 需要预先知道拓扑结构的复杂度，以分配足够的内存空间。

    • 这两行命令共同定义了模拟的物理空间（长方体盒子）和化学空间（原子类型、拓扑结构）。后续步骤通常会在这个盒子中放置原子（如create_atoms或read_data），并设置初始条件（如温度、速度）。

  • 使用 parameters.inc 文件来设置所有参数（质量、相互作用能、键平衡距离等）

  • molecule和create_atoms导入一个名为 water.mol 的分子模板，然后随机创建 700 个分子。

  • 将水分子的 OW 和 HW 类型的原子组织在一个名为 H2O 的组中，并执行一个小的能量最小化。

  • NPT平衡

  • 输出为图像

  • 每 500 步提取体积和密度

  • 运行 15 ps 的模拟, 将最终状态保存在名为 water.restart 的二进制文件中。

1.4.1.2. 运行命令#

使用lamms运行上述指令 lmp -in water.lmp, 你将得到 water.restart二进制文件(包含模拟的完整状态), log.lammps日志文件(终端的输出信息)以及每隔一定时间输出的图像(如下图).

1.4.2. 将 PEG 溶于水中#

下面将PEG聚合物 \(C_{16}H_{34}O_{9}\) 添加到上述的水盒子中.

创建一个名为 merge.lmp 的文件, 包含以下内容:

read_restart water.restart
kspace_style pppm 1e-5

molecule pegmol peg.mol
create_atoms 0 single 0 0 0 mol pegmol 454756

group PEG type C CPos H HC OAlc OE

delete_atoms overlap 2.0 H2O PEG mol yes

fix mynpt all npt temp 300 300 100 x 1 1 1000

fix myrct PEG recenter 0 0 0 shift all

dump viz all image 250 myimage-*.ppm type type size 1100 600 box no 0.1 shiny 0.1 view 0 90 zoom 3.3 fsaa yes bond atom 0.8
dump_modify viz backcolor white acolor OW red adiam OW 0.2 acolor OE darkred adiam OE 2.6 acolor HC white adiam HC 1.4 acolor H white adiam H 1.4 acolor CPos gray adiam CPos 2.8 acolor HW white adiam HW 0.2 acolor C gray  adiam C 2.8 acolor OAlc darkred adiam OAlc 2.6
thermo 500

timestep 1.0
run 10000

write_restart merge.restart

然后运行 lmp -in merge.lmp

1.4.3. 拉伸PEG分子#

接下来, 施加一个恒定的力到 PEG 分子的两端，直到它伸展.

创建一个名为 pull.lmp 的文件，其中包含以下内容：

kspace_style pppm 1e-5
read_restart merge.restart

group ends type OAlc
variable xcm equal xcm(ends,x)
variable oxies atom type==label2type(atom,OAlc)
variable end1 atom v_oxies*(x>v_xcm)
variable end2 atom v_oxies*(x<v_xcm)
group topull1 variable end1
group topull2 variable end2

dump viz all image 250 myimage-*.ppm type type shiny 0.1 box no 0.01 view 0 90 zoom 3.3 fsaa yes bond atom 0.8 size 1100 600
dump_modify viz backcolor white acolor OW red acolor HW white acolor OE darkred acolor OAlc darkred acolor C gray acolor CPos gray acolor H white acolor HC white adiam OW 0.2 adiam HW 0.2 adiam C 2.8 adiam CPos 2.8 adiam OAlc 2.6 adiam H 1.4 adiam HC 1.4 adiam OE 2.6

timestep 1.0
fix mynvt all nvt temp 300 300 100
fix myrct PEG recenter 0 0 0 shift all

compute rgyr PEG gyration
compute prop PEG property/local dtype
compute dphi PEG dihedral/local phi

thermo_style custom step temp etotal c_rgyr
thermo 250
dump mydmp all local 100 pull.dat index c_dphi c_prop

run 15000

fix myaf1 topull1 addforce 10 0 0
fix myaf2 topull2 addforce -10 0 0
run 15000

然后运行lmp -in pull.lmp

1.4.4. 结果分析#

https://2301032.xyz/peginwater-4.jpg

图. 1.4.1 PEG 分子的回转半径变化及二面角直方图#

https://2301032.xyz/peginwater-5.jpg

图. 1.4.2 水分子的氧原子和 PEG 分子的氧原子之间的径向分布图#

目录