Coverage for gpaw/xc/gllb/c

1import warnings

2from math import sqrt, pi

3import numpy as np

5from ase.units import Ha

6from gpaw import BadParallelization

7from gpaw.mpi import world

8from gpaw.density import redistribute_array, redistribute_atomic_matrices

9from gpaw.sphere.lebedev import weight_n

10from gpaw.utilities import (pack_atomic_matrices, pack_density,

11 unpack_atomic_matrices)

12from gpaw.xc.gllb import safe_sqr

13from gpaw.xc.gllb.contribution import Contribution

15# XXX Work in process

16debug = False

19def d(*args):

20 if debug:

21 print(args)

24class ResponsePotential:

25 """Container for response potential"""

26 def __init__(self, response, vt_sG, vt_sg, D_asp, Dresp_asp):

27 self.response = response

28 self.vt_sG = vt_sG

29 self.vt_sg = vt_sg

30 self.D_asp = D_asp

31 self.Dresp_asp = Dresp_asp

33 def redistribute(self, new_response):

34 old_response = self.response

35 new_wfs = new_response.wfs

36 new_density = new_response.density

37 self.vt_sG = redistribute_array(self.vt_sG,

38 old_response.density.gd,

39 new_density.gd,

40 new_wfs.nspins,

41 new_wfs.kptband_comm)

42 if self.vt_sg is not None:

43 self.vt_sg = redistribute_array(self.vt_sg,

44 old_response.density.finegd,

45 new_density.finegd,

46 new_wfs.nspins,

47 new_wfs.kptband_comm)

49 def redist_asp(D_asp):

50 return redistribute_atomic_matrices(D_asp,

51 new_density.gd,

52 new_wfs.nspins,

53 new_density.setups,

54 new_density.atom_partition,

55 new_wfs.kptband_comm)

57 self.D_asp = redist_asp(self.D_asp)

58 self.Dresp_asp = redist_asp(self.Dresp_asp)

59 self.response = new_response

62class C_Response(Contribution):

63 """Response contribution for GLLB functionals.

65 Parameters

66 ----------

67 weight

68 Weight of the contribution

69 coefficients

70 Coefficient calculator object

71 metallic

72 If True, then Fermi level is used as the reference

73 energy for coefficients instead of the HOMO energy.

74 This is necessary to get sensible results in metallic systems.

75 damp

76 Small value to damp divisions by zero

77 """

78 def __init__(self,

79 weight: float,

80 coefficients, *,

81 metallic: bool = False,

82 damp: float = 1e-10):

83 Contribution.__init__(self, weight)

84 d('In c_Response __init__', self)

85 self.coefficients = coefficients

86 self.vt_sg = None

87 self.vt_sG = None

88 self.D_asp = None

89 self.Dresp_asp = None

90 self.Ddist_asp = None

91 self.Drespdist_asp = None

92 self.damp = damp

93 self.fix_potential = False

94 self.metallic = metallic

96 # For logging reference energy

97 self.eref_s = None

98 self.eref_source_s = None

100 def set_damp(self, damp):

101 self.damp = damp

102

103 def get_name(self):

104 return 'RESPONSE'

105

106 def get_desc(self):

107 desc = []

108 if self.metallic:

109 desc += ['metallic']

110 desc += [self.coefficients.get_description()]

111 return ', '.join(desc)

112

113 def initialize(self, density, hamiltonian, wfs):

114 Contribution.initialize(self, density, hamiltonian, wfs)

115 self.coefficients.initialize(wfs)

116 if self.Dresp_asp is None:

117 assert self.density.D_asp is None

118 # XXX With the deprecated `fixdensity=True` option this

119 # initialize() function is called both before AND after read()!

120 # Thus, the second call of initialize() would discard the read

121 # potential unless we check if it was already allocated.

122 # Remove this if statement once `fixdensity=True` option has

123 # been removed from the calculator.

124 if self.vt_sg is None:

125 self.vt_sG = self.gd.empty(self.nspins)

