Coverage for gpaw/lrtddft2/lr_transitions.py: 73%

1import numpy as np

2import ase.units

3from scipy.linalg import eigh

5from gpaw.mpi import world

6from gpaw.lrtddft2.lr_layouts import LrDiagonalizeLayout

9class LrtddftTransitions:

10 def __init__(self, ks_singles, K_matrix, sl_lrtddft=None):

11 self.ks_singles = ks_singles

12 self.K_matrix = K_matrix

13 self.lr_comms = self.ks_singles.lr_comms

14 self.sl_lrtddft = sl_lrtddft

15 self.eigenvalues = None

16 self.eigenvectors = None

17 self.trans_prop_ready = False

19 # only for pros

20 self.custom_axes = None

22 def initialize(self):

23 self.trans_prop_ready = False

25 def read(self):

26 pass

28 def calculate(self):

29 self.diagonalize()

30 self.calculate_properties()

32 def diagonalize(self):

33 if self.sl_lrtddft is None:

34 self.diagonalize_serial()

35 else:

36 self.diagonalize_scalapack()

38 # now self.eigenvectors[iploc,kloc]

39 # is the "iploc"th local element of

40 # the "kloc"th local eigenvector

42 def diagonalize_serial(self):

43 nrows = len(self.ks_singles.kss_list)

45 # build local part of the full omega matrix

46 omega_matrix = np.zeros((nrows, nrows), dtype=float)

47 for (ip, kss_ip) in enumerate(self.ks_singles.kss_list):

48 for (jq, kss_jq) in enumerate(self.ks_singles.kss_list):

49 lip = self.lr_comms.get_local_eh_index(ip)

50 ljq = self.lr_comms.get_local_dd_index(jq)

51 if lip is None or ljq is None:

52 continue

54 # K-matrix

55 # if self.K_matrix.file_format == 1:

56 omega_matrix[ip, jq] = 2. * np.sqrt(

57 kss_ip.energy_diff * kss_jq.energy_diff * kss_ip.pop_diff *

58 kss_jq.pop_diff) * self.K_matrix.values[lip, ljq]

59 # old Casida format ... doing this already in k_matrix.py

60 # elif self.K_matrix.file_format == 0:

61 # omega_matrix[ip,jq] = self.K_matrix.values[lip,ljq]

62 # invalid

63 # else:

64 # raise RuntimeError('Invalid K-matrix file format')

66 if (ip == jq):

67 omega_matrix[ip,

68 jq] += kss_ip.energy_diff * kss_jq.energy_diff

70 # full omega matrix to each process

71 self.lr_comms.parent_comm.sum(omega_matrix)

73 # if self.lr_comms.parent_comm.rank == 0:

74 # print omega_matrix

76 # solve eigensystem

78 # debug

79 # omega_matrix[:] = 0

80 # for i in range(omega_matrix.shape[0]):

81 # for j in range(omega_matrix.shape[1]):

82 # if i == j: omega_matrix[i,j] = -2.

83 # if i+1 == j: omega_matrix[i,j] = omega_matrix[j,i] = 1.

85 self.eigenvalues, omega_matrix = eigh(omega_matrix)

87 # convert eigenvector back to local version

88 nlrows = self.K_matrix.values.shape[0]

89 nlcols = self.K_matrix.values.shape[1]

91 # debug

92 # (eh=2,dd=3) (0,0) (0,1) (0,2) (1,0) (1,1) (1,2)

93 # 00 01 02 03 00 03 01 02 10 13 11 12

94 # 10 11 12 13 => 20 23 21 22 30 33 31 32

95 # 20 21 22 23

96 # 30 31 32 33

98 self.eigenvectors = np.zeros((nlrows, nlcols), dtype=float)

99 for (ip, kss) in enumerate(self.ks_singles.kss_list):

100 for (jq, kss) in enumerate(self.ks_singles.kss_list):

101 lip = self.lr_comms.get_local_eh_index(ip)

102 ljq = self.lr_comms.get_local_dd_index(jq)

103 if lip is None or ljq is None:

104 continue

105 self.eigenvectors[lip,

106 ljq] = omega_matrix[ip,

107 jq] # debug: ip*10 + jq

108

109 def diagonalize_scalapack(self):

110 # total rows

111 nrows = len(self.ks_singles.kss_list)

112 # local

113 nlrows = self.K_matrix.values.shape[0]

114 nlcols = self.K_matrix.values.shape[1]

115

116 # create BLACS layout

117 layout = LrDiagonalizeLayout(self.sl_lrtddft, nrows, self.lr_comms)

118 if (nrows <

119 np.max([layout.mprocs, layout.nprocs]) * layout.block_size):

