1""" 

2IO functions for DMol3 file formats. 

3 

4read/write functionality for car, incoor and arc file formats 

5only car format is added to known ase file extensions 

6use format='dmol-arc' or 'dmol-incoor' for others 

7 

8car structure file - Angstrom and cellpar description of cell. 

9incoor structure file - Bohr and cellvector describption of cell. 

10 Note: incoor file not used if car file present. 

11arc multiple-structure file - Angstrom and cellpar description of cell. 

12 

13 

14The formats follow strict formatting 

15 

16car 

17---- 

18col: 1-5 atom name 

19col: 7-20 x Cartesian coordinate of atom in A 

20col: 22-35 y Cartesian coordinate of atom in A 

21col: 37-50 z Cartesian coordinate of atom in A 

22col: 52-55 type of residue containing atom 

23col: 57-63 residue sequence name relative to beginning of current molecule, 

24 left justified 

25col: 64-70 potential type of atom left justified 

26col: 72-73 element symbol 

27col: 75-80 partial charge on atom 

28 

29 

30incoor 

31------- 

32$cell vectors 

33 37.83609647462165 0.00000000000000 0.00000000000000 

34 0.00000000000000 37.60366016124745 0.00000000000000 

35 0.00000000000000 0.00000000000000 25.29020473078921 

36$coordinates 

37Si 15.94182672614820 1.85274838936809 16.01426481346124 

38Si 4.45559370448989 2.68957177851318 -0.05326937257442 

39$end 

40 

41 

42arc 

43---- 

44multiple images of car format separated with $end 

45 

46 

47""" 

48 

49from datetime import datetime 

50import numpy as np 

51 

52from ase import Atom, Atoms 

53from ase.geometry.cell import cell_to_cellpar, cellpar_to_cell 

54from ase.units import Bohr 

55from ase.utils import reader, writer 

56 

57 

58@writer 

59def write_dmol_car(fd, atoms): 

60 """ Write a dmol car-file from an Atoms object 

61 

62 Notes 

63 ----- 

64 The positions written to file are rotated as to align with the cell when 

65 reading (due to cellpar information) 

66 Can not handle multiple images. 

67 Only allows for pbc 111 or 000. 

68 """ 

69 

70 fd.write('!BIOSYM archive 3\n') 

71 dt = datetime.now() 

72 

73 symbols = atoms.get_chemical_symbols() 

74 if np.all(atoms.pbc): 

75 # Rotate positions so they will align with cellpar cell 

76 cellpar = cell_to_cellpar(atoms.cell) 

77 new_cell = cellpar_to_cell(cellpar) 

78 lstsq_fit = np.linalg.lstsq(atoms.cell, new_cell, rcond=-1) 

79 # rcond=-1 silences FutureWarning in numpy 1.14 

80 R = lstsq_fit[0] 

81 positions = np.dot(atoms.positions, R) 

82 

83 fd.write('PBC=ON\n\n') 

84 fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y')) 

85 fd.write('PBC %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n' % tuple(cellpar)) 

86 elif not np.any(atoms.pbc): # [False,False,False] 

87 fd.write('PBC=OFF\n\n') 

88 fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y')) 

89 positions = atoms.positions 

90 else: 

91 raise ValueError('PBC must be all true or all false for .car format') 

92 

93 for i, (sym, pos) in enumerate(zip(symbols, positions)): 

94 fd.write('%-6s %12.8f %12.8f %12.8f XXXX 1 xx %-2s ' 

95 '0.000\n' % (sym + str(i + 1), pos[0], pos[1], pos[2], sym)) 

96 fd.write('end\nend\n') 

97 

98 

99@reader 

100def read_dmol_car(fd): 

101 """ Read a dmol car-file and return an Atoms object. 

102 

103 Notes 

104 ----- 

105 Cell is constructed from cellpar so orientation of cell might be off. 

106 """ 

107 

108 lines = fd.readlines() 

109 atoms = Atoms() 

110 

111 start_line = 4 

112 

113 if lines[1][4:6] == 'ON': 

114 start_line += 1 

115 cell_dat = np.array([float(fld) for fld in lines[4].split()[1:7]]) 

116 cell = cellpar_to_cell(cell_dat) 

117 pbc = [True, True, True] 

118 else: 

119 cell = np.zeros((3, 3)) 

120 pbc = [False, False, False] 

121 

122 symbols = [] 

123 positions = [] 

124 for line in lines[start_line:]: 

125 if line.startswith('end'): 

126 break 

127 flds = line.split() 

128 symbols.append(flds[7]) 

129 positions.append(flds[1:4]) 

130 atoms.append(Atom(flds[7], flds[1:4])) 

131 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=pbc) 

132 return atoms 

133 

134 

135@writer 

136def write_dmol_incoor(fd, atoms, bohr=True): 

137 """ Write a dmol incoor-file from an Atoms object 

138 

139 Notes 

140 ----- 

141 Only used for pbc 111. 

142 Can not handle multiple images. 

143 DMol3 expect data in .incoor files to be in bohr, if bohr is false however 

144 the data is written in Angstroms. 

145 """ 

146 

147 if not np.all(atoms.pbc): 

148 raise ValueError('PBC must be all true for .incoor format') 

149 

150 if bohr: 

151 cell = atoms.cell / Bohr 

152 positions = atoms.positions / Bohr 

153 else: 

154 cell = atoms.cell 

155 positions = atoms.positions 

