Coverage for gpaw/lrtddft2/lr

1import numpy as np

3import ase.units

5from gpaw.lrtddft2.ks_singles import KohnShamSingleExcitation

6from gpaw.lrtddft2.lr_layouts import LrTDDFPTSolveLayout

9class LrResponse:

10 def __init__(self,

11 lrtddft2,

12 excitation_energy,

13 field_vector,

14 lorentzian_width,

15 sl_lrtddft=None):

16 self.lrtddft2 = lrtddft2

17 self.omega = excitation_energy

18 self.eta = lorentzian_width

19 self.field_vector = np.array(field_vector)

20 self.sl_lrtddft = sl_lrtddft

22 self.C_re = None

23 self.C_im = None

25 self.response_wf_ready = False

27 def read(self):

28 raise RuntimeError('Error in LrtddftResponseWF read: Not implemented')

30 # TD-DFPT:

31 #

32 # A C = -Vext

33 #

34 # ( E+KN+w KN -eta 0 ) ( C_+w,Re ) ( -dm_laser )

35 # ( KN E+KN-w 0 -eta ) ( C_-w,Re ) = ( -dm_laser )

36 # ( eta 0 E+KN+w -KN ) ( C_+w,Im ) ( 0 )

37 # ( 0 eta -KN E+KN-w ) ( C_-w,Im ) ( 0 )

38 #

39 # N is occupation difference matrix (just like in Casida)

40 # (at least K = 0 gives density diffs for both real and imaginary part)

41 def solve(self):

42 if self.sl_lrtddft is None:

43 self.solve_serial()

44 else:

45 self.solve_parallel()

47 def get_response_coefficients(self, units='eVang'):

48 C_re = self.C_re * 1.

49 C_im = self.C_im * 1.

51 if units == 'au':

52 pass

53 elif units == 'eVang':

54 # FIXME: where these units come from???

55 # S = 2 omega detla n C_im dm

56 # units: 1/eV = eV * e ang * C

57 # => units of C: 1/(eV * eV * e ang)

58 # maybe

59 # psi(w) = psi0(w) + E(w).mu psi1(w) + ...

60 # E(w).mu units = V/ang e ang = eV ?

61 # nope...

62 C_re *= 1. / (ase.units.Hartree**2 * ase.units.Bohr)

63 C_im *= 1. / (ase.units.Hartree**2 * ase.units.Bohr)

64 else:

65 raise RuntimeError(

66 'Error in get_response_coefficients: Invalid units.')

68 return (C_re, C_im)

70 def get_response_data(self, units='eVangcgs'):

71 """Get response data.

73 Returns matrix where transitions are in rows and columns are

74 occupied index, unoccupied index, occupied KS eigenvalue, unoccupied

75 KS eigenvalue, occupied occupation number, unoccupied occupation

76 number,

77 real part of response coefficent, imaginary part of response

78 coefficient, dipole moment x, y, and z-components, magnetic

79 moment x, y,

80 and z-components.

81 """

83 data = np.zeros((len(self.lrtddft2.ks_singles.kss_list),

84 2 + 2 + 2 + 2 + 3 + 3 + 3))

85 kpt_ind = self.lrtddft2.kpt_ind

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

88 i = kss_ip.occ_ind

89 p = kss_ip.unocc_ind

91 eps_i = self.lrtddft2.calc.wfs.kpt_u[kpt_ind].eps_n[i]

92 eps_p = self.lrtddft2.calc.wfs.kpt_u[kpt_ind].eps_n[p]

94 f_i = self.lrtddft2.calc.wfs.kpt_u[kpt_ind].f_n[i]

95 f_p = self.lrtddft2.calc.wfs.kpt_u[kpt_ind].f_n[p]

97 data[ip, 0] = i

98 data[ip, 1] = p

100 data[ip, 2] = eps_i

101 data[ip, 3] = eps_p

102

103 data[ip, 4] = f_i

104 data[ip, 5] = f_p

105

106 data[ip, 6] = self.C_re[ip]

107 data[ip, 7] = self.C_im[ip]

108

109 data[ip, 8] = kss_ip.dip_mom_r[0]

110 data[ip, 9] = kss_ip.dip_mom_r[1]

