1import numpy as np 

2 

3from ase.io.fortranfile import FortranFile 

4 

5 

6def read_rho(fname): 

7 "Read unformatted Siesta charge density file" 

8 

9 # TODO: 

10 # 

11 # Handle formatted and NetCDF files. 

12 # 

13 # Siesta source code (at least 2.0.2) can possibly also 

14 # save RHO as a _formatted_ file (the source code seems 

15 # prepared, but there seems to be no fdf-options for it though). 

16 # Siesta >= 3 has support for saving RHO as a NetCDF file 

17 # (according to manual) 

18 

19 fh = FortranFile(fname) 

20 

21 # Read (but ignore) unit cell vectors 

22 x = fh.readReals('d') 

23 if len(x) != 3 * 3: 

24 raise IOError('Failed to read cell vectors') 

25 

26 # Read number of grid points and spin components 

27 x = fh.readInts() 

28 if len(x) != 4: 

29 raise IOError('Failed to read grid size') 

30 gpts = x # number of 'X', 'Y', 'Z', 'spin' gridpoints 

31 

32 rho = np.zeros(gpts) 

33 for ispin in range(gpts[3]): 

34 for n3 in range(gpts[2]): 

35 for n2 in range(gpts[1]): 

36 x = fh.readReals('f') 

37 if len(x) != gpts[0]: 

38 raise IOError('Failed to read RHO[:,%i,%i,%i]' % 

39 (n2, n3, ispin)) 

40 rho[:, n2, n3, ispin] = x 

41 

42 fh.close() 

43 

44 return rho 

45 

46 

47def get_valence_charge(filename): 

48 """ Read the valence charge from '.psf'-file.""" 

49 with open(filename, 'r') as fd: 

50 fd.readline() 

51 fd.readline() 

52 fd.readline() 

53 valence = -float(fd.readline().split()[-1]) 

54 

55 return valence 

56 

57 

58def read_vca_synth_block(filename, species_number=None): 

59 """ Read the SyntheticAtoms block from the output of the 

60 'fractional' siesta utility. 

61 

62 Parameters: 

63 - filename: String with '.synth' output from fractional. 

64 - species_number: Optional argument to replace override the 

65 species number in the text block. 

66 

67 Returns: A string that can be inserted into the main '.fdf-file'. 

68 """ 

69 with open(filename, 'r') as fd: 

70 lines = fd.readlines() 

71 lines = lines[1:-1] 

72 

73 if species_number is not None: 

74 lines[0] = '%d\n' % species_number 

75 

76 block = ''.join(lines).strip() 

77 

78 return block 

79 

80 

81def readHSX(fname): 

82 """ 

83 Read unformatted siesta HSX file 

84 """ 

85 import collections 

86 

87 HSX_tuple = collections.namedtuple('HSX', 

88 ['norbitals', 'norbitals_sc', 'nspin', 

89 'nonzero', 'is_gamma', 'sc_orb2uc_orb', 

90 'row2nnzero', 'sparse_ind2column', 

91 'H_sparse', 'S_sparse', 

92 'aB2RaB_sparse', 'total_elec_charge', 

93 'temp']) 

94 

95 fh = FortranFile(fname) 

96 norbitals, norbitals_sc, nspin, nonzero = fh.readInts('i') 

97 is_gamma = fh.readInts('i')[0] 

98 

99 sc_orb2uc_orb = 0 

100 if is_gamma == 0: 

101 sc_orb2uc_orb = fh.readInts('i') 

102 

103 row2nnzero = fh.readInts('i') 

104 

105 sum_row2nnzero = np.sum(row2nnzero) 

106 if (sum_row2nnzero != nonzero): 

107 raise ValueError('sum_row2nnzero != nonzero: {0} != {1}' 

108 .format(sum_row2nnzero, nonzero)) 

109 

110 row2displ = np.zeros((norbitals), dtype=int) 

111 

112 for i in range(1, norbitals): 