120 raise RuntimeError(

121 'Linear response matrix is not large enough for the given '

122 'number of processes (or block size) in sl_lrtddft. Please, '

123 'use less processes (or smaller block size).'

124 )

125

126 # build local part

127 omega_matrix = np.zeros((nlrows, nlcols), dtype=float)

128 for (ip, kss_ip) in enumerate(self.ks_singles.kss_list):

129 for (jq, kss_jq) in enumerate(self.ks_singles.kss_list):

130 lip = self.lr_comms.get_local_eh_index(ip)

131 ljq = self.lr_comms.get_local_dd_index(jq)

132 if lip is None or ljq is None:

133 continue

134

135 # if self.K_matrix.file_format == 1:

136 omega_matrix[lip, ljq] = 2. * np.sqrt(

137 kss_ip.energy_diff * kss_jq.energy_diff * kss_ip.pop_diff *

138 kss_jq.pop_diff) * self.K_matrix.values[lip, ljq]

139 # old Casida format ... doing this already in k_matrix.py

140 # elif self.K_matrix.file_format == 0:

141 # omega_matrix[lip,ljq] = self.K_matrix.values[lip,ljq]

142 # else:

143 # raise RuntimeError('Invalid K-matrix file format')

144

145 if (ip == jq):

146 omega_matrix[

147 lip, ljq] += kss_ip.energy_diff * kss_jq.energy_diff

148

149 # solve eigen system (transpose for scalapack)

150 self.eigenvalues = np.zeros(nrows, dtype=float)

151 self.eigenvectors = np.zeros(

152 (omega_matrix.shape[1], omega_matrix.shape[0]), dtype=float)

153 self.eigenvectors[:, :] = np.transpose(omega_matrix)

154 layout.diagonalize(self.eigenvectors, self.eigenvalues)

155 omega_matrix[:, :] = np.transpose(self.eigenvectors)

156 self.eigenvectors = omega_matrix

157

158 def calculate_properties(self):

159 if self.custom_axes is not None:

160 self.custom_axes = np.array(self.custom_axes)

161 # else:

162 # self.custom_axes = np.array([[1.0,0.0,0.0],

163 # [0.0,1.0,0.0],

164 # [0.0,0.0,1.0]])

165

166 self.lr_comms.parent_comm.barrier()

167

168 # nlrows = self.eigenvectors.shape[0]

169 nleigs = self.eigenvectors.shape[1]

170

171 sqrtwloc = np.zeros(nleigs)

172 dmxloc = np.zeros(nleigs)

173 dmyloc = np.zeros(nleigs)

174 dmzloc = np.zeros(nleigs)

175 magnxloc = np.zeros(nleigs)

176 magnyloc = np.zeros(nleigs)

177 magnzloc = np.zeros(nleigs)

178 cloc_dm = np.zeros(nleigs)

179 cloc_magn = np.zeros(nleigs)

180

181 # see Autschbach et al., J. Chem. Phys., 116, 6930 (2002)

182 # see also Varsano et al., Phys.Chem.Chem.Phys. 11, 4481 (2009)

183 #

184 # mu_k = sum_jq sqrt(ediff_jq / omega_k)

185 # sqrt(population_jq) * F^(k)_jq * mu_jq

186 # m_k = sum_jq sqrt(omega_k / ediff_jq)

187 # sqrt(population_jq) * F^(k)_jq * m_jq

188 #

189 # FIXME: check once more

190 # mu_k = sum_jq c1_k * c2_k * mu_jq

191 # = sum_jq sqrt(ediff_jq / omega_k)

192 # sqrt(fdiff_jq) * F^(k)_jq mu_jq

193 # m_k = sum_jq c2_k / c1_k * m_jq

194 # = sum_jq sqrt(fdiff_jq) sqrt(omega_k / ediff_jq) * F^(k)_jq m_jq

195 #

196 # c1 = sqrt(ediff_jq / omega_k)

197 # = np.sqrt(kss_jq.energy_diff / self.get_excitation_energy(k))

198 #

199 # c2 = sqrt(fdiff_jq) * F^(k)_jq

200 # = np.sqrt(kss_jq.pop_diff) *

201 # self.get_local_eig_coeff(k,self.index_map[(i,p)])

202 #

203

204 # sqrt(omega_kloc)

205 for kloc in range(nleigs):

206 sqrtwloc[kloc] = np.sqrt(

207 np.sqrt(

208 self.eigenvalues[self.lr_comms.get_global_dd_index(kloc)]))

209

210 # loop over ks singles

211 for (ip, kss_ip) in enumerate(self.ks_singles.kss_list):

