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
|
from ase.units import Bohr
def read_turbomole(fd):
"""Method to read turbomole coord file
coords in bohr, atom types in lowercase, format:
$coord
x y z atomtype
x y z atomtype f
$end
Above 'f' means a fixed atom.
"""
from ase import Atoms
from ase.constraints import FixAtoms
lines = fd.readlines()
atoms_pos = []
atom_symbols = []
myconstraints = []
# find $coord section;
# does not necessarily have to be the first $<something> in file...
for i, l in enumerate(lines):
if l.strip().startswith('$coord'):
start = i
break
for line in lines[start + 1:]:
if line.startswith('$'): # start of new section
break
else:
x, y, z, symbolraw = line.split()[:4]
symbolshort = symbolraw.strip()
symbol = symbolshort[0].upper() + symbolshort[1:].lower()
# print(symbol)
atom_symbols.append(symbol)
atoms_pos.append(
[float(x) * Bohr, float(y) * Bohr, float(z) * Bohr]
)
cols = line.split()
if (len(cols) == 5):
fixedstr = line.split()[4].strip()
if (fixedstr == "f"):
myconstraints.append(True)
else:
myconstraints.append(False)
else:
myconstraints.append(False)
# convert Turbomole ghost atom Q to X
atom_symbols = [element if element !=
'Q' else 'X' for element in atom_symbols]
atoms = Atoms(positions=atoms_pos, symbols=atom_symbols, pbc=False)
c = FixAtoms(mask=myconstraints)
atoms.set_constraint(c)
return atoms
class TurbomoleFormatError(ValueError):
default_message = ('Data format in file does not correspond to known '
'Turbomole gradient format')
def __init__(self, *args, **kwargs):
if args or kwargs:
ValueError.__init__(self, *args, **kwargs)
else:
ValueError.__init__(self, self.default_message)
def read_turbomole_gradient(fd, index=-1):
""" Method to read turbomole gradient file """
# read entire file
lines = [x.strip() for x in fd.readlines()]
# find $grad section
start = end = -1
for i, line in enumerate(lines):
if not line.startswith('$'):
continue
if line.split()[0] == '$grad':
start = i
elif start >= 0:
end = i
break
if end <= start:
raise RuntimeError('File does not contain a valid \'$grad\' section')
# trim lines to $grad
del lines[:start + 1]
del lines[end - 1 - start:]
# Interpret $grad section
from ase import Atom, Atoms
from ase.calculators.singlepoint import SinglePointCalculator
from ase.units import Bohr, Hartree
images = []
while lines: # loop over optimization cycles
# header line
# cycle = 1 SCF energy = -267.6666811409 |dE/dxyz| = 0.157112 # noqa: E501
fields = lines[0].split('=')
try:
# cycle = int(fields[1].split()[0])
energy = float(fields[2].split()[0]) * Hartree
# gradient = float(fields[3].split()[0])
except (IndexError, ValueError) as e:
raise TurbomoleFormatError from e
# coordinates/gradient
atoms = Atoms()
forces = []
for line in lines[1:]:
fields = line.split()
if len(fields) == 4: # coordinates
# 0.00000000000000 0.00000000000000 0.00000000000000 c # noqa: E501
try:
symbol = fields[3].lower().capitalize()
# if dummy atom specified, substitute 'Q' with 'X'
if symbol == 'Q':
symbol = 'X'
position = tuple(Bohr * float(x) for x in fields[0:3])
except ValueError as e:
raise TurbomoleFormatError from e
atoms.append(Atom(symbol, position))
elif len(fields) == 3: # gradients
# -.51654903354681D-07 -.51654903206651D-07 0.51654903169644D-07 # noqa: E501
grad = []
for val in fields[:3]:
try:
grad.append(
-float(val.replace('D', 'E')) * Hartree / Bohr
)
except ValueError as e:
raise TurbomoleFormatError from e
forces.append(grad)
else: # next cycle
break
# calculator
calc = SinglePointCalculator(atoms, energy=energy, forces=forces)
atoms.calc = calc
# save frame
images.append(atoms)
# delete this frame from data to be handled
del lines[:2 * len(atoms) + 1]
return images[index]
def write_turbomole(fd, atoms):
""" Method to write turbomole coord file
"""
from ase.constraints import FixAtoms
coord = atoms.get_positions()
symbols = atoms.get_chemical_symbols()
# convert X to Q for Turbomole ghost atoms
symbols = [element if element != 'X' else 'Q' for element in symbols]
fix_indices = set()
if atoms.constraints:
for constr in atoms.constraints:
if isinstance(constr, FixAtoms):
fix_indices.update(constr.get_indices())
fix_str = []
for i in range(len(atoms)):
if i in fix_indices:
fix_str.append('f')
else:
fix_str.append('')
fd.write('$coord\n')
for (x, y, z), s, fix in zip(coord, symbols, fix_str):
fd.write('%20.14f %20.14f %20.14f %2s %2s \n'
% (x / Bohr, y / Bohr, z / Bohr, s.lower(), fix))
fd.write('$end\n')
|
