Coverage for gpaw/directmin/functional/fdpw/pz.py: 88%

1"""

2Potentials for orbital density dependent energy functionals

3"""

5import numpy as np

7from gpaw.directmin.tools import get_n_occ, d_matrix

8from gpaw.lfc import LFC

9from gpaw.poisson import PoissonSolver

10from gpaw.transformers import Transformer

11from gpaw.utilities import pack_density, unpack_hermitian

12from gpaw.utilities.ewald import madelung

13from gpaw.utilities.partition import AtomPartition

14import gpaw.cgpaw as cgpaw

17class PZSICFDPW:

19 """

20 Perdew-Zunger self-interaction corrections

22 """

23 def __init__(self, wfs, dens, ham, scaling_factor=(1.0, 1.0),

24 sic_coarse_grid=True, store_potentials=False,

25 poisson_solver='FPS'):

27 self.name = 'PZ-SIC'

28 # what we need from wfs

29 self.setups = wfs.setups

30 spos_ac = wfs.spos_ac

31 self.cgd = wfs.gd

33 # what we need from dens

34 self.finegd = dens.finegd

35 self.sic_coarse_grid = sic_coarse_grid

36 self.pd2 = None

37 self.pd3 = None

38 self.corr = None

39 if wfs.mode == 'pw':

40 from gpaw.wavefunctions.pw import PWLFC

42 assert self.sic_coarse_grid

43 self.pd2 = dens.pd2

44 self.pd3 = dens.pd3

45 self.ghat = PWLFC(

46 [setup.ghat_l for setup in self.setups], self.pd2)

47 rank_a = wfs.gd.get_ranks_from_positions(spos_ac)

48 atom_partition = AtomPartition(wfs.gd.comm, rank_a, name='gd')

49 self.ghat.set_positions(spos_ac, atom_partition)

50 self.corr = madelung(wfs.gd.cell_cv)

51 self.corr_q = 1.0

52 for nc in self.pd2.gd.N_c:

53 self.corr_q *= nc

54 self.corr_q *= self.corr

56 self.G2 = self.pd2.G2_qG[0].copy()

57 if self.pd2.gd.comm.rank == 0:

58 self.G2[0] = 1.0

59 else:

60 if self.sic_coarse_grid:

61 self.ghat = LFC(

62 self.cgd, [setup.ghat_l for setup in self.setups],

63 integral=np.sqrt(4 * np.pi), forces=True)

64 self.ghat.set_positions(spos_ac)

65 else:

66 self.ghat = dens.ghat # we usually solve poiss. on finegd

68 # what we need from ham

69 self.xc = ham.xc

71 if poisson_solver == 'FPS':

72 self.poiss = PoissonSolver(use_charge_center=True,

73 use_charged_periodic_corrections=True)

74 elif poisson_solver == 'GS':

75 self.poiss = PoissonSolver(

76 name='fd', relax=poisson_solver, eps=1.0e-16,

77 use_charge_center=True, use_charged_periodic_corrections=True)

79 if self.sic_coarse_grid is True:

80 self.poiss.set_grid_descriptor(self.cgd)

81 else:

82 self.poiss.set_grid_descriptor(self.finegd)

84 self.interpolator = Transformer(self.cgd, self.finegd, 3)

85 self.restrictor = Transformer(self.finegd, self.cgd, 3)

86 self.dtype = wfs.dtype

87 self.eigv_s = {}

88 self.lagr_diag_s = {}

89 self.e_sic_by_orbitals = {}

90 self.counter = 0 # number of calls of this class

91 # Scaling factor:

92 self.beta_c = scaling_factor[0]

93 self.beta_x = scaling_factor[1]

95 self.n_kps = wfs.kd.nibzkpts

96 self.store_potentials = store_potentials

97 self.grad = {}

98 self.total_sic = 0.0

100 if store_potentials:

101 self.old_pot = {}

102 for kpt in wfs.kpt_u:

103 k = self.n_kps * kpt.s + kpt.q

104 n_occ = get_n_occ(kpt)[0]

105 self.old_pot[k] = self.cgd.zeros(n_occ, dtype=float)

