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

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.

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.

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.

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 os

21import subprocess

22import warnings

24import numpy as np

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

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

30 input_files = tuple(input_files)

31 command = executable_name + ' ' + argument_format % input_files

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

35class DOS:

37 def __init__(self, calc):

38 self.calc = calc

39 self.dos_dict = {}

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

42 orbital='', spin_polarization=False):

43 """

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

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

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

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

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

49 False if only the total density of states has been

50 calculated

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

52 cell as specified in the OpenMX input file

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

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

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

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

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

58 up and down spin states.

59 :return: None

60 """

61 add = False

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

63 add = True

64 p = ''

65 if pdos:

66 p = 'P'

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

68 if pdos:

69 period = ''

70 if orbital != '':

71 period = '.'

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

74 with open(filename) as fd:

75 line = '\n'

76 number_of_lines = -1

77 while line != '':

78 line = fd.readline()

79 number_of_lines += 1

81 key = ''

82 atom_and_orbital = ''

83 if pdos:

84 key = 'p'

85 atom_and_orbital = str(atom_index) + orbital

86 key += 'dos'

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

88 number_of_lines)

89 if spin_polarization:

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

91 np.ndarray(number_of_lines)

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

93 np.ndarray(number_of_lines)

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

95 np.ndarray(number_of_lines)

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

97 np.ndarray(number_of_lines)

98 else:

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

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

101 np.ndarray(number_of_lines)

102 fd = open(filename)

103 if spin_polarization:

104 for i in range(number_of_lines):

105 line = fd.readline()

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

107 read_nth_to_last_value(line, 5)

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

109 read_nth_to_last_value(line, 4)

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

111 -float(read_nth_to_last_value(line, 3))

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

113 read_nth_to_last_value(line, 2)

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

115 read_nth_to_last_value(line)

116 elif add:

117 for i in range(number_of_lines):

118 line = fd.readline()

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

120 read_nth_to_last_value(line, 5)

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

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

123 float(read_nth_to_last_value(line, 3))

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

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

126 float(read_nth_to_last_value(line))

127 else:

128 for i in range(number_of_lines):

129 line = fd.readline()

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

131 read_nth_to_last_value(line, 3)

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

133 read_nth_to_last_value(line, 2)

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

135 read_nth_to_last_value(line)

136 fd.close()

137

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

139 atom_index=1, orbital='', spin='',

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

141 """

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

143 axis of a subplot.

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

145 plotted on any desired axis of a plot.

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

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

148 will be plotted

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

150 a given atom and orbital

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

152 PDOS to plot.

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

154 to plot.

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

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

157 state's PDOS will be plotted.

158 :return: None

159 """

160 p = ''

161 bottom_index = 0

162 atom_orbital = atom_orbital_spin = ''

163 if pdos:

164 p = 'p'

165 atom_orbital += str(atom_index) + orbital

166 atom_orbital_spin += atom_orbital + spin

167 key = p + 'dos'

168 density_color = 'r'

169 cum_color = 'b'

170 if spin == 'down':

171 density_color = 'c'

172 cum_color = 'm'

173 if density and cum:

174 axis_twin = axis.twinx()

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

176 self.dos_dict[key + atom_orbital_spin],

177 density_color)

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

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

180 cum_color)

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

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

183 if not max_density:

184 max_density = 1.

185 if not max_cum:

186 max_cum = 1

187 axis.set_ylim(ymax=max_density)

188 axis_twin.set_ylim(ymax=max_cum)

189 axis.set_ylim(ymin=0.)

190 axis_twin.set_ylim(ymin=0.)

191 label_index = 0

192 yticklabels = axis.get_yticklabels()

193 if spin == 'down':

194 bottom_index = len(yticklabels) - 1

195 for t in yticklabels:

196 if label_index == bottom_index or label_index == \

197 len(yticklabels) // 2:

198 t.set_color(density_color)

199 else:

200 t.set_visible(False)

201 label_index += 1

202 label_index = 0

203 yticklabels = axis_twin.get_yticklabels()

204 if spin == 'down':

205 bottom_index = len(yticklabels) - 1

206 for t in yticklabels:

207 if label_index == bottom_index or label_index == \

208 len(yticklabels) // 2:

209 t.set_color(cum_color)

210 else:

211 t.set_visible(False)

212 label_index += 1

213 if spin == 'down':

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

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

216 else:

217 color = density_color

218 if cum:

219 color = cum_color

220 key += '_cum_'

221 key += atom_orbital_spin

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

223 self.dos_dict[key], color)

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

225 if not maximum:

226 maximum = 1.

227 axis.set_ylim(ymax=maximum)

228 axis.set_ylim(ymin=0.)

229 label_index = 0

230 yticklabels = axis.get_yticklabels()

231 if spin == 'down':

232 bottom_index = len(yticklabels) - 1

233 for t in yticklabels:

234 if label_index == bottom_index or label_index == \

235 len(yticklabels) // 2:

236 t.set_color(color)

237 else:

238 t.set_visible(False)

239 label_index += 1

240 if spin == 'down':

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

242 if fermi_level:

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

244

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

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

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

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

249 """

250 Generates a graphical figure containing possible subplots of different

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

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

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

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

255 given atoms and orbitals

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

257 which atoms will have their PDOS plotted.

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

259 orbitals will have their PDOS plotted.

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

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

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

263 :param spin_polarization: If spin_polarization is False then spin

264 states will not be separated in different

265 PDOS's.

266 :param erange: range of energies to view DOS

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

268 """

269 import matplotlib.pyplot as plt

270 from matplotlib.lines import Line2D

271 if not spin_polarization:

272 spins = ['']

273 number_of_spins = len(spins)

274 if orbital_list is None:

275 orbital_list = ['']