126 self.vt_sg = self.finegd.empty(self.nspins)

127

128 def initialize_1d(self, ae):

129 Contribution.initialize_1d(self, ae)

130 self.coefficients.initialize_1d(ae)

131

132 def initialize_from_other_response(self, response):

133 pot = ResponsePotential(response,

134 response.vt_sG,

135 response.vt_sg,

136 response.D_asp,

137 response.Dresp_asp)

138 pot.redistribute(self)

139 self.vt_sG = pot.vt_sG

140 self.vt_sg = pot.vt_sg

141 self.D_asp = pot.D_asp

142 self.Dresp_asp = pot.Dresp_asp

143 self.Ddist_asp = self.distribute_D_asp(pot.D_asp)

144 self.Drespdist_asp = self.distribute_D_asp(pot.Dresp_asp)

145 self.eref_s = response.eref_s

146 self.eref_source_s = response.eref_source_s

147

148 # Calcualte the GLLB potential and energy 1d

149 def add_xc_potential_and_energy_1d(self, v_g):

150 w_i = self.coefficients.get_coefficients_1d()

151 u2_j = safe_sqr(self.ae.u_j)

152 v_g += self.weight * np.dot(w_i, u2_j) / (

153 np.dot(self.ae.f_j, u2_j) + self.damp)

154 return 0.0 # Response part does not contribute to energy

155

156 def update_potentials(self):

157 d('In update response potential')

158 if self.fix_potential:

159 # Skip the evaluation of the potential so that

160 # the existing potential is used. This is relevant for

161 # band structure calculation so that the potential

162 # does not get updated with the other k-point sampling.

163 return

164

165 if self.wfs.kpt_u[0].eps_n is None:

166 # This skips the evaluation of the potential so that

167 # the existing one is used.

168 # This happens typically after reading when occupations

169 # haven't been calculated yet and the potential read earlier

170 # is used.

171 return

172

173 # Calculate reference energy

174 self.eref_s = []

175 self.eref_source_s = []

176 if self.metallic:

177 # Use Fermi level as reference levels

178 fermilevel = self.wfs.fermi_level

179 assert isinstance(fermilevel, float), \

180 'GLLBSCM supports only a single Fermi level'

181 for s in range(self.wfs.nspins):

182 self.eref_source_s.append('Fermi level')

183 self.eref_s.append(fermilevel)

184 else:

185 # Find homo and lumo levels for each spin

186 for s in range(self.wfs.nspins):

187 homo, lumo = self.wfs.get_homo_lumo(s, _gllb=True)

188 # Check that homo and lumo are reasonable

189 if homo > lumo:

190 # HOMO higher than LUMO; set Fermi level as reference

191 fermilevel = self.wfs.fermi_level

192 self.eref_source_s.append('Fermi level')

193 self.eref_s.append(fermilevel)

194 else:

195 self.eref_source_s.append('HOMO')

196 self.eref_s.append(homo)

197

198 w_kn = self.coefficients.get_coefficients(self.wfs.kpt_u,

199 eref_s=self.eref_s)

200 f_kn = [kpt.f_n for kpt in self.wfs.kpt_u]

201 if w_kn is not None:

202 self.vt_sG[:] = 0.0

203 nt_sG = self.gd.zeros(self.nspins)

204

205 for kpt, w_n in zip(self.wfs.kpt_u, w_kn):

206 self.wfs.add_to_density_from_k_point_with_occupation(

207 self.vt_sG, kpt, w_n)

208 self.wfs.add_to_density_from_k_point(nt_sG, kpt)

209

210 self.wfs.kptband_comm.sum(nt_sG)

211 self.wfs.kptband_comm.sum(self.vt_sG)

212

213 if self.wfs.kd.symmetry:

214 for nt_G, vt_G in zip(nt_sG, self.vt_sG):

215 self.wfs.kd.symmetry.symmetrize(nt_G, self.gd)

216 self.wfs.kd.symmetry.symmetrize(vt_G, self.gd)

217

