Coverage for /builds/ase/ase/ase/io/dftb.py: 92.91%

1# fmt: off 

2 

3from typing import Sequence, Union 

4 

5import numpy as np 

6 

7from ase.atoms import Atoms 

8from ase.utils import reader, writer 

9 

10 

11@reader 

12def read_dftb(fd): 

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

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

15 

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

17 in a number of different formats. This reader supports 

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

19 

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

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

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

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

24 """ 

25 lines = fd.readlines() 

26 

27 atoms_pos = [] 

28 atom_symbols = [] 

29 type_names = [] 

30 my_pbc = False 

31 fractional = False 

32 mycell = [] 

33 

34 for iline, line in enumerate(lines): 

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

36 pass 

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

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

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

40 my_pbc = True 

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

42 my_pbc = True 

43 fractional = True 

44 

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

46 

47 for i in range(natoms): 

48 index = iline + 3 + i 

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

50 atom_symbols.append(symbols[aindex]) 

51 

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

53 atoms_pos.append(position) 

54 

55 if my_pbc: 

56 for i in range(3): 

57 index = iline + 4 + natoms + i 

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

59 mycell.append(cell) 

60 else: 

61 if 'TypeNames' in line: 

62 col = line.split() 

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

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

65 elif 'Periodic' in line: 

66 if 'Yes' in line: 

67 my_pbc = True 

68 elif 'LatticeVectors' in line: 

69 for imycell in range(3): 

70 extraline = lines[iline + imycell + 1] 

71 cols = extraline.split() 

72 mycell.append( 

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

74 else: 

75 pass 

76 

77 if not my_pbc: 

78 mycell = [0.] * 3 

79 

80 start_reading_coords = False 

81 stop_reading_coords = False 

82 for line in lines: 

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

84 pass 

85 else: 

86 if 'TypesAndCoordinates' in line: 

87 start_reading_coords = True 

88 if start_reading_coords: 

89 if '}' in line: 

90 stop_reading_coords = True 

91 if (start_reading_coords and not stop_reading_coords 

92 and 'TypesAndCoordinates' not in line): 

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

94 typeindex = int(typeindexstr) 

95 symbol = type_names[typeindex - 1] 

96 atom_symbols.append(symbol) 

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

98 

99 if fractional: 

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

101 cell=mycell, pbc=my_pbc) 

102 elif not fractional: 

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

104 cell=mycell, pbc=my_pbc) 

105 

106 return atoms 

107 

108 

109def read_dftb_velocities(atoms, filename): 

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

111 """ 

112 from ase.units import second 

113 

114 # AA/ps -> ase units 

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

116 

117 with open(filename) as fd: 

118 lines = fd.readlines() 

119 

120 # remove empty lines 

121 lines_ok = [] 

122 for line in lines: 

123 if line.rstrip(): 

124 lines_ok.append(line) 

125 

126 velocities = [] 

127 natoms = len(atoms) 

128 last_lines = lines_ok[-natoms:] 

129 for iline, line in enumerate(last_lines): 

130 inp = line.split() 

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

132 float(inp[6]) * AngdivPs2ASE, 

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

134 

135 atoms.set_velocities(velocities) 

136 return atoms 

137 

138 

139@reader 

140def read_dftb_lattice(fileobj, images=None): 

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

142 

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

144 if images is not None: 

145 append = True 

146 if hasattr(images, 'get_positions'): 

147 images = [images] 

148 else: 

149 append = False 

150 

151 fileobj.seek(0) 

152 lattices = [] 

153 for line in fileobj: 

154 if 'Lattice vectors' in line: 

155 vec = [] 

156 for _ in range(3): 

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

158 try: 

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

160 except ValueError: 

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

162 'type float.') 

163 vec.extend(line) 

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

165 

166 if append: 

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

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

169 'cell vectors found') 

170 

171 for i, atoms in enumerate(images): 

172 atoms.set_cell(lattices[i]) 

173 # DFTB+ only supports 3D PBC 

174 atoms.set_pbc(True) 

175 return None 

176 else: 

177 return lattices 

178 

179 

180@writer 

181def write_dftb( 

182 fileobj, 

183 images: Union[Atoms, Sequence[Atoms]], 

184 fractional: bool = False, 

185): 

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

187 Multiple snapshots are not allowed. """ 

188 from ase.io.gen import write_gen 

189 write_gen(fileobj, images, fractional=fractional) 

190 

191 

192def write_dftb_velocities(atoms, filename): 

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

194 to a file to be read by dftb+ 

195 """ 

196 from ase.units import AUT, Bohr 

197 

198 # ase units -> atomic units 

199 ASE2au = Bohr / AUT 

200 

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

202 velocities = atoms.get_velocities() 

203 for velocity in velocities: 

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

205 % (velocity[0] / ASE2au, 

206 velocity[1] / ASE2au, 

207 velocity[2] / ASE2au))