#!/usr/bin/env python3
import numpy as np
def _get_block (lines, keys) :
for idx in range(len(lines)) :
if keys in lines[idx] :
break
if idx == len(lines) - 1:
return None
idx_s = idx+2
idx = idx_s
ret = []
while True :
if len(lines[idx].split()) == 0 :
break
else :
ret.append(lines[idx])
idx += 1
return ret
[docs]def lmpbox2box(lohi, tilt) :
xy = tilt[0]
xz = tilt[1]
yz = tilt[2]
orig = np.array([lohi[0][0], lohi[1][0], lohi[2][0]])
lens = []
for dd in range(3) :
lens.append(lohi[dd][1] - lohi[dd][0])
xx = [lens[0], 0, 0]
yy = [xy, lens[1], 0]
zz= [xz, yz, lens[2]]
return orig, np.array([xx, yy, zz])
[docs]def box2lmpbox(orig, box) :
lohi = np.zeros([3,2])
for dd in range(3) :
lohi[dd][0] = orig[dd]
tilt = np.zeros(3)
tilt[0] = box[1][0]
tilt[1] = box[2][0]
tilt[2] = box[2][1]
lens = np.zeros(3)
lens[0] = box[0][0]
lens[1] = box[1][1]
lens[2] = box[2][2]
for dd in range(3) :
lohi[dd][1] = lohi[dd][0] + lens[dd]
return lohi, tilt
[docs]def get_atoms(lines) :
return _get_block(lines, 'Atoms')
[docs]def get_natoms(lines) :
for ii in lines :
if 'atoms' in ii :
return int(ii.split()[0])
return None
[docs]def get_natomtypes(lines) :
for ii in lines :
if 'atom types' in ii :
return int(ii.split()[0])
return None
def _atom_info_mol(line) :
vec = line.split()
# idx, mole_type, atom_type, charge, x, y, z
return int(vec[0]), int(vec[1]), int(vec[2]), float(vec[3]), float(vec[4]), float(vec[5]), float(vec[6])
def _atom_info_atom(line) :
vec = line.split()
# idx, atom_type, x, y, z
return int(vec[0]), int(vec[1]), float(vec[2]), float(vec[3]), float(vec[4])
[docs]def get_natoms_vec(lines) :
atype = get_atype(lines)
natoms_vec = []
natomtypes = get_natomtypes(lines)
for ii in range(natomtypes) :
natoms_vec.append(sum(atype == ii+1))
assert (sum(natoms_vec) == get_natoms(lines))
return natoms_vec
[docs]def get_atype(lines) :
alines = get_atoms(lines)
atype = []
for ii in alines :
# idx, mt, at, q, x, y, z = _atom_info_mol(ii)
idx, at, x, y, z = _atom_info_atom(ii)
atype.append(at)
return np.array(atype, dtype = int)
[docs]def get_posi(lines) :
atom_lines = get_atoms(lines)
posis = []
for ii in atom_lines :
# posis.append([float(jj) for jj in ii.split()[4:7]])
posis.append([float(jj) for jj in ii.split()[2:5]])
return np.array(posis)
[docs]def get_lmpbox(lines) :
box_info = []
tilt = np.zeros(3)
for ii in lines :
if 'xlo' in ii and 'xhi' in ii :
box_info.append([float(ii.split()[0]), float(ii.split()[1])])
break
for ii in lines :
if 'ylo' in ii and 'yhi' in ii :
box_info.append([float(ii.split()[0]), float(ii.split()[1])])
break
for ii in lines :
if 'zlo' in ii and 'zhi' in ii :
box_info.append([float(ii.split()[0]), float(ii.split()[1])])
break
for ii in lines :
if 'xy' in ii and 'xz' in ii and 'yz' in ii :
tilt = np.array([float(jj) for jj in ii.split()[0:3]])
return box_info, tilt
[docs]def system_data(lines) :
system = {}
system['atom_numbs'] = get_natoms_vec(lines)
system['atom_names'] = []
for ii in range(len(system['atom_numbs'])) :
system['atom_names'].append('Type_%d' % ii)
lohi, tilt = get_lmpbox(lines)
orig, cell = lmpbox2box(lohi, tilt)
system['orig'] = np.array(orig)
system['cell'] = np.array(cell)
natoms = sum(system['atom_numbs'])
system['atom_types'] = get_atype(lines)
system['coordinates'] = get_posi(lines)
return system
[docs]def to_system_data(lines) :
return system_data(lines)
[docs]def from_system_data(system) :
ret = ''
ret += '\n'
natoms = sum(system['atom_numbs'])
ntypes = len(system['atom_numbs'])
ret += '%d atoms\n' % natoms
ret += '%d atom types\n' % ntypes
ret += '0 %f xlo xhi\n' % system['cell'][0][0]
ret += '0 %f ylo yhi\n' % system['cell'][1][1]
ret += '0 %f zlo zhi\n' % system['cell'][2][2]
ret += '%f %f %f xy xz yz\n' % \
(system['cell'][1][0], system['cell'][2][0], system['cell'][2][1])
ret += '\n'
ret += 'Atoms # atomic\n'
ret += '\n'
for ii in range(natoms) :
ret += '%d %d %f %f %f\n' % \
(ii+1,
system['atom_types'][ii],
system['coordinates'][ii][0] - system['orig'][0],
system['coordinates'][ii][1] - system['orig'][1],
system['coordinates'][ii][2] - system['orig'][2]
)
return ret
if __name__ == '__main__' :
fname = 'water-SPCE.data'
lines = open(fname).read().split('\n')
bonds, tilt = get_lmpbox(lines)
# print(bonds, tilt)
orig, box = lmpbox2box(bonds, tilt)
# print(orig, box)
bonds1, tilt1 = box2lmpbox(orig, box)
# print(bonds1, tilt1)
print(bonds1 - bonds)
print(tilt1 - tilt)
print(box)
print(get_atype(lines))
print(get_posi(lines))