1""" 

2The ASE Calculator for OpenMX <http://www.openmx-square.org>: Python interface 

3to the software package for nano-scale material simulations based on density 

4functional theories. 

5 Copyright (C) 2017 Charles Thomas Johnson, JaeHwan Shim and JaeJun Yu 

6 

7 This program is free software: you can redistribute it and/or modify 

8 it under the terms of the GNU Lesser General Public License as published by 

9 the Free Software Foundation, either version 2.1 of the License, or 

10 (at your option) any later version. 

11 

12 This program is distributed in the hope that it will be useful, 

13 but WITHOUT ANY WARRANTY; without even the implied warranty of 

14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

15 GNU Lesser General Public License for more details. 

16 

17 You should have received a copy of the GNU Lesser General Public License 

18 along with ASE. If not, see <http://www.gnu.org/licenses/>. 

19""" 

20import numpy as np 

21import os 

22import subprocess 

23import warnings 

24 

25from ase.calculators.openmx.reader import rn as read_nth_to_last_value 

26 

27 

28def input_command(calc, executable_name, input_files, argument_format='%s'): 

29 input_files = tuple(input_files) 

30 command = executable_name + ' ' + argument_format % input_files 

31 subprocess.check_call(command, shell=True, cwd=calc.directory) 

32 

33 

34class DOS: 

35 

36 def __init__(self, calc): 

37 self.calc = calc 

38 self.dos_dict = {} 

39 

40 def read_dos(self, method='Tetrahedron', pdos=False, atom_index=1, 

41 orbital='', spin_polarization=False): 

42 """ 

43 function for reading DOS from the following OpenMX file extensions: 

44 ~.[DOS|PDOS].[Tetrahedron|Gaussian]<.atom(int).(orbital) 

45 :param method: the method which has been used to calculate the density 

46 of states ('Tetrahedron' or 'Gaussian') 

47 :param pdos: True if the pseudo-density of states have been calculated, 

48 False if only the total density of states has been 

49 calculated 

50 :param atom_index: positive integer, n. For the nth atom in the unit 

51 cell as specified in the OpenMX input file 

52 :param orbital: '' or 's1' or 'p1', 'p2', 'p3' or 'd1', 'd2', 'd3', 

53 'd4', 'd5' etc. If pdos is True then this specifies the 

54 pdos from a particular orbital to read from. If '' is 

55 given then the total pdos from the given atom is read. 

56 :param spin_polarization: if True this will read the separate pdos for 

57 up and down spin states. 

58 :return: None 

59 """ 

60 add = False 

61 if not spin_polarization and self.calc['initial_magnetic_moments']: 

62 add = True 

63 p = '' 

64 if pdos: 

65 p = 'P' 

66 filename = self.calc.label + '.' + p + 'DOS.' + method 

67 if pdos: 

68 period = '' 

69 if orbital != '': 

70 period = '.' 

71 filename += '.atom' + str(atom_index) + period + orbital 

72 

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

74 line = '\n' 

75 number_of_lines = -1 

76 while line != '': 

77 line = fd.readline() 

78 number_of_lines += 1 

79 

80 key = '' 

81 atom_and_orbital = '' 

82 if pdos: 

83 key = 'p' 

84 atom_and_orbital = str(atom_index) + orbital 

85 key += 'dos' 

86 self.dos_dict[key + '_energies_' + atom_and_orbital] = np.ndarray( 

87 number_of_lines) 

88 if spin_polarization: 

89 self.dos_dict[key + atom_and_orbital + 'up'] = \ 

90 np.ndarray(number_of_lines) 

91 self.dos_dict[key + atom_and_orbital + 'down'] = \ 

92 np.ndarray(number_of_lines) 

93 self.dos_dict[key + '_cum_' + atom_and_orbital + 'up'] = \ 

94 np.ndarray(number_of_lines) 

95 self.dos_dict[key + '_cum_' + atom_and_orbital + 'down'] = \ 

96 np.ndarray(number_of_lines) 

97 else: 

98 self.dos_dict[key + atom_and_orbital] = np.ndarray(number_of_lines) 

99 self.dos_dict[key + '_cum_' + atom_and_orbital] = \ 

100 np.ndarray(number_of_lines) 

101 fd = open(filename, 'r') 

102 if spin_polarization: 

103 for i in range(number_of_lines): 

104 line = fd.readline() 

105 self.dos_dict[key + '_energies_' + atom_and_orbital][i] = \ 

106 read_nth_to_last_value(line, 5) 

107 self.dos_dict[key + atom_and_orbital + 'up'][i] = \ 