218 d('response update D_asp', world.rank, self.Dresp_asp.keys(),

219 self.D_asp.keys())

220 self.wfs.calculate_atomic_density_matrices_with_occupation(

221 self.Dresp_asp, w_kn)

222 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp)

223 d('response update Drespdist_asp', world.rank,

224 self.Dresp_asp.keys(), self.D_asp.keys())

225 self.wfs.calculate_atomic_density_matrices_with_occupation(

226 self.D_asp, f_kn)

227 self.Ddist_asp = self.distribute_D_asp(self.D_asp)

228

229 self.vt_sG /= nt_sG + self.damp

230

231 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg)

232

233 def calculate(self, e_g, n_sg, v_sg):

234 self.update_potentials()

235 v_sg += self.weight * self.vt_sg

236

237 def distribute_D_asp(self, D_asp):

238 return self.density.atomdist.to_work(D_asp)

239

240 def calculate_energy_and_derivatives(self, setup, D_sp, H_sp, a,

241 addcoredensity=True):

242 # Get the XC-correction instance

243 c = setup.xc_correction

244 ncresp_g = setup.extra_xc_data['core_response'] / self.nspins

245 if not addcoredensity:

246 ncresp_g[:] = 0.0

247

248 for D_p, dEdD_p, Dresp_p in zip(self.Ddist_asp[a], H_sp,

249 self.Drespdist_asp[a]):

250 D_Lq = np.dot(c.B_pqL.T, D_p)

251 n_Lg = np.dot(D_Lq, c.n_qg) # Construct density

252 if addcoredensity:

253 n_Lg[0] += c.nc_g * sqrt(4 * pi) / self.nspins

254 nt_Lg = np.dot(

255 D_Lq, c.nt_qg

256 ) # Construct smooth density (_without smooth core density_)

257

258 Dresp_Lq = np.dot(c.B_pqL.T, Dresp_p)

259 nresp_Lg = np.dot(Dresp_Lq, c.n_qg) # Construct 'response density'

260 nrespt_Lg = np.dot(

261 Dresp_Lq, c.nt_qg

262 ) # Construct smooth 'response density' (w/o smooth core)

263

264 for w, Y_L in zip(weight_n, c.Y_nL):

265 nt_g = np.dot(Y_L, nt_Lg)

266 nrespt_g = np.dot(Y_L, nrespt_Lg)

267 x_g = nrespt_g / (nt_g + self.damp)

268 dEdD_p -= self.weight * w * np.dot(

269 np.dot(c.B_pqL, Y_L), np.dot(c.nt_qg, x_g * c.rgd.dv_g))

270

271 n_g = np.dot(Y_L, n_Lg)

272 nresp_g = np.dot(Y_L, nresp_Lg)

273 x_g = (nresp_g + ncresp_g) / (n_g + self.damp)

274

275 dEdD_p += self.weight * w * np.dot(

276 np.dot(c.B_pqL, Y_L), np.dot(c.n_qg, x_g * c.rgd.dv_g))

277

278 return 0.0

279

280 def integrate_sphere(self, a, Dresp_sp, D_sp, Dwf_p, spin=0):

281 c = self.wfs.setups[a].xc_correction

282 Dresp_p, D_p = Dresp_sp[spin], D_sp[spin]

283 D_Lq = np.dot(c.B_pqL.T, D_p)

284 n_Lg = np.dot(D_Lq, c.n_qg) # Construct density

285 n_Lg[0] += c.nc_g * sqrt(4 * pi) / len(D_sp)

286 nt_Lg = np.dot(D_Lq, c.nt_qg

287 ) # Construct smooth density (without smooth core)

288 Dresp_Lq = np.dot(c.B_pqL.T, Dresp_p) # Construct response

289 nresp_Lg = np.dot(Dresp_Lq, c.n_qg) # Construct 'response density'

290 nrespt_Lg = np.dot(

291 Dresp_Lq, c.nt_qg

292 ) # Construct smooth 'response density' (w/o smooth core)