276 number_of_atoms = 1

277 number_of_orbitals = 1

278 p = ''

279 if pdos:

280 p = 'P'

281 if atom_index_list is None:

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

283 number_of_atoms = len(atom_index_list)

284 number_of_orbitals = len(orbital_list)

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

286 number_of_atoms, sharex=True, sharey=False,

287 squeeze=False)

288 for i in range(number_of_orbitals):

289 for s in range(number_of_spins):

290 row_index = i * number_of_spins + s

291 for j in range(number_of_atoms):

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

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

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

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

296 cum=cum)

297 if j == 0 and pdos:

298 orbital = orbital_list[i]

299 if orbital == '':

300 orbital = 'All'

301 if spins[s]:

302 orbital += ' ' + spins[s]

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

304 if row_index == 0 and pdos:

305 atom_symbol = ''

306 if atoms:

307 atom_symbol = ' (' + \

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

309 axes[row_index][j].set_title(

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

311 if row_index == number_of_orbitals * number_of_spins - 1:

312 axes[row_index][j].set_xlabel(

313 'Energy above Fermi Level (eV)')

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

315 if density and cum:

316 figure.suptitle(self.calc.label)

317 xdata = (0., 1.)

318 ydata = (0., 0.)

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

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

321 if spin_polarization:

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

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

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

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

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

327 if spin_polarization:

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

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

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

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

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

333 elif density:

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

335 elif cum:

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

337 extra_margin = 0

338 if density and cum and spin_polarization:

339 extra_margin = 0.1

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

341 left=0.09, right=0.95)

342 if file_format:

343 orbitals = ''

344 if pdos:

345 atom_index_list = map(str, atom_index_list)

346 atoms = '&'.join(atom_index_list)

347 if '' in orbital_list:

348 all_index = orbital_list.index('')

349 orbital_list.remove('')

350 orbital_list.insert(all_index, 'all')

351 orbitals = ''.join(orbital_list)

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

353 method + '.atoms' + atoms + '.' + orbitals + '.' +

354 file_format)

355 if not file_format:

356 plt.show()

357 return figure, axes

358

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

360 atom_index_list=None):

361 """

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

363 Can automate the density of states

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

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

366 from eigenvalues and eigenvectors.

367 ('Tetrahedron' or 'Gaussian')

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

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

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

371 calculated.

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

373 is required (eV).

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

375 required to generate the pdos of each of those

376 specified atoms.

377 :return: None

378 """

379 method_code = '2\n'

380 if method == 'Tetrahedron':

381 method_code = '1\n'

382 pdos_code = '1\n'

383 if pdos:

384 pdos_code = '2\n'

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

386 fd.write(method_code)

387 if method == 'Gaussian':

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

389 fd.write(pdos_code)

390 if pdos:

391 atoms_code = ''

392 if atom_index_list is None:

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

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

395 else:

396 for i in atom_index_list:

397 atoms_code += str(i) + ' '

398 atoms_code += '\n'

399 fd.write(atoms_code)

400 fd.close()

401 executable_name = 'DosMain'

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

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

404 'std_dos.in'))

405 argument_format = '%s %s < %s'

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

407

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

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

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

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

412 fermi_level=True):

413 """

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

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

416 different PDOS plots against energy in one step.

417 :param atom_index_list:

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

419 from eigenvalues and eigenvectors.

420 ('Tetrahedron' or 'Gaussian')

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

422 is required (eV).

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

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

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

426 calculated.

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

428 required to generate the pdos of each of those

429 specified atoms.

430 :param spin_polarization: If spin_polarization is False then spin

431 states will not be separated in different

432 PDOS's.

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

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

435 be plotted

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

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

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

439 :return: matplotlib.figure.Figure object

440 """

441 if spin_polarization is None:

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

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

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

445 spin_polarization = False

446 if atom_index_list is None:

447 atom_index_list = [1]

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

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

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

451 method=method, gaussian_width=gaussian_width)

452 if pdos:

453 if orbital_list is None:

454 orbital_list = ['']

455 orbital_list = list(orbital_list)

456 if 's' in orbital_list:

457 s_index = orbital_list.index('s')

458 orbital_list.remove('s')

459 orbital_list.insert(s_index, 's1')

460 if 'p' in orbital_list:

461 p_index = orbital_list.index('p')

462 orbital_list.remove('p')

463 orbital_list.insert(p_index, 'p3')

464 orbital_list.insert(p_index, 'p2')

465 orbital_list.insert(p_index, 'p1')

466 if 'd' in orbital_list:

467 d_index = orbital_list.index('d')

468 orbital_list.remove('d')

469 orbital_list.insert(d_index, 'd5')

470 orbital_list.insert(d_index, 'd4')

471 orbital_list.insert(d_index, 'd3')

472 orbital_list.insert(d_index, 'd2')

473 orbital_list.insert(d_index, 'd1')

474 if 'f' in orbital_list:

475 f_index = orbital_list.index('f')

476 orbital_list.remove('f')

477 orbital_list.insert(f_index, 'f7')

478 orbital_list.insert(f_index, 'f6')

479 orbital_list.insert(f_index, 'f5')

480 orbital_list.insert(f_index, 'f4')

481 orbital_list.insert(f_index, 'f3')

482 orbital_list.insert(f_index, 'f2')

483 orbital_list.insert(f_index, 'f1')

484

485 for atom_index in atom_index_list:

486 for orbital in orbital_list:

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

488 pdos=pdos, orbital=orbital,

489 spin_polarization=spin_polarization)

490 else:

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

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

493 atom_index_list=atom_index_list, pdos=pdos,

494 orbital_list=orbital_list, erange=erange,

495 spin_polarization=spin_polarization,

496 file_format=file_format, method=method,

497 fermi_level=fermi_level)