212

213 # local ip index

214 lip = self.lr_comms.get_local_eh_index(ip)

215 # print self.lr_comms.parent_comm.rank, i,p,ip,lip

216 if lip is None:

217 continue

218

219 # now self.eigenvectors[iploc,kloc]

220 # is the "iploc"th local element of

221 # the "kloc"th local eigenvector

222 #

223 # init local mu and m sums with

224 # a continuous copy of ROW of eigenvectors

225 # i.e. a "ip"th element of each eigenvector

226 # (NOT the "jq"th eigenvector)

227 cloc_dm[:] = self.eigenvectors[lip, :]

228 cloc_magn[:] = self.eigenvectors[lip, :]

229

230 # c1 * c2 * F = sqrt(fdiff_ip * ediff_ip) / sqrt(omega_k) F^(k)_ip

231 cloc_dm *= np.sqrt(kss_ip.pop_diff * kss_ip.energy_diff)

232 cloc_dm /= sqrtwloc

233

234 # c2 / c1 * F = sqrt(fdiff_ip / ediff_ip) * sqrt(omega_k) F^(k)_ip

235 cloc_magn *= np.sqrt(kss_ip.pop_diff / kss_ip.energy_diff)

236 cloc_magn *= sqrtwloc

237

238 if self.custom_axes is None:

239 dmxloc += kss_ip.dip_mom_r[0] * cloc_dm

240 dmyloc += kss_ip.dip_mom_r[1] * cloc_dm

241 dmzloc += kss_ip.dip_mom_r[2] * cloc_dm

242

243 magnxloc += kss_ip.magn_mom[0] * cloc_magn

244 magnyloc += kss_ip.magn_mom[1] * cloc_magn

245 magnzloc += kss_ip.magn_mom[2] * cloc_magn

246 else:

247 dmxloc += np.dot(kss_ip.dip_mom_r,

248 self.custom_axes[0]) * cloc_dm

249 dmyloc += np.dot(kss_ip.dip_mom_r,

250 self.custom_axes[1]) * cloc_dm

251 dmzloc += np.dot(kss_ip.dip_mom_r,

252 self.custom_axes[2]) * cloc_dm

253

254 magnxloc += np.dot(kss_ip.magn_mom,

255 self.custom_axes[0]) * cloc_magn

256 magnyloc += np.dot(kss_ip.magn_mom,

257 self.custom_axes[1]) * cloc_magn

258 magnzloc += np.dot(kss_ip.magn_mom,

259 self.custom_axes[2]) * cloc_magn

260

261 # global properties, but local eigs

262 # sum over iploc dmx_iploc[kloc] to dmx[kloc]

263 props = np.array(

264 [dmxloc, dmyloc, dmzloc, magnxloc, magnyloc, magnzloc])

265 self.lr_comms.eh_comm.sum(props.ravel())

266 [dmxloc, dmyloc, dmzloc, magnxloc, magnyloc, magnzloc] = props

267

268 # global properties and eigs

269 # create local part of the global eigs

270 nrows = len(self.ks_singles.kss_list)

271 self.transition_properties = np.zeros([nrows, 1 + 3 + 3])

272 for k in range(nrows):

273 kloc = self.lr_comms.get_local_dd_index(k)

274 if kloc is None:

275 continue

276 self.transition_properties[k, 0] = sqrtwloc[kloc] * sqrtwloc[kloc]

277 self.transition_properties[k, 1] = dmxloc[kloc]

278 self.transition_properties[k, 2] = dmyloc[kloc]

279 self.transition_properties[k, 3] = dmzloc[kloc]

280 self.transition_properties[k, 4] = magnxloc[kloc]

281 self.transition_properties[k, 5] = magnyloc[kloc]

282 self.transition_properties[k, 6] = magnzloc[kloc]

283

284 # sum over kloc dmx[...kloc...] to dmx[k]

285 self.lr_comms.dd_comm.sum(self.transition_properties.ravel())

286

287 # print self.lr_comms.parent_comm.rank, self.transition_properties[:,1]

288 # return

289

290 self.trans_prop_ready = True

291

292 # omega_k = sqrt(lambda_k)

293 def get_excitation_energy(self, k, units='au'):

294 """Get excitation energy for kth interacting transition

295

296 Input parameters:

297

298 k

299 Transition index (indexing starts from zero).

300

301 units

302 Units for excitation energy: 'au' (Hartree) or 'eV'.

303 """

304 if not self.trans_prop_ready:

305 self.calculate()

306

307 if units == 'au':

308 return np.sqrt(self.eigenvalues[k])

309 elif units == 'eV' or units == 'eVcgs':