111 data[ip, 10] = kss_ip.dip_mom_r[2]

112

113 data[ip, 11] = kss_ip.magn_mom[0]

114 data[ip, 12] = kss_ip.magn_mom[1]

115 data[ip, 13] = kss_ip.magn_mom[2]

116

117 df = f_i - f_p

118

119 C_im_ip = self.C_im[ip]

120 data[ip, 14] = df * C_im_ip**2

121

122 dm = np.vdot(kss_ip.dip_mom_r, self.field_vector)

123 data[ip, 15] = 2. * self.omega * df * C_im_ip * dm

124 # data[ip,15] = 2. * (eps_p-eps_i) * df * C_im_ip * dm

125

126 # FIXME: IS THIS CORRECT ???

127 # maybe answer is here

128 # Varsano et al., Phys.Chem.Chem.Phys. 11, 4481 (2009)

129 mm = np.vdot(kss_ip.magn_mom, self.field_vector)

130 data[ip, 16] = -df * C_im_ip * mm

131

132 if units == 'au':

133 pass

134 elif units == 'eVangcgs':

135 data[:, 2] *= ase.units.Hartree

136 data[:, 3] *= ase.units.Hartree

137 data[:, 6] *= 1. / (ase.units.Hartree**2 * ase.units.Bohr)

138 data[:, 7] *= 1. / (ase.units.Hartree**2 * ase.units.Bohr)

139 data[:, 8] *= ase.units.Bohr

140 data[:, 9] *= ase.units.Bohr

141 data[:, 10] *= ase.units.Bohr

142 data[:, 11] *= float('NaN') # FIXME: units?

143 data[:, 12] *= float('NaN') # FIXME: units?

144 data[:, 13] *= float('NaN') # FIXME: units?

145 data[:, 14] *= float('NaN') # FIXME: units?

146 data[:, 15] *= float('NaN') # FIXME: units?

147 data[:, 16] *= float('NaN') # FIXME: units?

148 raise RuntimeError(

149 'Error in get_response_data: Unit conversion from '

150 'atomic units to eV ang cgs not implemented yet.'

151 )

152 else:

153 raise RuntimeError('Error in get_response_data: Invalid units.')

154

155 return data

156

157 # amplitude, dipole, rotatory

158 # dipole should integrate to absorption spectrum

159 # rotatory should integrate to CD spectrum

160 def get_transition_contribution_maps(self,

161 occ_min_energy,

162 occ_max_energy,

163 unocc_min_energy,

164 unocc_max_energy,

165 width=0.05,

166 energy_step=0.05,

167 units='eVangcgs'):

168 data = self.get_response_data(units='au')

169

170 eps_i = data[:, 2]

171 eps_p = data[:, 3]

172 a_ip = data[:, 14]

173 s_ip = data[:, 15]

174 r_ip = data[:, 16]

175

176 # occupied energy range

177 wx = np.arange(occ_min_energy, occ_max_energy, energy_step)

178 # unoccupied energy range

179 wy = np.arange(unocc_min_energy, unocc_max_energy, energy_step)

180

181 A = np.outer(wx, wy) * 0.

182 S = np.outer(wx, wy) * 0.

183 R = np.outer(wx, wy) * 0.

184

185 # A(omega_x, omega_y) =

186 # sum_ip delta n_ip |C^(im)_ip|**2

187 # * g(omega_x - eps_i) g(omega_y - eps_p)

188 #

189 # D(omega_x, omega_y) =

190 # sum_ip 2 * delta n_ip C^(im)_ip * mu_ip

191 # * g(omega_x - eps_i) g(omega_y - eps_p)

192 #

193 # where g(x) is gaussian

194 for (ei, ep, a, s, r) in zip(eps_i, eps_p, a_ip, s_ip, r_ip):

195 gx = np.exp((-.5 / width / width) *

196 np.power(wx - ei, 2)) / width / np.sqrt(2 * np.pi)

197 gy = np.exp((-.5 / width / width) *

198 np.power(wy - ep, 2)) / width / np.sqrt(2 * np.pi)

199

200 A += a * np.outer(gx, gy)

201 S += s * np.outer(gx, gy)

202 R += r * np.outer(gx, gy)