293 Dwf_Lq = np.dot(c.B_pqL.T, Dwf_p) # Construct lumo wf

294 nwf_Lg = np.dot(Dwf_Lq, c.n_qg)

295 nwft_Lg = np.dot(Dwf_Lq, c.nt_qg)

296 E = 0.0

297 for w, Y_L in zip(weight_n, c.Y_nL):

298 v = np.dot(Y_L, nwft_Lg) * np.dot(Y_L, nrespt_Lg) / (

299 np.dot(Y_L, nt_Lg) + self.damp)

300 E -= self.weight * w * np.dot(v, c.rgd.dv_g)

301 v = np.dot(Y_L, nwf_Lg) * np.dot(Y_L, nresp_Lg) / (

302 np.dot(Y_L, n_Lg) + self.damp)

303 E += self.weight * w * np.dot(v, c.rgd.dv_g)

304 return E

305

306 def add_smooth_xc_potential_and_energy_1d(self, vt_g):

307 w_ln = self.coefficients.get_coefficients_1d(smooth=True)

308 v_g = np.zeros(self.ae.N)

309 n_g = np.zeros(self.ae.N)

310 for w_n, f_n, u_n in zip(w_ln, self.ae.f_ln,

311 self.ae.s_ln): # For each angular momentum

312 u2_n = safe_sqr(u_n)

313 v_g += np.dot(w_n, u2_n)

314 n_g += np.dot(f_n, u2_n)

315

316 vt_g += self.weight * v_g / (n_g + self.damp)

317 return 0.0 # Response part does not contribute to energy

318

319 def calculate_delta_xc(self, homolumo=None):

320 warnings.warn(

321 'The function calculate_delta_xc() is deprecated. '

322 'Please use calculate_discontinuity_potential() instead. '

323 'See documentation on calculating band gap with GLLBSC.',

324 DeprecationWarning)

325

326 if homolumo is None:

327 # Calculate band gap

328

329 # This always happens, so we don't warn.

330 # We should perhaps print it as an ordinary message,

331 # but we do not have a file here to which to print.

332 # import warnings

333 # warnings.warn('Calculating KS-gap directly from the k-points, '

334 # 'can be inaccurate.')

335 homolumo = self.wfs.get_homo_lumo(_gllb=True)

336 homo, lumo = homolumo

337 dxc_pot = self.calculate_discontinuity_potential(homo * Ha, lumo * Ha)

338 self.Dxc_vt_sG = dxc_pot.vt_sG

339 self.Dxc_D_asp = dxc_pot.D_asp

340 self.Dxc_Dresp_asp = dxc_pot.Dresp_asp

341

342 def calculate_discontinuity_potential(self, homo, lumo):

343 r"""Calculate the derivative discontinuity potential.

344

345 This function calculates $`\Delta_{xc}(r)`$ as given by

346 Eq. (24) of https://doi.org/10.1103/PhysRevB.82.115106

347

348 Parameters:

349

350 homo:

351 homo energy (or energies in spin-polarized case) in eV

352 lumo:

353 lumo energy (or energies in spin-polarized case) in eV

354

355 Returns:

356

357 dxc_pot: Discontinuity potential

358 """

359 if self.nspins == 2:

360 eref_s = np.asarray(homo)

361 eref_lumo_s = np.asarray(lumo)

362 else:

363 eref_s = np.asarray([homo])

364 eref_lumo_s = np.asarray([lumo])

365 assert len(eref_s) == len(eref_lumo_s) == self.nspins

366 eref_s /= Ha

367 eref_lumo_s /= Ha

368

369 dxc_Dresp_asp = self.empty_atomic_matrix()

370 dxc_D_asp = self.empty_atomic_matrix()

371 for a in self.density.D_asp:

372 ni = self.wfs.setups[a].ni

373 dxc_Dresp_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2))

374 dxc_D_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2))

375

376 # Calculate new response potential with LUMO reference