108 read_nth_to_last_value(line, 4) 

109 self.dos_dict[key + atom_and_orbital + 'down'][i] = \ 

110 -float(read_nth_to_last_value(line, 3)) 

111 self.dos_dict[key + '_cum_' + atom_and_orbital + 'up'][i] = \ 

112 read_nth_to_last_value(line, 2) 

113 self.dos_dict[key + '_cum_' + atom_and_orbital + 'down'][i] = \ 

114 read_nth_to_last_value(line) 

115 elif add: 

116 for i in range(number_of_lines): 

117 line = fd.readline() 

118 self.dos_dict[key + '_energies_' + atom_and_orbital][i] = \ 

119 read_nth_to_last_value(line, 5) 

120 self.dos_dict[key + atom_and_orbital][i] = \ 

121 float(read_nth_to_last_value(line, 4)) - \ 

122 float(read_nth_to_last_value(line, 3)) 

123 self.dos_dict[key + '_cum_' + atom_and_orbital][i] = \ 

124 float(read_nth_to_last_value(line, 2)) + \ 

125 float(read_nth_to_last_value(line)) 

126 else: 

127 for i in range(number_of_lines): 

128 line = fd.readline() 

129 self.dos_dict[key + '_energies_' + atom_and_orbital][i] = \ 

130 read_nth_to_last_value(line, 3) 

131 self.dos_dict[key + atom_and_orbital][i] = \ 

132 read_nth_to_last_value(line, 2) 

133 self.dos_dict[key + '_cum_' + atom_and_orbital][i] = \ 

134 read_nth_to_last_value(line) 

135 fd.close() 

136 

137 def subplot_dos(self, axis, density=True, cum=False, pdos=False, 

138 atom_index=1, orbital='', spin='', 

139 erange=(-25, 20), fermi_level=True): 

140 """ 

141 Plots a graph of (pseudo-)density of states against energy onto a given 

142 axis of a subplot. 

143 :param axis: matplotlib.pyplot.Axes object. This allows the graph to 

144 plotted on any desired axis of a plot. 

145 :param density: If True, the density of states will be plotted 

146 :param cum: If True, the cumulative (or integrated) density of states 

147 will be plotted 

148 :param pdos: If True, the pseudo-density of states will be plotted for 

149 a given atom and orbital 

150 :param atom_index: If pdos is True, atom_index specifies which atom's 

151 PDOS to plot. 

152 :param orbital: If pdos is True, orbital specifies which orbital's PDOS 

153 to plot. 

154 :param spin: If '', density of states for both spin states will be 

155 combined into one plot. If 'up' or 'down', a given spin 

156 state's PDOS will be plotted. 

157 :return: None 

158 """ 

159 p = '' 

160 bottom_index = 0 

161 atom_orbital = atom_orbital_spin = '' 

162 if pdos: 

163 p = 'p' 

164 atom_orbital += str(atom_index) + orbital 

165 atom_orbital_spin += atom_orbital + spin 

166 key = p + 'dos' 

167 density_color = 'r' 

168 cum_color = 'b' 

169 if spin == 'down': 

170 density_color = 'c' 

171 cum_color = 'm' 

172 if density and cum: 

173 axis_twin = axis.twinx() 

174 axis.plot(self.dos_dict[key + '_energies_' + atom_orbital], 

175 self.dos_dict[key + atom_orbital_spin], 

176 density_color) 

177 axis_twin.plot(self.dos_dict[key + '_energies_' + atom_orbital], 

178 self.dos_dict[key + '_cum_' + atom_orbital_spin], 

179 cum_color) 

180 max_density = max(self.dos_dict[key + atom_orbital_spin]) 

181 max_cum = max(self.dos_dict[key + '_cum_' + atom_orbital_spin]) 

182 if not max_density: 

183 max_density = 1. 

184 if not max_cum: 

185 max_cum = 1 

186 axis.set_ylim(ymax=max_density) 

187 axis_twin.set_ylim(ymax=max_cum) 

188 axis.set_ylim(ymin=0.) 

189 axis_twin.set_ylim(ymin=0.) 

190 label_index = 0 

191 yticklabels = axis.get_yticklabels() 

192 if spin == 'down': 

193 bottom_index = len(yticklabels) - 1 

194 for t in yticklabels: 

195 if label_index == bottom_index or label_index == \ 

196 len(yticklabels) // 2: 

197 t.set_color(density_color) 

198 else: 

199 t.set_visible(False) 

200 label_index += 1 

201 label_index = 0 

202 yticklabels = axis_twin.get_yticklabels() 

203 if spin == 'down': 

