【分子模拟】使用lammps模拟水在多孔二氧化硅中的吸附

水在二氧化硅中的吸附**

简介

本教程的目的是模拟水在多孔二氧化硅中的吸附。首先构建一个包含多孔二氧化硅和水的体系，然后使用 LAMMPS 进行分子动力学模拟。最后，分析模拟结果，以了解水分子在二氧化硅孔隙中的行为。

1. 准备结构文件

在开始模拟之前，需要准备两个结构文件：一个用于多孔二氧化硅，另一个用于水。

1.1. 多孔二氧化硅结构

首先，需要创建一个多孔二氧化硅结构。有多种方法可以做到这一点，这里使用 VMD 及其内置的无机物构建器。

1. 打开 VMD。
2. 转到 Extensions -> Modeling -> Inorganic Builder。
3. 在 Material 下拉菜单中选择 Cristobalite (方石英)。
4. 设置晶胞参数，例如 a=5.0, b=5.0, c=5.0。
5. 点击 Create Unit Cell (创建晶胞)。
6. 现在，为了创建一个孔隙，可以删除一些原子。选择一些原子并删除它们。也可以使用 VMD 的 Tk Console (Tk 控制台) 来更精确地完成此操作。例如，要删除一个球体内的所有原子，可以使用以下命令：

   set sel [atomselect top "within 2.0 (x y z)"]
   $sel delete

   这将删除以坐标 (x, y, z) 为中心、半径为 2.0 Å 的球体内的所有原子。
7. 对结构满意后，将其保存为 data.sio2 文件。转到 File -> Save Data...，然后选择 LAMMPS Data 格式。
   注意: 也可以使用其他软件（如 Avogadro 或 Materials Studio）来创建多孔二氧化硅结构。确保最终导出为 LAMMPS 可以读取的格式，例如 .data 文件。

1.2. 水盒子

接下来，需要创建一个水盒子。这里将使用 Packmol 软件来构建它。

1. 创建一个名为 water.inp 的输入文件，内容如下：

   tolerance 2.0
   filetype xyz
   output water.xyz
   structure spc216.pdb
     number 1000
     inside box 0. 0. 0. 40. 40. 40.
   end structure

   • tolerance 2.0: 设置分子间的最小距离。
   • filetype xyz: 指定输出文件格式。
   • output water.xyz: 指定输出文件名。
   • structure spc216.pdb: 指定水分子模板文件（一个包含单个 SPC/E 水分子的 PDB 文件）。
   • number 1000: 指定要放置的水分子数量。
   • inside box ...: 指定放置水分子的盒子尺寸。
2. 运行 Packmol：

   packmol < water.inp

   这将生成一个名为 water.xyz 的文件。
3. 现在，需要将 water.xyz 文件转换为 LAMMPS 的 .data 文件。可以使用 VMD 来完成此操作。
   • 在 VMD 中打开 water.xyz 文件。
   • 打开 Tk Console (Extensions -> Tk Console)。
   • 输入以下命令：

     topo readlammpsdata water.data
     topo writelammpsdata water.data

     这将创建一个名为 water.data 的 LAMMPS 数据文件。

1.3. 合并结构

现在有了 data.sio2 和 water.data 两个文件，需要将它们合并成一个单一的 LAMMPS 数据文件。这可以通过手动编辑或使用脚本来完成。关键点在于：

• 将两个文件的原子坐标部分合并。
• 正确更新原子总数、原子类型数、键数等头部信息。
• 确保原子类型编号不冲突（例如，二氧化硅中的 Si 原子类型为 1，O 原子类型为 2，那么水中的 O 和 H 原子类型应从 3 和 4 开始）。
  这个过程可能有些繁琐。一个更简单的方法是使用 LAMMPS 本身来组合它们，但为了本教程的简洁性，假设你已经手动或通过自定义脚本成功创建了一个名为 system.data 的最终数据文件。

2. LAMMPS 输入脚本

现在，可以编写 LAMMPS 输入脚本了。脚本分成几个部分。

2.1. 初始化

# --- 初始化 ---
units           real
atom_style      full
boundary        p p p

• units real: 使用 real 单位制（Å, kcal/mol, fs 等）。
• atom_style full: 使用 full 原子样式，因为它支持分子、键、角和电荷。
• boundary p p p: 在所有三个方向上使用周期性边界条件。

2.2. 读取数据文件

# --- 读取数据 ---
read_data       system.data

这会读取之前准备的 system.data 文件。

2.3. 力场设置

使用 CLAYFF 力场来描述二氧化硅，使用 SPC/E 力场来描述水。

# --- 力场设置 ---
pair_style      lj/cut/coul/long 12.0
pair_modify     mix arithmetic
kspace_style    pppm 1e-4
# 二氧化硅 (CLAYFF)
pair_coeff      1 1 0.0      0.0      # Si-Si
pair_coeff      1 2 103033.252 3.302   # Si-O (二氧化硅)
pair_coeff      2 2 230400.0  3.166    # O-O (二氧化硅)
# 水 (SPC/E)
pair_coeff      3 3 0.1553   3.166    # O-O (水)
pair_coeff      4 4 0.0      0.0      # H-H
pair_coeff      3 4 0.0      0.0      # O-H (水)
# 二氧化硅-水相互作用
pair_coeff      1 3 103033.252 3.302   # Si-O (水)
pair_coeff      1 4 0.0      0.0      # Si-H
pair_coeff      2 3 230400.0  3.166    # O (二氧化硅)-O (水)
pair_coeff      2 4 0.0      0.0      # O (二氧化硅)-H
# 键和角参数 (水)
bond_style      harmonic
angle_style     harmonic
bond_coeff      1 1000.0 1.0
angle_coeff     1 100.0 109.47