310 return np.sqrt(self.eigenvalues[k]) * ase.units.Hartree

311 else:

312 raise RuntimeError('Unknown units.')

313

314 def get_oscillator_strength(self, k):

315 """Get oscillator strength for an interacting transition

316

317 Returns oscillator strength of kth interacting transition.

318

319 Input parameters:

320

321 k

322 Transition index (indexing starts from zero).

323 """

324 if not self.trans_prop_ready:

325 self.calculate()

326

327 omega = self.transition_properties[k][0]

328 dmx = self.transition_properties[k][1]

329 dmy = self.transition_properties[k][2]

330 dmz = self.transition_properties[k][3]

331

332 oscx = 2. * omega * dmx * dmx

333 oscy = 2. * omega * dmy * dmy

334 oscz = 2. * omega * dmz * dmz

335 osc = (oscx + oscy + oscz) / 3.

336 return osc, np.array([oscx, oscy, oscz])

337

338 def get_rotatory_strength(self, k, units='au'):

339 """Get rotatory strength for an interacting transition

340

341 Returns rotatory strength of kth interacting transition.

342

343 Input parameters:

344

345 k

346 Transition index (indexing starts from zero).

347

348 units

349 Units for rotatory strength: 'au' or 'cgs'.

350 """

351

352 if not self.trans_prop_ready:

353 self.calculate()

354

355 # omega = self.transition_properties[k][0]

356 dmx = self.transition_properties[k][1]

357 dmy = self.transition_properties[k][2]

358 dmz = self.transition_properties[k][3]

359 magnx = self.transition_properties[k][4]

360 magny = self.transition_properties[k][5]

361 magnz = self.transition_properties[k][6]

362

363 if units == 'au':

364 return -(dmx * magnx + dmy * magny + dmz * magnz)

365 elif units == 'cgs' or units == 'eVcgs':

366 # 10^-40 esu cm erg / G

367 # = 3.33564095 * 10^-15 A^2 m^3 s

368 # conversion factor 471.43 from

369 # T. B. Pedersen and A. E. Hansen,

370 # Chem. Phys. Lett. 246 (1995) 1

371 # From turbomole 64604.8164

372 #

373 # FIX ME: why?

374 # 64604.8164/471.43 = 137.040

375 # is the inverse of fine-structure constant

376 # OR it must be speed of light in atomic units

377 return -64604.8164 * (dmx * magnx + dmy * magny + dmz * magnz)

378 else:

379 raise RuntimeError('Unknown units.')

380

381 def get_transitions(self,

382 filename=None,

383 min_energy=0.0,

384 max_energy=30.0,

385 units='eVcgs'):

386 """Get transitions: energy, dipole strength and rotatory strength.

387

388 Returns transitions as (w,S,R, Sx,Sy,Sz) where w is an array of

389 frequencies,

390 S is an array of corresponding dipole strengths, and R is an array of

391 corresponding rotatory strengths.

392

393 Input parameters:

394

395 min_energy

396 Minimum energy

397

398 min_energy

399 Maximum energy

400

401 units

402 Units for spectrum: 'au' or 'eVcgs'

403 """

404 if not self.trans_prop_ready:

405 self.calculate()

406

407 w = []

408 S = []

409 R = []

410 Sx = []

411 Sy = []

412 Sz = []

413

414 for (k, omega) in enumerate(self.eigenvalues):

415 ww = self.get_excitation_energy(k, units)

416 if ww < min_energy or ww > max_energy:

417 continue

418

419 w.append(ww)

420 St, Sc = self.get_oscillator_strength(k)

421 S.append(St)

422 Sx.append(Sc[0])

423 Sy.append(Sc[1])

424 Sz.append(Sc[2])

425 R.append(self.get_rotatory_strength(k, units))

426

427 w = np.array(w)

428 S = np.array(S)

429 R = np.array(R)

430 Sx = np.array(Sx)

431 Sy = np.array(Sy)

432 Sz = np.array(Sz)

433

434 if filename is not None and world.rank == 0:

435 sfile = open(filename, 'w')

436 sfile.write("# %12s %12s %12s %12s %12s %12s %s\n" %

437 ('energy', 'osc str', 'rot str', 'osc str x',

438 'osc str y', 'osc str z', 'units: ' + units))

439 for (ww, SS, RR, SSx, SSy, SSz) in zip(w, S, R, Sx, Sy, Sz):

440 sfile.write(

441 (" %12.8lf %12.8lf %12.8lf "

442 "%12.8lf %12.8lf %12.8lf\n")

443 % (ww, SS, RR, SSx, SSy, SSz))