204 bottom_index = len(yticklabels) - 1 

205 for t in yticklabels: 

206 if label_index == bottom_index or label_index == \ 

207 len(yticklabels) // 2: 

208 t.set_color(cum_color) 

209 else: 

210 t.set_visible(False) 

211 label_index += 1 

212 if spin == 'down': 

213 axis.set_ylim(axis.get_ylim()[::-1]) 

214 axis_twin.set_ylim(axis_twin.get_ylim()[::-1]) 

215 else: 

216 color = density_color 

217 if cum: 

218 color = cum_color 

219 key += '_cum_' 

220 key += atom_orbital_spin 

221 axis.plot(self.dos_dict[p + 'dos_energies_' + atom_orbital], 

222 self.dos_dict[key], color) 

223 maximum = max(self.dos_dict[key]) 

224 if not maximum: 

225 maximum = 1. 

226 axis.set_ylim(ymax=maximum) 

227 axis.set_ylim(ymin=0.) 

228 label_index = 0 

229 yticklabels = axis.get_yticklabels() 

230 if spin == 'down': 

231 bottom_index = len(yticklabels) - 1 

232 for t in yticklabels: 

233 if label_index == bottom_index or label_index == \ 

234 len(yticklabels) // 2: 

235 t.set_color(color) 

236 else: 

237 t.set_visible(False) 

238 label_index += 1 

239 if spin == 'down': 

240 axis.set_ylim(axis.get_ylim()[::-1]) 

241 if fermi_level: 

242 axis.axvspan(erange[0], 0., color='y', alpha=0.5) 

243 

244 def plot_dos(self, density=True, cum=False, pdos=False, orbital_list=None, 

245 atom_index_list=None, spins=('up', 'down'), fermi_level=True, 

246 spin_polarization=False, erange=(-25, 20), atoms=None, 

247 method='Tetrahedron', file_format=None): 

248 """ 

249 Generates a graphical figure containing possible subplots of different 

250 PDOSs of different atoms, orbitals and spin state combinations. 

251 :param density: If True, density of states will be plotted 

252 :param cum: If True, cumulative density of states will be plotted 

253 :param pdos: If True, pseudo-density of states will be plotted for 

254 given atoms and orbitals 

255 :param atom_index_list: If pdos is True, atom_index_list specifies 

256 which atoms will have their PDOS plotted. 

257 :param orbital_list: If pdos is True, orbital_list specifies which 

258 orbitals will have their PDOS plotted. 

259 :param spins: If '' in spins, density of states for both spin states 

260 will be combined into one graph. If 'up' or 

261 'down' in spins, a given spin state's PDOS graph will be plotted. 

262 :param spin_polarization: If spin_polarization is False then spin 

263 states will not be separated in different 

264 PDOS's. 

265 :param erange: range of energies to view DOS 

266 :return: matplotlib.figure.Figure and matplotlib.axes.Axes object 

267 """ 

268 import matplotlib.pyplot as plt 

269 from matplotlib.lines import Line2D 

270 if not spin_polarization: 

271 spins = [''] 

272 number_of_spins = len(spins) 

273 if orbital_list is None: 

274 orbital_list = [''] 

275 number_of_atoms = 1 

276 number_of_orbitals = 1 

277 p = '' 

278 if pdos: 

279 p = 'P' 

280 if atom_index_list is None: 

281 atom_index_list = [i + 1 for i in range(len(atoms))] 

282 number_of_atoms = len(atom_index_list) 

283 number_of_orbitals = len(orbital_list) 

284 figure, axes = plt.subplots(number_of_orbitals * number_of_spins, 

285 number_of_atoms, sharex=True, sharey=False, 

286 squeeze=False) 

287 for i in range(number_of_orbitals): 

288 for s in range(number_of_spins): 

289 row_index = i * number_of_spins + s 

290 for j in range(number_of_atoms): 

291 self.subplot_dos(fermi_level=fermi_level, density=density, 

292 axis=axes[row_index][j], erange=erange, 

293 atom_index=atom_index_list[j], pdos=pdos, 

294 orbital=orbital_list[i], spin=spins[s], 

295 cum=cum) 

296 if j == 0 and pdos: 

297 orbital = orbital_list[i] 

298 if orbital == '': 

299 orbital = 'All' 

300 if spins[s]: 

301 orbital += ' ' + spins[s] 

302 axes[row_index][j].set_ylabel(orbital) 

303 if row_index == 0 and pdos: 

304 atom_symbol = '' 

305 if atoms: 

306 atom_symbol = ' (' + \ 

307 atoms[atom_index_list[j]].symbol + ')' 