106

107 def get_orbdens_compcharge_dm_kpt(self, wfs, kpt, n):

108

109 if wfs.mode == 'pw':

110 nt_G = wfs.pd.gd.zeros(global_array=True)

111 psit_G = wfs.pd.alltoall1(kpt.psit.array[n: n + 1], kpt.q)

112 if psit_G is not None:

113 psit_R = wfs.pd.ifft(psit_G, kpt.q, local=True, safe=False)

114 cgpaw.add_to_density(1.0, psit_R, nt_G)

115 wfs.pd.gd.comm.sum(nt_G)

116 nt_G = wfs.pd.gd.distribute(nt_G)

117 else:

118 nt_G = np.absolute(kpt.psit_nG[n]**2)

119

120 # paw

121 Q_aL = {}

122 D_ap = {}

123 for a, P_ni in kpt.P_ani.items():

124 P_i = P_ni[n]

125 D_ii = np.outer(P_i, P_i.conj()).real

126 D_ap[a] = D_p = pack_density(D_ii)

127 Q_aL[a] = np.dot(D_p, self.setups[a].Delta_pL)

128

129 return nt_G, Q_aL, D_ap

130

131 def get_energy_and_gradients_kpt(

132 self, wfs, kpt, grad_knG=None, U_k=None, add_grad=False,

133 ham=None, scalewithocc=True, exstate=False):

134

135 k = self.n_kps * kpt.s + kpt.q

136 n_occ = get_n_occ(kpt)[0]

137 self.grad[k] = np.zeros_like(kpt.psit_nG) if exstate \

138 else np.zeros_like(kpt.psit_nG[:n_occ])

139

140 if exstate:

141 self.get_gradient_ks_kpt(wfs, kpt, ham=ham)

142 esic = self.get_esic_add_sic_gradient_kpt(

143 wfs, kpt, grad_knG, U_k, add_grad, scalewithocc,

144 exstate=exstate)

145

146 return esic

147

148 def get_gradient_ks_kpt(self, wfs, kpt, ham=None):

149

150 k = self.n_kps * kpt.s + kpt.q

151 wfs.timer.start('KS e/g grid calculations')

152 wfs.apply_pseudo_hamiltonian(kpt, ham, kpt.psit_nG, self.grad[k])

153

154 c_axi = {}

155 for a, P_xi in kpt.P_ani.items():

156 dH_ii = unpack_hermitian(ham.dH_asp[a][kpt.s])

157 c_xi = np.dot(P_xi, dH_ii)

158 c_axi[a] = c_xi

159

160 # not sure about this:

161 ham.xc.add_correction(kpt, kpt.psit_nG, self.grad[k],

162 kpt.P_ani, c_axi, n_x=None,

163 calculate_change=False)

164 # add projectors to the H|psi_i>

165 wfs.pt.add(self.grad[k], c_axi, kpt.q)

166 # scale with occupation numbers

167 for i, f in enumerate(kpt.f_n):

168 self.grad[k][i] *= f

169

170 wfs.timer.stop('KS e/g grid calculations')

171

172 return 0.0

173

174 def get_esic_add_sic_gradient_kpt(

175 self, wfs, kpt, grad_knG=None, U_k=None, add_grad=False,

176 scalewithocc=True, exstate=False):

177

178 wfs.timer.start('SIC e/g grid calculations')

179 k = self.n_kps * kpt.s + kpt.q

180 n_occ = get_n_occ(kpt)[0]

181 e_total_sic = np.array([])

182

183 for i in range(n_occ):

184 if wfs.mode == 'pw':

185 e_sic, vt_G, dH_ap = self.get_si_pot_dh_pw(

186 wfs, kpt, i, exstate=exstate)

187 else:

188 nt_G, Q_aL, D_ap = \

189 self.get_orbdens_compcharge_dm_kpt(wfs, kpt, i)

190 e_sic, vt_G, dH_ap = \

191 self.get_pz_sic_ith_kpt(

192 nt_G, Q_aL, D_ap, i, k, wfs.timer)

193

194 e_total_sic = np.append(e_total_sic,

195 kpt.f_n[i] * e_sic, axis=0)

