Coverage for gpaw/lrtddft/__init_

1"""This module defines a linear response TDDFT-class."""

2import numbers

3import sys

4from math import sqrt

5from typing import Dict, Any

6import numpy as np

8from ase.units import Hartree

9from ase.utils.timing import Timer

11import gpaw.mpi as mpi

12from gpaw.xc import XC

13from gpaw.lrtddft.excitation import Excitation, ExcitationList, get_filehandle

14from gpaw.lrtddft.kssingle import KSSingles

15from gpaw.lrtddft.omega_matrix import OmegaMatrix

16from gpaw.lrtddft.apmb import ApmB

17from gpaw.lrtddft.spectrum import spectrum

19__all__ = ['LrTDDFT', 'photoabsorption_spectrum', 'spectrum']

22class LrTDDFT(ExcitationList):

23 """Linear Response TDDFT excitation list class

25 Input parameters:

27 calculator:

28 the calculator object after a ground state calculation

30 nspins:

31 number of spins considered in the calculation

32 Note: Valid only for unpolarised ground state calculation

34 eps:

35 Minimal occupation difference for a transition (default 0.001)

37 istart:

38 First occupied state to consider

39 jend:

40 Last unoccupied state to consider

42 xc:

43 Exchange-Correlation approximation in the Kernel

44 derivative_level:

45 0: use Exc, 1: use vxc, 2: use fxc if available

47 filename:

48 read from a file

49 """

51 default_parameters: Dict[str, Any] = {

52 'nspins': None,

53 'restrict': {},

54 'xc': None,

55 'derivative_level': 1,

56 'numscale': 0.00001,

57 'filename': None,

58 'finegrid': 2,

59 'force_ApmB': False, # for tests

60 'eh_comm': None, # parallelization over eh-pairs

61 'poisson': None} # use calculator's Poisson

63 def __init__(self, calculator=None, log=None, txt='-', **kwargs):

65 self.energy_to_eV_scale = Hartree

66 self.timer = Timer()

67 self.diagonalized = False

69 self.set(**kwargs)

70 self.calculator = calculator

72 ExcitationList.__init__(self, log=log, txt=txt)

74 if self.eh_comm is None:

75 self.eh_comm = mpi.serial_comm

76 elif isinstance(self.eh_comm, (mpi.world.__class__,

77 mpi.serial_comm.__class__)):

78 # Correct type already.

79 pass

80 else:

81 # world should be a list of ranks:

82 self.eh_comm = mpi.world.new_communicator(

83 np.asarray(self.eh_comm))

85 if calculator is not None and calculator.initialized:

86 # XXXX not ready for k-points

87 assert len(calculator.wfs.kd.bzk_kc) == 1

88 if calculator.wfs.mode not in ['fd', 'lcao']:

89 raise RuntimeError('LrTDDFT supports only fd and lcao modes')

90 if calculator.wfs.mode != 'lcao':

91 calculator.converge_wave_functions()

92 if calculator.density.nct_G is None:

93 spos_ac = calculator.initialize_positions()

94 calculator.wfs.initialize(calculator.density,

95 calculator.hamiltonian, spos_ac)

97 self.forced_update()

99 @property

100 def calculator(self):

101 return self._calc

102

103 @calculator.setter

104 def calculator(self, calc):

105 self._calc = calc

106

107 if self.xc is None and calc is not None:

108 if calc.initialized:

109 self.xc = calc.hamiltonian.xc

110 else:

111 self.xc = calc.parameters['xc']

112

113 def set(self, **kwargs):

114 """Change parameters."""

115 changed = []

116 for key, value in LrTDDFT.default_parameters.items():

117 if hasattr(self, key):

118 value = getattr(self, key) # do not overwrite

119 setattr(self, key, kwargs.pop(key, value))

120 if value != getattr(self, key):

121 changed.append(key)

122

123 for key in kwargs:

124 raise KeyError('Unknown key ' + key)

125

126 return changed

127

128 def analyse(self, what=None, out=None, min=0.1):

129 """Print info about the transitions.

130

131 Parameters:

132 1. what: I list of excitation indicees, None means all

133 2. out : I where to send the output, None means sys.stdout

134 3. min : I minimal contribution to list (0<min<1)

135 """

136 if what is None:

137 what = range(len(self))

138 elif isinstance(what, numbers.Integral):

139 what = [what]

140

141 if out is None:

142 out = sys.stdout

143

144 for i in what:

145 print(str(i) + ':', self[i].analyse(min=min), file=out)

146

147 def calculate(self, atoms):

148 self.calculator = atoms.calc

149 if not hasattr(self, 'Om') or self.calculator.check_state(atoms):