156 

157 fd.write('$cell vectors\n') 

158 fd.write(' %18.14f %18.14f %18.14f\n' % ( 

159 cell[0, 0], cell[0, 1], cell[0, 2])) 

160 fd.write(' %18.14f %18.14f %18.14f\n' % ( 

161 cell[1, 0], cell[1, 1], cell[1, 2])) 

162 fd.write(' %18.14f %18.14f %18.14f\n' % ( 

163 cell[2, 0], cell[2, 1], cell[2, 2])) 

164 

165 fd.write('$coordinates\n') 

166 for a, pos in zip(atoms, positions): 

167 fd.write('%-12s%18.14f %18.14f %18.14f \n' % ( 

168 a.symbol, pos[0], pos[1], pos[2])) 

169 fd.write('$end\n') 

170 

171 

172@reader 

173def read_dmol_incoor(fd, bohr=True): 

174 """ Reads an incoor file and returns an atoms object. 

175 

176 Notes 

177 ----- 

178 If bohr is True then incoor is assumed to be in bohr and the data 

179 is rescaled to Angstrom. 

180 """ 

181 

182 lines = fd.readlines() 

183 symbols = [] 

184 positions = [] 

185 for i, line in enumerate(lines): 

186 if line.startswith('$cell vectors'): 

187 cell = np.zeros((3, 3)) 

188 for j, line in enumerate(lines[i + 1:i + 4]): 

189 cell[j, :] = [float(fld) for fld in line.split()] 

190 if line.startswith('$coordinates'): 

191 j = i + 1 

192 while True: 

193 if lines[j].startswith('$end'): 

194 break 

195 flds = lines[j].split() 

196 symbols.append(flds[0]) 

197 positions.append(flds[1:4]) 

198 j += 1 

199 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True) 

200 if bohr: 

201 atoms.cell = atoms.cell * Bohr 

202 atoms.positions = atoms.positions * Bohr 

203 return atoms 

204 

205 

206@writer 

207def write_dmol_arc(fd, images): 

208 """ Writes all images to file filename in arc format. 

209 

210 Similar to the .car format only pbc 111 or 000 is supported. 

211 """ 

212 

213 fd.write('!BIOSYM archive 3\n') 

214 if np.all(images[0].pbc): 

215 fd.write('PBC=ON\n\n') 

216 # Rotate positions so they will align with cellpar cell 

217 elif not np.any(images[0].pbc): 

218 fd.write('PBC=OFF\n\n') 

219 else: 

220 raise ValueError('PBC must be all true or all false for .arc format') 

221 for atoms in images: 

222 dt = datetime.now() 

223 symbols = atoms.get_chemical_symbols() 

224 if np.all(atoms.pbc): 

225 cellpar = cell_to_cellpar(atoms.cell) 

226 new_cell = cellpar_to_cell(cellpar) 

227 lstsq_fit = np.linalg.lstsq(atoms.cell, new_cell, rcond=-1) 

228 R = lstsq_fit[0] 

229 fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y')) 

230 fd.write('PBC %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n' 

231 % tuple(cellpar)) 

232 positions = np.dot(atoms.positions, R) 

233 elif not np.any(atoms.pbc): # [False,False,False] 

234 fd.write('!DATE %s\n' % dt.strftime('%b %d %H:%m:%S %Y')) 

235 positions = atoms.positions 

236 else: 

237 raise ValueError( 

238 'PBC must be all true or all false for .arc format') 

239 for i, (sym, pos) in enumerate(zip(symbols, positions)): 

240 fd.write('%-6s %12.8f %12.8f %12.8f XXXX 1 xx %-2s ' 

241 '0.000\n' % (sym + str(i + 1), pos[0], pos[1], pos[2], sym)) 

242 fd.write('end\nend\n') 

243 fd.write('\n') 

244 

245 

246@reader 

247def read_dmol_arc(fd, index=-1): 

248 """ Read a dmol arc-file and return a series of Atoms objects (images). """ 

249 

250 lines = fd.readlines() 

251 images = [] 

252 

253 if lines[1].startswith('PBC=ON'): 

254 pbc = True 

255 elif lines[1].startswith('PBC=OFF'): 

256 pbc = False 

257 else: 

258 raise RuntimeError('Could not read pbc from second line in file') 

259 

260 i = 0 

261 while i < len(lines): 

262 cell = np.zeros((3, 3)) 

263 symbols = [] 

264 positions = [] 

265 # parse single image 

266 if lines[i].startswith('!DATE'): 

267 # read cell 

268 if pbc: 

269 cell_dat = np.array([float(fld) 

270 for fld in lines[i + 1].split()[1:7]]) 

271 cell = cellpar_to_cell(cell_dat) 

272 i += 1 

273 i += 1 

274 # read atoms 

275 while not lines[i].startswith('end'): 

276 flds = lines[i].split() 

277 symbols.append(flds[7]) 

278 positions.append(flds[1:4]) 

279 i += 1 

280 image = Atoms(symbols=symbols, positions=positions, cell=cell, 

281 pbc=pbc) 

282 images.append(image) 

283 if len(images) == index: 

284 return images[-1] 

285 i += 1 

286 

287 # return requested images, code borrowed from ase/io/trajectory.py 

288 if isinstance(index, int): 

289 return images[index] 

290 else: 

291 from ase.io.formats import index2range 

292 indices = index2range(index, len(images)) 

293 return [images[j] for j in indices]