308 axes[row_index][j].set_title( 

309 'Atom ' + str(atom_index_list[j]) + atom_symbol) 

310 if row_index == number_of_orbitals * number_of_spins - 1: 

311 axes[row_index][j].set_xlabel( 

312 'Energy above Fermi Level (eV)') 

313 plt.xlim(xmin=erange[0], xmax=erange[1]) 

314 if density and cum: 

315 figure.suptitle(self.calc.label) 

316 xdata = (0., 1.) 

317 ydata = (0., 0.) 

318 key_tuple = (Line2D(color='r', xdata=xdata, ydata=ydata), 

319 Line2D(color='b', xdata=xdata, ydata=ydata)) 

320 if spin_polarization: 

321 key_tuple = (Line2D(color='r', xdata=xdata, ydata=ydata), 

322 Line2D(color='b', xdata=xdata, ydata=ydata), 

323 Line2D(color='c', xdata=xdata, ydata=ydata), 

324 Line2D(color='m', xdata=xdata, ydata=ydata)) 

325 title_tuple = (p + 'DOS (eV^-1)', 'Number of States per Unit Cell') 

326 if spin_polarization: 

327 title_tuple = (p + 'DOS (eV^-1), spin up', 

328 'Number of States per Unit Cell, spin up', 

329 p + 'DOS (eV^-1), spin down', 

330 'Number of States per Unit Cell, spin down') 

331 figure.legend(key_tuple, title_tuple, 'lower center') 

332 elif density: 

333 figure.suptitle(self.calc.prefix + ': ' + p + 'DOS (eV^-1)') 

334 elif cum: 

335 figure.suptitle(self.calc.prefix + ': Number of States') 

336 extra_margin = 0 

337 if density and cum and spin_polarization: 

338 extra_margin = 0.1 

339 plt.subplots_adjust(hspace=0., bottom=0.2 + extra_margin, wspace=0.29, 

340 left=0.09, right=0.95) 

341 if file_format: 

342 orbitals = '' 

343 if pdos: 

344 atom_index_list = map(str, atom_index_list) 

345 atoms = '&'.join(atom_index_list) 

346 if '' in orbital_list: 

347 all_index = orbital_list.index('') 

348 orbital_list.remove('') 

349 orbital_list.insert(all_index, 'all') 

350 orbitals = ''.join(orbital_list) 

351 plt.savefig(filename=self.calc.label + '.' + p + 'DOS.' + 

352 method + '.atoms' + atoms + '.' + orbitals + '.' + 

353 file_format) 

354 if not file_format: 

355 plt.show() 

356 return figure, axes 

357 

358 def calc_dos(self, method='Tetrahedron', pdos=False, gaussian_width=0.1, 

359 atom_index_list=None): 

360 """ 

361 Python interface for DosMain (OpenMX's density of states calculator). 

362 Can automate the density of states 

363 calculations used in OpenMX by processing .Dos.val and .Dos.vec files. 

364 :param method: method to be used to calculate the density of states 

365 from eigenvalues and eigenvectors. 

366 ('Tetrahedron' or 'Gaussian') 

367 :param pdos: If True, the pseudo-density of states is calculated for a 

368 given list of atoms for each orbital. If the system is 

369 spin polarized, then each up and down state is also 

370 calculated. 

371 :param gaussian_width: If the method is 'Gaussian' then gaussian_width 

372 is required (eV). 

373 :param atom_index_list: If pdos is True, a list of atom indices are 

374 required to generate the pdos of each of those 

375 specified atoms. 

376 :return: None 

377 """ 

378 method_code = '2\n' 

379 if method == 'Tetrahedron': 

380 method_code = '1\n' 

381 pdos_code = '1\n' 

382 if pdos: 

383 pdos_code = '2\n' 

384 with open(os.path.join(self.calc.directory, 'std_dos.in'), 'w') as fd: 

385 fd.write(method_code) 

386 if method == 'Gaussian': 

387 fd.write(str(gaussian_width) + '\n') 

388 fd.write(pdos_code) 

389 if pdos: 

390 atoms_code = '' 

391 if atom_index_list is None: 

392 for i in range(len(self.calc.atoms)): 

393 atoms_code += str(i + 1) + ' ' 

394 else: 

395 for i in atom_index_list: 

396 atoms_code += str(i) + ' ' 

397 atoms_code += '\n' 

398 fd.write(atoms_code) 

399 fd.close() 

400 executable_name = 'DosMain' 

401 input_files = (self.calc.label + '.Dos.val', self.calc.label + 

402 '.Dos.vec', os.path.join(self.calc.directory, 

403 'std_dos.in')) 

