Coverage for /builds/ase/ase/ase/calculators/openmx/dos.py: 4.58%

1# fmt: off

3"""

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

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

6functional theories.

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

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

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

12 (at your option) any later version.

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

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

16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

17 GNU Lesser General Public License for more details.

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

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

21"""

22import os

23import subprocess

24import warnings

26import numpy as np

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

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

32 input_files = tuple(input_files)

33 command = executable_name + ' ' + argument_format % input_files

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

37class DOS:

39 def __init__(self, calc):

40 self.calc = calc

41 self.dos_dict = {}

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

44 orbital='', spin_polarization=False):

45 """

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

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

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

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

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

51 False if only the total density of states has been

52 calculated

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

54 cell as specified in the OpenMX input file

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

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

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

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

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

60 up and down spin states.

61 :return: None

62 """

63 add = False

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

65 add = True

66 p = ''

67 if pdos:

68 p = 'P'

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

70 if pdos:

71 period = ''

72 if orbital != '':

73 period = '.'

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

76 with open(filename) as fd:

77 line = '\n'

78 number_of_lines = -1

79 while line != '':

80 line = fd.readline()

81 number_of_lines += 1

83 key = ''

84 atom_and_orbital = ''

85 if pdos:

86 key = 'p'

87 atom_and_orbital = str(atom_index) + orbital

88 key += 'dos'

89 self.dos_dict[f'{key}_energies_{atom_and_orbital}'] = np.ndarray(

90 number_of_lines)

91 if spin_polarization:

92 self.dos_dict[f'{key}{atom_and_orbital}up'] = \

93 np.ndarray(number_of_lines)

94 self.dos_dict[f'{key}{atom_and_orbital}down'] = \

95 np.ndarray(number_of_lines)

96 self.dos_dict[f'{key}_cum_{atom_and_orbital}up'] = \

97 np.ndarray(number_of_lines)

98 self.dos_dict[f'{key}_cum_{atom_and_orbital}down'] = \

99 np.ndarray(number_of_lines)

100 else:

101 self.dos_dict[f'{key}{atom_and_orbital}'] = \

102 np.ndarray(number_of_lines)

103 self.dos_dict[f'{key}_cum_{atom_and_orbital}'] = \

104 np.ndarray(number_of_lines)

105 with open(filename) as fd:

106 if spin_polarization:

107 for i in range(number_of_lines):

108 line = fd.readline()

109 self.dos_dict[f'{key}_energies_{atom_and_orbital}'][i] = \

110 read_nth_to_last_value(line, 5)

111 self.dos_dict[f'{key}{atom_and_orbital}up'][i] = \

112 read_nth_to_last_value(line, 4)

113 self.dos_dict[f'{key}{atom_and_orbital}down'][i] = \

114 -float(read_nth_to_last_value(line, 3))

115 self.dos_dict[f'{key}_cum_{atom_and_orbital}up'][i] = \

116 read_nth_to_last_value(line, 2)

117 self.dos_dict[f'{key}_cum_{atom_and_orbital}down'][i] = \

118 read_nth_to_last_value(line)

119 elif add:

120 for i in range(number_of_lines):

121 line = fd.readline()

122 self.dos_dict[f'{key}_energies_{atom_and_orbital}'][i] = \

123 read_nth_to_last_value(line, 5)

124 self.dos_dict[f'{key}{atom_and_orbital}'][i] = \

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

126 float(read_nth_to_last_value(line, 3))

127 self.dos_dict[f'{key}_cum_{atom_and_orbital}'][i] = \

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

129 float(read_nth_to_last_value(line))

130 else:

131 for i in range(number_of_lines):

132 line = fd.readline()

133 self.dos_dict[f'{key}_energies_{atom_and_orbital}'][i] = \

134 read_nth_to_last_value(line, 3)

135 self.dos_dict[f'{key}{atom_and_orbital}'][i] = \

136 read_nth_to_last_value(line, 2)

137 self.dos_dict[f'{key}_cum_{atom_and_orbital}'][i] = \

138 read_nth_to_last_value(line)

139

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

141 atom_index=1, orbital='', spin='',

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

143 """

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

145 axis of a subplot.

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

147 plotted on any desired axis of a plot.

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

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

150 will be plotted

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

152 a given atom and orbital

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

154 PDOS to plot.

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

156 to plot.

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

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

159 state's PDOS will be plotted.

160 :return: None

161 """

162 p = ''

163 bottom_index = 0

