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

1from typing import Sequence, Union 

2 

3import numpy as np 

4 

5from ase.atoms import Atoms 

6from ase.utils import reader, writer 

7 

8 

9@reader 

10def read_dftb(fd): 

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

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

13 

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

15 in a number of different formats. This reader supports 

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

17 

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

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

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

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

22 """ 

23 lines = fd.readlines() 

24 

25 atoms_pos = [] 

26 atom_symbols = [] 

27 type_names = [] 

28 my_pbc = False 

29 fractional = False 

30 mycell = [] 

31 

32 for iline, line in enumerate(lines): 

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

34 pass 

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

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

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

38 my_pbc = True 

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

40 my_pbc = True 

41 fractional = True 

42 

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

44 

45 for i in range(natoms): 

46 index = iline + 3 + i 

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

48 atom_symbols.append(symbols[aindex]) 

49 

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

51 atoms_pos.append(position) 

52 

53 if my_pbc: 

54 for i in range(3): 

55 index = iline + 4 + natoms + i 

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

57 mycell.append(cell) 

58 else: 

59 if 'TypeNames' in line: 

60 col = line.split() 

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

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

63 elif 'Periodic' in line: 

64 if 'Yes' in line: 

65 my_pbc = True 

66 elif 'LatticeVectors' in line: 

67 for imycell in range(3): 

68 extraline = lines[iline + imycell + 1] 

69 cols = extraline.split() 

70 mycell.append( 

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

72 else: 

73 pass 

74 

75 if not my_pbc: 

76 mycell = [0.] * 3 

77 

78 start_reading_coords = False 

79 stop_reading_coords = False 

80 for line in lines: 

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

82 pass 

83 else: 

84 if 'TypesAndCoordinates' in line: 

85 start_reading_coords = True 

86 if start_reading_coords: 

87 if '}' in line: 

88 stop_reading_coords = True 

89 if (start_reading_coords and not stop_reading_coords 

90 and 'TypesAndCoordinates' not in line): 

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

92 typeindex = int(typeindexstr) 

93 symbol = type_names[typeindex - 1] 

94 atom_symbols.append(symbol) 

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

96 

97 if fractional: 

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

99 cell=mycell, pbc=my_pbc) 

100 elif not fractional: 

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

102 cell=mycell, pbc=my_pbc) 

103 

104 return atoms 

105 

106 

107def read_dftb_velocities(atoms, filename): 

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

109 """ 

110 from ase.units import second 

111 

112 # AA/ps -> ase units 

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

114 

115 with open(filename) as fd: 

116 lines = fd.readlines() 

117 

118 # remove empty lines 

119 lines_ok = [] 

120 for line in lines: 

121 if line.rstrip(): 

122 lines_ok.append(line) 

123 

124 velocities = [] 

125 natoms = len(atoms) 

126 last_lines = lines_ok[-natoms:] 

127 for iline, line in enumerate(last_lines): 

128 inp = line.split() 

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

130 float(inp[6]) * AngdivPs2ASE, 

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

132 

133 atoms.set_velocities(velocities) 

134 return atoms 

135 

136 

137@reader 

138def read_dftb_lattice(fileobj, images=None): 

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

140 

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

142 if images is not None: 

143 append = True 

144 if hasattr(images, 'get_positions'): 

145 images = [images] 

146 else: 

147 append = False 

148 

149 fileobj.seek(0) 

150 lattices = [] 

151 for line in fileobj: 

152 if 'Lattice vectors' in line: 

153 vec = [] 

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

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

156 try: 

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

158 except ValueError: 

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

160 'type float.') 

161 vec.extend(line) 

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

163 

164 if append: 

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

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

167 'cell vectors found') 

168 

169 for i, atoms in enumerate(images): 

170 atoms.set_cell(lattices[i]) 

171 # DFTB+ only supports 3D PBC 

172 atoms.set_pbc(True) 

173 return 

174 else: 

175 return lattices 

176 

177 

178@writer 

179def write_dftb( 

180 fileobj, 

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

182 fractional: bool = False, 

183): 

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

185 Multiple snapshots are not allowed. """ 

186 from ase.io.gen import write_gen 

187 write_gen(fileobj, images, fractional=fractional) 

188 

189 

190def write_dftb_velocities(atoms, filename): 

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

192 to a file to be read by dftb+ 

193 """ 

194 from ase.units import AUT, Bohr 

195 

196 # ase units -> atomic units 

197 ASE2au = Bohr / AUT 

198 

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

200 velocities = atoms.get_velocities() 

201 for velocity in velocities: 

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

203 % (velocity[0] / ASE2au, 

204 velocity[1] / ASE2au, 

205 velocity[2] / ASE2au))