事情是这样的,跑了几个小时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文件格式如下: