- 科研
- 科学计算
- 高性能计算
【免费下载链接】lammps
Public development project of the LAMMPS MD software package
导读
本文基于 LAMMPS 官方手册中的Howto bulk2slab指南（源文档位于 doc/src/Howto_bulk2slab.rst），讲解如何把已经运行一段时间、达到平衡的体相（bulk）体系安全地转换为带有真空层与自由表面的slab 体系。文中以 LAMMPS 自带的 Rhodopsin 基准算例（bench/in.rhodo）为演示对象，给出可直接复制的完整输入脚本，并深入解析reset_atoms、write_dump、change_box、read_dump、kspace_modify slab auto等核心命令的底层原理。读完本文，你将掌握一套通用的 bulk-to-slab 转换工作流，并理解 image flag（镜像标志）在盒子变形中的关键作用。
问题背景：为什么不能直接改盒子尺寸？
体相体系经过长时间 MD 模拟后，原子会发生扩散，很多原子带有非零的 image flag。这些标志记录了原子相对盒子原点经历了多少次周期平移，用于还原其 "unwrapped"（解卷绕）坐标：真实坐标 = 局部坐标 + image×盒长。
如果直接使用 change_box 命令 改变盒子的 x、y、z 尺寸，或直接修改 data 文件中的盒子边界，那么这些带有非零 image flag 的原子会发生非物理的位移——因为盒子长度变化后，image flag 所代表的平移距离也随之改变，使它们相对零 image flag 的邻居原子发生错位。
而如果强制把所有 image flag 清零，同样会出问题：对于跨越周期边界的分子（例如一条键的其中一个原子在体系顶部、另一个在底部），强行清零会把它们"拆散"到盒子的两端，破坏分子完整性。
因此，正确做法是：先"解卷绕"原子坐标并保存，再安全地扩大盒子，最后把解卷绕坐标读回。这正是本文工作流的思路。
适用范围与限制（原文档声明）
必须注意，该工作流适用于大多数体相体系，但并非全部：
- 无法转换的体系：例如分子通过周期边界与自身成键的"无限长聚合物"（infinite polymers），此类体系原则上无法切分成 slab；
- 真空层厚度的选择：需要添加多少真空取决于体系在切分处的分子长度（本文示例在含短分子（水）的区域切分）；
- 可能需要先重新居中：某些体系需要先用 displace_atoms 命令 将体系重新居中到盒子中心再操作；
- 热化与平衡时间：强热化与平衡所需时间取决于具体体系及其热力学条件，不能一概而论。
完整工作流：基于 Rhodopsin 基准算例
原文档给出的流程以 Speed_bench 中的 Rhodopsin 基准输入 bench/in.rhodo 为演示对象。该算例是"溶剂化脂质双层中的视紫红质蛋白"体系，采用 CHARMM 力场（lj/charmm/coul/long对势 +pppm长程库仑求解），原始输入片段如下：
# Rhodopsin model units real neigh_modify delay 5 every 1 atom_style full bond_style harmonic angle_style charmm dihedral_style charmm improper_style harmonic pair_style lj/charmm/coul/long 8.0 10.0 pair_modify mix arithmetic kspace_style pppm 1e-4 read_data data.rhodo fix 1 all shake 0.0001 5 0 m 1.0 a 232 fix 2 all npt temp 300.0 300.0 100.0 & z 0.0 0.0 1000.0 mtk no pchain 0 tchain 1 special_bonds charmm thermo 50 thermo_style multi timestep 2.0 run 100整个转换流程的实现方式是：复制一份in.rhodo输入文件，保留从开头到read_data命令为止的所有内容不变，然后在read_data之后插入下述命令，向 z 方向的上、下两侧添加真空：
variable delta index 10.0 reset_atoms image all write_dump all custom rhodo-unwrap.lammpstrj id xu yu zu change_box all z final $(zlo-2.0*v_delta) $(zhi+2.0*v_delta) & boundary p p f read_dump rhodo-unwrap.lammpstrj 0 x y z box no replace yes kspace_modify slab auto下面逐条拆解这六条命令的作用与底层原理。
1. 用variable delta控制真空层厚度
variable delta index 10.0variable命令 定义了一个值为 10.0 的标量变量。它表示向 z 方向每一侧添加的真空距离：上下各加2.0*delta（即各 20 Å），因此体系总共新增40 Å的真空。读者可根据自己的体系自由调整该值——真空越多，slab 之间通过长程库仑的相互作用越弱。
2. 用reset_atoms image规范化镜像标志
reset_atoms image allreset_atoms image命令 把所有原子的 image flag 重置为 0 或 ±1，使每个原子距离盒子中心最近，同时保持成键原子之间的正确相对距离。这一步是后续坐标读回的基石：它清除了长时模拟累积的"多圈"平移信息，只保留最小周期像。
3. 用write_dump保存解卷绕坐标
write_dump all custom rhodo-unwrap.lammpstrj id xu yu zuwrite_dump命令 把体系中每个原子的unwrapped 坐标（xu yu zu）写入 dump 文件。解卷绕坐标不受盒子尺寸影响，因此能如实记录原子在真实空间中的位置。之后扩大盒子时，原本在盒子外的原子（解卷绕后位于盒子外）将被"收进"盒内，越出周期边界的原子会被卷绕（wrap）回盒子，其在对应方向的 image flag 也随之恢复。
4. 用change_box扩大盒子并切换边界条件
change_box all z final $(zlo-2.0*v_delta) $(zhi+2.0*v_delta) & boundary p p fchange_box命令 将 z 方向的低、高边界分别向外扩展2.0*delta，同时通过 boundary 命令 的关键字把边界条件改为p p f（x、y 方向保持周期性，z 方向改为固定非周期边界）。
如背景所述，如果此刻体系里还有非零 image flag 的原子，改边界会迫使 z 方向 image flag 归零，从而产生不希望的位移。所以这一行必须放在reset_atoms image与后续坐标读回之间执行。
5. 用read_dump读回解卷绕坐标
read_dump rhodo-unwrap.lammpstrj 0 x y z box no replace yesread_dump命令 把刚才保存的解卷绕坐标读回并替换当前体系中部分不正确的坐标。两个关键子参数：
box no：忽略 dump 文件中记录的盒子尺寸——我们必须保留刚被change_box扩大后的盒子，而不是 dump 里的旧盒子；replace yes：用读入坐标替换现有原子坐标。
此时，原本因 image flag 清零而可能错位的原子，其位置被解卷绕坐标精确恢复，而 z 方向边界已是非周期f，新盒子的上下两侧自然形成真空与自由表面。
6. 用kspace_modify slab auto开启 slab 长程校正
kspace_modify slab autoz 方向改为非周期边界后，体系沿 z 方向实际上是一个"2D 周期"结构。若继续使用 PPPM（或 Ewald）这类 3D 长程库仑求解器，就必须开启 slab 校正（基于 Yeh & Berkowitz 的 EW3DC 方法），否则周期像之间会产生虚假的库仑相互作用。
slab auto（推荐）：让求解器根据当前长程精度目标和盒子几何自动选择近似最优的虚拟真空体积因子（volfactor），无需用户手动调参。该选项由 PPPM 的auto_slab_volfactor()实现（见 src/KSPACE/pppm.cpp）：它从归一化力容差accuracy / two_charge_force出发，取横向尺寸项max(xprd,yprd)*log(1/tol)/(2π)与倒空间衰减项sqrt(log(1/tol))/alpha中的较大者，再结合盒子 z 长度估算slab_volfactor = (zprd + max(lateral, reciprocal)) / zprd；若该值与上一轮估计差异超过阈值，会迭代最多 5 次直至收敛（src/KSPACE/pppm.cpp）。Ewald 求解器也提供同样的自动选项（src/KSPACE/ewald.cpp）。- 手动指定（备用）：如果希望人工覆盖，仍可使用固定值，例如
kspace_modify slab 3.0。kspace_modify中slab关键字还可接受nozforce（关闭 z 方向 kspace 力）等变体，但需要注意其支持范围与限制（详见 kspace_modify 文档）。值得注意的是，slab auto需要正的 gewald 值，且要求归一化力容差落在 (0,1) 区间，否则求解器会直接报错（src/KSPACE/pppm.cpp）。
从源码角度看，slab 校正在倒空间求和时把 z 方向长度按slab_volfactor拉伸为zprd_slab = zprd * slab_volfactor（见 src/KSPACE/ewald.cpp 与 src/KSPACE/pppm.cpp），并用slabcorr()加入偶极矩校正项来消除相邻 slab 间的相互作用（src/KSPACE/ewald.cpp）。
后续平衡：更换系综、加软墙、恢复恒温器
由于体系现在是开放体系，z 方向不再需要调整晶胞尺寸，因此要把原输入中的fix npt替换为fix nvt：
fix 2 all nvt temp 300.0 300.0 10.0fix nvt命令 使用 Nose–Hoover 恒温器控制温度。两个原因促成了这次替换：
- z 方向盒子不再变化，
npt中的 z 方向控压项失去意义； - 在真空插入处（体相水区域）切开体系会暴露高势能的表面原子。为尽快清除由此产生的多余动能，恒温器时间常数从原来的 100.0 减小到10.0，使控温更"强力"。
添加软谐波墙防止原子逃逸
高势能的表面原子有可能被"弹射"出 slab。为抑制这一现象，需要先求出体系原子在 z 方向的极值，再在上下两侧各放置一个软谐波墙：
compute zmin all reduce min z compute zmax all reduce max z thermo_style custom zlo c_zmin zhi c_zmax run 0 post no fix 3 all wall/harmonic zhi $(c_zmax+v_delta) 10.0 0.0 ${delta} & zlo $(c_zmin-v_delta) 10.0 0.0 ${delta}逐行说明：
- 两条 compute reduce 命令 分别求所有原子 z 坐标的最小值（
c_zmin）与最大值（c_zmax）； - compute 是惰性求值的，必须被"消费"才会执行。这里用 thermo_style custom 把
zlo c_zmin zhi c_zmax加入热力学输出，紧接着用 run 0 post no 触发一次零步运行，从而拿到 min/max 值，避免在尚未求值时就用c_zmin/c_zmax计算墙位置而报错； - 两条
fix wall/harmonic分别在zmax + delta（上）与zmin - delta（下）处放置软墙：力常数为10.0 kcal/mol/Ų（经验取值，可在 fix wall/harmonic 文档 中查到完整参数含义），截止距离为 0.0，并复用之前定义的delta变量决定墙与原子极值间的距离； - 只有超出 z 方向 min/max 范围的原子才会感受到回复力并被"推回" slab，普通原子完全不受影响。
添加这些回复性软墙还能让上下自由表面更加平整，避免出现粗糙或起伏的表面。墙的效果可通过调整力常数、截止距离和墙的位置来改变。
分两阶段运行并写出新数据文件
最后，把原输入中的run 100替换为：
run 1000 post no unfix 3 fix 2 all nvt temp 300.0 300.0 100.0 run 1000 post no write_data data.rhodo-slab阶段说明：
- 第一阶段：带软墙 + 强恒温器（时间常数 10.0）跑 1000 步，让新生成的表面迅速弛豫、多余动能被快速带走；
- 用 unfix 3 移除软墙，并把恒温器时间常数恢复为 100.0（回归常规温和控温）；
- 第二阶段：再跑 1000 步，让体系在无墙、常规恒温下继续弛豫；
- 最后用 write_data 命令 把最终 slab 几何写入新的 data 文件
data.rhodo-slab，供后续模拟直接read_data使用。
注意：这里 1000 步（对应 2 ps，因为timestep 2.0fs）只是为演示方便而选取的小步数，真实体系达到平衡几乎肯定需要更多 MD 步数。
小结：bulk-to-slab 转换五步法
| 步骤 | 命令 | 目的 |
|---|---|---|
| 1 | variable delta index 10.0 | 定义真空层厚度变量（上下各 2×delta，共 4×delta 真空） |
| 2 | reset_atoms image all | 规范化 image flag 为 0/±1，保持键合原子相对距离 |
| 3 | write_dump ... id xu yu zu | 保存所有原子的解卷绕坐标 |
| 4 | change_box ... boundary p p f | 扩展 z 边界并切换为非周期固定边界 |
| 5 | read_dump ... box no replace yes | 读回解卷绕坐标，保留扩大后的盒子 |
最后辅以kspace_modify slab auto处理非周期 z 方向的长程库仑，再通过"nvt 强控温 + 软谐波墙 → 恢复恒温器 → 再平衡"两阶段运行得到稳定的 slab，并write_data落盘。这套流程的核心思想是借助 unwrapped 坐标在盒子变形前后架起一座桥梁：先用reset_atoms清理 image flag，再保存真实坐标，扩大盒子，最后用read_dump把真实坐标"搬"进新盒子。理解 image flag 与 unwrapped 坐标的语义，是掌握此工作流乃至后续自定义变体的关键。
想更深入理解底层实现，可继续阅读：
- 手册原文档：Howto_bulk2slab.rst
- 演示输入：bench/in.rhodo 与数据文件 bench/data.rhodo
- 相关命令手册：change_box、reset_atoms、read_dump、kspace_modify
- slab 自动体积因子源码：src/KSPACE/pppm.cpp、src/KSPACE/ewald.cpp
- 科研
- 科学计算
- 高性能计算
【免费下载链接】lammps
Public development project of the LAMMPS MD software package
相关推荐
LAMMPS 平台抽象层解析：`LAMMPS_NS::platform` 命名空间的函数体系与源码实现
LAMMPS 平台抽象层解析： LAMMPS_NS::platform 命名空间的函数体系与源码实现 本篇技术指南围绕 LAMMPS 开发者文档中的 Devel
科研科学计算高性能计算Spring 格式化体系源码解析：Printer 接口如何将对象转换为字符串
Spring 格式化体系源码解析：Printer 接口如何将对象转换为字符串 本文围绕 Spring 框架 org.springframework.format
文档教程知识库如何快速将TTF转换为WOFF？Node.js字体转换神器ttf2woff的完整指南 🚀
如何快速将TTF转换为WOFF？Node.js字体转换神器ttf2woff的完整指南 🚀 ttf2woff是一款专为Node.js开发的 TTF转WOFF字体
创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
