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
|
import numpy as np
from ase.atoms import Atoms
from ase.utils import reader, writer
@reader
def read_dftb(fd):
"""Method to read coordinates from the Geometry section
of a DFTB+ input file (typically called "dftb_in.hsd").
As described in the DFTB+ manual, this section can be
in a number of different formats. This reader supports
the GEN format and the so-called "explicit" format.
The "explicit" format is unique to DFTB+ input files.
The GEN format can also be used in a stand-alone fashion,
as coordinate files with a `.gen` extension. Reading and
writing such files is implemented in `ase.io.gen`.
"""
lines = fd.readlines()
atoms_pos = []
atom_symbols = []
type_names = []
my_pbc = False
fractional = False
mycell = []
for iline, line in enumerate(lines):
if line.strip().startswith('#'):
pass
elif 'genformat' in line.lower():
natoms = int(lines[iline + 1].split()[0])
if lines[iline + 1].split()[1].lower() == 's':
my_pbc = True
elif lines[iline + 1].split()[1].lower() == 'f':
my_pbc = True
fractional = True
symbols = lines[iline + 2].split()
for i in range(natoms):
index = iline + 3 + i
aindex = int(lines[index].split()[1]) - 1
atom_symbols.append(symbols[aindex])
position = [float(p) for p in lines[index].split()[2:]]
atoms_pos.append(position)
if my_pbc:
for i in range(3):
index = iline + 4 + natoms + i
cell = [float(c) for c in lines[index].split()]
mycell.append(cell)
else:
if 'TypeNames' in line:
col = line.split()
for i in range(3, len(col) - 1):
type_names.append(col[i].strip("\""))
elif 'Periodic' in line:
if 'Yes' in line:
my_pbc = True
elif 'LatticeVectors' in line:
for imycell in range(3):
extraline = lines[iline + imycell + 1]
cols = extraline.split()
mycell.append(
[float(cols[0]), float(cols[1]), float(cols[2])])
else:
pass
if not my_pbc:
mycell = [0.] * 3
start_reading_coords = False
stop_reading_coords = False
for line in lines:
if line.strip().startswith('#'):
pass
else:
if 'TypesAndCoordinates' in line:
start_reading_coords = True
if start_reading_coords:
if '}' in line:
stop_reading_coords = True
if (start_reading_coords and not stop_reading_coords
and 'TypesAndCoordinates' not in line):
typeindexstr, xxx, yyy, zzz = line.split()[:4]
typeindex = int(typeindexstr)
symbol = type_names[typeindex - 1]
atom_symbols.append(symbol)
atoms_pos.append([float(xxx), float(yyy), float(zzz)])
if fractional:
atoms = Atoms(scaled_positions=atoms_pos, symbols=atom_symbols,
cell=mycell, pbc=my_pbc)
elif not fractional:
atoms = Atoms(positions=atoms_pos, symbols=atom_symbols,
cell=mycell, pbc=my_pbc)
return atoms
def read_dftb_velocities(atoms, filename='geo_end.xyz'):
"""Method to read velocities (AA/ps) from DFTB+ output file geo_end.xyz
"""
from ase.units import second
# AA/ps -> ase units
AngdivPs2ASE = 1.0 / (1e-12 * second)
with open(filename) as fd:
lines = fd.readlines()
# remove empty lines
lines_ok = []
for line in lines:
if line.rstrip():
lines_ok.append(line)
velocities = []
natoms = len(atoms)
last_lines = lines_ok[-natoms:]
for iline, line in enumerate(last_lines):
inp = line.split()
velocities.append([float(inp[5]) * AngdivPs2ASE,
float(inp[6]) * AngdivPs2ASE,
float(inp[7]) * AngdivPs2ASE])
atoms.set_velocities(velocities)
return atoms
@reader
def read_dftb_lattice(fileobj, images=None):
"""Read lattice vectors from MD and return them as a list.
If a molecules are parsed add them there."""
if images is not None:
append = True
if hasattr(images, 'get_positions'):
images = [images]
else:
append = False
fileobj.seek(0)
lattices = []
for line in fileobj:
if 'Lattice vectors' in line:
vec = []
for i in range(3): # DFTB+ only supports 3D PBC
line = fileobj.readline().split()
try:
line = [float(x) for x in line]
except ValueError:
raise ValueError('Lattice vector elements should be of '
'type float.')
vec.extend(line)
lattices.append(np.array(vec).reshape((3, 3)))
if append:
if len(images) != len(lattices):
raise ValueError('Length of images given does not match number of '
'cell vectors found')
for i, atoms in enumerate(images):
atoms.set_cell(lattices[i])
# DFTB+ only supports 3D PBC
atoms.set_pbc(True)
return
else:
return lattices
@writer
def write_dftb(fileobj, images):
"""Write structure in GEN format (refer to DFTB+ manual).
Multiple snapshots are not allowed. """
from ase.io.gen import write_gen
write_gen(fileobj, images)
def write_dftb_velocities(atoms, filename='velocities.txt'):
"""Method to write velocities (in atomic units) from ASE
to a file to be read by dftb+
"""
from ase.units import AUT, Bohr
# ase units -> atomic units
ASE2au = Bohr / AUT
with open(filename, 'w') as fd:
velocities = atoms.get_velocities()
for velocity in velocities:
fd.write(' %19.16f %19.16f %19.16f \n'
% (velocity[0] / ASE2au,
velocity[1] / ASE2au,
velocity[2] / ASE2au))
|