164 atom_orbital = atom_orbital_spin = ''

165 if pdos:

166 p = 'p'

167 atom_orbital += str(atom_index) + orbital

168 atom_orbital_spin += atom_orbital + spin

169 key = p + 'dos'

170 density_color = 'r'

171 cum_color = 'b'

172 if spin == 'down':

173 density_color = 'c'

174 cum_color = 'm'

175 if density and cum:

176 axis_twin = axis.twinx()

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

178 self.dos_dict[key + atom_orbital_spin],

179 density_color)

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

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

182 cum_color)

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

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

185 if not max_density:

186 max_density = 1.

187 if not max_cum:

188 max_cum = 1

189 axis.set_ylim(ymax=max_density)

190 axis_twin.set_ylim(ymax=max_cum)

191 axis.set_ylim(ymin=0.)

192 axis_twin.set_ylim(ymin=0.)

193 label_index = 0

194 yticklabels = axis.get_yticklabels()

195 if spin == 'down':

196 bottom_index = len(yticklabels) - 1

197 for t in yticklabels:

198 if label_index == bottom_index or label_index == \

199 len(yticklabels) // 2:

200 t.set_color(density_color)

201 else:

202 t.set_visible(False)

203 label_index += 1

204 label_index = 0

205 yticklabels = axis_twin.get_yticklabels()

206 if spin == 'down':

207 bottom_index = len(yticklabels) - 1

208 for t in yticklabels:

209 if label_index == bottom_index or label_index == \

210 len(yticklabels) // 2:

211 t.set_color(cum_color)

212 else:

213 t.set_visible(False)

214 label_index += 1

215 if spin == 'down':

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

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

218 else:

219 color = density_color

220 if cum:

221 color = cum_color

222 key += '_cum_'

223 key += atom_orbital_spin

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

225 self.dos_dict[key], color)

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

227 if not maximum:

228 maximum = 1.

229 axis.set_ylim(ymax=maximum)

230 axis.set_ylim(ymin=0.)

231 label_index = 0

232 yticklabels = axis.get_yticklabels()

233 if spin == 'down':

234 bottom_index = len(yticklabels) - 1

235 for t in yticklabels:

236 if label_index == bottom_index or label_index == \

237 len(yticklabels) // 2:

238 t.set_color(color)

239 else:

240 t.set_visible(False)

241 label_index += 1

242 if spin == 'down':

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

244 if fermi_level:

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

246

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

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

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

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

251 """

252 Generates a graphical figure containing possible subplots of different

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

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

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

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

257 given atoms and orbitals

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

259 which atoms will have their PDOS plotted.

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

261 orbitals will have their PDOS plotted.

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

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

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

265 :param spin_polarization: If spin_polarization is False then spin

266 states will not be separated in different

267 PDOS's.

268 :param erange: range of energies to view DOS

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

270 """

271 import matplotlib.pyplot as plt

272 from matplotlib.lines import Line2D

273 if not spin_polarization:

274 spins = ['']

275 number_of_spins = len(spins)

276 if orbital_list is None:

277 orbital_list = ['']

278 number_of_atoms = 1

279 number_of_orbitals = 1

280 p = ''

281 if pdos:

282 p = 'P'

283 if atom_index_list is None:

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

285 number_of_atoms = len(atom_index_list)

286 number_of_orbitals = len(orbital_list)

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

288 number_of_atoms, sharex=True, sharey=False,

289 squeeze=False)

290 for i in range(number_of_orbitals):

291 for s in range(number_of_spins):

292 row_index = i * number_of_spins + s

293 for j in range(number_of_atoms):

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

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

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

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

298 cum=cum)

299 if j == 0 and pdos:

300 orbital = orbital_list[i]

301 if orbital == '':

302 orbital = 'All'

303 if spins[s]:

304 orbital += ' ' + spins[s]

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

306 if row_index == 0 and pdos:

307 atom_symbol = ''

308 if atoms:

309 atom_symbol = ' (' + \

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

311 axes[row_index][j].set_title(

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

313 if row_index == number_of_orbitals * number_of_spins - 1:

314 axes[row_index][j].set_xlabel(

315 'Energy above Fermi Level (eV)')

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

317 if density and cum:

318 figure.suptitle(self.calc.label)

319 xdata = (0., 1.)

320 ydata = (0., 0.)

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

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

323 if spin_polarization:

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

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

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

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

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

329 if spin_polarization:

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

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

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

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

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

335 elif density:

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

337 elif cum:

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

339 extra_margin = 0

340 if density and cum and spin_polarization:

341 extra_margin = 0.1

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

343 left=0.09, right=0.95)