150 self.calculator.get_potential_energy(atoms)

151 self.forced_update()

152 return self

153

154 def forced_update(self):

155 """Recalc yourself."""

156 Om = OmegaMatrix

157 name = 'LrTDDFT'

158

159 hybrid = False

160 if self.xc and self.xc != 'RPA':

161 xc = self.xc

162 if isinstance(self.xc, str):

163 xc = XC(self.xc)

164 hybrid = hasattr(xc, 'hybrid') and xc.hybrid > 0.0

165

166 if hybrid or self.force_ApmB:

167 Om = ApmB

168 name = 'LrTDDFThyb'

169

170 kss = KSSingles(restrict=self.restrict,

171 log=self.log)

172 atoms = self.calculator.get_atoms()

173 kss.calculate(atoms, self.nspins)

174

175 self.Om = Om(self.calculator, kss,

176 self.xc, self.derivative_level, self.numscale,

177 finegrid=self.finegrid, eh_comm=self.eh_comm,

178 poisson=self.poisson, log=self.log)

179 self.name = name

180

181 def diagonalize(self, **kwargs):

182 """Diagonalize and save new Eigenvalues and Eigenvectors"""

183 self.set(**kwargs)

184 self.timer.start('diagonalize')

185 self.timer.start('omega')

186 self.Om.diagonalize(kwargs.pop('restrict', {}))

187 self.timer.stop('omega')

188 self.diagonalized = True

189

190 # remove old stuff

191 self.timer.start('clean')

192 while len(self):

193 self.pop()

194 self.timer.stop('clean')

195

196 self.log('LrTDDFT digonalized:')

197 self.timer.start('build')

198 for j in range(len(self.Om.kss)):

199 self.append(LrTDDFTExcitation(self.Om, j))

200 self.log(' ', str(self[-1]))

201 self.timer.stop('build')

202 self.timer.stop('diagonalize')

203

204 @classmethod

205 def read(cls, filename=None, fh=None, restrict={}, log=None, txt=None):

206 """Read myself from a file"""

207 lr = cls(log=log, txt=txt)

208 timer = lr.timer

209 timer.start('name')

210 if fh is None:

211 f = get_filehandle(lr, filename)

212 else:

213 f = fh

214 timer.stop('name')

215

216 timer.start('header')

217 # get my name

218 s = f.readline().strip()

219 lr.name = s.split()[1]

220

221 xc = f.readline().strip().split()[0]

222 if xc == 'RPA':

223 lr.xc = xc

224 else:

225 lr.xc = XC(xc)

226 values = f.readline().split()

227 eps = float(values[0])

228 if len(values) > 1:

229 lr.derivative_level = int(values[1])

230 lr.numscale = float(values[2])

231 lr.finegrid = int(values[3])

232 else:

233 # old writing style, use old defaults

234 lr.numscale = 0.001

235 timer.stop('header')

236

237 timer.start('init_kss')

238 kss = KSSingles.read(fh=f, log=log)

239 assert eps == kss.restrict['eps']

240 lr.restrict = kss.restrict.values

241 timer.stop('init_kss')

242 timer.start('init_obj')

243 if lr.name == 'LrTDDFT':

244 lr.Om = OmegaMatrix(kss=kss, filehandle=f, log=lr.log)

245 else:

246 lr.Om = ApmB(kss=kss, filehandle=f, log=lr.log)

247 timer.stop('init_obj')

248

249 if not len(restrict):

250 timer.start('read diagonalized')

251 # check if already diagonalized

252 p = f.tell()

253 s = f.readline()

254 if s != '# Eigenvalues\n' or len(restrict):

255 # no further info or selection of

256 # Kohn-Sham states changed

257 # go back to previous position

258 f.seek(p)

259 else:

260 lr.diagonalized = True

261 # load the eigenvalues

262 n = int(f.readline().split()[0])

263 for i in range(n):

264 lr.append(LrTDDFTExcitation(string=f.readline()))

265 # load the eigenvectors

266 timer.start('read eigenvectors')

267 f.readline()

268 for i in range(n):

269 lr[i].f = np.array([float(x) for x in

270 f.readline().split()])

271 lr[i].kss = lr.kss

272 timer.stop('read eigenvectors')

273 timer.stop('read diagonalized')

274 else:

275 timer.start('diagonalize')

276 lr.diagonalize(restrict=restrict)

277 timer.stop('diagonalize')

278

279 if fh is None:

280 f.close()

281

282 return lr

283

284 @property

285 def kss(self):

286 return self.Om.kss

287

288 def singlets_triplets(self):