196 if wfs.mode == 'pw':

197 vt_G = wfs.pd.gd.collect(vt_G, broadcast=True)

198 Q_G = wfs.pd.Q_qG[kpt.q]

199 psit_G = wfs.pd.alltoall1(kpt.psit_nG[i: i + 1], kpt.q)

200 if psit_G is not None:

201 psit_R = wfs.pd.ifft(

202 psit_G, kpt.q, local=True, safe=False)

203 psit_R *= vt_G

204 wfs.pd.fftplan.execute()

205 vtpsit_G = wfs.pd.tmp_Q.ravel()[Q_G]

206 else:

207 vtpsit_G = wfs.pd.tmp_G

208 tmp = np.zeros_like(self.grad[k][i: i + 1])

209 wfs.pd.alltoall2(vtpsit_G, kpt.q, tmp)

210 self.grad[k][i] += tmp[0]

211 if scalewithocc:

212 self.grad[k][i] *= kpt.f_n[i]

213 else:

214 if scalewithocc:

215 self.grad[k][i] += kpt.psit_nG[i] * vt_G * kpt.f_n[i]

216 else:

217 self.grad[k][i] += kpt.psit_nG[i] * vt_G

218 c_axi = {}

219 for a in kpt.P_ani.keys():

220 dH_ii = unpack_hermitian(dH_ap[a])

221 c_xi = np.dot(kpt.P_ani[a][i], dH_ii)

222 c_axi[a] = c_xi * kpt.f_n[i]

223 # add projectors to

224 wfs.pt.add(self.grad[k][i], c_axi, kpt.q)

225

226 if exstate:

227 if U_k is not None:

228 grad_knG[k][:] += np.tensordot(

229 U_k[k].conj(), self.grad[k], axes=1)

230 if add_grad:

231 grad_knG[k][:] += self.grad[k]

232 else:

233 if add_grad:

234 if U_k is not None:

235 grad_knG[k][:n_occ] += np.tensordot(

236 U_k[k][:n_occ, :n_occ].conj(), self.grad[k][:n_occ],

237 axes=1)

238 else:

239 grad_knG[k][:n_occ] += self.grad[k][:n_occ]

240 else:

241 if U_k is not None:

242 self.grad[k][:] = np.tensordot(

243 U_k[k][:n_occ, :n_occ].conj(), self.grad[k][:n_occ],

244 axes=1)

245

246 self.e_sic_by_orbitals[k] = \

