lammps输出轨迹转标准格式xyz

事情是这样的，跑了几个小时lammps的到结果后想用脚本对dump的xyz轨迹进行分析，结果发现居然不支持。原来lammps的轨迹太过自由了，并不是标准的xyz文件。很多软件和脚本都只支持标准的xyz文件（如MS，VMD，MDAnalysis）等。

 标准xyz文件格式如下：

12
benzene example
  C        0.00000        1.40272        0.00000
  H        0.00000        2.49029        0.00000
  C       -1.21479        0.70136        0.00000
  H       -2.15666        1.24515        0.00000
  C       -1.21479       -0.70136        0.00000
  H       -2.15666       -1.24515        0.00000
  C        0.00000       -1.40272        0.00000
  H        0.00000       -2.49029        0.00000
  C        1.21479       -0.70136        0.00000
  H        2.15666       -1.24515        0.00000
  C        1.21479        0.70136        0.00000
  H        2.15666        1.24515        0.00000

第一列为元素符号，后面三列为xyz坐标。第一行为原子数，第二行随意。

在跑lammps之前，可以按下面连接的方法设置dump格式为标准的xyz：https://www.cnblogs.com/sysu/p/17006091.html

我在这次跑仿真前忘记设置输出格式了，使用的以下命令进行dump：

dump 1 all custom 1000 shear.xyz type x y z

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
3128
ITEM: BOX BOUNDS xy xz yz pp pp pp
-3.7243842623499729e-01 3.2234124596235063e+01 0.0000000000000000e+00
4.1885043476486233e-01 3.3025413457234471e+01 0.0000000000000000e+00
8.5273312137650272e+00 4.1133894236234468e+01 0.0000000000000000e+00
ITEM: ATOMS type x y z
2 6.27403 7.45496 12.7015
9 2.02983 8.40239 14.0077
9 3.38754 3.88625 10.3987
3 2.91287 4.68163 14.4446
9 5.12803 5.28437 11.3195
9 3.91993 6.71038 8.8808
......
ITEM: TIMESTEP
1000
ITEM: NUMBER OF ATOMS
3128
ITEM: BOX BOUNDS xy xz yz pp pp pp
-3.3698973156626344e-01 3.2329102153652570e+01 0.0000000000000000e+00
4.6178284437435391e-01 3.2982481047625178e+01 1.3042625208987602e-01
8.6006685257616358e+00 4.1060556924237886e+01 0.0000000000000000e+00
ITEM: ATOMS type x y z
2 6.29072 7.49291 12.7049
9 2.18091 8.43777 14.0266
9 3.35011 3.92172 10.4029

需要处理的是将第一列的序号换成元素符号，在原子坐标之前的许多多余行进行替换。

不想花钱按之前连接写的方法重跑仿真，于是写了个简单的python进行修改：

import re

# 字母和数字的对应关系
atom_dict = {
    '1': 'O', '2': 'C', '3': 'N', '4': 'C', '5': 'C', '6': 'O', '7': 'C',
    '8': 'O', '9': 'H', '10': 'O', '11': 'C', '12': 'H', '13': 'C', '14': 'H'
}

# 读取数据文件
with open('shear.xyz', 'r') as f:
    data = f.readlines()

# 处理数据
for i in range(len(data)):
    if re.match(r'^ITEM: BOX BOUNDS xy xz yz pp pp pp', data[i]):
        # 替换第一列数字为字母
        for j in range(i+1, len(data)):
            if re.match(r'^ITEM: TIMESTEP', data[j]):
                break
            else:
                line = data[j].strip().split()
                if len(line) == 4:
                    line[0] = atom_dict[line[0]]
                    data[j] = ' '.join(line) + '\n'

i = 0
output = []
while i < len(data):
    if data[i].startswith("ITEM: TIMESTEP"):
        output.extend(["3128\n", f"Atoms. Timestep: {i//9}\n"])
        i += 9
    else:
        output.append(data[i])
        i += 1

# 将处理后的数据写入新文件
with open('new.xyz', 'w') as f:
    f.writelines(output)

最终处理后的data文件格式如下：