289 """Split yourself into a singlet and triplet object"""

290

291 slr = LrTDDFT(nspins=self.nspins, xc=self.xc,

292 restrict=self.kss.restrict.values,

293 derivative_level=self.derivative_level,

294 numscale=self.numscale)

295 tlr = LrTDDFT(nspins=self.nspins, xc=self.xc,

296 restrict=self.kss.restrict.values,

297 derivative_level=self.derivative_level,

298 numscale=self.numscale)

299 slr.Om, tlr.Om = self.Om.singlets_triplets()

300

301 return slr, tlr

302

303 def single_pole_approximation(self, i, j):

304 """Return the excitation according to the

305 single pole approximation. See e.g.:

306 Grabo et al, Theochem 501 (2000) 353-367

307 """

308 for ij, kss in enumerate(self.kss):

309 if kss.i == i and kss.j == j:

310 return sqrt(self.Om.full[ij][ij]) * Hartree

311 return self.Om.full[ij][ij] / kss.energy * Hartree

312

313 def __str__(self):

314 string = ExcitationList.__str__(self)

315 string += '# derived from:\n'

316 string += self.Om.kss.__str__()

317 return string

318

319 def write(self, filename=None, fh=None):

320 """Write current state to a file.

321

322 'filename' is the filename. If the filename ends in .gz,

323 the file is automatically saved in compressed gzip format.

324

325 'fh' is a filehandle. This can be used to write into already

326 opened files.

327 """

328

329 if self.calculator is None:

330 rank = mpi.world.rank

331 else:

332 rank = self.calculator.wfs.world.rank

333

334 if rank == 0:

335 if fh is None:

336 f = get_filehandle(self, filename, mode='w')

337 else:

338 f = fh

339

340 f.write('# ' + self.name + '\n')

341 if isinstance(self.xc, str):

342 xc = self.xc

343 else:

344 xc = self.xc.tostring()

345 if xc is None:

346 xc = 'RPA'

347 if self.calculator is not None:

348 xc += ' ' + self.calculator.get_xc_functional()

349 f.write(xc + '\n')

350 f.write('%g %d %g %d' % (self.kss.restrict['eps'],

351 int(self.derivative_level),

352 self.numscale, int(self.finegrid)) + '\n')

353 self.kss.write(fh=f)

354 self.Om.write(fh=f)

355

356 if len(self):

357 f.write('# Eigenvalues\n')

358 f.write(f'{len(self)}\n')

359 for ex in self:

360 f.write(ex.outstring())

361 f.write('# Eigenvectors\n')

362 for ex in self:

363 for w in ex.f:

364 f.write('%g ' % w)

365 f.write('\n')

366

367 if fh is None:

368 f.close()

369 mpi.world.barrier()

370

371 def overlap(self, ov_nn, other):

372 """Matrix element overlap determined from pair density overlaps.

373

374 Parameters

375 ----------

376 ov_nn: array

377 Wave function overlap factors from a displaced calculator.

378

379 Index 0 corresponds to our own wavefunctions and

380 index 1 to the others wavefunctions

381

382 Returns

383 -------

384 ov_pp: array

385 Overlap

386 """

387 # XXX ov_pp = self.kss.overlap(ov_nn, other.kss)

388 ov_pp = self.Om.kss.overlap(ov_nn, other.Om.kss)

389 self.diagonalize()

390 other.diagonalize()

391 # ov[pLm, pLo] = Om[pLm, :pKm]* ov[:pKm, pLo]

392 return np.dot(self.Om.eigenvectors.conj(),

393 # ov[pKm, pLo] = ov[pKm, :pKo] Om[pLo, :pKo].T

394 np.dot(ov_pp, other.Om.eigenvectors.T))

395

396 def __getitem__(self, i):

397 if not self.diagonalized:

398 self.diagonalize()

399 return list.__getitem__(self, i)

400

401 def __iter__(self):

402 if not self.diagonalized:

403 self.diagonalize()

404 return list.__iter__(self)

405

406 def __len__(self):

407 if not self.diagonalized:

408 self.diagonalize()

409 return list.__len__(self)

410

411 def __del__(self):

412 try:

413 self.timer.write(self.log.fd)

414 except AttributeError:

415 pass

416

417

418class LrTDDFTExcitation(Excitation):

419

420 def __init__(self, Om=None, i=None,

421 e=None, m=None, string=None):

422

423 # multiplicity comes from Kohn-Sham contributions

424 self.fij = 1

425

426 if string is not None:

427 self.fromstring(string)

428 return None

429

430 # define from the diagonalized Omega matrix

431 if Om is not None:

432 if i is None:

433 raise RuntimeError