113 row2displ[i] = row2displ[i - 1] + row2nnzero[i - 1] 

114 

115 max_nonzero = np.max(row2nnzero) 

116 int_buff = np.zeros((max_nonzero), dtype=int) 

117 sparse_ind2column = np.zeros((nonzero)) 

118 

119 # Fill the rows for each index in *_sparse arrays 

120 for irow in range(norbitals): 

121 f = row2nnzero[irow] 

122 int_buff[0:f] = fh.readInts('i') 

123 # read set of rows where nonzero elements reside 

124 d = row2displ[irow] 

125 sparse_ind2column[d:d + f] = int_buff[0:f] 

126 # END of Fill the rows for each index in *_sparse arrays 

127 

128 # allocate H, S and X matrices 

129 sp_buff = np.zeros((max_nonzero), dtype=float) 

130 

131 H_sparse = np.zeros((nonzero, nspin), dtype=float) 

132 S_sparse = np.zeros((nonzero), dtype=float) 

133 aB2RaB_sparse = np.zeros((3, nonzero), dtype=float) 

134 

135 # Read the data to H_sparse array 

136 for ispin in range(nspin): 

137 for irow in range(norbitals): 

138 d = row2displ[irow] 

139 f = row2nnzero[irow] 

140 sp_buff[0:f] = fh.readReals('f') 

141 H_sparse[d:d + f, ispin] = sp_buff[0:f] 

142 

143 # Read the data to S_sparse array 

144 for irow in range(norbitals): 

145 f = row2nnzero[irow] 

146 d = row2displ[irow] 

147 sp_buff[0:f] = fh.readReals('f') 

148 S_sparse[d:d + f] = sp_buff[0:f] 

149 

150 total_elec_charge, temp = fh.readReals('d') 

151 

152 sp_buff = np.zeros((3 * max_nonzero), dtype=float) 

153 # Read the data to S_sparse array 

154 for irow in range(norbitals): 

155 f = row2nnzero[irow] 

156 d = row2displ[irow] 

157 sp_buff[0: 3 * f] = fh.readReals('f') 

158 aB2RaB_sparse[0, d:d + f] = sp_buff[0:f] 

159 aB2RaB_sparse[1, d:d + f] = sp_buff[f:2 * f] 

160 aB2RaB_sparse[2, d:d + f] = sp_buff[2 * f:3 * f] 

161 

162 fh.close() 

163 

164 return HSX_tuple(norbitals, norbitals_sc, nspin, nonzero, is_gamma, 

165 sc_orb2uc_orb, row2nnzero, sparse_ind2column, H_sparse, 

166 S_sparse, aB2RaB_sparse, total_elec_charge, temp) 

167 

168 

169def readDIM(fname): 

170 """ 

171 Read unformatted siesta DIM file 

172 """ 

173 import collections 

174 

175 DIM_tuple = collections.namedtuple('DIM', ['natoms_sc', 'norbitals_sc', 

176 'norbitals', 'nspin', 

177 'nnonzero', 

178 'natoms_interacting']) 

179 

180 fh = FortranFile(fname) 

181 

182 natoms_sc = fh.readInts('i')[0] 

183 norbitals_sc = fh.readInts('i')[0] 

184 norbitals = fh.readInts('i')[0] 

185 nspin = fh.readInts('i')[0] 

186 nnonzero = fh.readInts('i')[0] 

187 natoms_interacting = fh.readInts('i')[0] 

188 fh.close() 

189 

190 return DIM_tuple(natoms_sc, norbitals_sc, norbitals, nspin, 

191 nnonzero, natoms_interacting) 

192 

193 

194def readPLD(fname, norbitals, natoms): 

195 """ 

196 Read unformatted siesta PLD file 

197 """ 

198 import collections 

199 # use struct library to read mixed data type from binary 

200 import struct 

201 

