5.转换lammps的轨迹文件为XDATCAR
from ase.io import read, write
from sys import argv
if __name__ == '__main__':
# lammps strj
strj_file = argv[1]
write_dir = argv[2]
print(f"From {strj_file} to {write_dir}")
data = read(strj_file, index=":")
if not os.path.exists(write_dir):
os.mkdir(write_dir)
for i in range(len(data)):
poscar_name = os.path.join(write_dir, str(i) + ".vasp")
write(poscar_name, data[i])
另一种方法：需要搞到晶格参数（ovito或者上面的脚本）
import numpy as np
import time
def timer(func):
def wrapper(*args, **kwargs):
start_time = time.time()
result = func(*args, **kwargs)
end_time = time.time()
print(f"Time elapsed: {end_time - start_time:.4f} seconds")
return result
return wrapper
@timer
def read_lmp_trj(file_name: str):
f = open(file_name, 'r')
"""data segment"""
n_atoms_list = []
lat = []
pos = []
num1 = 0
#跳过文件的行数，需要修改
num_lines_to_skip = 0
current_line_num = 0
"""temp"""
# n_atoms = None
while True:
line = f.readline()
if not line:
break
current_line_num += 1
if current_line_num <= num_lines_to_skip:
continue
if "NUMBER OF ATOMS" in line:
line = f.readline()
n_atoms = int(line)
n_atoms_list.append(n_atoms)
f.readline()
"""lat"""
for i in range(3):
lines = f.readline()
lat.append(lines)
f.readline()
"""pos"""
for i in range(n_atoms):
pos.append(f.readline())
#每隔1000个时间间隔读取一次 需要改
num1 +=1
"""post treatment"""
#固定一个三维数组的两个维度，第三个维度也固定了，numpy的逻辑应该是3-2-1，第三个维度先排布，依次向前，第一个是深度
lat = np.array([i.split() for i in lat]).astype(float).reshape(-1, 3, 3)
#需要看这里的8需不需要改
pos = np.array([i.split() for i in pos]).astype(float).reshape(len(lat), -1, 8)
# idx = pos[:, :, 0].argsort(1)
#lexsort的顺序是从最后一个传入的参数开始，依次到第一个传入的参数
idx =np.lexsort((pos[:,:,0],pos[:,:,1]))
sort_pos = pos[np.arange(len(pos))[:, None], idx, :]
start3 = sort_pos[:,:,0:2].reshape(-1,2)
n_atoms_list = np.array(n_atoms_list)
return n_atoms_list, lat, sort_pos[:, :, [2, 3, 4]],start3
if __name__ == '__main__':
file_name = "A.lammpstrj"
print("Loading lmp trajectory")
n_atoms, _, x,stc = read_lmp_trj(file_name=file_name)
#np.savetxt('first_element.txt', x[0])
#np.savetxt('end_element.txt',x[-1])
#30这里需要修改，是移动的原子个数
with open('XDATCAR', 'w') as f:
for i, element in enumerate(x):
f.write(f'Direct configuration = {i + 1}\n') # 注意这里的 i + 1
np.savetxt(f, element, fmt='%.10f', delimiter='\t')
转载请注明来源 有问题可通过github提交issue