From bbacdfb068f82ac70c7c3b8494c5ab9e49824f14 Mon Sep 17 00:00:00 2001 From: zezhong-zhang Date: Tue, 9 Aug 2022 10:05:45 +0200 Subject: [PATCH 1/6] add unwarp support for lammps dump --- dpdata/lammps/dump.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index 7136d3e63..bdc84dd7d 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -59,16 +59,17 @@ def get_coordtype_and_scalefactor(keys): key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter lmp_coor_type = [key_pc,key_uc,key_s,key_su] sf = [0,0,1,1] + unwarp = [0, 1, 0, 1] for k in range(4): if all(i in keys for i in lmp_coor_type[k]): - return lmp_coor_type[k],sf[k] + return lmp_coor_type[k], sf[k], unwarp[k] def safe_get_posi(lines,cell,orig=np.zeros(3)) : blk, head = _get_block(lines, 'ATOMS') keys = head.split() coord_tp_and_sf = get_coordtype_and_scalefactor(keys) assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' - coordtype, sf = coord_tp_and_sf + coordtype, sf, unwarp = coord_tp_and_sf id_idx = keys.index('id') - 2 xidx = keys.index(coordtype[0])-2 yidx = keys.index(coordtype[1])-2 @@ -81,8 +82,11 @@ def safe_get_posi(lines,cell,orig=np.zeros(3)) : posis.sort() posis = np.array(posis)[:,1:4] if not sf: - posis = (posis-orig)@np.linalg.inv(cell)# Convert to scaled coordinates for unscaled coordinates - return (posis%1)@cell # Convert scaled coordinates back to Cartesien coordinates + posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates + if not unwarp: + return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with warping at periodic boundary conditions + else: + return posis @ cell # without warping at periodic boundary conditions def get_dumpbox(lines) : blk, h = _get_block(lines, 'BOX BOUNDS') @@ -145,7 +149,7 @@ def load_file(fname, begin = 0, step = 1) : cc += 1 if cc >= begin and (cc - begin) % step == 0 : buff.append(line) - + def system_data(lines, type_map = None, type_idx_zero = True) : array_lines = split_traj(lines) @@ -181,7 +185,7 @@ def split_traj(dump_lines) : marks = [] for idx,ii in enumerate(dump_lines) : if 'ITEM: TIMESTEP' in ii : - marks.append(idx) + marks.append(idx) if len(marks) == 0 : return None elif len(marks) == 1 : @@ -191,9 +195,9 @@ def split_traj(dump_lines) : ret = [] for ii in marks : ret.append(dump_lines[ii:ii+block_size]) - # for ii in range(len(marks)-1): + # for ii in range(len(marks)-1): # assert(marks[ii+1] - marks[ii] == block_size) - return ret + return ret return None From 87473dbd326fc9d5ea39b13334016e7a36cfac3a Mon Sep 17 00:00:00 2001 From: Zezhong Zhang Date: Thu, 11 Aug 2022 00:30:41 +0200 Subject: [PATCH 2/6] let user choose to unwrap and add unittest case --- dpdata/lammps/dump.py | 236 +++++++++++++++++-------------- dpdata/plugins/lammps.py | 17 +-- tests/poscars/poscar_ref_oh.py | 62 ++++---- tests/test_lammps_dump_unfold.py | 42 ++++-- 4 files changed, 201 insertions(+), 156 deletions(-) diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index bdc84dd7d..e39e72aa7 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -2,206 +2,230 @@ import os, sys import numpy as np + lib_path = os.path.dirname(os.path.realpath(__file__)) sys.path.append(lib_path) import lmp -def _get_block (lines, key) : - for idx in range(len(lines)) : - if ('ITEM: ' + key) in lines[idx] : + +def _get_block(lines, key): + for idx in range(len(lines)): + if ("ITEM: " + key) in lines[idx]: break idx_s = idx + 1 - for idx in range(idx_s, len(lines)) : - if ('ITEM: ') in lines[idx] : + for idx in range(idx_s, len(lines)): + if ("ITEM: ") in lines[idx]: break idx_e = idx - if idx_e == len(lines)-1 : + if idx_e == len(lines) - 1: idx_e += 1 - return lines[idx_s:idx_e], lines[idx_s-1] + return lines[idx_s:idx_e], lines[idx_s - 1] -def get_atype(lines, type_idx_zero = False) : - blk, head = _get_block(lines, 'ATOMS') + +def get_atype(lines, type_idx_zero=False): + blk, head = _get_block(lines, "ATOMS") keys = head.split() - id_idx = keys.index('id') - 2 - tidx = keys.index('type') - 2 + id_idx = keys.index("id") - 2 + tidx = keys.index("type") - 2 atype = [] - for ii in blk : + for ii in blk: atype.append([int(ii.split()[id_idx]), int(ii.split()[tidx])]) atype.sort() - atype = np.array(atype, dtype = int) - if type_idx_zero : - return atype[:,1] - 1 - else : - return atype[:,1] - -def get_natoms(lines) : - blk, head = _get_block(lines, 'NUMBER OF ATOMS') + atype = np.array(atype, dtype=int) + if type_idx_zero: + return atype[:, 1] - 1 + else: + return atype[:, 1] + + +def get_natoms(lines): + blk, head = _get_block(lines, "NUMBER OF ATOMS") return int(blk[0]) -def get_natomtypes(lines) : + +def get_natomtypes(lines): atype = get_atype(lines) return max(atype) -def get_natoms_vec(lines) : + +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)) + for ii in range(natomtypes): + natoms_vec.append(sum(atype == ii + 1)) + assert sum(natoms_vec) == get_natoms(lines) return natoms_vec + def get_coordtype_and_scalefactor(keys): # 4 types in total,with different scaling factor - key_pc=['x','y','z'] # plain cartesian, sf = 1 - key_uc=['xu','yu','zu'] # unwraped cartesian, sf = 1 - key_s=['xs','ys','zs'] # scaled by lattice parameter, sf = lattice parameter - key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter - lmp_coor_type = [key_pc,key_uc,key_s,key_su] - sf = [0,0,1,1] - unwarp = [0, 1, 0, 1] + key_pc = ["x", "y", "z"] # plain cartesian, sf = 1 + key_uc = ["xu", "yu", "zu"] # unwraped cartesian, sf = 1 + key_s = ["xs", "ys", "zs"] # scaled by lattice parameter, sf = lattice parameter + key_su = ["xsu", "ysu", "zsu"] # scaled and unfolded,sf = lattice parameter + lmp_coor_type = [key_pc, key_uc, key_s, key_su] + sf = [0, 0, 1, 1] + uw = [0, 1, 0, 1] # unwraped or not for k in range(4): if all(i in keys for i in lmp_coor_type[k]): - return lmp_coor_type[k], sf[k], unwarp[k] + return lmp_coor_type[k], sf[k], uw[k] + -def safe_get_posi(lines,cell,orig=np.zeros(3)) : - blk, head = _get_block(lines, 'ATOMS') +def safe_get_posi(lines, cell, orig=np.zeros(3), unwrap=False): + blk, head = _get_block(lines, "ATOMS") keys = head.split() coord_tp_and_sf = get_coordtype_and_scalefactor(keys) - assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' - coordtype, sf, unwarp = coord_tp_and_sf - id_idx = keys.index('id') - 2 - xidx = keys.index(coordtype[0])-2 - yidx = keys.index(coordtype[1])-2 - zidx = keys.index(coordtype[2])-2 + assert coord_tp_and_sf is not None, "Dump file does not contain atomic coordinates!" + coordtype, sf, uw = coord_tp_and_sf + id_idx = keys.index("id") - 2 + xidx = keys.index(coordtype[0]) - 2 + yidx = keys.index(coordtype[1]) - 2 + zidx = keys.index(coordtype[2]) - 2 sel = (xidx, yidx, zidx) posis = [] - for ii in blk : + for ii in blk: words = ii.split() - posis.append([float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])]) + posis.append( + [ + float(words[id_idx]), + float(words[xidx]), + float(words[yidx]), + float(words[zidx]), + ] + ) posis.sort() - posis = np.array(posis)[:,1:4] + posis = np.array(posis)[:, 1:4] if not sf: - posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates - if not unwarp: - return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with warping at periodic boundary conditions + posis = (posis - orig) @ np.linalg.inv( + cell + ) # Convert to scaled coordinates for unscaled coordinates + if uw and unwrap: + return posis @ cell # unwrap at the periodic boundaries else: - return posis @ cell # without warping at periodic boundary conditions + return ( + posis % 1 + ) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions -def get_dumpbox(lines) : - blk, h = _get_block(lines, 'BOX BOUNDS') - bounds = np.zeros([3,2]) + +def get_dumpbox(lines): + blk, h = _get_block(lines, "BOX BOUNDS") + bounds = np.zeros([3, 2]) tilt = np.zeros([3]) - load_tilt = 'xy xz yz' in h - for dd in range(3) : + load_tilt = "xy xz yz" in h + for dd in range(3): info = [float(jj) for jj in blk[dd].split()] bounds[dd][0] = info[0] bounds[dd][1] = info[1] - if load_tilt : + if load_tilt: tilt[dd] = info[2] return bounds, tilt -def dumpbox2box(bounds, tilt) : + +def dumpbox2box(bounds, tilt): xy = tilt[0] xz = tilt[1] yz = tilt[2] - xlo = bounds[0][0] - min(0.0,xy,xz,xy+xz) - xhi = bounds[0][1] - max(0.0,xy,xz,xy+xz) - ylo = bounds[1][0] - min(0.0,yz) - yhi = bounds[1][1] - max(0.0,yz) + xlo = bounds[0][0] - min(0.0, xy, xz, xy + xz) + xhi = bounds[0][1] - max(0.0, xy, xz, xy + xz) + ylo = bounds[1][0] - min(0.0, yz) + yhi = bounds[1][1] - max(0.0, yz) zlo = bounds[2][0] zhi = bounds[2][1] info = [[xlo, xhi], [ylo, yhi], [zlo, zhi]] return lmp.lmpbox2box(info, tilt) -def box2dumpbox(orig, box) : + +def box2dumpbox(orig, box): lohi, tilt = lmp.box2lmpbox(orig, box) xy = tilt[0] xz = tilt[1] yz = tilt[2] - bounds = np.zeros([3,2]) - bounds[0][0] = lohi[0][0] + min(0.0,xy,xz,xy+xz) - bounds[0][1] = lohi[0][1] + max(0.0,xy,xz,xy+xz) - bounds[1][0] = lohi[1][0] + min(0.0,yz) - bounds[1][1] = lohi[1][1] + max(0.0,yz) + bounds = np.zeros([3, 2]) + bounds[0][0] = lohi[0][0] + min(0.0, xy, xz, xy + xz) + bounds[0][1] = lohi[0][1] + max(0.0, xy, xz, xy + xz) + bounds[1][0] = lohi[1][0] + min(0.0, yz) + bounds[1][1] = lohi[1][1] + max(0.0, yz) bounds[2][0] = lohi[2][0] bounds[2][1] = lohi[2][1] return bounds, tilt -def load_file(fname, begin = 0, step = 1) : +def load_file(fname, begin=0, step=1): lines = [] buff = [] cc = -1 with open(fname) as fp: while True: - line = fp.readline().rstrip('\n') - if not line : - if cc >= begin and (cc - begin) % step == 0 : + line = fp.readline().rstrip("\n") + if not line: + if cc >= begin and (cc - begin) % step == 0: lines += buff buff = [] cc += 1 return lines - if 'ITEM: TIMESTEP' in line : - if cc >= begin and (cc - begin) % step == 0 : + if "ITEM: TIMESTEP" in line: + if cc >= begin and (cc - begin) % step == 0: lines += buff buff = [] cc += 1 - if cc >= begin and (cc - begin) % step == 0 : + if cc >= begin and (cc - begin) % step == 0: buff.append(line) -def system_data(lines, type_map = None, type_idx_zero = True) : +def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): array_lines = split_traj(lines) lines = array_lines[0] system = {} - system['atom_numbs'] = get_natoms_vec(lines) - system['atom_names'] = [] - if type_map == None : - for ii in range(len(system['atom_numbs'])) : - system['atom_names'].append('TYPE_%d' % ii) - else : - assert(len(type_map) >= len(system['atom_numbs'])) - for ii in range(len(system['atom_numbs'])) : - system['atom_names'].append(type_map[ii]) + system["atom_numbs"] = get_natoms_vec(lines) + system["atom_names"] = [] + if type_map == None: + for ii in range(len(system["atom_numbs"])): + system["atom_names"].append("TYPE_%d" % ii) + else: + assert len(type_map) >= len(system["atom_numbs"]) + for ii in range(len(system["atom_numbs"])): + system["atom_names"].append(type_map[ii]) bounds, tilt = get_dumpbox(lines) orig, cell = dumpbox2box(bounds, tilt) - system['orig'] = np.array(orig) - np.array(orig) - system['cells'] = [np.array(cell)] - natoms = sum(system['atom_numbs']) - system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero) - system['coords'] = [safe_get_posi(lines, cell, np.array(orig))] - for ii in range(1, len(array_lines)) : + system["orig"] = np.array(orig) - np.array(orig) + system["cells"] = [np.array(cell)] + natoms = sum(system["atom_numbs"]) + system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero) + system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] + for ii in range(1, len(array_lines)): bounds, tilt = get_dumpbox(array_lines[ii]) orig, cell = dumpbox2box(bounds, tilt) - system['cells'].append(cell) - system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig))) - system['cells'] = np.array(system['cells']) - system['coords'] = np.array(system['coords']) + system["cells"].append(cell) + system["coords"].append( + safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap) + ) + system["cells"] = np.array(system["cells"]) + system["coords"] = np.array(system["coords"]) return system -def split_traj(dump_lines) : +def split_traj(dump_lines): marks = [] - for idx,ii in enumerate(dump_lines) : - if 'ITEM: TIMESTEP' in ii : + for idx, ii in enumerate(dump_lines): + if "ITEM: TIMESTEP" in ii: marks.append(idx) - if len(marks) == 0 : + if len(marks) == 0: return None - elif len(marks) == 1 : + elif len(marks) == 1: return [dump_lines] - else : + else: block_size = marks[1] - marks[0] ret = [] - for ii in marks : - ret.append(dump_lines[ii:ii+block_size]) + for ii in marks: + ret.append(dump_lines[ii : ii + block_size]) # for ii in range(len(marks)-1): # assert(marks[ii+1] - marks[ii] == block_size) return ret return None -if __name__ == '__main__' : +if __name__ == "__main__": # fname = 'dump.hti' # lines = open(fname).read().split('\n') # # print(get_natoms(lines)) @@ -215,9 +239,9 @@ def split_traj(dump_lines) : # print(box) # np.savetxt('tmp.out', posi - orig, fmt='%.6f') # print(system_data(lines)) - lines = load_file('conf_unfold.dump', begin = 0, step = 1) + lines = load_file("conf_unfold.dump", begin=0, step=1) al = split_traj(lines) - s = system_data(lines,['O','H']) - #l = np.linalg.norm(s['cells'][1],axis=1) - #p = s['coords'][0] + l - #np.savetxt('p',p,fmt='%1.10f') + s = system_data(lines, ["O", "H"]) + # l = np.linalg.norm(s['cells'][1],axis=1) + # p = s['coords'][0] + l + # np.savetxt('p',p,fmt='%1.10f') diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py index f1356bb1d..d4bce01b4 100644 --- a/dpdata/plugins/lammps.py +++ b/dpdata/plugins/lammps.py @@ -9,7 +9,7 @@ class LAMMPSLmpFormat(Format): @Format.post("shift_orig_zero") def from_system(self, file_name, type_map=None, **kwargs): with open(file_name) as fp: - lines = [line.rstrip('\n') for line in fp] + lines = [line.rstrip("\n") for line in fp] return dpdata.lammps.lmp.to_system_data(lines, type_map) def to_system(self, data, file_name, frame_idx=0, **kwargs): @@ -25,9 +25,9 @@ def to_system(self, data, file_name, frame_idx=0, **kwargs): frame_idx : int The index of the frame to dump """ - assert(frame_idx < len(data['coords'])) + assert frame_idx < len(data["coords"]) w_str = dpdata.lammps.lmp.from_system_data(data, frame_idx) - with open(file_name, 'w') as fp: + with open(file_name, "w") as fp: fp.write(w_str) @@ -35,11 +35,8 @@ def to_system(self, data, file_name, frame_idx=0, **kwargs): @Format.register("lammps/dump") class LAMMPSDumpFormat(Format): @Format.post("shift_orig_zero") - def from_system(self, - file_name, - type_map=None, - begin=0, - step=1, - **kwargs): + def from_system( + self, file_name, type_map=None, begin=0, step=1, unwrap=False, **kwargs + ): lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step) - return dpdata.lammps.dump.system_data(lines, type_map) + return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py index 23bab644f..f120183ed 100644 --- a/tests/poscars/poscar_ref_oh.py +++ b/tests/poscars/poscar_ref_oh.py @@ -1,38 +1,50 @@ import numpy as np -class TestPOSCARoh : +class TestPOSCARoh: def test_atom_numbs(self): - self.assertEqual(self.system.data['atom_numbs'], [1,1]) + self.assertEqual(self.system.data["atom_numbs"], [1, 1]) def test_atom_names(self): - self.assertEqual(self.system.data['atom_names'], ['O','H']) + self.assertEqual(self.system.data["atom_names"], ["O", "H"]) def test_atom_types(self): - self.assertEqual(self.system.data['atom_types'][0], 0) - self.assertEqual(self.system.data['atom_types'][1], 1) + self.assertEqual(self.system.data["atom_types"][0], 0) + self.assertEqual(self.system.data["atom_types"][1], 1) def test_orig(self): - for d0 in range(3) : - self.assertEqual(self.system.data['orig'][d0], 0) + for d0 in range(3): + self.assertEqual(self.system.data["orig"][d0], 0) def test_cell(self): - ovito_cell = np.array([[2.5243712, 0.0000000, 0.0000000], - [1.2621856, 2.0430257, 0.0000000], - [1.2874292, 0.7485898, 2.2254033]]) - for ii in range(3) : - for jj in range(3) : - self.assertAlmostEqual(self.system.data['cells'][0][ii][jj], - ovito_cell[ii][jj], - places = 6, - msg = 'cell[%d][%d] failed' % (ii,jj)) + ovito_cell = np.array( + [ + [2.5243712, 0.0000000, 0.0000000], + [1.2621856, 2.0430257, 0.0000000], + [1.2874292, 0.7485898, 2.2254033], + ] + ) + for ii in range(3): + for jj in range(3): + self.assertAlmostEqual( + self.system.data["cells"][0][ii][jj], + ovito_cell[ii][jj], + places=6, + msg="cell[%d][%d] failed" % (ii, jj), + ) - def test_frame(self): - ovito_posis = np.array([[0, 0, 0], - [1.2621856, 0.7018028, 0.5513885]]) - for ii in range(2) : - for jj in range(3) : - self.assertAlmostEqual(self.system.data['coords'][0][ii][jj], - ovito_posis[ii][jj], - places = 6, - msg = 'posis[%d][%d] failed' % (ii,jj)) + def test_frame(self): + if hasattr(self, "unwrap") and self.unwrap is True: + ovito_posis = np.array( + [[5.0739861, 2.7916155, 2.2254033], [6.3361717, 3.4934183, 2.7767918]] + ) + else: + ovito_posis = np.array([[0, 0, 0], [1.2621856, 0.7018028, 0.5513885]]) + for ii in range(2): + for jj in range(3): + self.assertAlmostEqual( + self.system.data["coords"][0][ii][jj], + ovito_posis[ii][jj], + places=6, + msg="posis[%d][%d] failed" % (ii, jj), + ) diff --git a/tests/test_lammps_dump_unfold.py b/tests/test_lammps_dump_unfold.py index 4ffd5f730..b4c61195d 100644 --- a/tests/test_lammps_dump_unfold.py +++ b/tests/test_lammps_dump_unfold.py @@ -1,26 +1,38 @@ +from inspect import unwrap import os import numpy as np import unittest from context import dpdata -from poscars.poscar_ref_oh import TestPOSCARoh +from poscars.poscar_ref_oh import TestPOSCARoh + class TestDump(unittest.TestCase, TestPOSCARoh): - - def setUp(self): - self.system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), - type_map = ['O', 'H']) - + def setUp(self): + self.system = dpdata.System( + os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] + ) + + class TestDump2(unittest.TestCase, TestPOSCARoh): - - def setUp(self): - self.tmp_system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), - type_map = ['O', 'H']) + def setUp(self): + self.tmp_system = dpdata.System( + os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] + ) self.system = self.tmp_system.sub_system([1]) - def test_nframes (self) : + def test_nframes(self): self.assertEqual(self.tmp_system.get_nframes(), 2) - - -if __name__ == '__main__': + + +class TestDumpUnwrap(unittest.TestCase, TestPOSCARoh): + def setUp(self): + self.unwrap = True + self.system = dpdata.System( + os.path.join("poscars", "conf_unfold.dump"), + type_map=["O", "H"], + unwrap=self.unwrap, + ) + + +if __name__ == "__main__": unittest.main() - From bd571260e3470743cfd78755238c94b30b15127a Mon Sep 17 00:00:00 2001 From: zezhong-zhang Date: Thu, 11 Aug 2022 12:04:26 +0200 Subject: [PATCH 3/6] temp to revert format changes --- dpdata/lammps/dump.py | 228 ++++++------- dpdata/plugins/lammps.py | 16 +- my.patch | 531 +++++++++++++++++++++++++++++++ tests/poscars/poscar_ref_oh.py | 60 ++-- tests/test_lammps_dump_unfold.py | 36 +-- 5 files changed, 686 insertions(+), 185 deletions(-) create mode 100644 my.patch diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index e39e72aa7..9de734e26 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -2,230 +2,206 @@ import os, sys import numpy as np - lib_path = os.path.dirname(os.path.realpath(__file__)) sys.path.append(lib_path) import lmp - -def _get_block(lines, key): - for idx in range(len(lines)): - if ("ITEM: " + key) in lines[idx]: +def _get_block (lines, key) : + for idx in range(len(lines)) : + if ('ITEM: ' + key) in lines[idx] : break idx_s = idx + 1 - for idx in range(idx_s, len(lines)): - if ("ITEM: ") in lines[idx]: + for idx in range(idx_s, len(lines)) : + if ('ITEM: ') in lines[idx] : break idx_e = idx - if idx_e == len(lines) - 1: + if idx_e == len(lines)-1 : idx_e += 1 - return lines[idx_s:idx_e], lines[idx_s - 1] + return lines[idx_s:idx_e], lines[idx_s-1] - -def get_atype(lines, type_idx_zero=False): - blk, head = _get_block(lines, "ATOMS") +def get_atype(lines, type_idx_zero = False) : + blk, head = _get_block(lines, 'ATOMS') keys = head.split() - id_idx = keys.index("id") - 2 - tidx = keys.index("type") - 2 + id_idx = keys.index('id') - 2 + tidx = keys.index('type') - 2 atype = [] - for ii in blk: + for ii in blk : atype.append([int(ii.split()[id_idx]), int(ii.split()[tidx])]) atype.sort() - atype = np.array(atype, dtype=int) - if type_idx_zero: - return atype[:, 1] - 1 - else: - return atype[:, 1] - - -def get_natoms(lines): - blk, head = _get_block(lines, "NUMBER OF ATOMS") + atype = np.array(atype, dtype = int) + if type_idx_zero : + return atype[:,1] - 1 + else : + return atype[:,1] + +def get_natoms(lines) : + blk, head = _get_block(lines, 'NUMBER OF ATOMS') return int(blk[0]) - -def get_natomtypes(lines): +def get_natomtypes(lines) : atype = get_atype(lines) return max(atype) - -def get_natoms_vec(lines): +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) + for ii in range(natomtypes) : + natoms_vec.append(sum(atype == ii+1)) + assert (sum(natoms_vec) == get_natoms(lines)) return natoms_vec - def get_coordtype_and_scalefactor(keys): # 4 types in total,with different scaling factor - key_pc = ["x", "y", "z"] # plain cartesian, sf = 1 - key_uc = ["xu", "yu", "zu"] # unwraped cartesian, sf = 1 - key_s = ["xs", "ys", "zs"] # scaled by lattice parameter, sf = lattice parameter - key_su = ["xsu", "ysu", "zsu"] # scaled and unfolded,sf = lattice parameter - lmp_coor_type = [key_pc, key_uc, key_s, key_su] - sf = [0, 0, 1, 1] + key_pc=['x','y','z'] # plain cartesian, sf = 1 + key_uc=['xu','yu','zu'] # unwraped cartesian, sf = 1 + key_s=['xs','ys','zs'] # scaled by lattice parameter, sf = lattice parameter + key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter + lmp_coor_type = [key_pc,key_uc,key_s,key_su] + sf = [0,0,1,1] uw = [0, 1, 0, 1] # unwraped or not for k in range(4): if all(i in keys for i in lmp_coor_type[k]): return lmp_coor_type[k], sf[k], uw[k] - -def safe_get_posi(lines, cell, orig=np.zeros(3), unwrap=False): - blk, head = _get_block(lines, "ATOMS") +def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : + blk, head = _get_block(lines, 'ATOMS') keys = head.split() coord_tp_and_sf = get_coordtype_and_scalefactor(keys) - assert coord_tp_and_sf is not None, "Dump file does not contain atomic coordinates!" + assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' coordtype, sf, uw = coord_tp_and_sf - id_idx = keys.index("id") - 2 - xidx = keys.index(coordtype[0]) - 2 - yidx = keys.index(coordtype[1]) - 2 - zidx = keys.index(coordtype[2]) - 2 + id_idx = keys.index('id') - 2 + xidx = keys.index(coordtype[0])-2 + yidx = keys.index(coordtype[1])-2 + zidx = keys.index(coordtype[2])-2 sel = (xidx, yidx, zidx) posis = [] - for ii in blk: + for ii in blk : words = ii.split() - posis.append( - [ - float(words[id_idx]), - float(words[xidx]), - float(words[yidx]), - float(words[zidx]), - ] - ) + posis.append([float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])]) posis.sort() - posis = np.array(posis)[:, 1:4] + posis = np.array(posis)[:,1:4] if not sf: - posis = (posis - orig) @ np.linalg.inv( - cell - ) # Convert to scaled coordinates for unscaled coordinates + posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates if uw and unwrap: - return posis @ cell # unwrap at the periodic boundaries + return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries else: - return ( - posis % 1 - ) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions + return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions - -def get_dumpbox(lines): - blk, h = _get_block(lines, "BOX BOUNDS") - bounds = np.zeros([3, 2]) +def get_dumpbox(lines) : + blk, h = _get_block(lines, 'BOX BOUNDS') + bounds = np.zeros([3,2]) tilt = np.zeros([3]) - load_tilt = "xy xz yz" in h - for dd in range(3): + load_tilt = 'xy xz yz' in h + for dd in range(3) : info = [float(jj) for jj in blk[dd].split()] bounds[dd][0] = info[0] bounds[dd][1] = info[1] - if load_tilt: + if load_tilt : tilt[dd] = info[2] return bounds, tilt - -def dumpbox2box(bounds, tilt): +def dumpbox2box(bounds, tilt) : xy = tilt[0] xz = tilt[1] yz = tilt[2] - xlo = bounds[0][0] - min(0.0, xy, xz, xy + xz) - xhi = bounds[0][1] - max(0.0, xy, xz, xy + xz) - ylo = bounds[1][0] - min(0.0, yz) - yhi = bounds[1][1] - max(0.0, yz) + xlo = bounds[0][0] - min(0.0,xy,xz,xy+xz) + xhi = bounds[0][1] - max(0.0,xy,xz,xy+xz) + ylo = bounds[1][0] - min(0.0,yz) + yhi = bounds[1][1] - max(0.0,yz) zlo = bounds[2][0] zhi = bounds[2][1] info = [[xlo, xhi], [ylo, yhi], [zlo, zhi]] return lmp.lmpbox2box(info, tilt) - -def box2dumpbox(orig, box): +def box2dumpbox(orig, box) : lohi, tilt = lmp.box2lmpbox(orig, box) xy = tilt[0] xz = tilt[1] yz = tilt[2] - bounds = np.zeros([3, 2]) - bounds[0][0] = lohi[0][0] + min(0.0, xy, xz, xy + xz) - bounds[0][1] = lohi[0][1] + max(0.0, xy, xz, xy + xz) - bounds[1][0] = lohi[1][0] + min(0.0, yz) - bounds[1][1] = lohi[1][1] + max(0.0, yz) + bounds = np.zeros([3,2]) + bounds[0][0] = lohi[0][0] + min(0.0,xy,xz,xy+xz) + bounds[0][1] = lohi[0][1] + max(0.0,xy,xz,xy+xz) + bounds[1][0] = lohi[1][0] + min(0.0,yz) + bounds[1][1] = lohi[1][1] + max(0.0,yz) bounds[2][0] = lohi[2][0] bounds[2][1] = lohi[2][1] return bounds, tilt -def load_file(fname, begin=0, step=1): +def load_file(fname, begin = 0, step = 1) : lines = [] buff = [] cc = -1 with open(fname) as fp: while True: - line = fp.readline().rstrip("\n") - if not line: - if cc >= begin and (cc - begin) % step == 0: + line = fp.readline().rstrip('\n') + if not line : + if cc >= begin and (cc - begin) % step == 0 : lines += buff buff = [] cc += 1 return lines - if "ITEM: TIMESTEP" in line: - if cc >= begin and (cc - begin) % step == 0: + if 'ITEM: TIMESTEP' in line : + if cc >= begin and (cc - begin) % step == 0 : lines += buff buff = [] cc += 1 - if cc >= begin and (cc - begin) % step == 0: + if cc >= begin and (cc - begin) % step == 0 : buff.append(line) -def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): +def system_data(lines, type_map = None, type_idx_zero = True, unwrap=False) : array_lines = split_traj(lines) lines = array_lines[0] system = {} - system["atom_numbs"] = get_natoms_vec(lines) - system["atom_names"] = [] - if type_map == None: - for ii in range(len(system["atom_numbs"])): - system["atom_names"].append("TYPE_%d" % ii) - else: - assert len(type_map) >= len(system["atom_numbs"]) - for ii in range(len(system["atom_numbs"])): - system["atom_names"].append(type_map[ii]) + system['atom_numbs'] = get_natoms_vec(lines) + system['atom_names'] = [] + if type_map == None : + for ii in range(len(system['atom_numbs'])) : + system['atom_names'].append('TYPE_%d' % ii) + else : + assert(len(type_map) >= len(system['atom_numbs'])) + for ii in range(len(system['atom_numbs'])) : + system['atom_names'].append(type_map[ii]) bounds, tilt = get_dumpbox(lines) orig, cell = dumpbox2box(bounds, tilt) - system["orig"] = np.array(orig) - np.array(orig) - system["cells"] = [np.array(cell)] - natoms = sum(system["atom_numbs"]) - system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero) - system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] - for ii in range(1, len(array_lines)): + system['orig'] = np.array(orig) - np.array(orig) + system['cells'] = [np.array(cell)] + natoms = sum(system['atom_numbs']) + system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero) + system['coords'] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] + for ii in range(1, len(array_lines)) : bounds, tilt = get_dumpbox(array_lines[ii]) orig, cell = dumpbox2box(bounds, tilt) - system["cells"].append(cell) - system["coords"].append( - safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap) - ) - system["cells"] = np.array(system["cells"]) - system["coords"] = np.array(system["coords"]) + system['cells'].append(cell) + system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)) + system['cells'] = np.array(system['cells']) + system['coords'] = np.array(system['coords']) return system -def split_traj(dump_lines): +def split_traj(dump_lines) : marks = [] - for idx, ii in enumerate(dump_lines): - if "ITEM: TIMESTEP" in ii: + for idx,ii in enumerate(dump_lines) : + if 'ITEM: TIMESTEP' in ii : marks.append(idx) - if len(marks) == 0: + if len(marks) == 0 : return None - elif len(marks) == 1: + elif len(marks) == 1 : return [dump_lines] - else: + else : block_size = marks[1] - marks[0] ret = [] - for ii in marks: - ret.append(dump_lines[ii : ii + block_size]) + for ii in marks : + ret.append(dump_lines[ii:ii+block_size]) # for ii in range(len(marks)-1): # assert(marks[ii+1] - marks[ii] == block_size) return ret return None -if __name__ == "__main__": +if __name__ == '__main__' : # fname = 'dump.hti' # lines = open(fname).read().split('\n') # # print(get_natoms(lines)) @@ -239,9 +215,9 @@ def split_traj(dump_lines): # print(box) # np.savetxt('tmp.out', posi - orig, fmt='%.6f') # print(system_data(lines)) - lines = load_file("conf_unfold.dump", begin=0, step=1) + lines = load_file('conf_unfold.dump', begin = 0, step = 1) al = split_traj(lines) - s = system_data(lines, ["O", "H"]) - # l = np.linalg.norm(s['cells'][1],axis=1) - # p = s['coords'][0] + l - # np.savetxt('p',p,fmt='%1.10f') + s = system_data(lines,['O','H']) + #l = np.linalg.norm(s['cells'][1],axis=1) + #p = s['coords'][0] + l + #np.savetxt('p',p,fmt='%1.10f') diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py index d4bce01b4..8f9962967 100644 --- a/dpdata/plugins/lammps.py +++ b/dpdata/plugins/lammps.py @@ -9,7 +9,7 @@ class LAMMPSLmpFormat(Format): @Format.post("shift_orig_zero") def from_system(self, file_name, type_map=None, **kwargs): with open(file_name) as fp: - lines = [line.rstrip("\n") for line in fp] + lines = [line.rstrip('\n') for line in fp] return dpdata.lammps.lmp.to_system_data(lines, type_map) def to_system(self, data, file_name, frame_idx=0, **kwargs): @@ -25,9 +25,9 @@ def to_system(self, data, file_name, frame_idx=0, **kwargs): frame_idx : int The index of the frame to dump """ - assert frame_idx < len(data["coords"]) + assert(frame_idx < len(data['coords'])) w_str = dpdata.lammps.lmp.from_system_data(data, frame_idx) - with open(file_name, "w") as fp: + with open(file_name, 'w') as fp: fp.write(w_str) @@ -35,8 +35,12 @@ def to_system(self, data, file_name, frame_idx=0, **kwargs): @Format.register("lammps/dump") class LAMMPSDumpFormat(Format): @Format.post("shift_orig_zero") - def from_system( - self, file_name, type_map=None, begin=0, step=1, unwrap=False, **kwargs - ): + def from_system(self, + file_name, + type_map=None, + begin=0, + step=1, + unwrap=False, + **kwargs): lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step) return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) diff --git a/my.patch b/my.patch new file mode 100644 index 000000000..b203a16d1 --- /dev/null +++ b/my.patch @@ -0,0 +1,531 @@ +diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py +index e39e72a..bdc84dd 100644 +--- a/dpdata/lammps/dump.py ++++ b/dpdata/lammps/dump.py +@@ -2,230 +2,206 @@ + + import os, sys + import numpy as np +- + lib_path = os.path.dirname(os.path.realpath(__file__)) + sys.path.append(lib_path) + import lmp + +- +-def _get_block(lines, key): +- for idx in range(len(lines)): +- if ("ITEM: " + key) in lines[idx]: ++def _get_block (lines, key) : ++ for idx in range(len(lines)) : ++ if ('ITEM: ' + key) in lines[idx] : + break + idx_s = idx + 1 +- for idx in range(idx_s, len(lines)): +- if ("ITEM: ") in lines[idx]: ++ for idx in range(idx_s, len(lines)) : ++ if ('ITEM: ') in lines[idx] : + break + idx_e = idx +- if idx_e == len(lines) - 1: ++ if idx_e == len(lines)-1 : + idx_e += 1 +- return lines[idx_s:idx_e], lines[idx_s - 1] ++ return lines[idx_s:idx_e], lines[idx_s-1] + +- +-def get_atype(lines, type_idx_zero=False): +- blk, head = _get_block(lines, "ATOMS") ++def get_atype(lines, type_idx_zero = False) : ++ blk, head = _get_block(lines, 'ATOMS') + keys = head.split() +- id_idx = keys.index("id") - 2 +- tidx = keys.index("type") - 2 ++ id_idx = keys.index('id') - 2 ++ tidx = keys.index('type') - 2 + atype = [] +- for ii in blk: ++ for ii in blk : + atype.append([int(ii.split()[id_idx]), int(ii.split()[tidx])]) + atype.sort() +- atype = np.array(atype, dtype=int) +- if type_idx_zero: +- return atype[:, 1] - 1 +- else: +- return atype[:, 1] +- +- +-def get_natoms(lines): +- blk, head = _get_block(lines, "NUMBER OF ATOMS") ++ atype = np.array(atype, dtype = int) ++ if type_idx_zero : ++ return atype[:,1] - 1 ++ else : ++ return atype[:,1] ++ ++def get_natoms(lines) : ++ blk, head = _get_block(lines, 'NUMBER OF ATOMS') + return int(blk[0]) + +- +-def get_natomtypes(lines): ++def get_natomtypes(lines) : + atype = get_atype(lines) + return max(atype) + +- +-def get_natoms_vec(lines): ++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) ++ for ii in range(natomtypes) : ++ natoms_vec.append(sum(atype == ii+1)) ++ assert (sum(natoms_vec) == get_natoms(lines)) + return natoms_vec + +- + def get_coordtype_and_scalefactor(keys): + # 4 types in total,with different scaling factor +- key_pc = ["x", "y", "z"] # plain cartesian, sf = 1 +- key_uc = ["xu", "yu", "zu"] # unwraped cartesian, sf = 1 +- key_s = ["xs", "ys", "zs"] # scaled by lattice parameter, sf = lattice parameter +- key_su = ["xsu", "ysu", "zsu"] # scaled and unfolded,sf = lattice parameter +- lmp_coor_type = [key_pc, key_uc, key_s, key_su] +- sf = [0, 0, 1, 1] +- uw = [0, 1, 0, 1] # unwraped or not ++ key_pc=['x','y','z'] # plain cartesian, sf = 1 ++ key_uc=['xu','yu','zu'] # unwraped cartesian, sf = 1 ++ key_s=['xs','ys','zs'] # scaled by lattice parameter, sf = lattice parameter ++ key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter ++ lmp_coor_type = [key_pc,key_uc,key_s,key_su] ++ sf = [0,0,1,1] ++ uw = [0, 1, 0, 1] # unwraped or not + for k in range(4): + if all(i in keys for i in lmp_coor_type[k]): +- return lmp_coor_type[k], sf[k], uw[k] +- ++ return lmp_coor_type[k], sf[k], uw[k] + +-def safe_get_posi(lines, cell, orig=np.zeros(3), unwrap=False): +- blk, head = _get_block(lines, "ATOMS") ++def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : ++ blk, head = _get_block(lines, 'ATOMS') + keys = head.split() + coord_tp_and_sf = get_coordtype_and_scalefactor(keys) +- assert coord_tp_and_sf is not None, "Dump file does not contain atomic coordinates!" +- coordtype, sf, uw = coord_tp_and_sf +- id_idx = keys.index("id") - 2 +- xidx = keys.index(coordtype[0]) - 2 +- yidx = keys.index(coordtype[1]) - 2 +- zidx = keys.index(coordtype[2]) - 2 ++ assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' ++ coordtype, sf, uw = coord_tp_and_sf ++ id_idx = keys.index('id') - 2 ++ xidx = keys.index(coordtype[0])-2 ++ yidx = keys.index(coordtype[1])-2 ++ zidx = keys.index(coordtype[2])-2 + sel = (xidx, yidx, zidx) + posis = [] +- for ii in blk: ++ for ii in blk : + words = ii.split() +- posis.append( +- [ +- float(words[id_idx]), +- float(words[xidx]), +- float(words[yidx]), +- float(words[zidx]), +- ] +- ) ++ posis.append([float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])]) + posis.sort() +- posis = np.array(posis)[:, 1:4] ++ posis = np.array(posis)[:,1:4] + if not sf: +- posis = (posis - orig) @ np.linalg.inv( +- cell +- ) # Convert to scaled coordinates for unscaled coordinates +- if uw and unwrap: +- return posis @ cell # unwrap at the periodic boundaries ++ posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates ++ if uw and unwrap: ++ return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries + else: +- return ( +- posis % 1 +- ) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions ++ return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions + +- +-def get_dumpbox(lines): +- blk, h = _get_block(lines, "BOX BOUNDS") +- bounds = np.zeros([3, 2]) ++def get_dumpbox(lines) : ++ blk, h = _get_block(lines, 'BOX BOUNDS') ++ bounds = np.zeros([3,2]) + tilt = np.zeros([3]) +- load_tilt = "xy xz yz" in h +- for dd in range(3): ++ load_tilt = 'xy xz yz' in h ++ for dd in range(3) : + info = [float(jj) for jj in blk[dd].split()] + bounds[dd][0] = info[0] + bounds[dd][1] = info[1] +- if load_tilt: ++ if load_tilt : + tilt[dd] = info[2] + return bounds, tilt + +- +-def dumpbox2box(bounds, tilt): ++def dumpbox2box(bounds, tilt) : + xy = tilt[0] + xz = tilt[1] + yz = tilt[2] +- xlo = bounds[0][0] - min(0.0, xy, xz, xy + xz) +- xhi = bounds[0][1] - max(0.0, xy, xz, xy + xz) +- ylo = bounds[1][0] - min(0.0, yz) +- yhi = bounds[1][1] - max(0.0, yz) ++ xlo = bounds[0][0] - min(0.0,xy,xz,xy+xz) ++ xhi = bounds[0][1] - max(0.0,xy,xz,xy+xz) ++ ylo = bounds[1][0] - min(0.0,yz) ++ yhi = bounds[1][1] - max(0.0,yz) + zlo = bounds[2][0] + zhi = bounds[2][1] + info = [[xlo, xhi], [ylo, yhi], [zlo, zhi]] + return lmp.lmpbox2box(info, tilt) + +- +-def box2dumpbox(orig, box): ++def box2dumpbox(orig, box) : + lohi, tilt = lmp.box2lmpbox(orig, box) + xy = tilt[0] + xz = tilt[1] + yz = tilt[2] +- bounds = np.zeros([3, 2]) +- bounds[0][0] = lohi[0][0] + min(0.0, xy, xz, xy + xz) +- bounds[0][1] = lohi[0][1] + max(0.0, xy, xz, xy + xz) +- bounds[1][0] = lohi[1][0] + min(0.0, yz) +- bounds[1][1] = lohi[1][1] + max(0.0, yz) ++ bounds = np.zeros([3,2]) ++ bounds[0][0] = lohi[0][0] + min(0.0,xy,xz,xy+xz) ++ bounds[0][1] = lohi[0][1] + max(0.0,xy,xz,xy+xz) ++ bounds[1][0] = lohi[1][0] + min(0.0,yz) ++ bounds[1][1] = lohi[1][1] + max(0.0,yz) + bounds[2][0] = lohi[2][0] + bounds[2][1] = lohi[2][1] + return bounds, tilt + + +-def load_file(fname, begin=0, step=1): ++def load_file(fname, begin = 0, step = 1) : + lines = [] + buff = [] + cc = -1 + with open(fname) as fp: + while True: +- line = fp.readline().rstrip("\n") +- if not line: +- if cc >= begin and (cc - begin) % step == 0: ++ line = fp.readline().rstrip('\n') ++ if not line : ++ if cc >= begin and (cc - begin) % step == 0 : + lines += buff + buff = [] + cc += 1 + return lines +- if "ITEM: TIMESTEP" in line: +- if cc >= begin and (cc - begin) % step == 0: ++ if 'ITEM: TIMESTEP' in line : ++ if cc >= begin and (cc - begin) % step == 0 : + lines += buff + buff = [] + cc += 1 +- if cc >= begin and (cc - begin) % step == 0: ++ if cc >= begin and (cc - begin) % step == 0 : + buff.append(line) + + +-def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): ++def system_data(lines, type_map = None, type_idx_zero = True, unwrap=False) : + array_lines = split_traj(lines) + lines = array_lines[0] + system = {} +- system["atom_numbs"] = get_natoms_vec(lines) +- system["atom_names"] = [] +- if type_map == None: +- for ii in range(len(system["atom_numbs"])): +- system["atom_names"].append("TYPE_%d" % ii) +- else: +- assert len(type_map) >= len(system["atom_numbs"]) +- for ii in range(len(system["atom_numbs"])): +- system["atom_names"].append(type_map[ii]) ++ system['atom_numbs'] = get_natoms_vec(lines) ++ system['atom_names'] = [] ++ if type_map == None : ++ for ii in range(len(system['atom_numbs'])) : ++ system['atom_names'].append('TYPE_%d' % ii) ++ else : ++ assert(len(type_map) >= len(system['atom_numbs'])) ++ for ii in range(len(system['atom_numbs'])) : ++ system['atom_names'].append(type_map[ii]) + bounds, tilt = get_dumpbox(lines) + orig, cell = dumpbox2box(bounds, tilt) +- system["orig"] = np.array(orig) - np.array(orig) +- system["cells"] = [np.array(cell)] +- natoms = sum(system["atom_numbs"]) +- system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero) +- system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] +- for ii in range(1, len(array_lines)): ++ system['orig'] = np.array(orig) - np.array(orig) ++ system['cells'] = [np.array(cell)] ++ natoms = sum(system['atom_numbs']) ++ system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero) ++ system['coords'] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] ++ for ii in range(1, len(array_lines)) : + bounds, tilt = get_dumpbox(array_lines[ii]) + orig, cell = dumpbox2box(bounds, tilt) +- system["cells"].append(cell) +- system["coords"].append( +- safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap) +- ) +- system["cells"] = np.array(system["cells"]) +- system["coords"] = np.array(system["coords"]) ++ system['cells'].append(cell) ++ system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)) ++ system['cells'] = np.array(system['cells']) ++ system['coords'] = np.array(system['coords']) + return system + + +-def split_traj(dump_lines): ++def split_traj(dump_lines) : + marks = [] +- for idx, ii in enumerate(dump_lines): +- if "ITEM: TIMESTEP" in ii: ++ for idx,ii in enumerate(dump_lines) : ++ if 'ITEM: TIMESTEP' in ii : + marks.append(idx) +- if len(marks) == 0: ++ if len(marks) == 0 : + return None +- elif len(marks) == 1: ++ elif len(marks) == 1 : + return [dump_lines] +- else: ++ else : + block_size = marks[1] - marks[0] + ret = [] +- for ii in marks: +- ret.append(dump_lines[ii : ii + block_size]) ++ for ii in marks : ++ ret.append(dump_lines[ii:ii+block_size]) + # for ii in range(len(marks)-1): + # assert(marks[ii+1] - marks[ii] == block_size) + return ret + return None + + +-if __name__ == "__main__": ++if __name__ == '__main__' : + # fname = 'dump.hti' + # lines = open(fname).read().split('\n') + # # print(get_natoms(lines)) +@@ -239,9 +215,9 @@ if __name__ == "__main__": + # print(box) + # np.savetxt('tmp.out', posi - orig, fmt='%.6f') + # print(system_data(lines)) +- lines = load_file("conf_unfold.dump", begin=0, step=1) ++ lines = load_file('conf_unfold.dump', begin = 0, step = 1) + al = split_traj(lines) +- s = system_data(lines, ["O", "H"]) +- # l = np.linalg.norm(s['cells'][1],axis=1) +- # p = s['coords'][0] + l +- # np.savetxt('p',p,fmt='%1.10f') ++ s = system_data(lines,['O','H']) ++ #l = np.linalg.norm(s['cells'][1],axis=1) ++ #p = s['coords'][0] + l ++ #np.savetxt('p',p,fmt='%1.10f') +diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py +index d4bce01..f1356bb 100644 +--- a/dpdata/plugins/lammps.py ++++ b/dpdata/plugins/lammps.py +@@ -9,7 +9,7 @@ class LAMMPSLmpFormat(Format): + @Format.post("shift_orig_zero") + def from_system(self, file_name, type_map=None, **kwargs): + with open(file_name) as fp: +- lines = [line.rstrip("\n") for line in fp] ++ lines = [line.rstrip('\n') for line in fp] + return dpdata.lammps.lmp.to_system_data(lines, type_map) + + def to_system(self, data, file_name, frame_idx=0, **kwargs): +@@ -25,9 +25,9 @@ class LAMMPSLmpFormat(Format): + frame_idx : int + The index of the frame to dump + """ +- assert frame_idx < len(data["coords"]) ++ assert(frame_idx < len(data['coords'])) + w_str = dpdata.lammps.lmp.from_system_data(data, frame_idx) +- with open(file_name, "w") as fp: ++ with open(file_name, 'w') as fp: + fp.write(w_str) + + +@@ -35,8 +35,11 @@ class LAMMPSLmpFormat(Format): + @Format.register("lammps/dump") + class LAMMPSDumpFormat(Format): + @Format.post("shift_orig_zero") +- def from_system( +- self, file_name, type_map=None, begin=0, step=1, unwrap=False, **kwargs +- ): ++ def from_system(self, ++ file_name, ++ type_map=None, ++ begin=0, ++ step=1, ++ unwrap=False, ++ **kwargs): + lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step) +- return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) ++ return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) +diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py +index f120183..23bab64 100644 +--- a/tests/poscars/poscar_ref_oh.py ++++ b/tests/poscars/poscar_ref_oh.py +@@ -1,50 +1,38 @@ + import numpy as np + ++class TestPOSCARoh : + +-class TestPOSCARoh: + def test_atom_numbs(self): +- self.assertEqual(self.system.data["atom_numbs"], [1, 1]) ++ self.assertEqual(self.system.data['atom_numbs'], [1,1]) + + def test_atom_names(self): +- self.assertEqual(self.system.data["atom_names"], ["O", "H"]) ++ self.assertEqual(self.system.data['atom_names'], ['O','H']) + + def test_atom_types(self): +- self.assertEqual(self.system.data["atom_types"][0], 0) +- self.assertEqual(self.system.data["atom_types"][1], 1) ++ self.assertEqual(self.system.data['atom_types'][0], 0) ++ self.assertEqual(self.system.data['atom_types'][1], 1) + + def test_orig(self): +- for d0 in range(3): +- self.assertEqual(self.system.data["orig"][d0], 0) ++ for d0 in range(3) : ++ self.assertEqual(self.system.data['orig'][d0], 0) + + def test_cell(self): +- ovito_cell = np.array( +- [ +- [2.5243712, 0.0000000, 0.0000000], +- [1.2621856, 2.0430257, 0.0000000], +- [1.2874292, 0.7485898, 2.2254033], +- ] +- ) +- for ii in range(3): +- for jj in range(3): +- self.assertAlmostEqual( +- self.system.data["cells"][0][ii][jj], +- ovito_cell[ii][jj], +- places=6, +- msg="cell[%d][%d] failed" % (ii, jj), +- ) ++ ovito_cell = np.array([[2.5243712, 0.0000000, 0.0000000], ++ [1.2621856, 2.0430257, 0.0000000], ++ [1.2874292, 0.7485898, 2.2254033]]) ++ for ii in range(3) : ++ for jj in range(3) : ++ self.assertAlmostEqual(self.system.data['cells'][0][ii][jj], ++ ovito_cell[ii][jj], ++ places = 6, ++ msg = 'cell[%d][%d] failed' % (ii,jj)) + +- def test_frame(self): +- if hasattr(self, "unwrap") and self.unwrap is True: +- ovito_posis = np.array( +- [[5.0739861, 2.7916155, 2.2254033], [6.3361717, 3.4934183, 2.7767918]] +- ) +- else: +- ovito_posis = np.array([[0, 0, 0], [1.2621856, 0.7018028, 0.5513885]]) +- for ii in range(2): +- for jj in range(3): +- self.assertAlmostEqual( +- self.system.data["coords"][0][ii][jj], +- ovito_posis[ii][jj], +- places=6, +- msg="posis[%d][%d] failed" % (ii, jj), +- ) ++ def test_frame(self): ++ ovito_posis = np.array([[0, 0, 0], ++ [1.2621856, 0.7018028, 0.5513885]]) ++ for ii in range(2) : ++ for jj in range(3) : ++ self.assertAlmostEqual(self.system.data['coords'][0][ii][jj], ++ ovito_posis[ii][jj], ++ places = 6, ++ msg = 'posis[%d][%d] failed' % (ii,jj)) +diff --git a/tests/test_lammps_dump_unfold.py b/tests/test_lammps_dump_unfold.py +index b4c6119..4ffd5f7 100644 +--- a/tests/test_lammps_dump_unfold.py ++++ b/tests/test_lammps_dump_unfold.py +@@ -1,38 +1,26 @@ +-from inspect import unwrap + import os + import numpy as np + import unittest + from context import dpdata +-from poscars.poscar_ref_oh import TestPOSCARoh +- ++from poscars.poscar_ref_oh import TestPOSCARoh + + class TestDump(unittest.TestCase, TestPOSCARoh): +- def setUp(self): +- self.system = dpdata.System( +- os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] +- ) +- +- ++ ++ def setUp(self): ++ self.system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), ++ type_map = ['O', 'H']) ++ + class TestDump2(unittest.TestCase, TestPOSCARoh): +- def setUp(self): +- self.tmp_system = dpdata.System( +- os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] +- ) ++ ++ def setUp(self): ++ self.tmp_system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), ++ type_map = ['O', 'H']) + self.system = self.tmp_system.sub_system([1]) + +- def test_nframes(self): ++ def test_nframes (self) : + self.assertEqual(self.tmp_system.get_nframes(), 2) +- +- +-class TestDumpUnwrap(unittest.TestCase, TestPOSCARoh): +- def setUp(self): +- self.unwrap = True +- self.system = dpdata.System( +- os.path.join("poscars", "conf_unfold.dump"), +- type_map=["O", "H"], +- unwrap=self.unwrap, +- ) +- +- +-if __name__ == "__main__": ++ ++ ++if __name__ == '__main__': + unittest.main() ++ diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py index f120183ed..9b12c1511 100644 --- a/tests/poscars/poscar_ref_oh.py +++ b/tests/poscars/poscar_ref_oh.py @@ -1,50 +1,42 @@ import numpy as np +class TestPOSCARoh : -class TestPOSCARoh: def test_atom_numbs(self): - self.assertEqual(self.system.data["atom_numbs"], [1, 1]) + self.assertEqual(self.system.data['atom_numbs'], [1,1]) def test_atom_names(self): - self.assertEqual(self.system.data["atom_names"], ["O", "H"]) + self.assertEqual(self.system.data['atom_names'], ['O','H']) def test_atom_types(self): - self.assertEqual(self.system.data["atom_types"][0], 0) - self.assertEqual(self.system.data["atom_types"][1], 1) + self.assertEqual(self.system.data['atom_types'][0], 0) + self.assertEqual(self.system.data['atom_types'][1], 1) def test_orig(self): - for d0 in range(3): - self.assertEqual(self.system.data["orig"][d0], 0) + for d0 in range(3) : + self.assertEqual(self.system.data['orig'][d0], 0) def test_cell(self): - ovito_cell = np.array( - [ - [2.5243712, 0.0000000, 0.0000000], - [1.2621856, 2.0430257, 0.0000000], - [1.2874292, 0.7485898, 2.2254033], - ] - ) - for ii in range(3): - for jj in range(3): - self.assertAlmostEqual( - self.system.data["cells"][0][ii][jj], - ovito_cell[ii][jj], - places=6, - msg="cell[%d][%d] failed" % (ii, jj), - ) + ovito_cell = np.array([[2.5243712, 0.0000000, 0.0000000], + [1.2621856, 2.0430257, 0.0000000], + [1.2874292, 0.7485898, 2.2254033]]) + for ii in range(3) : + for jj in range(3) : + self.assertAlmostEqual(self.system.data['cells'][0][ii][jj], + ovito_cell[ii][jj], + places = 6, + msg = 'cell[%d][%d] failed' % (ii,jj)) def test_frame(self): if hasattr(self, "unwrap") and self.unwrap is True: - ovito_posis = np.array( - [[5.0739861, 2.7916155, 2.2254033], [6.3361717, 3.4934183, 2.7767918]] - ) + ovito_posis = np.array([[5.0739861, 2.7916155, 2.2254033], + [6.3361717, 3.4934183, 2.7767918]]) else: - ovito_posis = np.array([[0, 0, 0], [1.2621856, 0.7018028, 0.5513885]]) - for ii in range(2): - for jj in range(3): - self.assertAlmostEqual( - self.system.data["coords"][0][ii][jj], - ovito_posis[ii][jj], - places=6, - msg="posis[%d][%d] failed" % (ii, jj), - ) + ovito_posis = np.array([[0, 0, 0], + [1.2621856, 0.7018028, 0.5513885]]) + for ii in range(2) : + for jj in range(3) : + self.assertAlmostEqual(self.system.data['coords'][0][ii][jj], + ovito_posis[ii][jj], + places = 6, + msg = 'posis[%d][%d] failed' % (ii,jj)) diff --git a/tests/test_lammps_dump_unfold.py b/tests/test_lammps_dump_unfold.py index b4c61195d..68aa2c8b7 100644 --- a/tests/test_lammps_dump_unfold.py +++ b/tests/test_lammps_dump_unfold.py @@ -1,38 +1,36 @@ -from inspect import unwrap import os import numpy as np import unittest from context import dpdata -from poscars.poscar_ref_oh import TestPOSCARoh - +from poscars.poscar_ref_oh import TestPOSCARoh class TestDump(unittest.TestCase, TestPOSCARoh): - def setUp(self): - self.system = dpdata.System( - os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] - ) - - + + def setUp(self): + self.system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), + type_map = ['O', 'H']) + class TestDump2(unittest.TestCase, TestPOSCARoh): - def setUp(self): - self.tmp_system = dpdata.System( - os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] - ) + + def setUp(self): + self.tmp_system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), + type_map = ['O', 'H']) self.system = self.tmp_system.sub_system([1]) - def test_nframes(self): + def test_nframes (self) : self.assertEqual(self.tmp_system.get_nframes(), 2) - - + + class TestDumpUnwrap(unittest.TestCase, TestPOSCARoh): def setUp(self): self.unwrap = True self.system = dpdata.System( - os.path.join("poscars", "conf_unfold.dump"), - type_map=["O", "H"], + os.path.join('poscars', 'conf_unfold.dump'), + type_map=['O', 'H'], unwrap=self.unwrap, ) -if __name__ == "__main__": +if __name__ == '__main__': unittest.main() + From affe5be8b2593b4879fcb9202d84d1cfd2526891 Mon Sep 17 00:00:00 2001 From: zezhong-zhang Date: Thu, 11 Aug 2022 14:07:25 +0200 Subject: [PATCH 4/6] clean up --- my.patch | 531 ------------------------------------------------------- 1 file changed, 531 deletions(-) delete mode 100644 my.patch diff --git a/my.patch b/my.patch deleted file mode 100644 index b203a16d1..000000000 --- a/my.patch +++ /dev/null @@ -1,531 +0,0 @@ -diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py -index e39e72a..bdc84dd 100644 ---- a/dpdata/lammps/dump.py -+++ b/dpdata/lammps/dump.py -@@ -2,230 +2,206 @@ - - import os, sys - import numpy as np -- - lib_path = os.path.dirname(os.path.realpath(__file__)) - sys.path.append(lib_path) - import lmp - -- --def _get_block(lines, key): -- for idx in range(len(lines)): -- if ("ITEM: " + key) in lines[idx]: -+def _get_block (lines, key) : -+ for idx in range(len(lines)) : -+ if ('ITEM: ' + key) in lines[idx] : - break - idx_s = idx + 1 -- for idx in range(idx_s, len(lines)): -- if ("ITEM: ") in lines[idx]: -+ for idx in range(idx_s, len(lines)) : -+ if ('ITEM: ') in lines[idx] : - break - idx_e = idx -- if idx_e == len(lines) - 1: -+ if idx_e == len(lines)-1 : - idx_e += 1 -- return lines[idx_s:idx_e], lines[idx_s - 1] -+ return lines[idx_s:idx_e], lines[idx_s-1] - -- --def get_atype(lines, type_idx_zero=False): -- blk, head = _get_block(lines, "ATOMS") -+def get_atype(lines, type_idx_zero = False) : -+ blk, head = _get_block(lines, 'ATOMS') - keys = head.split() -- id_idx = keys.index("id") - 2 -- tidx = keys.index("type") - 2 -+ id_idx = keys.index('id') - 2 -+ tidx = keys.index('type') - 2 - atype = [] -- for ii in blk: -+ for ii in blk : - atype.append([int(ii.split()[id_idx]), int(ii.split()[tidx])]) - atype.sort() -- atype = np.array(atype, dtype=int) -- if type_idx_zero: -- return atype[:, 1] - 1 -- else: -- return atype[:, 1] -- -- --def get_natoms(lines): -- blk, head = _get_block(lines, "NUMBER OF ATOMS") -+ atype = np.array(atype, dtype = int) -+ if type_idx_zero : -+ return atype[:,1] - 1 -+ else : -+ return atype[:,1] -+ -+def get_natoms(lines) : -+ blk, head = _get_block(lines, 'NUMBER OF ATOMS') - return int(blk[0]) - -- --def get_natomtypes(lines): -+def get_natomtypes(lines) : - atype = get_atype(lines) - return max(atype) - -- --def get_natoms_vec(lines): -+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) -+ for ii in range(natomtypes) : -+ natoms_vec.append(sum(atype == ii+1)) -+ assert (sum(natoms_vec) == get_natoms(lines)) - return natoms_vec - -- - def get_coordtype_and_scalefactor(keys): - # 4 types in total,with different scaling factor -- key_pc = ["x", "y", "z"] # plain cartesian, sf = 1 -- key_uc = ["xu", "yu", "zu"] # unwraped cartesian, sf = 1 -- key_s = ["xs", "ys", "zs"] # scaled by lattice parameter, sf = lattice parameter -- key_su = ["xsu", "ysu", "zsu"] # scaled and unfolded,sf = lattice parameter -- lmp_coor_type = [key_pc, key_uc, key_s, key_su] -- sf = [0, 0, 1, 1] -- uw = [0, 1, 0, 1] # unwraped or not -+ key_pc=['x','y','z'] # plain cartesian, sf = 1 -+ key_uc=['xu','yu','zu'] # unwraped cartesian, sf = 1 -+ key_s=['xs','ys','zs'] # scaled by lattice parameter, sf = lattice parameter -+ key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter -+ lmp_coor_type = [key_pc,key_uc,key_s,key_su] -+ sf = [0,0,1,1] -+ uw = [0, 1, 0, 1] # unwraped or not - for k in range(4): - if all(i in keys for i in lmp_coor_type[k]): -- return lmp_coor_type[k], sf[k], uw[k] -- -+ return lmp_coor_type[k], sf[k], uw[k] - --def safe_get_posi(lines, cell, orig=np.zeros(3), unwrap=False): -- blk, head = _get_block(lines, "ATOMS") -+def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : -+ blk, head = _get_block(lines, 'ATOMS') - keys = head.split() - coord_tp_and_sf = get_coordtype_and_scalefactor(keys) -- assert coord_tp_and_sf is not None, "Dump file does not contain atomic coordinates!" -- coordtype, sf, uw = coord_tp_and_sf -- id_idx = keys.index("id") - 2 -- xidx = keys.index(coordtype[0]) - 2 -- yidx = keys.index(coordtype[1]) - 2 -- zidx = keys.index(coordtype[2]) - 2 -+ assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!' -+ coordtype, sf, uw = coord_tp_and_sf -+ id_idx = keys.index('id') - 2 -+ xidx = keys.index(coordtype[0])-2 -+ yidx = keys.index(coordtype[1])-2 -+ zidx = keys.index(coordtype[2])-2 - sel = (xidx, yidx, zidx) - posis = [] -- for ii in blk: -+ for ii in blk : - words = ii.split() -- posis.append( -- [ -- float(words[id_idx]), -- float(words[xidx]), -- float(words[yidx]), -- float(words[zidx]), -- ] -- ) -+ posis.append([float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])]) - posis.sort() -- posis = np.array(posis)[:, 1:4] -+ posis = np.array(posis)[:,1:4] - if not sf: -- posis = (posis - orig) @ np.linalg.inv( -- cell -- ) # Convert to scaled coordinates for unscaled coordinates -- if uw and unwrap: -- return posis @ cell # unwrap at the periodic boundaries -+ posis = (posis - orig) @ np.linalg.inv(cell) # Convert to scaled coordinates for unscaled coordinates -+ if uw and unwrap: -+ return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries - else: -- return ( -- posis % 1 -- ) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions -+ return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions - -- --def get_dumpbox(lines): -- blk, h = _get_block(lines, "BOX BOUNDS") -- bounds = np.zeros([3, 2]) -+def get_dumpbox(lines) : -+ blk, h = _get_block(lines, 'BOX BOUNDS') -+ bounds = np.zeros([3,2]) - tilt = np.zeros([3]) -- load_tilt = "xy xz yz" in h -- for dd in range(3): -+ load_tilt = 'xy xz yz' in h -+ for dd in range(3) : - info = [float(jj) for jj in blk[dd].split()] - bounds[dd][0] = info[0] - bounds[dd][1] = info[1] -- if load_tilt: -+ if load_tilt : - tilt[dd] = info[2] - return bounds, tilt - -- --def dumpbox2box(bounds, tilt): -+def dumpbox2box(bounds, tilt) : - xy = tilt[0] - xz = tilt[1] - yz = tilt[2] -- xlo = bounds[0][0] - min(0.0, xy, xz, xy + xz) -- xhi = bounds[0][1] - max(0.0, xy, xz, xy + xz) -- ylo = bounds[1][0] - min(0.0, yz) -- yhi = bounds[1][1] - max(0.0, yz) -+ xlo = bounds[0][0] - min(0.0,xy,xz,xy+xz) -+ xhi = bounds[0][1] - max(0.0,xy,xz,xy+xz) -+ ylo = bounds[1][0] - min(0.0,yz) -+ yhi = bounds[1][1] - max(0.0,yz) - zlo = bounds[2][0] - zhi = bounds[2][1] - info = [[xlo, xhi], [ylo, yhi], [zlo, zhi]] - return lmp.lmpbox2box(info, tilt) - -- --def box2dumpbox(orig, box): -+def box2dumpbox(orig, box) : - lohi, tilt = lmp.box2lmpbox(orig, box) - xy = tilt[0] - xz = tilt[1] - yz = tilt[2] -- bounds = np.zeros([3, 2]) -- bounds[0][0] = lohi[0][0] + min(0.0, xy, xz, xy + xz) -- bounds[0][1] = lohi[0][1] + max(0.0, xy, xz, xy + xz) -- bounds[1][0] = lohi[1][0] + min(0.0, yz) -- bounds[1][1] = lohi[1][1] + max(0.0, yz) -+ bounds = np.zeros([3,2]) -+ bounds[0][0] = lohi[0][0] + min(0.0,xy,xz,xy+xz) -+ bounds[0][1] = lohi[0][1] + max(0.0,xy,xz,xy+xz) -+ bounds[1][0] = lohi[1][0] + min(0.0,yz) -+ bounds[1][1] = lohi[1][1] + max(0.0,yz) - bounds[2][0] = lohi[2][0] - bounds[2][1] = lohi[2][1] - return bounds, tilt - - --def load_file(fname, begin=0, step=1): -+def load_file(fname, begin = 0, step = 1) : - lines = [] - buff = [] - cc = -1 - with open(fname) as fp: - while True: -- line = fp.readline().rstrip("\n") -- if not line: -- if cc >= begin and (cc - begin) % step == 0: -+ line = fp.readline().rstrip('\n') -+ if not line : -+ if cc >= begin and (cc - begin) % step == 0 : - lines += buff - buff = [] - cc += 1 - return lines -- if "ITEM: TIMESTEP" in line: -- if cc >= begin and (cc - begin) % step == 0: -+ if 'ITEM: TIMESTEP' in line : -+ if cc >= begin and (cc - begin) % step == 0 : - lines += buff - buff = [] - cc += 1 -- if cc >= begin and (cc - begin) % step == 0: -+ if cc >= begin and (cc - begin) % step == 0 : - buff.append(line) - - --def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False): -+def system_data(lines, type_map = None, type_idx_zero = True, unwrap=False) : - array_lines = split_traj(lines) - lines = array_lines[0] - system = {} -- system["atom_numbs"] = get_natoms_vec(lines) -- system["atom_names"] = [] -- if type_map == None: -- for ii in range(len(system["atom_numbs"])): -- system["atom_names"].append("TYPE_%d" % ii) -- else: -- assert len(type_map) >= len(system["atom_numbs"]) -- for ii in range(len(system["atom_numbs"])): -- system["atom_names"].append(type_map[ii]) -+ system['atom_numbs'] = get_natoms_vec(lines) -+ system['atom_names'] = [] -+ if type_map == None : -+ for ii in range(len(system['atom_numbs'])) : -+ system['atom_names'].append('TYPE_%d' % ii) -+ else : -+ assert(len(type_map) >= len(system['atom_numbs'])) -+ for ii in range(len(system['atom_numbs'])) : -+ system['atom_names'].append(type_map[ii]) - bounds, tilt = get_dumpbox(lines) - orig, cell = dumpbox2box(bounds, tilt) -- system["orig"] = np.array(orig) - np.array(orig) -- system["cells"] = [np.array(cell)] -- natoms = sum(system["atom_numbs"]) -- system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero) -- system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] -- for ii in range(1, len(array_lines)): -+ system['orig'] = np.array(orig) - np.array(orig) -+ system['cells'] = [np.array(cell)] -+ natoms = sum(system['atom_numbs']) -+ system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero) -+ system['coords'] = [safe_get_posi(lines, cell, np.array(orig), unwrap)] -+ for ii in range(1, len(array_lines)) : - bounds, tilt = get_dumpbox(array_lines[ii]) - orig, cell = dumpbox2box(bounds, tilt) -- system["cells"].append(cell) -- system["coords"].append( -- safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap) -- ) -- system["cells"] = np.array(system["cells"]) -- system["coords"] = np.array(system["coords"]) -+ system['cells'].append(cell) -+ system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)) -+ system['cells'] = np.array(system['cells']) -+ system['coords'] = np.array(system['coords']) - return system - - --def split_traj(dump_lines): -+def split_traj(dump_lines) : - marks = [] -- for idx, ii in enumerate(dump_lines): -- if "ITEM: TIMESTEP" in ii: -+ for idx,ii in enumerate(dump_lines) : -+ if 'ITEM: TIMESTEP' in ii : - marks.append(idx) -- if len(marks) == 0: -+ if len(marks) == 0 : - return None -- elif len(marks) == 1: -+ elif len(marks) == 1 : - return [dump_lines] -- else: -+ else : - block_size = marks[1] - marks[0] - ret = [] -- for ii in marks: -- ret.append(dump_lines[ii : ii + block_size]) -+ for ii in marks : -+ ret.append(dump_lines[ii:ii+block_size]) - # for ii in range(len(marks)-1): - # assert(marks[ii+1] - marks[ii] == block_size) - return ret - return None - - --if __name__ == "__main__": -+if __name__ == '__main__' : - # fname = 'dump.hti' - # lines = open(fname).read().split('\n') - # # print(get_natoms(lines)) -@@ -239,9 +215,9 @@ if __name__ == "__main__": - # print(box) - # np.savetxt('tmp.out', posi - orig, fmt='%.6f') - # print(system_data(lines)) -- lines = load_file("conf_unfold.dump", begin=0, step=1) -+ lines = load_file('conf_unfold.dump', begin = 0, step = 1) - al = split_traj(lines) -- s = system_data(lines, ["O", "H"]) -- # l = np.linalg.norm(s['cells'][1],axis=1) -- # p = s['coords'][0] + l -- # np.savetxt('p',p,fmt='%1.10f') -+ s = system_data(lines,['O','H']) -+ #l = np.linalg.norm(s['cells'][1],axis=1) -+ #p = s['coords'][0] + l -+ #np.savetxt('p',p,fmt='%1.10f') -diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py -index d4bce01..f1356bb 100644 ---- a/dpdata/plugins/lammps.py -+++ b/dpdata/plugins/lammps.py -@@ -9,7 +9,7 @@ class LAMMPSLmpFormat(Format): - @Format.post("shift_orig_zero") - def from_system(self, file_name, type_map=None, **kwargs): - with open(file_name) as fp: -- lines = [line.rstrip("\n") for line in fp] -+ lines = [line.rstrip('\n') for line in fp] - return dpdata.lammps.lmp.to_system_data(lines, type_map) - - def to_system(self, data, file_name, frame_idx=0, **kwargs): -@@ -25,9 +25,9 @@ class LAMMPSLmpFormat(Format): - frame_idx : int - The index of the frame to dump - """ -- assert frame_idx < len(data["coords"]) -+ assert(frame_idx < len(data['coords'])) - w_str = dpdata.lammps.lmp.from_system_data(data, frame_idx) -- with open(file_name, "w") as fp: -+ with open(file_name, 'w') as fp: - fp.write(w_str) - - -@@ -35,8 +35,11 @@ class LAMMPSLmpFormat(Format): - @Format.register("lammps/dump") - class LAMMPSDumpFormat(Format): - @Format.post("shift_orig_zero") -- def from_system( -- self, file_name, type_map=None, begin=0, step=1, unwrap=False, **kwargs -- ): -+ def from_system(self, -+ file_name, -+ type_map=None, -+ begin=0, -+ step=1, -+ unwrap=False, -+ **kwargs): - lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step) -- return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) -+ return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap) -diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py -index f120183..23bab64 100644 ---- a/tests/poscars/poscar_ref_oh.py -+++ b/tests/poscars/poscar_ref_oh.py -@@ -1,50 +1,38 @@ - import numpy as np - -+class TestPOSCARoh : - --class TestPOSCARoh: - def test_atom_numbs(self): -- self.assertEqual(self.system.data["atom_numbs"], [1, 1]) -+ self.assertEqual(self.system.data['atom_numbs'], [1,1]) - - def test_atom_names(self): -- self.assertEqual(self.system.data["atom_names"], ["O", "H"]) -+ self.assertEqual(self.system.data['atom_names'], ['O','H']) - - def test_atom_types(self): -- self.assertEqual(self.system.data["atom_types"][0], 0) -- self.assertEqual(self.system.data["atom_types"][1], 1) -+ self.assertEqual(self.system.data['atom_types'][0], 0) -+ self.assertEqual(self.system.data['atom_types'][1], 1) - - def test_orig(self): -- for d0 in range(3): -- self.assertEqual(self.system.data["orig"][d0], 0) -+ for d0 in range(3) : -+ self.assertEqual(self.system.data['orig'][d0], 0) - - def test_cell(self): -- ovito_cell = np.array( -- [ -- [2.5243712, 0.0000000, 0.0000000], -- [1.2621856, 2.0430257, 0.0000000], -- [1.2874292, 0.7485898, 2.2254033], -- ] -- ) -- for ii in range(3): -- for jj in range(3): -- self.assertAlmostEqual( -- self.system.data["cells"][0][ii][jj], -- ovito_cell[ii][jj], -- places=6, -- msg="cell[%d][%d] failed" % (ii, jj), -- ) -+ ovito_cell = np.array([[2.5243712, 0.0000000, 0.0000000], -+ [1.2621856, 2.0430257, 0.0000000], -+ [1.2874292, 0.7485898, 2.2254033]]) -+ for ii in range(3) : -+ for jj in range(3) : -+ self.assertAlmostEqual(self.system.data['cells'][0][ii][jj], -+ ovito_cell[ii][jj], -+ places = 6, -+ msg = 'cell[%d][%d] failed' % (ii,jj)) - -- def test_frame(self): -- if hasattr(self, "unwrap") and self.unwrap is True: -- ovito_posis = np.array( -- [[5.0739861, 2.7916155, 2.2254033], [6.3361717, 3.4934183, 2.7767918]] -- ) -- else: -- ovito_posis = np.array([[0, 0, 0], [1.2621856, 0.7018028, 0.5513885]]) -- for ii in range(2): -- for jj in range(3): -- self.assertAlmostEqual( -- self.system.data["coords"][0][ii][jj], -- ovito_posis[ii][jj], -- places=6, -- msg="posis[%d][%d] failed" % (ii, jj), -- ) -+ def test_frame(self): -+ ovito_posis = np.array([[0, 0, 0], -+ [1.2621856, 0.7018028, 0.5513885]]) -+ for ii in range(2) : -+ for jj in range(3) : -+ self.assertAlmostEqual(self.system.data['coords'][0][ii][jj], -+ ovito_posis[ii][jj], -+ places = 6, -+ msg = 'posis[%d][%d] failed' % (ii,jj)) -diff --git a/tests/test_lammps_dump_unfold.py b/tests/test_lammps_dump_unfold.py -index b4c6119..4ffd5f7 100644 ---- a/tests/test_lammps_dump_unfold.py -+++ b/tests/test_lammps_dump_unfold.py -@@ -1,38 +1,26 @@ --from inspect import unwrap - import os - import numpy as np - import unittest - from context import dpdata --from poscars.poscar_ref_oh import TestPOSCARoh -- -+from poscars.poscar_ref_oh import TestPOSCARoh - - class TestDump(unittest.TestCase, TestPOSCARoh): -- def setUp(self): -- self.system = dpdata.System( -- os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] -- ) -- -- -+ -+ def setUp(self): -+ self.system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), -+ type_map = ['O', 'H']) -+ - class TestDump2(unittest.TestCase, TestPOSCARoh): -- def setUp(self): -- self.tmp_system = dpdata.System( -- os.path.join("poscars", "conf_unfold.dump"), type_map=["O", "H"] -- ) -+ -+ def setUp(self): -+ self.tmp_system = dpdata.System(os.path.join('poscars', 'conf_unfold.dump'), -+ type_map = ['O', 'H']) - self.system = self.tmp_system.sub_system([1]) - -- def test_nframes(self): -+ def test_nframes (self) : - self.assertEqual(self.tmp_system.get_nframes(), 2) -- -- --class TestDumpUnwrap(unittest.TestCase, TestPOSCARoh): -- def setUp(self): -- self.unwrap = True -- self.system = dpdata.System( -- os.path.join("poscars", "conf_unfold.dump"), -- type_map=["O", "H"], -- unwrap=self.unwrap, -- ) -- -- --if __name__ == "__main__": -+ -+ -+if __name__ == '__main__': - unittest.main() -+ From afbf7d7a6e4fd548d0837b46ce99535857b6603e Mon Sep 17 00:00:00 2001 From: zezhong-zhang Date: Sat, 13 Aug 2022 12:11:42 +0200 Subject: [PATCH 5/6] add warning for unwrap --- dpdata/lammps/dump.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index 9de734e26..50b5a54fe 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -5,6 +5,9 @@ lib_path = os.path.dirname(os.path.realpath(__file__)) sys.path.append(lib_path) import lmp +import warnings +warnings.simplefilter('once', UserWarning) + def _get_block (lines, key) : for idx in range(len(lines)) : @@ -59,7 +62,7 @@ def get_coordtype_and_scalefactor(keys): key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter lmp_coor_type = [key_pc,key_uc,key_s,key_su] sf = [0,0,1,1] - uw = [0, 1, 0, 1] # unwraped or not + uw = [0,1,0,1] # unwraped or not for k in range(4): if all(i in keys for i in lmp_coor_type[k]): return lmp_coor_type[k], sf[k], uw[k] @@ -86,6 +89,8 @@ def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : if uw and unwrap: return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries else: + if uw and not unwrap: + warnings.warn('Your dump file contains unwrapped coordinates, but you did not specify unwrapping (unwrap = True). The default is wrapping at periodic boundaries (unwrap = False).\n') return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions def get_dumpbox(lines) : From 4a0c73d11501028674fa858896e6d507fd1ff320 Mon Sep 17 00:00:00 2001 From: zezhong-zhang Date: Mon, 15 Aug 2022 10:14:09 +0200 Subject: [PATCH 6/6] add UnwrapWarning class --- dpdata/lammps/dump.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py index 50b5a54fe..e40dc400b 100644 --- a/dpdata/lammps/dump.py +++ b/dpdata/lammps/dump.py @@ -6,7 +6,9 @@ sys.path.append(lib_path) import lmp import warnings -warnings.simplefilter('once', UserWarning) +class UnwrapWarning(UserWarning): + pass +warnings.simplefilter('once', UnwrapWarning) def _get_block (lines, key) : @@ -90,7 +92,7 @@ def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) : return posis @ cell # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries else: if uw and not unwrap: - warnings.warn('Your dump file contains unwrapped coordinates, but you did not specify unwrapping (unwrap = True). The default is wrapping at periodic boundaries (unwrap = False).\n') + warnings.warn(message='Your dump file contains unwrapped coordinates, but you did not specify unwrapping (unwrap = True). The default is wrapping at periodic boundaries (unwrap = False).\n',category=UnwrapWarning) return (posis % 1) @ cell # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions def get_dumpbox(lines) :