Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

# coding: utf-8 

# Copyright (c) Pymatgen Development Team. 

# Distributed under the terms of the MIT License. 

 

from __future__ import division, unicode_literals 

import numpy as np 

 

import six 

 

from pymatgen.electronic_structure.core import Spin, Orbital 

from pymatgen.core.periodic_table import get_el_sp 

from pymatgen.core.structure import Structure 

from pymatgen.util.coord_utils import get_linear_interpolated_value 

from monty.json import MSONable 

 

""" 

This module defines classes to represent the density of states, etc. 

""" 

 

__author__ = "Shyue Ping Ong" 

__copyright__ = "Copyright 2012, The Materials Project" 

__version__ = "2.0" 

__maintainer__ = "Shyue Ping Ong" 

__email__ = "shyuep@gmail.com" 

__date__ = "Mar 20, 2012" 

 

 

class Dos(MSONable): 

""" 

Basic DOS object. All other DOS objects are extended versions of this 

object. 

 

Args: 

efermi: Fermi level energy 

energies: A sequences of energies 

densities ({Spin: np.array}): representing the density of states 

for each Spin. 

 

.. attribute: energies 

 

The sequence of energies 

 

.. attribute: densities 

 

A dict of spin densities, e.g., {Spin.up: [...], Spin.down: [...]} 

 

.. attribute: efermi 

 

Fermi level 

""" 

 

def __init__(self, efermi, energies, densities): 

self.efermi = efermi 

self.energies = np.array(energies) 

self.densities = {k: np.array(d) for k, d in densities.items()} 

 

def get_densities(self, spin=None): 

""" 

Returns the density of states for a particular spin. 

 

Args: 

spin: Spin 

 

Returns: 

Returns the density of states for a particular spin. If Spin is 

None, the sum of all spins is returned. 

""" 

if self.densities is None: 

result = None 

elif spin is None: 

if Spin.down in self.densities: 

result = self.densities[Spin.up] + self.densities[Spin.down] 

else: 

result = self.densities[Spin.up] 

else: 

result = self.densities[spin] 

return result 

 

def get_smeared_densities(self, sigma): 

""" 

Returns the Dict representation of the densities, {Spin: densities}, 

but with a Gaussian smearing of std dev sigma applied about the fermi 

level. 

 

Args: 

sigma: Std dev of Gaussian smearing function. 

 

Returns: 

Dict of Gaussian-smeared densities. 

""" 

from scipy.ndimage.filters import gaussian_filter1d 

smeared_dens = {} 

diff = [self.energies[i + 1] - self.energies[i] 

for i in range(len(self.energies) - 1)] 

avgdiff = sum(diff) / len(diff) 

for spin, dens in self.densities.items(): 

smeared_dens[spin] = gaussian_filter1d(dens, sigma / avgdiff) 

return smeared_dens 

 

def __add__(self, other): 

""" 

Adds two DOS together. Checks that energy scales are the same. 

Otherwise, a ValueError is thrown. 

 

Args: 

other: Another DOS object. 

 

Returns: 

Sum of the two DOSs. 

""" 

if not all(np.equal(self.energies, other.energies)): 

raise ValueError("Energies of both DOS are not compatible!") 

densities = {spin: self.densities[spin] + other.densities[spin] 

for spin in self.densities.keys()} 

return Dos(self.efermi, self.energies, densities) 

 

def get_interpolated_value(self, energy): 

""" 

Returns interpolated density for a particular energy. 

 

Args: 

energy: Energy to return the density for. 

""" 

f = {} 

for spin in self.densities.keys(): 

f[spin] = get_linear_interpolated_value(self.energies, 

self.densities[spin], 

energy) 

return f 

 

def get_interpolated_gap(self, tol=0.001, abs_tol=False, spin=None): 

""" 

Expects a DOS object and finds the gap 

 

Args: 

tol: tolerance in occupations for determining the gap 

abs_tol: Set to True for an absolute tolerance and False for a 

relative one. 

spin: Possible values are None - finds the gap in the summed 

densities, Up - finds the gap in the up spin channel, 

Down - finds the gap in the down spin channel. 

 

Returns: 

(gap, cbm, vbm): 

Tuple of floats in eV corresponding to the gap, cbm and vbm. 

""" 

 

