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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
import numpy as np
import xml.etree.ElementTree as ET
from xml.dom import minidom
from ase.io.xsd import SetChild, _write_xsd_html
from ase import Atoms
_image_header = ' ' * 74 + '0.0000\n!DATE Jan 01 00:00:00 2000\n'
_image_footer = 'end\nend\n'
def _get_atom_str(an, xyz):
s = '{:<5}'.format(an)
s += '{:>15.9f}{:>15.9f}{:>15.9f}'.format(xyz[0], xyz[1], xyz[2])
s += ' XXXX 1 xx '
s += '{:<2}'.format(an)
s += ' 0.000\n'
return s
def write_xtd(filename, images, connectivity=None, moviespeed=10):
"""Takes Atoms object, and write materials studio file
atoms: Atoms object
filename: path of the output file
moviespeed: speed of animation. between 0 and 10
note: material studio file cannot use a partial periodic system. If partial
perodic system was inputted, full periodicity was assumed.
"""
if moviespeed < 0 or moviespeed > 10:
raise ValueError('moviespeed only between 0 and 10 allowed')
if hasattr(images, 'get_positions'):
images = [images]
XSD, ATR = _write_xsd_html(images, connectivity)
ATR.attrib['NumChildren'] = '2'
natoms = len(images[0])
bonds = list()
if connectivity is not None:
for i in range(connectivity.shape[0]):
for j in range(i + 1, connectivity.shape[0]):
if connectivity[i, j]:
bonds.append([i, j])
# non-periodic system
s = '!BIOSYM archive 3\n'
if not images[0].pbc.all():
# Write trajectory
SetChild(ATR, 'Trajectory', dict(
ID=str(natoms + 3 + len(bonds)),
Increment='-1',
End=str(len(images)),
Type='arc',
Speed=str(moviespeed),
FrameClassType='Atom',
))
# write frame information file
s += 'PBC=OFF\n'
for image in images:
s += _image_header
s += '\n'
an = image.get_chemical_symbols()
xyz = image.get_positions()
for i in range(natoms):
s += _get_atom_str(an[i], xyz[i, :])
s += _image_footer
# periodic system
else:
SetChild(ATR, 'Trajectory', dict(
ID=str(natoms + 9 + len(bonds)),
Increment='-1',
End=str(len(images)),
Type='arc',
Speed=str(moviespeed),
FrameClassType='Atom',
))
# write frame information file
s += 'PBC=ON\n'
for image in images:
s += _image_header
s += 'PBC'
vec = image.cell.lengths()
s += '{:>10.4f}{:>10.4f}{:>10.4f}'.format(vec[0], vec[1], vec[2])
angles = image.cell.angles()
s += '{:>10.4f}{:>10.4f}{:>10.4f}'.format(*angles)
s += '\n'
an = image.get_chemical_symbols()
angrad = np.deg2rad(angles)
cell = np.zeros((3, 3))
cell[0, :] = [vec[0], 0, 0]
cell[1, :] = (np.array([np.cos(angrad[2]), np.sin(angrad[2]), 0])
* vec[1])
cell[2, 0] = vec[2] * np.cos(angrad[1])
cell[2, 1] = ((vec[1] * vec[2] * np.cos(angrad[0])
- cell[1, 0] * cell[2, 0]) / cell[1, 1])
cell[2, 2] = np.sqrt(vec[2]**2 - cell[2, 0]**2 - cell[2, 1]**2)
xyz = np.dot(image.get_scaled_positions(), cell)
for i in range(natoms):
s += _get_atom_str(an[i], xyz[i, :])
s += _image_footer
# print arc file
if isinstance(filename, str):
farcname = filename[:-3] + 'arc'
else:
farcname = filename.name[:-3] + 'arc'
with open(farcname, 'w') as farc:
farc.write(s)
# check if file is an object or not.
openandclose = False
try:
if isinstance(filename, str):
fd = open(filename, 'w')
openandclose = True
else: # Assume it's a 'file-like object'
fd = filename
# Return a pretty-printed XML string for the Element.
rough_string = ET.tostring(XSD, 'utf-8')
reparsed = minidom.parseString(rough_string)
Document = reparsed.toprettyxml(indent='\t')
# write
fd.write(Document)
finally:
if openandclose:
fd.close()
def read_xtd(filename, index=-1):
"""Import xtd file (Materials Studio)
Xtd files always come with arc file, and arc file
contains all the relevant information to make atoms
so only Arc file needs to be read
"""
if isinstance(filename, str):
arcfilename = filename[:-3] + 'arc'
else:
arcfilename = filename.name[:-3] + 'arc'
# This trick with opening a totally different file is a gross violation of
# common sense.
with open(arcfilename, 'r') as fd:
return read_arcfile(fd, index)
def read_arcfile(fd, index):
images = []
# the first line is comment
fd.readline()
pbc = 'ON' in fd.readline()
L = fd.readline()
while L != '':
if '!' not in L: # flag for the start of an image
L = fd.readline()
continue
if pbc:
L = fd.readline()
cell = [float(d) for d in L.split()[1:]]
else:
fd.readline()
symbols = []
coords = []
while True:
line = fd.readline()
L = line.split()
if not line or 'end' in L:
break
symbols.append(L[0])
coords.append([float(x) for x in L[1:4]])
if pbc:
image = Atoms(symbols, positions=coords, cell=cell, pbc=pbc)
else:
image = Atoms(symbols, positions=coords, pbc=pbc)
images.append(image)
L = fd.readline()
if not index:
return images
else:
return images[index]
|