404 argument_format = '%s %s < %s' 

405 input_command(self.calc, executable_name, input_files, argument_format) 

406 

407 def get_dos(self, atom_index_list=None, method='Tetrahedron', 

408 gaussian_width=0.1, pdos=False, orbital_list=None, 

409 spin_polarization=None, density=True, cum=False, 

410 erange=(-25, 20), file_format=None, atoms=None, 

411 fermi_level=True): 

412 """ 

413 Wraps all the density of states processing functions. Can go from 

414 .Dos.val and .Dos.vec files to a graphical figure showing many 

415 different PDOS plots against energy in one step. 

416 :param atom_index_list: 

417 :param method: method to be used to calculate the density of states 

418 from eigenvalues and eigenvectors. 

419 ('Tetrahedron' or 'Gaussian') 

420 :param gaussian_width: If the method is 'Gaussian' then gaussian_width 

421 is required (eV). 

422 :param pdos: If True, the pseudo-density of states is calculated for a 

423 given list of atoms for each orbital. If the system is 

424 spin polarized, then each up and down state is also 

425 calculated. 

426 :param orbital_list: If pdos is True, a list of atom indices are 

427 required to generate the pdos of each of those 

428 specified atoms. 

429 :param spin_polarization: If spin_polarization is False then spin 

430 states will not be separated in different 

431 PDOS's. 

432 :param density: If True, density of states will be plotted 

433 :param cum: If True, cumulative (or integrated) density of states will 

434 be plotted 

435 :param erange: range of energies to view the DOS 

436 :param file_format: If not None, a file will be saved automatically in 

437 that format ('pdf', 'png', 'jpeg' etc.) 

438 :return: matplotlib.figure.Figure object 

439 """ 

440 if spin_polarization is None: 

441 spin_polarization = bool(self.calc['initial_magnetic_moments']) 

442 if spin_polarization and not self.calc['initial_magnetic_moments']: 

443 warnings.warn('No spin polarization calculations provided') 

444 spin_polarization = False 

445 if atom_index_list is None: 

446 atom_index_list = [1] 

447 if method == 'Tetrahedron' and self.calc['dos_kgrid'] == (1, 1, 1): 

448 raise ValueError('Not enough k-space grid points.') 

449 self.calc_dos(atom_index_list=atom_index_list, pdos=pdos, 

450 method=method, gaussian_width=gaussian_width) 

451 if pdos: 

452 if orbital_list is None: 

453 orbital_list = [''] 

454 orbital_list = list(orbital_list) 

455 if 's' in orbital_list: 

456 s_index = orbital_list.index('s') 

457 orbital_list.remove('s') 

458 orbital_list.insert(s_index, 's1') 

459 if 'p' in orbital_list: 

460 p_index = orbital_list.index('p') 

461 orbital_list.remove('p') 

462 orbital_list.insert(p_index, 'p3') 

463 orbital_list.insert(p_index, 'p2') 

464 orbital_list.insert(p_index, 'p1') 

465 if 'd' in orbital_list: 

466 d_index = orbital_list.index('d') 

467 orbital_list.remove('d') 

468 orbital_list.insert(d_index, 'd5') 

469 orbital_list.insert(d_index, 'd4') 

470 orbital_list.insert(d_index, 'd3') 

471 orbital_list.insert(d_index, 'd2') 

472 orbital_list.insert(d_index, 'd1') 

473 if 'f' in orbital_list: 

474 f_index = orbital_list.index('f') 

475 orbital_list.remove('f') 

476 orbital_list.insert(f_index, 'f7') 

477 orbital_list.insert(f_index, 'f6') 

478 orbital_list.insert(f_index, 'f5') 

479 orbital_list.insert(f_index, 'f4') 

480 orbital_list.insert(f_index, 'f3') 

481 orbital_list.insert(f_index, 'f2') 

482 orbital_list.insert(f_index, 'f1') 

483 

484 for atom_index in atom_index_list: 

485 for orbital in orbital_list: 

486 self.read_dos(method=method, atom_index=atom_index, 

487 pdos=pdos, orbital=orbital, 

488 spin_polarization=spin_polarization) 

489 else: 

490 self.read_dos(method=method, spin_polarization=spin_polarization) 

491 return self.plot_dos(density=density, cum=cum, atoms=atoms, 

492 atom_index_list=atom_index_list, pdos=pdos, 

493 orbital_list=orbital_list, erange=erange, 

494 spin_polarization=spin_polarization, 

495 file_format=file_format, method=method, 

496 fermi_level=fermi_level)