203

204 # FIXME: units of A, S, and R

205

206 if units == 'au':

207 pass

208 elif units == 'eVangcgs':

209 raise RuntimeError(

210 'Error in get_transition_contribution_maps: Unit conversion '

211 'from atomic units to eV ang cgs not implemented yet.'

212 )

213 else:

214 raise RuntimeError(

215 'Error in get_transition_contribution_maps: Invalid units.')

216

217 return (wx, wy, A, S, R)

218

219 def get_induced_density(self,

220 units='au',

221 collect=False,

222 amplitude_filter=1e-5):

223 # Init pair densities

224 dnt_Gip = self.lrtddft2.calc.wfs.gd.empty()

225 dnt_gip = self.lrtddft2.calc.density.finegd.empty()

226 drhot_gip = self.lrtddft2.calc.density.finegd.empty()

227

228 drhot_g = self.lrtddft2.calc.density.finegd.empty()

229 drhot_g[:] = 0.0

230

231 C_im = self.C_im

232

233 maxC = np.max(abs(C_im))

234

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

236 # distribute ips over eh comm

237 if self.lrtddft2.lr_comms.get_local_eh_index(ip) is None:

238 continue

239 if abs(C_im[ip]) < amplitude_filter * maxC:

240 continue

241 dnt_Gip[:] = 0.0

242 dnt_gip[:] = 0.0

243 drhot_gip[:] = 0.0

244 kss_ip.calculate_pair_density(dnt_Gip, dnt_gip, drhot_gip)

245 # ?FIXME?: SHOULD HERE BE OCCUPATION DIFFERENCE? yes ####

246 drhot_g += kss_ip.pop_diff * C_im[ip] * drhot_gip

247

248 # sum over eh comm

249 self.lrtddft2.lr_comms.eh_comm.sum(drhot_g)

250

251 if collect:

252 drhot_g = self.lrtddft2.calc.density.finegd.collect(drhot_g)

253

254 if units == 'au':

255 pass

256 elif units == 'ang':

257 drhot_g *= 1. / ase.units.Bohr**3

258 else:

259 raise RuntimeError('Error in get_induced_density: Invalid units.')

260

261 return drhot_g

262

263 def get_approximate_electron_and_hole_densities(self,

264 units='au',

265 collect=False,

266 amplitude_filter=1e-4):

267 # Init pair densities

268 dnt_Gip = self.lrtddft2.calc.wfs.gd.empty()

269 dnt_gip = self.lrtddft2.calc.density.finegd.empty()

270 drhot_gip = self.lrtddft2.calc.density.finegd.empty()

271

272 drhot_ge = self.lrtddft2.calc.density.finegd.empty()

273 drhot_ge[:] = 0.0

274 drhot_gh = self.lrtddft2.calc.density.finegd.empty()

275 drhot_gh[:] = 0.0

276

277 #

278 # Sum of Slater determinants:

279 # diag ip ip => |c_ip|**2 ( |psi_p|**2 - |psi_i]**2 )

280 # offdiag ip iq => c_ip c_iq psi_p psi_q

281 # offdiag ip jp => - c_ip c_jp psi_i psi_j

282 # offdiag ip jq => 0

283 #

284 # remember both ip,jq and jq,ip

285 #

286

287 # occupations for electron and hole

288 f_n = self.lrtddft2.calc.wfs.kpt_u[self.lrtddft2.kpt_ind].f_n

289 fe_n = np.zeros(len(f_n))

290 fh_n = np.zeros(len(f_n))

291

292 C_im = self.C_im

293

294 maxC = np.max(abs(C_im))

295

296 # diagonal ip,ip

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

298 if self.lrtddft2.lr_comms.get_local_eh_index(ip) is None:

299 continue

300 if abs(C_im[ip]) < amplitude_filter * maxC:

301 continue

302

303 # decrease in density

304 fh_n[kss_ip.occ_ind] -= C_im[ip] * C_im[ip] * kss_ip.pop_diff

305 # increase in density

306 fe_n[kss_ip.unocc_ind] += C_im[ip] * C_im[ip] * kss_ip.pop_diff

307

308 for k in range(len(f_n)):

309 if (abs(fh_n[k]) < amplitude_filter * maxC * maxC

310 and abs(fe_n[k]) < amplitude_filter * maxC * maxC):

311 continue

312

313 dnt_Gip[:] = 0.0

314 dnt_gip[:] = 0.0

315 drhot_gip[:] = 0.0

316

317 kss_ip = KohnShamSingleExcitation(self.lrtddft2.calc,

318 self.lrtddft2.kpt_ind, k, k)

319 kss_ip.calculate_pair_density(dnt_Gip, dnt_gip, drhot_gip)

320

321 # fh_n < 0

322 drhot_gh += fh_n[k] * drhot_gip

323 # fe_n > 0

324 drhot_ge += fe_n[k] * drhot_gip

325

326 # offdiagonal

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

328 if self.lrtddft2.lr_comms.get_local_eh_index(ip) is None:

329 continue

330 if abs(C_im[ip]) < amplitude_filter * maxC:

331 continue

332

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

334 # if diagonal, skip because already done

335 if (kss_ip.occ_ind == kss_jq.occ_ind

336 and kss_ip.unocc_ind == kss_jq.unocc_ind):

337 continue

338

339 # if both are different, then orthogonal

340 if (kss_ip.occ_ind != kss_jq.occ_ind

341 and kss_ip.unocc_ind != kss_jq.unocc_ind):

342 continue

343

344 if abs(C_im[ip] * C_im[jq]) < amplitude_filter * maxC * maxC:

345 continue

346

347 # if occupied state is the same, then only electron density

348 # changes (ie use kss_??.unocc_ind)

349 # this gives direction, and should integrate to zero

350 if (kss_ip.occ_ind == kss_jq.occ_ind

351 and kss_ip.unocc_ind != kss_jq.unocc_ind):

352 dnt_Gip[:] = 0.0

353 dnt_gip[:] = 0.0

354 drhot_gip[:] = 0.0

355

356 kss_pq = KohnShamSingleExcitation(self.lrtddft2.calc,

357 self.lrtddft2.kpt_ind,

358 kss_ip.unocc_ind,

359 kss_jq.unocc_ind)

360 kss_pq.calculate_pair_density(dnt_Gip, dnt_gip, drhot_gip)

361

362 drhot_ge += (C_im[ip] * C_im[jq] * np.sqrt(

363 kss_ip.pop_diff * kss_jq.pop_diff)) * drhot_gip

364

365 # if occupied state is the same, then only hole density

366 # changes (ie use kss_??.occ_ind)

367 # this gives direction, and should integrate to zero

368 if (kss_ip.occ_ind != kss_jq.occ_ind

369 and kss_ip.unocc_ind == kss_jq.unocc_ind):

370 dnt_Gip[:] = 0.0

371 dnt_gip[:] = 0.0

372 drhot_gip[:] = 0.0

373 kss_ij = KohnShamSingleExcitation(self.lrtddft2.calc,

374 self.lrtddft2.kpt_ind,

375 kss_ip.occ_ind,

376 kss_jq.occ_ind)

377 kss_ij.calculate_pair_density(dnt_Gip, dnt_gip, drhot_gip)

378

379 drhot_gh -= (C_im[ip] * C_im[jq] * np.sqrt(

380 kss_ip.pop_diff * kss_jq.pop_diff)) * drhot_gip

381

382 self.lrtddft2.lr_comms.eh_comm.sum(drhot_ge)

383 self.lrtddft2.lr_comms.eh_comm.sum(drhot_gh)

384

385 drhot_geh = drhot_ge + drhot_gh

386

387 # Ige = self.lrtddft2.calc.density.finegd.integrate(drhot_ge)

388 # Igh = self.lrtddft2.calc.density.finegd.integrate(drhot_gh)

389 # Igeh = self.lrtddft2.calc.density.finegd.integrate(drhot_geh)

390

391 # if self.lrtddft2.lr_comms.parent_comm.rank == 0:

392 # print 'drho_ge ', Ige

393 # print 'drho_gh ', Igh

394 # print 'drho_geh', Igeh

395

396 if collect:

397 drhot_ge = self.lrtddft2.calc.density.finegd.collect(drhot_ge)