tdos = self.get_densities(spin) 

if not abs_tol: 

tol = tol * tdos.sum() / tdos.shape[0] 

energies = self.energies 

below_fermi = [i for i in range(len(energies)) 

if energies[i] < self.efermi and tdos[i] > tol] 

above_fermi = [i for i in range(len(energies)) 

if energies[i] > self.efermi and tdos[i] > tol] 

vbm_start = max(below_fermi) 

cbm_start = min(above_fermi) 

if vbm_start == cbm_start: 

return 0.0, self.efermi, self.efermi 

else: 

# Interpolate between adjacent values 

terminal_dens = tdos[vbm_start:vbm_start + 2][::-1] 

terminal_energies = energies[vbm_start:vbm_start + 2][::-1] 

start = get_linear_interpolated_value(terminal_dens, 

terminal_energies, tol) 

terminal_dens = tdos[cbm_start - 1:cbm_start + 1] 

terminal_energies = energies[cbm_start - 1:cbm_start + 1] 

end = get_linear_interpolated_value(terminal_dens, 

terminal_energies, tol) 

return end - start, end, start 

 

def get_cbm_vbm(self, tol=0.001, abs_tol=False, spin=None): 

""" 

Expects a DOS object and finds the cbm and vbm. 

 

Args: 

tol: tolerance in occupations for determining the gap 

abs_tol: An absolute tolerance (True) and a relative one (False) 

spin: Possible values are None - finds the gap in the summed 

densities, Up - finds the gap in the up spin channel, 

Down - finds the gap in the down spin channel. 

 

Returns: 

(cbm, vbm): float in eV corresponding to the gap 

""" 

#determine tolerance 

tdos = self.get_densities(spin) 

if not abs_tol: 

tol = tol * tdos.sum() / tdos.shape[0] 

 

# find index of fermi energy 

i_fermi = 0 

while self.energies[i_fermi] <= self.efermi: 

i_fermi += 1 

 

# work backwards until tolerance is reached 

i_gap_start = i_fermi 

while i_gap_start - 1 >= 0 and tdos[i_gap_start - 1] <= tol: 

i_gap_start -= 1 

 

# work forwards until tolerance is reached 

i_gap_end = i_gap_start 

while i_gap_end < tdos.shape[0] and tdos[i_gap_end] <= tol: 

i_gap_end += 1 

i_gap_end -= 1 

return self.energies[i_gap_end], self.energies[i_gap_start] 

 

def get_gap(self, tol=0.001, abs_tol=False, spin=None): 

""" 

Expects a DOS object and finds the gap. 

 

Args: 

tol: tolerance in occupations for determining the gap 

abs_tol: An absolute tolerance (True) and a relative one (False) 

spin: Possible values are None - finds the gap in the summed 

densities, Up - finds the gap in the up spin channel, 

Down - finds the gap in the down spin channel. 

 

Returns: 

gap in eV 

""" 

(cbm, vbm) = self.get_cbm_vbm(tol, abs_tol, spin) 

return max(cbm - vbm, 0.0) 

 

def __str__(self): 

""" 

Returns a string which can be easily plotted (using gnuplot). 

""" 

if Spin.down in self.densities: 

stringarray = ["#{:30s} {:30s} {:30s}".format("Energy", 

"DensityUp", 

"DensityDown")] 

for i, energy in enumerate(self.energies): 

stringarray.append("{:.5f} {:.5f} {:.5f}" 

.format(energy, self.densities[Spin.up][i], 

self.densities[Spin.down][i])) 

else: 

stringarray = ["#{:30s} {:30s}".format("Energy", "DensityUp")] 

for i, energy in enumerate(self.energies): 

stringarray.append("{:.5f} {:.5f}" 

.format(energy, self.densities[Spin.up][i])) 

