1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
from ase.atoms import Atoms
from ase.quaternions import Quaternions
from ase.calculators.singlepoint import SinglePointCalculator
from ase.parallel import paropen
def read_lammps_dump(fileobj, index=-1, order=True, atomsobj=Atoms):
"""Method which reads a LAMMPS dump file.
order: Order the particles according to their id. Might be faster to
switch it off.
"""
if isinstance(fileobj, str):
f = paropen(fileobj)
else:
f = fileobj
# load everything into memory
lines = f.readlines()
natoms = 0
images = []
while len(lines) > natoms:
line = lines.pop(0)
if 'ITEM: TIMESTEP' in line:
lo = []
hi = []
tilt = []
id = []
types = []
positions = []
scaled_positions = []
velocities = []
forces = []
quaternions = []
if 'ITEM: NUMBER OF ATOMS' in line:
line = lines.pop(0)
natoms = int(line.split()[0])
if 'ITEM: BOX BOUNDS' in line:
# save labels behind "ITEM: BOX BOUNDS" in
# triclinic case (>=lammps-7Jul09)
tilt_items = line.split()[3:]
for i in range(3):
line = lines.pop(0)
fields = line.split()
lo.append(float(fields[0]))
hi.append(float(fields[1]))
if (len(fields) >= 3):
tilt.append(float(fields[2]))
# determine cell tilt (triclinic case!)
if (len(tilt) >= 3):
# for >=lammps-7Jul09 use labels behind
# "ITEM: BOX BOUNDS" to assign tilt (vector) elements ...
if (len(tilt_items) >= 3):
xy = tilt[tilt_items.index('xy')]
xz = tilt[tilt_items.index('xz')]
yz = tilt[tilt_items.index('yz')]
# ... otherwise assume default order in 3rd column
# (if the latter was present)
else:
xy = tilt[0]
xz = tilt[1]
yz = tilt[2]
else:
xy = xz = yz = 0
xhilo = (hi[0] - lo[0]) - (xy**2)**0.5 - (xz**2)**0.5
yhilo = (hi[1] - lo[1]) - (yz**2)**0.5
zhilo = (hi[2] - lo[2])
if xy < 0:
if xz < 0:
celldispx = lo[0] - xy - xz
else:
celldispx = lo[0] - xy
else:
celldispx = lo[0]
celldispy = lo[1]
celldispz = lo[2]
cell = [[xhilo, 0, 0], [xy, yhilo, 0], [xz, yz, zhilo]]
celldisp = [[celldispx, celldispy, celldispz]]
def add_quantity(fields, var, labels):
for label in labels:
if label not in atom_attributes:
return
var.append([float(fields[atom_attributes[label]])
for label in labels])
if 'ITEM: ATOMS' in line:
# (reliably) identify values by labels behind
# "ITEM: ATOMS" - requires >=lammps-7Jul09
# create corresponding index dictionary before
# iterating over atoms to (hopefully) speed up lookups...
atom_attributes = {}
for (i, x) in enumerate(line.split()[2:]):
atom_attributes[x] = i
for n in range(natoms):
line = lines.pop(0)
fields = line.split()
id.append(int(fields[atom_attributes['id']]))
types.append(int(fields[atom_attributes['type']]))
add_quantity(fields, positions, ['x', 'y', 'z'])
add_quantity(fields, scaled_positions, ['xs', 'ys', 'zs'])
add_quantity(fields, velocities, ['vx', 'vy', 'vz'])
add_quantity(fields, forces, ['fx', 'fy', 'fz'])
add_quantity(fields, quaternions, ['c_q[1]', 'c_q[2]',
'c_q[3]', 'c_q[4]'])
if order:
def reorder(inlist):
if not len(inlist):
return inlist
outlist = [None] * len(id)
for i, v in zip(id, inlist):
outlist[i - 1] = v
return outlist
types = reorder(types)
positions = reorder(positions)
scaled_positions = reorder(scaled_positions)
velocities = reorder(velocities)
forces = reorder(forces)
quaternions = reorder(quaternions)
if len(quaternions):
images.append(Quaternions(symbols=types,
positions=positions,
cell=cell, celldisp=celldisp,
quaternions=quaternions))
elif len(positions):
images.append(atomsobj(
symbols=types, positions=positions,
celldisp=celldisp, cell=cell))
elif len(scaled_positions):
images.append(atomsobj(
symbols=types, scaled_positions=scaled_positions,
celldisp=celldisp, cell=cell))
if len(velocities):
images[-1].set_velocities(velocities)
if len(forces):
calculator = SinglePointCalculator(images[-1],
energy=0.0, forces=forces)
images[-1].set_calculator(calculator)
return images[index]
|