398 drhot_gh = self.lrtddft2.calc.density.finegd.collect(drhot_gh)

399 drhot_geh = self.lrtddft2.calc.density.finegd.collect(drhot_geh)

400

401 if units == 'au':

402 pass

403 elif units == 'ang':

404 drhot_geh *= 1. / ase.units.Bohr**3

405 drhot_gh *= 1. / ase.units.Bohr**3

406 drhot_ge *= 1. / ase.units.Bohr**3

407 else:

408 raise RuntimeError(

409 'Error in get_approximate_electron_and_hole_densities: '

410 'Invalid units.'

411 )

412

413 return (drhot_ge, drhot_gh, drhot_geh)

414

415 def solve_serial(self):

416 nrows = len(self.lrtddft2.ks_singles.kss_list)

417

418 A_matrix = np.zeros((nrows * 4, nrows * 4))

419 K_matrix = self.lrtddft2.K_matrix.values

420

421 # add K N

422 # slow but if you're using serial it's ok anyway

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

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

425 lip = self.lrtddft2.lr_comms.get_local_eh_index(ip)

426 ljq = self.lrtddft2.lr_comms.get_local_dd_index(jq)

427 if lip is None or ljq is None:

428 continue

429

430 A_matrix[ip * 4 + 0,

431 jq * 4 + 0] = K_matrix[lip, ljq] * kss_jq.pop_diff

432 A_matrix[ip * 4 + 0,

433 jq * 4 + 1] = K_matrix[lip, ljq] * kss_jq.pop_diff

434 A_matrix[ip * 4 + 1,

435 jq * 4 + 0] = K_matrix[lip, ljq] * kss_jq.pop_diff

436 A_matrix[ip * 4 + 1,

437 jq * 4 + 1] = K_matrix[lip, ljq] * kss_jq.pop_diff

438

439 A_matrix[ip * 4 + 2,

440 jq * 4 + 2] = K_matrix[lip, ljq] * kss_jq.pop_diff

441 A_matrix[ip * 4 + 2,

442 jq * 4 + 3] = -K_matrix[lip, ljq] * kss_jq.pop_diff

443 A_matrix[ip * 4 + 3,

444 jq * 4 + 2] = -K_matrix[lip, ljq] * kss_jq.pop_diff

445 A_matrix[ip * 4 + 3,

446 jq * 4 + 3] = K_matrix[lip, ljq] * kss_jq.pop_diff

447

448 # diagonal stuff: E, w, and eta

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

450 lip = self.lrtddft2.lr_comms.get_local_eh_index(ip)

451 ljq = self.lrtddft2.lr_comms.get_local_dd_index(ip)

452 if lip is None or ljq is None:

453 continue

454

455 # energy difference and excitation energy

456 A_matrix[ip * 4 + 0, ip * 4 + 0] += kss_ip.energy_diff + self.omega

457 A_matrix[ip * 4 + 1, ip * 4 + 1] += kss_ip.energy_diff - self.omega

458 A_matrix[ip * 4 + 2, ip * 4 + 2] += kss_ip.energy_diff + self.omega

459 A_matrix[ip * 4 + 3, ip * 4 + 3] += kss_ip.energy_diff - self.omega

460

461 # life-time parameter

462 A_matrix[ip * 4 + 0, ip * 4 + 2] += -self.eta

463 A_matrix[ip * 4 + 1, ip * 4 + 3] += -self.eta

464 A_matrix[ip * 4 + 2, ip * 4 + 0] += self.eta

465 A_matrix[ip * 4 + 3, ip * 4 + 1] += self.eta

466

467 # RHS

468 V_rhs = np.zeros(nrows * 4)

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

470 ljq = self.lrtddft2.lr_comms.get_local_dd_index(jq)

471 if ljq is None:

472 continue

473

474 # 2 cos(wt) = exp(+iwt) + exp(iwt)

475 # fixme: no minus here???

476 # -dm_laser = - (-er) cos(wt) ????

477 V_rhs[jq * 4 + 0] = np.dot(kss_jq.dip_mom_r, self.field_vector)

478 V_rhs[jq * 4 + 1] = np.dot(kss_jq.dip_mom_r, self.field_vector)

479 V_rhs[jq * 4 + 2] = 0.