247 e_total_sic.reshape(e_total_sic.shape[0] // 2, 2)

248

249 wfs.timer.stop('SIC e/g grid calculations')

250

251 return e_total_sic.sum()

252

253 def get_pseudo_pot(self, nt, Q_aL, i, kpoint=None):

254

255 if self.sic_coarse_grid is False:

256 # fine grid

257 vt_sg = self.finegd.zeros(2)

258 v_ht_g = self.finegd.zeros()

259 nt_sg = self.finegd.zeros(2)

260 else:

261 # coarse grid

262 vt_sg = self.cgd.zeros(2)

263 v_ht_g = self.cgd.zeros()

264 nt_sg = self.cgd.zeros(2)

265

266 if self.sic_coarse_grid is False:

267 self.interpolator.apply(nt, nt_sg[0])

268 nt_sg[0] *= self.cgd.integrate(nt) / \

269 self.finegd.integrate(nt_sg[0])

270 e_xc = self.xc.calculate(self.finegd, nt_sg, vt_sg)

271 else:

272 nt_sg[0] = nt

273 e_xc = self.xc.calculate(self.cgd, nt_sg, vt_sg)

274

275 vt_sg[0] *= -self.beta_x

276

277 self.ghat.add(nt_sg[0], Q_aL)

278

279 if self.store_potentials:

280 if self.sic_coarse_grid:

281 v_ht_g = self.old_pot[kpoint][i]

282 else:

283 self.interpolator.apply(self.old_pot[kpoint][i],

284 v_ht_g)

285

286 self.poiss.solve(v_ht_g, nt_sg[0],

287 zero_initial_phi=False)

288

289 if self.store_potentials:

290 if self.sic_coarse_grid is True:

291 self.old_pot[kpoint][i] = v_ht_g.copy()

292 else:

293 self.restrictor.apply(v_ht_g, self.old_pot[kpoint][i])

294

295 if self.sic_coarse_grid is False:

296 ec = 0.5 * self.finegd.integrate(nt_sg[0] * v_ht_g)

297 else:

298 ec = 0.5 * self.cgd.integrate(nt_sg[0] * v_ht_g)

299

300 vt_sg[0] -= v_ht_g * self.beta_c

301

302 if self.sic_coarse_grid is False:

303 vt_G = self.cgd.zeros()

304 self.restrictor.apply(vt_sg[0], vt_G)

305 else:

306 vt_G = vt_sg[0]

307

308 return np.array([-ec * self.beta_c, -e_xc * self.beta_x]), \

309 vt_G, v_ht_g

310

311 def get_paw_corrections(self, D_ap, vHt_g):

312

313 # XC-PAW

314 dH_ap = {}

315

316 exc = 0.0

317 for a, D_p in D_ap.items():

318 setup = self.setups[a]

319 dH_sp = np.zeros((2, len(D_p)))

320 D_sp = np.array([D_p, np.zeros_like(D_p)])

321 exc += self.xc.calculate_paw_correction(

322 setup, D_sp, dH_sp, addcoredensity=False)

323 dH_ap[a] = -dH_sp[0] * self.beta_x

324

325 # Hartree-PAW

326 ec = 0.0

327 W_aL = self.ghat.dict()

328 self.ghat.integrate(vHt_g, W_aL)

329

330 for a, D_p in D_ap.items():

331 setup = self.setups[a]

332 M_p = np.dot(setup.M_pp, D_p)

333 ec += np.dot(D_p, M_p)

334

335 dH_ap[a] += -(2.0 * M_p + np.dot(setup.Delta_pL,

336 W_aL[a])) * self.beta_c

337

338 if self.sic_coarse_grid is False:

339 ec = self.finegd.comm.sum_scalar(ec)

340 exc = self.finegd.comm.sum_scalar(exc)

341 else:

342 ec = self.cgd.comm.sum_scalar(ec)

343 exc = self.cgd.comm.sum_scalar(exc)

344

345 return np.array([-ec * self.beta_c, -exc * self.beta_x]), dH_ap

346

347 def get_energy_and_gradients_inner_loop(

348 self, wfs, kpt, a_mat, evals, evec, ham=None,

349 exstate=False):

350

351 if exstate:

352 ndim = wfs.bd.nbands

353 else:

354 ndim = 0

355 for f in kpt.f_n:

356 if f > 1.0e-10:

357 ndim += 1

358

359 k = self.n_kps * kpt.s + kpt.q

360 self.grad[k] = np.zeros_like(kpt.psit_nG[:ndim])

361 e_sic = self.get_energy_and_gradients_kpt(

362 wfs, kpt, ham=ham, exstate=exstate)

363 wfs.timer.start('Unitary gradients')

364 l_odd = wfs.integrate(kpt.psit_nG[:ndim], self.grad[k][:ndim], True)

365 f = np.ones(ndim)

366 indz = np.absolute(l_odd) > 1.0e-4

367 l_c = 2.0 * l_odd[indz]

368 l_odd = f[:, np.newaxis] * l_odd.T.conj() - f * l_odd

369 kappa = np.max(np.absolute(l_odd[indz]) / np.absolute(l_c))

370

371 if a_mat is None:

372 wfs.timer.stop('Unitary gradients')

373 return 2.0 * l_odd.T.conj(), e_sic, kappa

374 else:

375 g_mat = evec.T.conj() @ l_odd.T.conj() @ evec

376 g_mat = g_mat * d_matrix(evals)

377 g_mat = evec @ g_mat @ evec.T.conj()

378 for i in range(g_mat.shape[0]):

379 g_mat[i][i] *= 0.5

380 wfs.timer.stop('Unitary gradients')

381 if a_mat.dtype == float:

382 g_mat = g_mat.real

383 return 2.0 * g_mat, e_sic, kappa

384

385 def get_odd_corrections_to_forces(self, F_av, wfs, kpt, exstate=False):

386

387 n_occ = get_n_occ(kpt)[0]

388 n_kps = self.n_kps

389

390 dP_amiv = wfs.pt.dict(n_occ, derivative=True)

391 wfs.pt.derivative(kpt.psit_nG[:n_occ], dP_amiv)

392 k = n_kps * kpt.s + kpt.q

393 for m in range(n_occ):

394 # calculate Hartree pot, compans. charge and PAW corrects

395 if wfs.mode == 'fd':

396 nt_G, Q_aL, D_ap = \

397 self.get_orbdens_compcharge_dm_kpt(wfs, kpt, m)

398 e_sic, vt_G, v_ht_g = \

399 self.get_pseudo_pot(nt_G, Q_aL, m, kpoint=k)

400 e_sic_paw_m, dH_ap = \

401 self.get_paw_corrections(D_ap, v_ht_g)

402 elif wfs.mode == 'pw':

403 e_sic, vt_G, dH_ap, Q_aL, v_ht_g = \

404 self.get_si_pot_dh_pw(

405 wfs, kpt, m, returnQalandVhq=True, exstate=exstate)

406 else:

407 raise NotImplementedError

408

409 # Force from compensation charges:

410 dF_aLv = self.ghat.dict(derivative=True)

411 self.ghat.derivative(v_ht_g, dF_aLv)

412 for a, dF_Lv in dF_aLv.items():

413 F_av[a] -= kpt.f_n[m] * self.beta_c * \

414 np.dot(Q_aL[a], dF_Lv)

415

416 # Force from projectors

417 for a, dP_miv in dP_amiv.items():

418 dP_vi = dP_miv[m].T.conj()

419 dH_ii = unpack_hermitian(dH_ap[a])

420 P_i = kpt.P_ani[a][m]

421 F_v = np.dot(np.dot(dP_vi, dH_ii), P_i)

422 F_av[a] += kpt.f_n[m] * 2.0 * F_v.real

423

424 def get_pz_sic_ith_kpt(self, nt_G, Q_aL, D_ap, i, k, timer):

425

426 """

427 :param nt_G: one-electron orbital density

428 :param Q_aL: its compensation charge

429 :param D_ap: its density matrix

430 :param i: number of orbital

431 :param k: k-point number, k = n_kperspin * kpt.s + kpt.q

432 :param timer:

433 :return: E, v and dH

434 E = -(beta_c * E_Hartree[n_i] + beta_x * E_xc[n_i])

435 v = dE / dn_i

436 dH - paw corrections

437

438 """

439

440 # calculate sic energy,

441 # sic pseudo-potential and Hartree

442 timer.start('Get Pseudo Potential')

443 # calculate sic energy, sic pseudo-potential and Hartree

444 e_pz, vt_G, v_ht_g = \

445 self.get_pseudo_pot(nt_G, Q_aL, i, kpoint=k)

446 timer.stop('Get Pseudo Potential')

447

448 # calculate PAW corrections

449 timer.start('PAW')

450 # calculate PAW corrections

451 e_pz_paw_m, dH_ap = self.get_paw_corrections(D_ap, v_ht_g)

452 timer.stop('PAW')

453

454 # total sic:

455 e_pz += e_pz_paw_m

456

457 return e_pz, vt_G, dH_ap

458

459 def get_si_pot_dh_pw(

460 self, wfs, kpt, n, returnQalandVhq=False, exstate=False):

461

462 wfs.timer.start("IFFT: Get density on real grid")

463 nt_G = wfs.pd.gd.zeros(global_array=True)

464 psit_G = wfs.pd.alltoall1(kpt.psit.array[n:n + 1], kpt.q)

465 if psit_G is not None:

466 # real space wfs:

467 psit_R = wfs.pd.ifft(psit_G, kpt.q,

468 local=True, safe=False)

469 cgpaw.add_to_density(1.0, psit_R, nt_G)

470 wfs.pd.gd.comm.sum(nt_G)

471 nt_G = wfs.pd.gd.distribute(nt_G) # this is real space grid

472 wfs.timer.stop("IFFT: Get density on real grid")

473

474 # multipole moments and atomic dm

475 wfs.timer.start("Multipole Moments and Atomic Dens. Mat.")

476 Q_aL = {}

477 D_ap = {}

478 for a, P_ni in kpt.P_ani.items():

479 P_i = P_ni[n]

480 D_ii = np.outer(P_i, P_i.conj()).real

481 D_ap[a] = D_p = pack_density(D_ii)

482 Q_aL[a] = np.dot(D_p, self.setups[a].Delta_pL)

483 wfs.timer.stop("Multipole Moments and Atomic Dens. Mat.")

484

485 # xc

486 wfs.timer.start("Calc. XC on pseudo density")

487 nt_sg = wfs.gd.zeros(2)

488 nt_sg[0] = nt_G

489 vt_sg = wfs.gd.zeros(2)

490 exc = self.xc.calculate(wfs.gd, nt_sg, vt_sg)

491 vt_G = - vt_sg[0] * self.beta_x

492 dH_ap = {}

493 wfs.timer.stop("Calc. XC on pseudo density")

494

495 wfs.timer.start("Calc. PAW-XC")

496 excpaw = 0.0

497 for a, D_p in D_ap.items():

498 setup = self.setups[a]

499 dH_sp = np.zeros((2, len(D_p)))

500 D_sp = np.array([D_p, np.zeros_like(D_p)])

501 excpaw += self.xc.calculate_paw_correction(

502 setup, D_sp, dH_sp, addcoredensity=False)

503 dH_ap[a] = -dH_sp[0] * self.beta_x

504 excpaw = wfs.gd.comm.sum_scalar(excpaw)

505 exc += excpaw

506 wfs.timer.stop("Calc. PAW-XC")

507

508 wfs.timer.start("FFT density")

509 nt_Q = self.pd2.fft(nt_G)

510 self.ghat.add(nt_Q, Q_aL)

511 wfs.timer.stop("FFT density")

512

513 wfs.timer.start("Calc. Hartree on pseudo density")

514 if exstate:

515 nt_G = self.pd2.ifft(nt_Q)

516 vHt = np.zeros_like(nt_G)

517 self.poiss.solve(vHt, nt_G, zero_initial_phi=False)

518 ehart = 0.5 * self.cgd.integrate(nt_G, vHt)

519 vHt_q = self.pd2.fft(vHt)

520 else:

521 if self.pd2.gd.comm.rank == 0:

522 nt_Q[0] = 0.0

523 vHt_q = 4 * np.pi * nt_Q / self.G2

524 ehart = 0.5 * self.pd2.integrate(vHt_q, nt_Q)

525 # correct for uniform background

526 if self.pd2.gd.comm.rank == 0:

527 vHt_q[0] += self.corr_q

528 vHt = self.pd2.ifft(vHt_q)

529 ehart += self.corr / 2.0

530 wfs.timer.stop("Calc. Hartree on pseudo density")

531

532 # PAW to Hartree

533 ehartpaw = 0.0

534 W_aL = self.ghat.dict()

535 self.ghat.integrate(vHt_q, W_aL)

536

537 wfs.timer.start("Calc. PAW-Hartree")

538 for a, D_p in D_ap.items():

539 setup = self.setups[a]

540 M_p = np.dot(setup.M_pp, D_p)

541 ehartpaw += np.dot(D_p, M_p)

542 dH_ap[a] += -(2.0 * M_p + np.dot(setup.Delta_pL,

543 W_aL[a])) * self.beta_c

544 ehartpaw = wfs.gd.comm.sum_scalar(ehartpaw)

545 ehart += ehartpaw

546 wfs.timer.stop("Calc. PAW-Hartree")

547

548 vt_G -= vHt * self.beta_c

549

550 if returnQalandVhq:

551 return np.array([-ehart * self.beta_c, -exc * self.beta_x]), \

552 vt_G, dH_ap, Q_aL, vHt_q

553 else:

554 return np.array([-ehart * self.beta_c, -exc * self.beta_x]), \

555 vt_G, dH_ap