444 sfile.close()

445

446 return (w, S, R, Sx, Sy, Sz)

447

448 def get_spectrum(self,

449 filename=None,

450 min_energy=0.0,

451 max_energy=30.0,

452 energy_step=0.01,

453 width=0.1,

454 units='eVcgs'):

455 """Get spectrum for dipole and rotatory strength.

456

457 Returns folded spectrum as (w,S,R) where w is an array of frequencies,

458 S is an array of corresponding dipole strengths, and R is an array of

459 corresponding rotatory strengths.

460

461 Input parameters:

462

463 min_energy

464 Minimum energy

465

466 min_energy

467 Maximum energy

468

469 energy_step

470 Spacing between calculated energies

471

472 width

473 Width of the Gaussian

474

475 units

476 Units for spectrum: 'au' or 'eVcgs'

477 """

478

479 if not self.trans_prop_ready:

480 self.calculate()

481

482 (ww, SS, RR, SSx, SSy,

483 SSz) = self.get_transitions(min_energy=min_energy - 5 * width,

484 max_energy=max_energy + 5 * width,

485 units=units)

486

487 if units == 'eVcgs':

488 pass # convf = 1/ase.units.Hartree

489 elif units == 'au':

490 assert 0 # convf = 1.

491 else:

492 raise RuntimeError('Invalid units')

493

494 w = np.arange(min_energy, max_energy, energy_step)

495 S = np.zeros_like(w)

496 Sx = np.zeros_like(w)

497 Sy = np.zeros_like(w)

498 Sz = np.zeros_like(w)

499 R = np.zeros_like(w)

500

501 for (k, www) in enumerate(ww):

502

503 c = SS[k] / width / np.sqrt(2 * np.pi)

504 S += c * np.exp((-.5 / width / width) * np.power(w - ww[k], 2))

505

506 c = SSx[k] / width / np.sqrt(2 * np.pi)

507 Sx += c * np.exp((-.5 / width / width) * np.power(w - ww[k], 2))

508

509 c = SSy[k] / width / np.sqrt(2 * np.pi)

510 Sy += c * np.exp((-.5 / width / width) * np.power(w - ww[k], 2))

511

512 c = SSz[k] / width / np.sqrt(2 * np.pi)

513 Sz += c * np.exp((-.5 / width / width) * np.power(w - ww[k], 2))

514

515 c = RR[k] / width / np.sqrt(2 * np.pi)

516 R += c * np.exp((-.5 / width / width) * np.power(w - ww[k], 2))

517

518 if filename is not None and world.rank == 0:

519 sfile = open(filename, 'w')

520 sfile.write("# %12s %12s %12s %12s %12s %12s %s\n" %

521 ('energy', 'osc str', 'rot str', 'osc str x',

522 'osc str y', 'osc str z', 'units: ' + units))

523 for (ww, SS, RR, SSx, SSy, SSz) in zip(w, S, R, Sx, Sy, Sz):

524 sfile.write(

525 (" %12.8lf %12.8lf %12.8lf "

526 "%12.8lf %12.8lf %12.8lf\n")

527 % (ww, SS, RR, SSx, SSy, SSz))

528 sfile.close()

529

530 return (w, S, R, Sx, Sy, Sz)

531

532 def get_transition_contributions(self, index_of_transition):

533 """Get contributions of Kohn-Sham singles to a given transition.

534

535 Includes population difference.

536

537 For large systems this is slow.

538

539 Input parameters:

540

541 index_of_transition:

542 index of transition starting from zero

543 """

544

545 neigs = len(self.ks_singles.kss_list)

546 f2 = np.zeros(neigs)

547

548 if index_of_transition < 0:

549 raise RuntimeError(

550 'Error in get_transition_contributions: Index < zero.')

551 if index_of_transition >= neigs:

552 raise RuntimeError(

553 'Error in get_transition_contributions: '

554 'Index >= number of transitions'

555 )

556

557 k = index_of_transition

558 # local k index

559 kloc = self.lr_comms.get_local_dd_index(k)

560

561 if kloc is not None:

562

563 for (ip, kss_ip) in enumerate(self.ks_singles.kss_list):

564 # local ip index

565 lip = self.lr_comms.get_local_eh_index(ip)

566 if lip is None:

567 continue

568

569 # self.eigenvectors[iploc,kloc]

570 # is the "iploc"th local element of

571 # the "kloc"th local eigenvector

572 f2[ip] = self.eigenvectors[lip, kloc] * self.eigenvectors[

573 lip, kloc] * kss_ip.pop_diff

574

575 self.lr_comms.parent_comm.sum(f2)

576

577 return f2