377 w_kn = self.coefficients.get_coefficients_for_lumo_perturbation(

378 self.wfs.kpt_u, eref_s=eref_s, eref_lumo_s=eref_lumo_s)

379 f_kn = [kpt.f_n for kpt in self.wfs.kpt_u]

380

381 dxc_vt_sG = self.gd.zeros(self.nspins)

382 nt_sG = self.gd.zeros(self.nspins)

383 for kpt, w_n in zip(self.wfs.kpt_u, w_kn):

384 self.wfs.add_to_density_from_k_point_with_occupation(dxc_vt_sG,

385 kpt, w_n)

386 self.wfs.add_to_density_from_k_point(nt_sG, kpt)

387

388 self.wfs.kptband_comm.sum(nt_sG)

389 self.wfs.kptband_comm.sum(dxc_vt_sG)

390

391 if self.wfs.kd.symmetry:

392 for nt_G, dxc_vt_G in zip(nt_sG, dxc_vt_sG):

393 self.wfs.kd.symmetry.symmetrize(nt_G, self.gd)

394 self.wfs.kd.symmetry.symmetrize(dxc_vt_G, self.gd)

395

396 dxc_vt_sG /= nt_sG + self.damp

397

398 self.wfs.calculate_atomic_density_matrices_with_occupation(

399 dxc_Dresp_asp, w_kn)

400 self.wfs.calculate_atomic_density_matrices_with_occupation(

401 dxc_D_asp, f_kn)

402 dxc_pot = ResponsePotential(self, dxc_vt_sG, None,

403 dxc_D_asp, dxc_Dresp_asp)

404 return dxc_pot

405

406 def calculate_delta_xc_perturbation(self):

407 warnings.warn(

408 'The function calculate_delta_xc_perturbation() is deprecated. '

409 'Please use calculate_discontinuity() instead. '

410 'See documentation on calculating band gap with GLLBSC.',

411 DeprecationWarning)

412 dxc_pot = ResponsePotential(self, self.Dxc_vt_sG, None, self.Dxc_D_asp,

413 self.Dxc_Dresp_asp)

414 if self.nspins != 1:

415 ret = []

416 for spin in range(self.nspins):

417 ret.append(self.calculate_discontinuity(dxc_pot, spin=spin))

418 return ret

419 return self.calculate_discontinuity(dxc_pot)

420

421 def calculate_delta_xc_perturbation_spin(self, s=0):

422 warnings.warn(

423 'The function calculate_delta_xc_perturbation_spin() is '

424 'deprecated. Please use calculate_discontinuity_spin() instead. '

425 'See documentation on calculating band gap with GLLBSC.',

426 DeprecationWarning)

427 dxc_pot = ResponsePotential(self, self.Dxc_vt_sG, None, self.Dxc_D_asp,

428 self.Dxc_Dresp_asp)

429 return self.calculate_discontinuity(dxc_pot, spin=s)

430

431 def calculate_discontinuity(self,

432 dxc_pot: ResponsePotential,

433 spin: int = None):

434 r"""Calculate the derivative discontinuity.

435

436 This function evaluates the perturbation theory expression

437 for the derivative discontinuity value as given by

438 Eq. (25) of https://doi.org/10.1103/PhysRevB.82.115106

439

440 Parameters:

441

442 dxc_pot:

443 Discontinuity potential.

444 spin:

445 Spin value.

446

447 Returns:

448

449 KSgap:

450 Kohn-Sham gap in eV

451 dxc:

452 Derivative discontinuity in eV

453 """

454 if spin is None:

455 if self.nspins != 1:

456 raise ValueError('Specify spin for the discontinuity.')

457 spin = 0

458

459 # Redistribute discontinuity potential

460 if dxc_pot.response is not self:

461 dxc_pot.redistribute(self)

462

463 homo, lumo = self.wfs.get_homo_lumo(spin, _gllb=True)

464 KSgap = lumo - homo

465

466 # Find the lumo-orbital of this spin

467 if self.wfs.bd.comm.size != 1:

468 raise BadParallelization(

469 'Band parallelization is not supported by '

470 'calculate_discontinuity().')

471 eps_n = self.wfs.kpt_qs[0][spin].eps_n

472 lumo_n = (eps_n < self.wfs.fermi_level).sum()

473

474 # Calculate the expectation value for all k points, and keep

475 # the smallest energy value

476 nt_G = self.gd.empty()

477 min_energy = np.inf

478 for u, kpt in enumerate(self.wfs.kpt_u):

479 if kpt.s == spin:

480 nt_G[:] = 0.0

481 self.wfs.add_orbital_density(nt_G, kpt, lumo_n)

482 E = 0.0

483 for a in dxc_pot.D_asp:

484 D_sp = dxc_pot.D_asp[a]

485 Dresp_sp = dxc_pot.Dresp_asp[a]

486 P_ni = kpt.P_ani[a]

487 Dwf_p = pack_density(

488 np.outer(P_ni[lumo_n].T.conj(), P_ni[lumo_n]).real)

489 E += self.integrate_sphere(a, Dresp_sp, D_sp, Dwf_p,

490 spin=spin)

491 E = self.gd.comm.sum_scalar(E)

492 E += self.gd.integrate(nt_G * dxc_pot.vt_sG[spin])

493 E += kpt.eps_n[lumo_n] - lumo

494 min_energy = min(min_energy, E)

495

496 # Take the smallest value over all distributed k points

497 dxc = self.wfs.kd.comm.min_scalar(min_energy)

498 return KSgap * Ha, dxc * Ha

499

500 def calculate_discontinuity_from_average(self,

501 dxc_pot: ResponsePotential,

502 spin: int,

503 testing: bool = False):

504 msg = ('This function estimates discontinuity by calculating '

505 'the average of the discontinuity potential. '

506 'Use only for testing and debugging.')

507 if not testing:

508 raise RuntimeError(msg)

509 else:

510 warnings.warn(msg)

511 assert self.wfs.world.size == 1

512

513 # Calculate average of lumo reference response potential

514 dxc = np.average(dxc_pot.vt_sG[spin])

515 return dxc * Ha

516

517 def initialize_from_atomic_orbitals(self, basis_functions):

518 # Initialize 'response-density' and density-matrices

519 # print('Initializing from atomic orbitals')

520 self.Dresp_asp = self.empty_atomic_matrix()

521 setups = self.wfs.setups

522

523 for a in self.density.D_asp.keys():

524 ni = setups[a].ni

525 self.Dresp_asp[a] = np.zeros((self.nspins, ni * (ni + 1) // 2))

526

527 self.D_asp = self.empty_atomic_matrix()

528 f_asi = {}

529 w_asi = {}

530

531 for a in basis_functions.atom_indices:

532 if setups[a].type == 'ghost':

533 w_j = [0.0]

534 else:

535 w_j = setups[a].extra_xc_data['w_j']

536

537 # Basis function coefficients based of response weights

538 w_si = setups[a].calculate_initial_occupation_numbers(

539 0, False,

540 charge=0,

541 nspins=self.nspins,

542 f_j=w_j)

543 # Basis function coefficients based on density

544 f_si = setups[a].calculate_initial_occupation_numbers(

545 0, False,

546 charge=0,

547 nspins=self.nspins)

548

549 if a in basis_functions.my_atom_indices:

550 self.Dresp_asp[a] = setups[a].initialize_density_matrix(w_si)

551 self.D_asp[a] = setups[a].initialize_density_matrix(f_si)

552

553 f_asi[a] = f_si

554 w_asi[a] = w_si

555

556 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp)

557 self.Ddist_asp = self.distribute_D_asp(self.D_asp)

558 nt_sG = self.gd.zeros(self.nspins)

559 basis_functions.add_to_density(nt_sG, f_asi)

560 self.vt_sG = self.gd.zeros(self.nspins)

561 basis_functions.add_to_density(self.vt_sG, w_asi)

562 # Update vt_sG to correspond atomic response potential. This will be

563 # used until occupations and eigenvalues are available.

564 self.vt_sG /= nt_sG + self.damp

565 self.vt_sg = self.finegd.zeros(self.nspins)

566 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg)

