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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
|
"""
IO functions for DMol3 file formats.
read/write functionality for car, incoor and arc file formats
only car format is added to known ase file extensions
use format='dmol-arc' or 'dmol-incoor' for others
car structure file - Angstrom and cellpar description of cell.
incoor structure file - Bohr and cellvector describption of cell.
Note: incoor file not used if car file present.
arc multiple-structure file - Angstrom and cellpar description of cell.
The formats follow strict formatting
car
----
col: 1-5 atom name
col: 7-20 x Cartesian coordinate of atom in A
col: 22-35 y Cartesian coordinate of atom in A
col: 37-50 z Cartesian coordinate of atom in A
col: 52-55 type of residue containing atom
col: 57-63 residue sequence name relative to beginning of current molecule,
left justified
col: 64-70 potential type of atom left justified
col: 72-73 element symbol
col: 75-80 partial charge on atom
incoor
-------
$cell vectors
37.83609647462165 0.00000000000000 0.00000000000000
0.00000000000000 37.60366016124745 0.00000000000000
0.00000000000000 0.00000000000000 25.29020473078921
$coordinates
Si 15.94182672614820 1.85274838936809 16.01426481346124
Si 4.45559370448989 2.68957177851318 -0.05326937257442
$end
arc
----
multiple images of car format separated with $end
"""
from datetime import datetime
import numpy as np
from ase import Atom, Atoms
from ase.geometry.cell import cell_to_cellpar, cellpar_to_cell
from ase.units import Bohr
from ase.utils import reader, writer
@writer
def write_dmol_car(fd, atoms):
""" Write a dmol car-file from an Atoms object
Notes
-----
The positions written to file are rotated as to align with the cell when
reading (due to cellpar information)
Can not handle multiple images.
Only allows for pbc 111 or 000.
"""
fd.write('!BIOSYM archive 3\n')
dt = datetime.now()
symbols = atoms.get_chemical_symbols()
if np.all(atoms.pbc):
# Rotate positions so they will align with cellpar cell
cellpar = cell_to_cellpar(atoms.cell)
new_cell = cellpar_to_cell(cellpar)
lstsq_fit = np.linalg.lstsq(atoms.cell, new_cell, rcond=-1)
# rcond=-1 silences FutureWarning in numpy 1.14
R = lstsq_fit[0]
positions = np.dot(atoms.positions, R)
fd.write('PBC=ON\n\n')
fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y'))
fd.write('PBC %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n' % tuple(cellpar))
elif not np.any(atoms.pbc): # [False,False,False]
fd.write('PBC=OFF\n\n')
fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y'))
positions = atoms.positions
else:
raise ValueError('PBC must be all true or all false for .car format')
for i, (sym, pos) in enumerate(zip(symbols, positions)):
fd.write('%-6s %12.8f %12.8f %12.8f XXXX 1 xx %-2s '
'0.000\n' % (sym + str(i + 1), pos[0], pos[1], pos[2], sym))
fd.write('end\nend\n')
@reader
def read_dmol_car(fd):
""" Read a dmol car-file and return an Atoms object.
Notes
-----
Cell is constructed from cellpar so orientation of cell might be off.
"""
lines = fd.readlines()
atoms = Atoms()
start_line = 4
if lines[1][4:6] == 'ON':
start_line += 1
cell_dat = np.array([float(fld) for fld in lines[4].split()[1:7]])
cell = cellpar_to_cell(cell_dat)
pbc = [True, True, True]
else:
cell = np.zeros((3, 3))
pbc = [False, False, False]
symbols = []
positions = []
for line in lines[start_line:]:
if line.startswith('end'):
break
flds = line.split()
symbols.append(flds[7])
positions.append(flds[1:4])
atoms.append(Atom(flds[7], flds[1:4]))
atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=pbc)
return atoms
@writer
def write_dmol_incoor(fd, atoms, bohr=True):
""" Write a dmol incoor-file from an Atoms object
Notes
-----
Only used for pbc 111.
Can not handle multiple images.
DMol3 expect data in .incoor files to be in bohr, if bohr is false however
the data is written in Angstroms.
"""
if not np.all(atoms.pbc):
raise ValueError('PBC must be all true for .incoor format')
if bohr:
cell = atoms.cell / Bohr
positions = atoms.positions / Bohr
else:
cell = atoms.cell
positions = atoms.positions
fd.write('$cell vectors\n')
fd.write(' %18.14f %18.14f %18.14f\n' % (
cell[0, 0], cell[0, 1], cell[0, 2]))
fd.write(' %18.14f %18.14f %18.14f\n' % (
cell[1, 0], cell[1, 1], cell[1, 2]))
fd.write(' %18.14f %18.14f %18.14f\n' % (
cell[2, 0], cell[2, 1], cell[2, 2]))
fd.write('$coordinates\n')
for a, pos in zip(atoms, positions):
fd.write('%-12s%18.14f %18.14f %18.14f \n' % (
a.symbol, pos[0], pos[1], pos[2]))
fd.write('$end\n')
@reader
def read_dmol_incoor(fd, bohr=True):
""" Reads an incoor file and returns an atoms object.
Notes
-----
If bohr is True then incoor is assumed to be in bohr and the data
is rescaled to Angstrom.
"""
lines = fd.readlines()
symbols = []
positions = []
for i, line in enumerate(lines):
if line.startswith('$cell vectors'):
cell = np.zeros((3, 3))
for j, line in enumerate(lines[i + 1:i + 4]):
cell[j, :] = [float(fld) for fld in line.split()]
if line.startswith('$coordinates'):
j = i + 1
while True:
if lines[j].startswith('$end'):
break
flds = lines[j].split()
symbols.append(flds[0])
positions.append(flds[1:4])
j += 1
atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True)
if bohr:
atoms.cell = atoms.cell * Bohr
atoms.positions = atoms.positions * Bohr
return atoms
@writer
def write_dmol_arc(fd, images):
""" Writes all images to file filename in arc format.
Similar to the .car format only pbc 111 or 000 is supported.
"""
fd.write('!BIOSYM archive 3\n')
if np.all(images[0].pbc):
fd.write('PBC=ON\n\n')
# Rotate positions so they will align with cellpar cell
elif not np.any(images[0].pbc):
fd.write('PBC=OFF\n\n')
else:
raise ValueError('PBC must be all true or all false for .arc format')
for atoms in images:
dt = datetime.now()
symbols = atoms.get_chemical_symbols()
if np.all(atoms.pbc):
cellpar = cell_to_cellpar(atoms.cell)
new_cell = cellpar_to_cell(cellpar)
lstsq_fit = np.linalg.lstsq(atoms.cell, new_cell, rcond=-1)
R = lstsq_fit[0]
fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y'))
fd.write('PBC %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'
% tuple(cellpar))
positions = np.dot(atoms.positions, R)
elif not np.any(atoms.pbc): # [False,False,False]
fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y'))
positions = atoms.positions
else:
raise ValueError(
'PBC must be all true or all false for .arc format')
for i, (sym, pos) in enumerate(zip(symbols, positions)):
fd.write('%-6s %12.8f %12.8f %12.8f XXXX 1 xx %-2s '
'0.000\n' % (sym + str(i + 1), pos[0], pos[1], pos[2], sym))
fd.write('end\nend\n')
fd.write('\n')
@reader
def read_dmol_arc(fd, index=-1):
""" Read a dmol arc-file and return a series of Atoms objects (images). """
lines = fd.readlines()
images = []
if lines[1].startswith('PBC=ON'):
pbc = True
elif lines[1].startswith('PBC=OFF'):
pbc = False
else:
raise RuntimeError('Could not read pbc from second line in file')
i = 0
while i < len(lines):
cell = np.zeros((3, 3))
symbols = []
positions = []
# parse single image
if lines[i].startswith('!DATE'):
# read cell
if pbc:
cell_dat = np.array([float(fld)
for fld in lines[i + 1].split()[1:7]])
cell = cellpar_to_cell(cell_dat)
i += 1
i += 1
# read atoms
while not lines[i].startswith('end'):
flds = lines[i].split()
symbols.append(flds[7])
positions.append(flds[1:4])
i += 1
image = Atoms(symbols=symbols, positions=positions, cell=cell,
pbc=pbc)
images.append(image)
if len(images) == index:
return images[-1]
i += 1
# return requested images, code borrowed from ase/io/trajectory.py
if isinstance(index, int):
return images[index]
else:
from ase.io.formats import index2range
indices = index2range(index, len(images))
return [images[j] for j in indices]
|