344 if file_format:

345 orbitals = ''

346 if pdos:

347 atom_index_list = map(str, atom_index_list)

348 atoms = '&'.join(atom_index_list)

349 if '' in orbital_list:

350 all_index = orbital_list.index('')

351 orbital_list.remove('')

352 orbital_list.insert(all_index, 'all')

353 orbitals = ''.join(orbital_list)

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

355 method + '.atoms' + atoms + '.' + orbitals + '.' +

356 file_format)

357 if not file_format:

358 plt.show()

359 return figure, axes

360

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

362 atom_index_list=None):

363 """

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

365 Can automate the density of states

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

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

368 from eigenvalues and eigenvectors.

369 ('Tetrahedron' or 'Gaussian')

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

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

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

373 calculated.

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

375 is required (eV).

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

377 required to generate the pdos of each of those

378 specified atoms.

379 :return: None

380 """

381 method_code = '2\n'

382 if method == 'Tetrahedron':

383 method_code = '1\n'

384 pdos_code = '1\n'

385 if pdos:

386 pdos_code = '2\n'

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

388 fd.write(method_code)

389 if method == 'Gaussian':

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

391 fd.write(pdos_code)

392 if pdos:

393 atoms_code = ''

394 if atom_index_list is None:

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

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

397 else:

398 for i in atom_index_list:

399 atoms_code += str(i) + ' '

400 atoms_code += '\n'

401 fd.write(atoms_code)

402 fd.close()

403 executable_name = 'DosMain'

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

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

406 'std_dos.in'))

407 argument_format = '%s %s < %s'

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

409

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

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

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

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

414 fermi_level=True):

415 """

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

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

418 different PDOS plots against energy in one step.

419 :param atom_index_list:

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

421 from eigenvalues and eigenvectors.

422 ('Tetrahedron' or 'Gaussian')

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

424 is required (eV).

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

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

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

428 calculated.

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

430 required to generate the pdos of each of those

431 specified atoms.

432 :param spin_polarization: If spin_polarization is False then spin

433 states will not be separated in different

434 PDOS's.

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

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

437 be plotted

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

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

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

441 :return: matplotlib.figure.Figure object

442 """

443 if spin_polarization is None:

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

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

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

447 spin_polarization = False

448 if atom_index_list is None:

449 atom_index_list = [1]

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

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

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

453 method=method, gaussian_width=gaussian_width)

454 if pdos:

455 if orbital_list is None:

456 orbital_list = ['']

457 orbital_list = list(orbital_list)

458 if 's' in orbital_list:

459 s_index = orbital_list.index('s')

460 orbital_list.remove('s')

461 orbital_list.insert(s_index, 's1')

462 if 'p' in orbital_list:

463 p_index = orbital_list.index('p')

464 orbital_list.remove('p')

465 orbital_list.insert(p_index, 'p3')

466 orbital_list.insert(p_index, 'p2')

467 orbital_list.insert(p_index, 'p1')

468 if 'd' in orbital_list:

469 d_index = orbital_list.index('d')

470 orbital_list.remove('d')

471 orbital_list.insert(d_index, 'd5')

472 orbital_list.insert(d_index, 'd4')

473 orbital_list.insert(d_index, 'd3')

474 orbital_list.insert(d_index, 'd2')

475 orbital_list.insert(d_index, 'd1')

476 if 'f' in orbital_list:

477 f_index = orbital_list.index('f')

478 orbital_list.remove('f')

479 orbital_list.insert(f_index, 'f7')

480 orbital_list.insert(f_index, 'f6')

481 orbital_list.insert(f_index, 'f5')

482 orbital_list.insert(f_index, 'f4')

483 orbital_list.insert(f_index, 'f3')

484 orbital_list.insert(f_index, 'f2')

485 orbital_list.insert(f_index, 'f1')

486

487 for atom_index in atom_index_list:

488 for orbital in orbital_list:

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

490 pdos=pdos, orbital=orbital,

491 spin_polarization=spin_polarization)

492 else:

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

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

495 atom_index_list=atom_index_list, pdos=pdos,

496 orbital_list=orbital_list, erange=erange,

497 spin_polarization=spin_polarization,

498 file_format=file_format, method=method,

499 fermi_level=fermi_level)