434

435 ev = Om.eigenvalues[i]

436 if ev < 0:

437 # we reached an instability, mark it with a negative value

438 self.energy = -sqrt(-ev)

439 else:

440 self.energy = sqrt(ev)

441 self.f = Om.eigenvectors[i]

442 self.kss = Om.kss

443

444 self.kss.set_arrays()

445 self.me = np.dot(self.f, self.kss.me)

446 erat_k = np.sqrt(self.kss.energies / self.energy)

447 wght_k = np.sqrt(self.kss.fij) * self.f

448 ew_k = erat_k * wght_k

449 self.mur = np.dot(ew_k, self.kss.mur)

450 if self.kss.muv is not None:

451 self.muv = np.dot(ew_k, self.kss.muv)

452 else:

453 self.muv = None

454 if self.kss.magn is not None:

455 self.magn = np.dot(wght_k / erat_k, self.kss.magn)

456 else:

457 self.magn = None

458

459 return

460

461 # define from energy and matrix element

462 if e is not None:

463 self.energy = e

464 if m is None:

465 raise RuntimeError

466 self.me = m

467 return

468

469 raise RuntimeError

470

471 def density_change(self, paw):

472 """get the density change associated with this transition"""

473 raise NotImplementedError

474

475 def fromstring(self, string):

476 l = string.split()

477 self.energy = float(l.pop(0))

478 if len(l) == 3: # old writing style

479 self.me = np.array([float(l.pop(0)) for i in range(3)])

480 else:

481 self.mur = np.array([float(l.pop(0)) for i in range(3)])

482 self.me = - self.mur * sqrt(self.energy)

483 self.muv = np.array([float(l.pop(0)) for i in range(3)])

484 self.magn = np.array([float(l.pop(0)) for i in range(3)])

485

486 def outstring(self):

487 str = '%g ' % self.energy

488 str += ' '

489 for m in self.mur:

490 str += '%12.4e' % m

491 str += ' '

492 for m in self.muv:

493 str += '%12.4e' % m

494 str += ' '

495 for m in self.magn:

496 str += '%12.4e' % m

497 str += '\n'

498 return str

499

500 def __str__(self):

501 m2 = np.sum(self.me * self.me)

502 m = sqrt(m2)

503 if m > 0:

504 me = self.me / m

505 else:

506 me = self.me

507 str = '<LrTDDFTExcitation> om=%g[eV] |me|=%g (%.2f,%.2f,%.2f)' % \

508 (self.energy * Hartree, m, me[0], me[1], me[2])

509 return str

510

511 def analyse(self, min=.1):

512 """Return an analysis string of the excitation"""

513 osc = self.get_oscillator_strength()

514 s = ('E=%.3f' % (self.energy * Hartree) + ' eV, ' +

515 'f=%.5g' % osc[0] + ', (%.5g,%.5g,%.5g) ' %

516 (osc[1], osc[2], osc[3]) + '\n')

517 # 'R=%.5g' % self.get_rotatory_strength() + ' cgs\n')

518

519 def sqr(x):

520 return x * x

521 spin = ['u', 'd']

522 min2 = sqr(min)

523 rest = np.sum(self.f**2)

524 for f, k in zip(self.f, self.kss):

525 f2 = sqr(f)

526 if f2 > min2:

527 s += ' %d->%d ' % (k.i, k.j) + spin[k.pspin] + ' '

528 s += '%.3g \n' % f2

529 rest -= f2

530 s += ' rest=%.3g' % rest

531 return s

532

533

534def photoabsorption_spectrum(excitation_list, spectrum_file=None,

535 e_min=None, e_max=None, delta_e=None,

536 energyunit='eV',

537 folding='Gauss', width=0.1, comment=None):

538 """Uniform absorption spectrum interface

539

540 Parameters:

541 ================= ===================================================

542 ``exlist`` ExcitationList

543 ``spectrum_file`` File name for the output file, STDOUT if not given

544 ``e_min`` min. energy, set to cover all energies if not given

545 ``e_max`` max. energy, set to cover all energies if not given

546 ``delta_e`` energy spacing

547 ``energyunit`` Energy unit, default 'eV'

548 ``folding`` Gauss (default) or Lorentz

549 ``width`` folding width in terms of the chosen energyunit

550 ================= ===================================================

551 all energies in [eV]

552 """

553

554 spectrum(exlist=excitation_list, filename=spectrum_file,

555 emin=e_min, emax=e_max,

556 de=delta_e, energyunit=energyunit,

557 folding=folding, width=width,

558 comment=comment)

Coverage for gpaw/lrtddft/init.py: 90%

324 statements