1import numpy as np 

2from ase.atoms import Atoms 

3from ase.utils import reader, writer 

4 

5 

6@reader 

7def read_dftb(fd): 

8 """Method to read coordinates from the Geometry section 

9 of a DFTB+ input file (typically called "dftb_in.hsd"). 

10 

11 As described in the DFTB+ manual, this section can be 

12 in a number of different formats. This reader supports 

13 the GEN format and the so-called "explicit" format. 

14 

15 The "explicit" format is unique to DFTB+ input files. 

16 The GEN format can also be used in a stand-alone fashion, 

17 as coordinate files with a `.gen` extension. Reading and 

18 writing such files is implemented in `ase.io.gen`. 

19 """ 

20 lines = fd.readlines() 

21 

22 atoms_pos = [] 

23 atom_symbols = [] 

24 type_names = [] 

25 my_pbc = False 

26 fractional = False 

27 mycell = [] 

28 

29 for iline, line in enumerate(lines): 

30 if line.strip().startswith('#'): 

31 pass 

32 elif 'genformat' in line.lower(): 

33 natoms = int(lines[iline + 1].split()[0]) 

34 if lines[iline + 1].split()[1].lower() == 's': 

35 my_pbc = True 

36 elif lines[iline + 1].split()[1].lower() == 'f': 

37 my_pbc = True 

38 fractional = True 

39 

40 symbols = lines[iline + 2].split() 

41 

42 for i in range(natoms): 

43 index = iline + 3 + i 

44 aindex = int(lines[index].split()[1]) - 1 

45 atom_symbols.append(symbols[aindex]) 

46 

47 position = [float(p) for p in lines[index].split()[2:]] 

48 atoms_pos.append(position) 

49 

50 if my_pbc: 

51 for i in range(3): 

52 index = iline + 4 + natoms + i 

53 cell = [float(c) for c in lines[index].split()] 

54 mycell.append(cell) 

55 else: 

56 if 'TypeNames' in line: 

57 col = line.split() 

58 for i in range(3, len(col) - 1): 

59 type_names.append(col[i].strip("\"")) 

60 elif 'Periodic' in line: 

61 if 'Yes' in line: 

62 my_pbc = True 

63 elif 'LatticeVectors' in line: 

64 for imycell in range(3): 

65 extraline = lines[iline + imycell + 1] 

66 cols = extraline.split() 

67 mycell.append( 

68 [float(cols[0]), float(cols[1]), float(cols[2])]) 

69 else: 

70 pass 

71 

72 if not my_pbc: 

73 mycell = [0.] * 3 

74 

75 start_reading_coords = False 

76 stop_reading_coords = False 

77 for line in lines: 

78 if line.strip().startswith('#'): 

79 pass 

80 else: 

81 if 'TypesAndCoordinates' in line: 

82 start_reading_coords = True 

83 if start_reading_coords: 

84 if '}' in line: 

85 stop_reading_coords = True 

86 if (start_reading_coords and not stop_reading_coords 

87 and 'TypesAndCoordinates' not in line): 

88 typeindexstr, xxx, yyy, zzz = line.split()[:4] 

89 typeindex = int(typeindexstr) 

90 symbol = type_names[typeindex - 1] 

91 atom_symbols.append(symbol) 

92 atoms_pos.append([float(xxx), float(yyy), float(zzz)]) 

93 

94 if fractional: 

95 atoms = Atoms(scaled_positions=atoms_pos, symbols=atom_symbols, 

96 cell=mycell, pbc=my_pbc) 

97 elif not fractional: 

98 atoms = Atoms(positions=atoms_pos, symbols=atom_symbols, 

99 cell=mycell, pbc=my_pbc) 

100 

101 return atoms 

102 

103 

104def read_dftb_velocities(atoms, filename): 

105 """Method to read velocities (AA/ps) from DFTB+ output file geo_end.xyz 

106 """ 

107 from ase.units import second 

108 # AA/ps -> ase units 

109 AngdivPs2ASE = 1.0 / (1e-12 * second) 

110 

111 with open(filename) as fd: 

112 lines = fd.readlines() 

113 

114 # remove empty lines 

115 lines_ok = [] 

116 for line in lines: 

117 if line.rstrip(): 

118 lines_ok.append(line) 

119 

120 velocities = [] 

121 natoms = len(atoms) 

122 last_lines = lines_ok[-natoms:] 

123 for iline, line in enumerate(last_lines): 

124 inp = line.split() 

125 velocities.append([float(inp[5]) * AngdivPs2ASE, 

126 float(inp[6]) * AngdivPs2ASE, 

127 float(inp[7]) * AngdivPs2ASE]) 

128 

129 atoms.set_velocities(velocities) 

130 return atoms 

131 

132 

133@reader 

134def read_dftb_lattice(fileobj, images=None): 

135 """Read lattice vectors from MD and return them as a list. 

136 

137 If a molecules are parsed add them there.""" 

138 if images is not None: 

139 append = True 

140 if hasattr(images, 'get_positions'): 

141 images = [images] 

142 else: 

143 append = False 

144 

145 fileobj.seek(0) 

146 lattices = [] 

147 for line in fileobj: 

148 if 'Lattice vectors' in line: 

149 vec = [] 

150 for i in range(3): # DFTB+ only supports 3D PBC 

151 line = fileobj.readline().split() 

152 try: 

153 line = [float(x) for x in line] 

154 except ValueError: 

155 raise ValueError('Lattice vector elements should be of ' 

156 'type float.') 

157 vec.extend(line) 

158 lattices.append(np.array(vec).reshape((3, 3))) 

159 

160 if append: 

161 if len(images) != len(lattices): 

162 raise ValueError('Length of images given does not match number of ' 

163 'cell vectors found') 

164 

165 for i, atoms in enumerate(images): 

166 atoms.set_cell(lattices[i]) 

167 # DFTB+ only supports 3D PBC 

168 atoms.set_pbc(True) 

169 return 

170 else: 

171 return lattices 

172 

173 

174@writer 

175def write_dftb(fileobj, images): 

176 """Write structure in GEN format (refer to DFTB+ manual). 

177 Multiple snapshots are not allowed. """ 

178 from ase.io.gen import write_gen 

179 write_gen(fileobj, images) 

180 

181 

182def write_dftb_velocities(atoms, filename): 

183 """Method to write velocities (in atomic units) from ASE 

184 to a file to be read by dftb+ 

185 """ 

186 from ase.units import AUT, Bohr 

187 # ase units -> atomic units 

188 ASE2au = Bohr / AUT 

189 

190 with open(filename, 'w') as fd: 

191 velocities = atoms.get_velocities() 

192 for velocity in velocities: 

193 fd.write(' %19.16f %19.16f %19.16f \n' 

194 % (velocity[0] / ASE2au, 

195 velocity[1] / ASE2au, 

196 velocity[2] / ASE2au))