return "\n".join(stringarray) 

 

@classmethod 

def from_dict(cls, d): 

""" 

Returns Dos object from dict representation of Dos. 

""" 

return Dos(d["efermi"], d["energies"], 

{Spin(int(k)): v 

for k, v in d["densities"].items()}) 

 

def as_dict(self): 

""" 

Json-serializable dict representation of Dos. 

""" 

return {"@module": self.__class__.__module__, 

"@class": self.__class__.__name__, "efermi": self.efermi, 

"energies": list(self.energies), 

"densities": {str(spin): list(dens) 

for spin, dens in self.densities.items()}} 

 

 

class CompleteDos(Dos): 

""" 

This wrapper class defines a total dos, and also provides a list of PDos. 

Mainly used by pymatgen.io.vasp.Vasprun to create a complete Dos from 

a vasprun.xml file. You are unlikely to try to generate this object 

manually. 

 

Args: 

structure: Structure associated with this particular DOS. 

total_dos: total Dos for structure 

pdoss: The pdoss are supplied as an {Site:{Orbital:{ 

Spin:Densities}}} 

 

.. attribute:: structure 

 

Structure associated with the CompleteDos. 

 

.. attribute:: pdos 

 

Dict of partial densities of the form {Site:{Orbital:{Spin:Densities}}} 

""" 

def __init__(self, structure, total_dos, pdoss): 

super(CompleteDos, self).__init__( 

total_dos.efermi, energies=total_dos.energies, 

densities={k: np.array(d) for k, d in total_dos.densities.items()}) 

self.pdos = pdoss 

self.structure = structure 

 

def get_site_orbital_dos(self, site, orbital): 

""" 

Get the Dos for a particular orbital of a particular site. 

 

Args: 

site: Site in Structure associated with CompleteDos. 

orbital: Orbital in the site. 

 

Returns: 

Dos containing densities for orbital of site. 

""" 

return Dos(self.efermi, self.energies, self.pdos[site][orbital]) 

 

def get_site_dos(self, site): 

""" 

Get the total Dos for a site (all orbitals). 

 

Args: 

site: Site in Structure associated with CompleteDos. 

 

Returns: 

Dos containing summed orbital densities for site. 

""" 

site_dos = six.moves.reduce(add_densities, self.pdos[site].values()) 

return Dos(self.efermi, self.energies, site_dos) 

 

def get_site_spd_dos(self, site): 

""" 

Get orbital projected Dos of a particular site 

 

Args: 

site: Site in Structure associated with CompleteDos. 

 

Returns: 

dict of {orbital: Dos}, e.g. {"s": Dos object, ...} 

""" 

spd_dos = dict() 

for orb, pdos in self.pdos[site].items(): 

orbital_type = _get_orb_type(orb) 

if orbital_type in spd_dos: 

spd_dos[orbital_type] = add_densities(spd_dos[orbital_type], pdos) 

else: 

spd_dos[orbital_type] = pdos 

return {orb: Dos(self.efermi, self.energies, densities) 

for orb, densities in spd_dos.items()} 

 

def get_site_t2g_eg_resolved_dos(self, site): 

""" 

Get the t2g, eg projected DOS for a particular site. 

 

Args: 

site: Site in Structure associated with CompleteDos. 

 

Returns: 

A dict {"e_g": Dos, "t2g": Dos} containing summed e_g and t2g DOS 

for the site. 

""" 

t2g_dos = [] 

eg_dos = [] 

for s, atom_dos in self.pdos.items(): 

if s == site: 

for orb, pdos in atom_dos.items(): 

if orb in (Orbital.dxy, Orbital.dxz, Orbital.dyz): 

t2g_dos.append(pdos) 

elif orb in (Orbital.dx2, Orbital.dz2): 

eg_dos.append(pdos) 

return {"t2g": Dos(self.efermi, self.energies, 

six.moves.reduce(add_densities, t2g_dos)), 

"e_g": Dos(self.efermi, self.energies, 

six.moves.reduce(add_densities, eg_dos))} 

 

