如何使用Python3从块体LAMMPS数据构建40埃半径γ-Fe₂O₃球形纳米颗粒

代码存在的核心问题 #

• 语法错误：条件判断abs(xmid-x) = 40使用单等号赋值，并非相等判断，会直接触发运行报错
• 几何逻辑错误：原代码分别限制x、y、z三个轴相对于中心的偏移范围在±40埃，最终切割出的是轴对齐正立方块，完全不符合球形结构的距离判定要求
• 文件读取错误：LAMMPS数据文件包含系统参数、模拟框尺寸等非原子信息的头部内容，直接用pandas按tab分隔全量读取，会把杂行混入原子数据，导致原子属性、坐标列完全错位
• 效率问题：逐行循环调用pandas的append方法拼接数据的写法，在新版本pandas中已废弃，处理16万原子规模的结构时运行效率极低

正确的切割逻辑 #

以大立方块的几何中心为球心，对所有原子计算三维欧氏距离，仅保留距离球心小于等于目标半径（40埃）的原子，即可得到中心对齐的球形纳米颗粒结构。球形距离计算公式为：

dist = sqrt( (x - xmid)² + (y - ymid)² + (z - zmid)² )
其中xmid、ymid、zmid为立方框的几何中心坐标。

可直接运行的修正代码 #

import numpy as np

# 配置参数
input_data_path = "data.supercell443.txt"
output_data_path = "data.spherical_np.txt"
target_radius = 40  # 目标球半径，单位：埃

# 读取原始LAMMPS文件
with open(input_data_path, 'r') as f:
    all_lines = f.readlines()

# 提取模拟框边界信息，定位原子块起始位置
box_bound = {}
atom_block_start = 0
for line_idx, line in enumerate(all_lines):
    line_strip = line.strip()
    # 解析x/y/z方向的框边界
    if "xlo xhi" in line_strip:
        xlo, xhi = map(float, line_strip.split()[:2])
        box_bound['x'] = (xlo, xhi)
    if "ylo yhi" in line_strip:
        ylo, yhi = map(float, line_strip.split()[:2])
        box_bound['y'] = (ylo, yhi)
    if "zlo zhi" in line_strip:
        zlo, zhi = map(float, line_strip.split()[:2])
        box_bound['z'] = (zlo, zhi)
    # 找到原子块的起始位置
    if line_strip == "Atoms":
        atom_block_start = line_idx + 2  # 跳过Atoms标题行与相邻空行
        break

# 计算球心坐标
xmid = (box_bound['x'][0] + box_bound['x'][1]) / 2
ymid = (box_bound['y'][0] + box_bound['y'][1]) / 2
zmid = (box_bound['z'][0] + box_bound['z'][1]) / 2

# 筛选符合半径要求的原子
kept_atoms = []
for line in all_lines[atom_block_start:]:
    line_strip = line.strip()
    if not line_strip:
        break
    line_parts = line_strip.split()
    # 原子属性：类型、x、y、z（坐标取每行最后三个值，兼容带电荷列的文件格式）
    atom_type = int(line_parts[1])
    x, y, z = map(float, line_parts[-3:])
    # 计算到球心的距离
    dist = np.sqrt( (x - xmid)**2 + (y - ymid)**2 + (z - zmid)**2 )
    if dist <= target_radius:
        kept_atoms.append( (atom_type, x, y, z) )

# 写出新的LAMMPS格式文件
with open(output_data_path, 'w') as f:
    f.write("Gamma-Fe2O3 spherical nanoparticle, radius 4nm\n\n")
    f.write(f"{len(kept_atoms)} atoms\n")
    f.write("2 atom types\n\n")
    f.write(f"{box_bound['x'][0]} {box_bound['x'][1]} xlo xhi\n")
    f.write(f"{box_bound['y'][0]} {box_bound['y'][1]} ylo yhi\n")
    f.write(f"{box_bound['z'][0]} {box_bound['z'][1]} zlo zhi\n\n")
    f.write("Atoms\n\n")
    # 重新分配连续的原子ID，避免不连续的ID触发LAMMPS报错
    for atom_id, (atype, x, y, z) in enumerate(kept_atoms, start=1):
        f.write(f"{atom_id} {atype} {x:.6f} {y:.6f} {z:.6f}\n")

使用注意事项 #

• 代码采用“取每行最后三个值为坐标”的兼容写法，可适配带电荷列、不带电荷列的不同LAMMPS原子块格式，不需要硬编码列索引
• 如果原始结构存在周期性边界导致的坐标偏移，可在计算距离前先做最小镜像校正，避免球心定位偏差
• 输出文件生成后可直接在VMD中打开，确认球形轮廓是否符合目标要求

内容的提问来源于stack exchange，提问作者Sowmiya R