• pair_style lj/cut/coul/long 12.0: 使用 Lennard-Jones (LJ) 和长程库仑相互作用，截断半径为 12.0 Å。
• pair_modify mix arithmetic: 使用算术平均规则来混合不同原子类型的 LJ 参数。
• kspace_style pppm 1e-4: 使用 PPPM 算法计算长程库仑相互作用，精度为 1e-4。
• pair_coeff: 为不同原子类型对定义 LJ 参数 (ε, σ)。注意: 你需要根据你的 system.data 文件中的原子类型编号来调整这些系数。这里假设 1=Si, 2=O_sio2, 3=O_water, 4=H_water。
• bond_style 和 angle_style: 定义键和角的势能函数为谐振子势。
• bond_coeff 和 angle_coeff: 为水分子定义键和角的力常数与平衡值。

2.4. 能量最小化

在开始动力学模拟之前，需要最小化体系的能量，以消除任何不合理的原子重叠。

# --- 能量最小化 ---
min_style       cg
minimize        1.0e-4 1.0e-6 1000 10000

• min_style cg: 使用共轭梯度算法进行最小化。
• minimize: 设置最小化收敛标准。

2.5. 平衡模拟

现在，将对体系进行平衡，使其达到目标温度和压力。

# --- 平衡 ---
reset_timestep  0
timestep        1.0
# 速度初始化
velocity        all create 300.0 12345
# NVT 平衡
fix             1 all nvt temp 300.0 300.0 100.0
thermo          100
thermo_style    custom step temp pe ke etotal press vol
run             10000
unfix           1
# NPT 平衡
fix             1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run             20000
unfix           1

• reset_timestep 0: 将时间步重置为 0。
• timestep 1.0: 设置时间步长为 1 fs。
• velocity all create 300.0 12345: 根据随机数种子 12345，为所有原子赋予 300 K 的初始速度。
• fix 1 all nvt ...: 在 NVT（正则）系综下运行模拟，使用 Nosé-Hoover 恒温器将温度维持在 300 K。
• thermo 100: 每 100 步输出一次热力学信息。
• thermo_style: 自定义热力学输出的内容。
• run 10000: 运行 10000 步（10 ps）的 NVT 模拟。
• fix 1 all npt ...: 在 NPT（等温等压）系综下运行模拟，使用 Nosé-Hoover 恒温器和恒压器，将温度维持在 300 K，压力维持在 1 atm。
• run 20000: 运行 20000 步（20 ps）的 NPT 模拟。

2.6. 生产模拟

体系平衡后，运行一个更长的生产模拟来收集数据。

# --- 生产模拟 ---
reset_timestep  0
fix             1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
# 输出轨迹
dump            1 all custom 1000 dump.lammpstrj id type x y z
# 运行模拟
run             100000
# 清理
undump          1
unfix           1
print           "Simulation finished"

• dump 1 all custom 1000 dump.lammpstrj ...: 每 1000 步将所有原子的坐标、类型和 ID 输出到一个名为 dump.lammpstrj 的轨迹文件中。
• run 100000: 运行 100000 步（100 ps）的生产模拟。
• undump 和 unfix: 在模拟结束时移除 dump 和 fix 命令，这是一个好习惯。

3. 运行模拟

将上述所有部分保存到一个名为 in.adsorption 的文件中。然后，在终端中运行 LAMMPS：

lmp -in in.adsorption

模拟完成后，你将得到一个名为 dump.lammpstrj 的轨迹文件和一个名为 log.lammps 的日志文件。

4. 分析结果

现在，可以分析模拟结果了。我们将使用 VMD 和一些自定义脚本来进行分析。

4.1. 可视化

使用 VMD 打开 dump.lammpstrj 文件，可视化水分子在二氧化硅孔隙中的运动和吸附行为。

4.2. 密度分布图

计算水分子沿 z 轴的密度分布，以观察水在孔隙中的富集情况。

1. 在 VMD 中加载轨迹文件。
2. 打开 Tk Console。
3. 运行以下 Tcl 脚本：

   set sel [atomselect top "name O and type 3"] # 选择水中的氧原子
   set nf [molinfo top get numframes]
   set outfile [open "density_z.dat" w]
   puts $outfile "# z density"
   for {set i 0} {$i < $nf} {incr i} {
       $sel frame $i
       set z [$sel get z]
       set z_bins [histogram $z -bins 100 -min 0 -max 50]
       puts $outfile "$z_bins"
   }
   close $outfile

   这将生成一个名为 density_z.dat 的文件，其中包含沿 z 轴的密度数据。可以使用 Gnuplot 或 Python/Matplotlib 对其进行绘图。

4.3. 径向分布函数 (RDF)

RDF 可以告诉我们水分子在二氧化硅表面的结构排布信息。

1. 在 VMD Tk Console 中运行：

   set sel1 [atomselect top "name O and type 2"] # 选择二氧化硅中的氧原子
   set sel2 [atomselect top "name O and type 3"] # 选择水中的氧原子
   set r [measure rdf $sel1 $sel2 -rmax 10.0 -bins 100 -useframe all -update 100]
   set outfile [open "rdf_oo.dat" w]
   foreach {r g_r} $r {
       puts $outfile "$r $g_r"
   }
   close $outfile

   这将计算并输出二氧化硅氧原子与水氧原子之间的 RDF 到 rdf_oo.dat 文件中。

结论

本教程向你展示了如何模拟水在多孔二氧化硅中的吸附。涵盖了从准备结构、编写 LAMMPS 输入脚本到分析模拟结果的整个过程。可以通过修改体系（例如，改变孔隙大小、使用不同的表面或不同的液体）来扩展这个教程，以研究其他吸附体系。

---