def get_spd_dos(self): 

""" 

Get orbital projected Dos. 

 

Returns: 

dict of {orbital: Dos}, e.g. {"s": Dos object, ...} 

""" 

spd_dos = {} 

for atom_dos in self.pdos.values(): 

for orb, pdos in atom_dos.items(): 

orbital_type = _get_orb_type(orb) 

if orbital_type not in spd_dos: 

spd_dos[orbital_type] = pdos 

else: 

spd_dos[orbital_type] = \ 

add_densities(spd_dos[orbital_type], pdos) 

return {orb: Dos(self.efermi, self.energies, densities) 

for orb, densities in spd_dos.items()} 

 

def get_element_dos(self): 

""" 

Get element projected Dos. 

 

Returns: 

dict of {Element: Dos} 

""" 

 

el_dos = {} 

for site, atom_dos in self.pdos.items(): 

el = site.specie 

for pdos in atom_dos.values(): 

if el not in el_dos: 

el_dos[el] = pdos 

else: 

el_dos[el] = add_densities(el_dos[el], pdos) 

return {el: Dos(self.efermi, self.energies, densities) 

for el, densities in el_dos.items()} 

 

def get_element_spd_dos(self, el): 

""" 

Get element and spd projected Dos 

 

Args: 

el: Element in Structure.composition associated with CompleteDos 

 

Returns: 

dict of {Element: {"S": densities, "P": densities, "D": densities}} 

""" 

el = get_el_sp(el) 

el_dos = {} 

for site, atom_dos in self.pdos.items(): 

if site.specie == el: 

for orb, pdos in atom_dos.items(): 

orbital_type = _get_orb_type(orb) 

if orbital_type not in el_dos: 

el_dos[orbital_type] = pdos 

else: 

el_dos[orbital_type] = \ 

add_densities(el_dos[orbital_type], pdos) 

 

return {orb: Dos(self.efermi, self.energies, densities) 

for orb, densities in el_dos.items()} 

 

@classmethod 

def from_dict(cls, d): 

""" 

Returns CompleteDos object from dict representation. 

""" 

tdos = Dos.from_dict(d) 

struct = Structure.from_dict(d["structure"]) 

pdoss = {} 

for i in range(len(d["pdos"])): 

at = struct[i] 

orb_dos = {} 

for orb_str, odos in d["pdos"][i].items(): 

orb = Orbital[orb_str] 

orb_dos[orb] = {Spin(int(k)): v 

for k, v in odos["densities"].items()} 

pdoss[at] = orb_dos 

return CompleteDos(struct, tdos, pdoss) 

 

def as_dict(self): 

""" 

Json-serializable dict representation of CompleteDos. 

""" 

d = {"@module": self.__class__.__module__, 

"@class": self.__class__.__name__, "efermi": self.efermi, 

"structure": self.structure.as_dict(), 

"energies": list(self.energies), 

"densities": {str(spin): list(dens) 

for spin, dens in self.densities.items()}, 

"pdos": []} 

if len(self.pdos) > 0: 

for at in self.structure: 

dd = {} 

for orb, pdos in self.pdos[at].items(): 

dd[str(orb)] = {"densities": {str(int(spin)): list(dens) 

for spin, 

dens in pdos.items()}} 

d["pdos"].append(dd) 

d["atom_dos"] = {str(at): dos.as_dict() for at, 

dos in self.get_element_dos().items()} 

d["spd_dos"] = {str(orb): dos.as_dict() for orb, 

dos in self.get_spd_dos().items()} 

return d 

 

def __str__(self): 

return "Complete DOS for " + str(self.structure) 

 

 

def add_densities(density1, density2): 

""" 

Method to sum two densities. 

 

Args: 

density1: First density. 

density2: Second density. 

 

Returns: 

Dict of {spin: density}. 

""" 

return {spin: np.array(density1[spin]) + np.array(density2[spin]) 

for spin in density1.keys()} 

 

 

def _get_orb_type(orb): 

try: 

return orb.orbital_type 

except AttributeError: 

return orb