202 PLD_tuple = collections.namedtuple('PLD', ['max_rcut', 'orb2ao', 

203 'orb2uorb', 'orb2occ', 

204 'atm2sp', 'atm2shift', 

205 'coord_sc', 'cell', 

206 'nunit_cells']) 

207 

208 fh = FortranFile(fname) 

209 

210 orb2ao = np.zeros((norbitals), dtype=int) 

211 orb2uorb = np.zeros((norbitals), dtype=int) 

212 orb2occ = np.zeros((norbitals), dtype=float) 

213 

214 max_rcut = fh.readReals('d') 

215 for iorb in range(norbitals): 

216 dat = fh.readRecord() 

217 dat_size = struct.calcsize('iid') 

218 val_list = struct.unpack('iid', dat[0:dat_size]) 

219 orb2ao[iorb] = val_list[0] 

220 orb2uorb[iorb] = val_list[1] 

221 orb2occ[iorb] = val_list[2] 

222 

223 atm2sp = np.zeros((natoms), dtype=int) 

224 atm2shift = np.zeros((natoms + 1), dtype=int) 

225 for iatm in range(natoms): 

226 atm2sp[iatm] = fh.readInts('i')[0] 

227 

228 for iatm in range(natoms + 1): 

229 atm2shift[iatm] = fh.readInts('i')[0] 

230 

231 cell = np.zeros((3, 3), dtype=float) 

232 nunit_cells = np.zeros((3), dtype=int) 

233 for i in range(3): 

234 cell[i, :] = fh.readReals('d') 

235 nunit_cells = fh.readInts('i') 

236 

237 coord_sc = np.zeros((natoms, 3), dtype=float) 

238 for iatm in range(natoms): 

239 coord_sc[iatm, :] = fh.readReals('d') 

240 

241 fh.close() 

242 return PLD_tuple(max_rcut, orb2ao, orb2uorb, orb2occ, atm2sp, atm2shift, 

243 coord_sc, cell, nunit_cells) 

244 

245 

246def readWFSX(fname): 

247 """ 

248 Read unformatted siesta WFSX file 

249 """ 

250 import collections 

251 # use struct library to read mixed data type from binary 

252 import struct 

253 

254 WFSX_tuple = collections.namedtuple('WFSX', 

255 ['nkpoints', 'nspin', 'norbitals', 

256 'gamma', 'orb2atm', 'orb2strspecies', 

257 'orb2ao', 'orb2n', 'orb2strsym', 

258 'kpoints', 'DFT_E', 'DFT_X', 

259 'mo_spin_kpoint_2_is_read']) 

260 

261 fh = FortranFile(fname) 

262 

263 nkpoints, gamma = fh.readInts('i') 

264 nspin = fh.readInts('i')[0] 

265 norbitals = fh.readInts('i')[0] 

266 

267 orb2atm = np.zeros((norbitals), dtype=int) 

268 orb2strspecies = [] 

269 orb2ao = np.zeros((norbitals), dtype=int) 

270 orb2n = np.zeros((norbitals), dtype=int) 

271 orb2strsym = [] 

272 # for string list are better to select all the string length 

273 

274 dat_size = struct.calcsize('i20sii20s') 

275 dat = fh.readRecord() 

276 

277 ind_st = 0 

278 ind_fn = dat_size 

279 for iorb in range(norbitals): 

280 val_list = struct.unpack('i20sii20s', dat[ind_st:ind_fn]) 

281 orb2atm[iorb] = val_list[0] 

282 orb2strspecies.append(val_list[1]) 

283 orb2ao[iorb] = val_list[2] 

284 orb2n[iorb] = val_list[3] 

285 orb2strsym.append(val_list[4]) 

286 ind_st = ind_st + dat_size 

287 ind_fn = ind_fn + dat_size 

288 orb2strspecies = np.array(orb2strspecies) 

289 orb2strsym = np.array(orb2strsym) 

290 

291 kpoints = np.zeros((3, nkpoints), dtype=np.float64) 