480 V_rhs[jq * 4 + 3] = 0.

481

482 # no transpose

483 # A_matrix[:] = A_matrix.transpose()

484

485 # solve A C = V

486 self.lrtddft2.lr_comms.parent_comm.sum(A_matrix)

487 self.lrtddft2.lr_comms.dd_comm.sum(V_rhs)

488

489 # np.set_printoptions(precision=5, suppress=True)

490 # print A_matrix

491 C = np.linalg.solve(A_matrix, V_rhs)

492 # print C

493

494 # response wavefunction

495 # perturbation was exp(iwt) + exp(-iwt) = 2 cos(wt)

496 # => real part of the polarizability is propto 2 cos(wt)

497 # => imag part of the polarizability is propto 2 sin(wt)

498 #

499 # exp(iwt)

500 # dn+ = phi-**h phi0 + phi0**h phi+ = phi0 (phi-**h + phi+)

501 # exp(-iwt)

502 # dn- = phi+**h phi0 + phi0**h phi- = phi0 (phi+**h + phi-) = dn+**h

503 #

504 # 2 cos(wt)

505 # dn_2cos = dn+ + dn- = dn+ + dn+**h = Re[dn+] + Re[dn-]

506 # 2 sin(wt)

507 # dn_2sin = -i(dn+ - dn-) = -i (dn+ - dn-**h) = Im[dn+] - Im[dn-]

508

509 self.C_re = np.zeros(len(self.lrtddft2.ks_singles.kss_list))

510 self.C_im = np.zeros(len(self.lrtddft2.ks_singles.kss_list))

511

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

513 # Re[phi-**h + phi_+]xo

514 self.C_re[jq] = C[jq * 4 + 0] + C[jq * 4 + 1]

515 # Im[phi-**h + phi_+]

516 self.C_im[jq] = C[jq * 4 + 2] - C[jq * 4 + 3]

517

518 # normalization... where it comes from? fourier (or lorentzian)...

519 self.C_re *= 1. / np.pi

520 self.C_im *= 1. / np.pi

521

522 # anyway, you can check that

523 #

524 # osc_z(omega) / (pi eta) =

525 # 2 * omega * sum_ip n_ip C^(im)_ip(omega) * mu^(z)_ip

526 #

527 # where osc_z, you can get from Casida

528 # and 1/(pi eta) is because of Lorentzian folding

529

530 # spectrum

531 # S_z(omega) = 2 * omega sum_ip n_ip C^(im)_ip(omega) * mu^(z)_ip

532

533 return (self.C_re, self.C_im)

534

535 def solve_parallel(self):

536 # total rows

537 nrows = len(self.lrtddft2.ks_singles.kss_list)

538 # local

539 nlrows = self.lrtddft2.K_matrix.values.shape[0]

540 nlcols = self.lrtddft2.K_matrix.values.shape[1]

541

542 # create BLACS layout

543 layout = LrTDDFPTSolveLayout(self.sl_lrtddft, nrows,

544 self.lrtddft2.lr_comms)