567

568 def get_extra_setup_data(self, extra_data):

569 ae = self.ae

570 njcore = ae.njcore

571 w_ln = self.coefficients.get_coefficients_1d(smooth=True)

572 extra_data['w_ln'] = w_ln

573

574 w_j = self.coefficients.get_coefficients_1d()

575 x_g = np.dot(w_j[:njcore], safe_sqr(ae.u_j[:njcore]))

576 x_g[1:] /= ae.r[1:] ** 2 * 4 * np.pi

577 x_g[0] = x_g[1]

578 extra_data['core_response'] = x_g

579

580 def write(self, writer):

581 """Writes response specific data."""

582 wfs = self.wfs

583 kpt_comm = wfs.kd.comm

584 gd = wfs.gd

585

586 # Write the pseudodensity on the coarse grid:

587 shape = (wfs.nspins,) + tuple(gd.get_size_of_global_array())

588 writer.add_array('gllb_pseudo_response_potential', shape)

589 if kpt_comm.rank == 0:

590 for vt_G in self.vt_sG:

591 writer.fill(gd.collect(vt_G) * Ha)

592

593 writer.write('gllb_atomic_density_matrices',

594 pack_atomic_matrices(self.D_asp))

595 writer.write('gllb_atomic_response_matrices',

596 pack_atomic_matrices(self.Dresp_asp))

597

598 writer.write(eref_s=self.eref_s)

599 writer.write(eref_source_s=self.eref_source_s)

600

601 def empty_atomic_matrix(self):

602 assert self.wfs.atom_partition is self.density.atom_partition

603 return self.wfs.setups.empty_atomic_matrix(self.wfs.nspins,

604 self.wfs.atom_partition)

605

606 def read(self, reader):

607 r = reader.hamiltonian.xc

608 wfs = self.wfs

609

610 d('Reading vt_sG')

611 self.gd.distribute(r.gllb_pseudo_response_potential / reader.ha,

612 self.vt_sG)

613 self.density.distribute_and_interpolate(self.vt_sG, self.vt_sg)

614

615 def unpack_density(D_sP):

616 return unpack_atomic_matrices(D_sP, wfs.setups)

617

618 # Read atomic density matrices and non-local part of hamiltonian:

619 D_asp = unpack_density(r.gllb_atomic_density_matrices)

620 Dresp_asp = unpack_density(r.gllb_atomic_response_matrices)

621

622 # All density matrices are loaded to all cores

623 # First distribute them to match density.D_asp

624 self.D_asp = self.empty_atomic_matrix()

625 self.Dresp_asp = self.empty_atomic_matrix()

626 for a in self.density.D_asp:

627 self.D_asp[a][:] = D_asp[a]

628 self.Dresp_asp[a][:] = Dresp_asp[a]

629 assert len(self.Dresp_asp) == len(self.density.D_asp)

630

631 # Further distributes the density matricies for xc-corrections

632 self.Ddist_asp = self.distribute_D_asp(self.D_asp)

633 self.Drespdist_asp = self.distribute_D_asp(self.Dresp_asp)

634

635 if 'eref_s' in r:

636 self.eref_s = r.eref_s

637 self.eref_source_s = r.eref_source_s

638

639 def heeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeelp(self, olddens):

640 # XXX This function should be removed once the deprecated

641 # `fixdensity=True` option is removed.

642 from gpaw.density import redistribute_array

643 self.vt_sg = redistribute_array(self.vt_sg,

644 olddens.finegd, self.finegd,

645 self.wfs.nspins, self.wfs.kptband_comm)

646 self.vt_sG = redistribute_array(self.vt_sG,

647 olddens.gd, self.gd,

648 self.wfs.nspins, self.wfs.kptband_comm)

Coverage for gpaw/xc/gllb/c_response.py: 96%

365 statements