292 DFT_E = np.zeros((norbitals, nspin, nkpoints), dtype=np.float64) 

293 

294 if (gamma == 1): 

295 DFT_X = np.zeros((1, norbitals, norbitals, nspin, nkpoints), 

296 dtype=np.float64) 

297 eigenvector = np.zeros((1, norbitals), dtype=float) 

298 else: 

299 DFT_X = np.zeros((2, norbitals, norbitals, nspin, nkpoints), 

300 dtype=np.float64) 

301 eigenvector = np.zeros((2, norbitals), dtype=float) 

302 

303 mo_spin_kpoint_2_is_read = np.zeros((norbitals, nspin, nkpoints), 

304 dtype=bool) 

305 mo_spin_kpoint_2_is_read[0:norbitals, 0:nspin, 0:nkpoints] = False 

306 

307 dat_size = struct.calcsize('iddd') 

308 for ikpoint in range(nkpoints): 

309 for ispin in range(nspin): 

310 dat = fh.readRecord() 

311 val_list = struct.unpack('iddd', dat[0:dat_size]) 

312 ikpoint_in = val_list[0] - 1 

313 kpoints[0:3, ikpoint] = val_list[1:4] 

314 if (ikpoint != ikpoint_in): 

315 raise ValueError('siesta_get_wfsx: ikpoint != ikpoint_in') 

316 ispin_in = fh.readInts('i')[0] - 1 

317 if (ispin_in > nspin - 1): 

318 msg = 'siesta_get_wfsx: err: ispin_in>nspin\n \ 

319 siesta_get_wfsx: ikpoint, ispin, ispin_in = \ 

320 {0} {1} {2}\n siesta_get_wfsx'.format(ikpoint, 

321 ispin, ispin_in) 

322 raise ValueError(msg) 

323 

324 norbitals_in = fh.readInts('i')[0] 

325 if (norbitals_in > norbitals): 

326 msg = 'siesta_get_wfsx: err: norbitals_in>norbitals\n \ 

327 siesta_get_wfsx: ikpoint, norbitals, norbitals_in = \ 

328 {0} {1} {2}\n siesta_get_wfsx'.format(ikpoint, 

329 norbitals, 

330 norbitals_in) 

331 raise ValueError(msg) 

332 

333 for imolecular_orb in range(norbitals_in): 

334 imolecular_orb_in = fh.readInts('i')[0] - 1 

335 if (imolecular_orb_in > norbitals - 1): 

336 msg = """ 

337 siesta_get_wfsx: err: imolecular_orb_in>norbitals\n 

338 siesta_get_wfsx: ikpoint, norbitals, 

339 imolecular_orb_in = {0} {1} {2}\n 

340 siesta_get_wfsx""".format(ikpoint, norbitals, 

341 imolecular_orb_in) 

342 raise ValueError(msg) 

343 

344 real_E_eV = fh.readReals('d')[0] 

345 eigenvector = fh.readReals('f') 

346 DFT_E[imolecular_orb_in, ispin_in, 

347 ikpoint] = real_E_eV / 13.60580 

348 DFT_X[:, :, imolecular_orb_in, ispin_in, 

349 ikpoint] = eigenvector 

350 mo_spin_kpoint_2_is_read[imolecular_orb_in, ispin_in, 

351 ikpoint] = True 

352 

353 if (not all(mo_spin_kpoint_2_is_read[:, ispin_in, ikpoint])): 

354 msg = 'siesta_get_wfsx: warn: .not. all(mo_spin_k_2_is_read)' 

355 print('mo_spin_kpoint_2_is_read = ', mo_spin_kpoint_2_is_read) 

356 raise ValueError(msg) 

357 

358 fh.close() 

359 return WFSX_tuple(nkpoints, nspin, norbitals, gamma, orb2atm, 

360 orb2strspecies, orb2ao, orb2n, orb2strsym, 

361 kpoints, DFT_E, DFT_X, mo_spin_kpoint_2_is_read)