545 if (nrows <

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

547 raise RuntimeError(

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

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

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

551 )

552

553 # build local part

554 # NOTE: scalapack needs TRANSPOSED matrix!!!

555 K_matrix = self.lrtddft2.K_matrix.values

556 KN_matrix_T = np.zeros([nlcols, nlrows])

557 A_matrix_T = np.zeros((nlcols * 4, nlrows * 4))

558

559 KN_matrix_T[:] = K_matrix.T

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

561 ljq = self.lrtddft2.lr_comms.get_local_dd_index(jq)

562 if ljq is None:

563 continue

564 KN_matrix_T[ljq, :] = K_matrix[:, ljq] * kss_jq.pop_diff

565

566 lip_list = []

567 ljq_list = []

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

569 lip = self.lrtddft2.lr_comms.get_local_eh_index(ip)

570 ljq = self.lrtddft2.lr_comms.get_local_dd_index(ip)

571 if lip is not None:

572 lip_list.append(lip)

573 if ljq is not None:

574 ljq_list.append(ljq)

575

576 lip4_list = np.array(lip_list) * 4

577

578 # add K N

579 for ljq in ljq_list:

580 A_matrix_T[ljq * 4 + 0, lip4_list + 0] = KN_matrix_T[ljq, :]

581 A_matrix_T[ljq * 4 + 0, lip4_list + 1] = KN_matrix_T[ljq, :]

582 A_matrix_T[ljq * 4 + 1, lip4_list + 0] = KN_matrix_T[ljq, :]

583 A_matrix_T[ljq * 4 + 1, lip4_list + 1] = KN_matrix_T[ljq, :]

584

585 A_matrix_T[ljq * 4 + 2, lip4_list + 2] = KN_matrix_T[ljq, :]

586 A_matrix_T[ljq * 4 + 2, lip4_list + 3] = -KN_matrix_T[ljq, :]

587 A_matrix_T[ljq * 4 + 3, lip4_list + 2] = -KN_matrix_T[ljq, :]

588 A_matrix_T[ljq * 4 + 3, lip4_list + 3] = KN_matrix_T[ljq, :]

589

590 # diagonal stuff: E, w, and eta

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

592 lip = self.lrtddft2.lr_comms.get_local_eh_index(ip)

593 ljq = self.lrtddft2.lr_comms.get_local_dd_index(ip)

594 if lip is None or ljq is None:

595 continue

596

597 # energy difference and excitation energy

598 A_matrix_T[ljq * 4 + 0,

599 lip * 4 + 0] += kss_ip.energy_diff + self.omega

600 A_matrix_T[ljq * 4 + 1,

601 lip * 4 + 1] += kss_ip.energy_diff - self.omega

602 A_matrix_T[ljq * 4 + 2,

603 lip * 4 + 2] += kss_ip.energy_diff + self.omega

604 A_matrix_T[ljq * 4 + 3,

605 lip * 4 + 3] += kss_ip.energy_diff - self.omega

606

607 # life-time parameter, transposed again!

608 A_matrix_T[ljq * 4 + 0, lip * 4 + 2] += self.eta

609 A_matrix_T[ljq * 4 + 1, lip * 4 + 3] += self.eta

610 A_matrix_T[ljq * 4 + 2, lip * 4 + 0] += -self.eta

611 A_matrix_T[ljq * 4 + 3, lip * 4 + 1] += -self.eta

612

613 # RHS

614 V_rhs = []

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

616 if (jq % self.lrtddft2.lr_comms.parent_comm.size !=

617 self.lrtddft2.lr_comms.parent_comm.rank):

618 continue

619

620 # 2 cos(wt) = exp(+iwt) + exp(iwt)

621 # fixme: no minus here???

622 # -dm_laser = - (-er) cos(wt) ????

623 V_rhs.append(np.dot(kss_jq.dip_mom_r, self.field_vector))

624 V_rhs.append(np.dot(kss_jq.dip_mom_r, self.field_vector))

625 V_rhs.append(0.)

626 V_rhs.append(0.)

627

628 V_rhs_T = np.array([V_rhs])

629

630 # solve A C = V

631 C = layout.solve(A_matrix_T, V_rhs_T)[0]

632

633 # see solve serial for comments

634 self.C_re = np.zeros(len(self.lrtddft2.ks_singles.kss_list))

635 self.C_im = np.zeros(len(self.lrtddft2.ks_singles.kss_list))

636

637 l = 0

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

639 if (jq % self.lrtddft2.lr_comms.parent_comm.size !=

640 self.lrtddft2.lr_comms.parent_comm.rank):

641 continue

642 # Re[phi-**h + phi_+]

643 self.C_re[jq] = C[l * 4 + 0] + C[l * 4 + 1]

644 # Im[phi-**h + phi_+]

645 self.C_im[jq] = C[l * 4 + 2] - C[l * 4 + 3]

646 l += 1

647

648 self.lrtddft2.lr_comms.parent_comm.sum(self.C_re)

649 self.lrtddft2.lr_comms.parent_comm.sum(self.C_im)

650

651 # normalization... where it comes from? fourier (or lorentzian)...

652 self.C_re *= 1. / np.pi

653 self.C_im *= 1. / np.pi

654

655 return (self.C_re, self.C_im)

Coverage for gpaw/lrtddft2/lr_response.py: 5%

325 statements