Coverage for gpaw/directmin/functional/lcao/pz.py: 87%

1"""

2Potentials for orbital density dependent energy functionals

3"""

4import numpy as np

6from gpaw.directmin.tools import d_matrix

7from gpaw.lfc import LFC

8from gpaw.poisson import PoissonSolver

9from gpaw.transformers import Transformer

10from gpaw.utilities import pack_density, unpack_hermitian

13class PZSICLCAO:

14 """

15 Perdew-Zunger self-interaction corrections

17 """

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

19 sic_coarse_grid=True, store_potentials=False,

20 poisson_solver='FPS'):

22 self.name = 'PZ-SIC'

23 # what we need from wfs

24 self.setups = wfs.setups

25 self.bfs = wfs.basis_functions

26 spos_ac = wfs.spos_ac

27 self.cgd = wfs.gd

29 # what we need from dens

30 self.finegd = dens.finegd

31 self.sic_coarse_grid = sic_coarse_grid

33 if self.sic_coarse_grid:

34 self.ghat_cg = LFC(self.cgd,

35 [setup.ghat_l for setup

36 in self.setups],

37 integral=np.sqrt(4 * np.pi),

38 forces=True)

39 self.ghat_cg.set_positions(spos_ac)

40 self.ghat = None

41 else:

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

43 self.ghat_cg = None

45 # what we need from ham

46 self.xc = ham.xc

48 if poisson_solver == 'FPS':

49 self.poiss = PoissonSolver(use_charge_center=True,

50 use_charged_periodic_corrections=True)

51 elif poisson_solver == 'GS':

52 self.poiss = PoissonSolver(name='fd',

53 relax=poisson_solver,

54 eps=1.0e-16,

55 use_charge_center=True,

56 use_charged_periodic_corrections=True)

58 if self.sic_coarse_grid is True:

59 self.poiss.set_grid_descriptor(self.cgd)

60 else:

61 self.poiss.set_grid_descriptor(self.finegd)

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

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

65 self.dtype = wfs.dtype

66 self.eigv_s = {}

67 self.lagr_diag_s = {}

68 self.e_sic_by_orbitals = {}

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

70 # Scaling factor:

71 self.beta_c = scaling_factor[0]

72 self.beta_x = scaling_factor[1]

74 self.n_kps = wfs.kd.nibzkpts

75 self.store_potentials = store_potentials

76 if store_potentials:

77 self.old_pot = {}

78 for kpt in wfs.kpt_u:

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

80 n_occ = 0

81 nbands = len(kpt.f_n)

82 while n_occ < nbands and kpt.f_n[n_occ] > 1e-10:

83 n_occ += 1

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

86 def get_gradients(self, h_mm, c_nm, f_n, evec, evals, kpt,

87 wfs, timer, matrix_exp, repr_name,

88 ind_up, constraints, occupied_only=False):

89 """

90 :param C_nM: coefficients of orbitals

91 :return: matrix G - gradients, and orbital SI energies

93 which is G_{ij} = (1-delta_{ij}/2)*(int_0^1 e^{tA} L e^{-tA} dt )_{ji}

95 Lambda_ij = (C_i, F_j C_j )

97 L_{ij} = Lambda_ji^{cc} - Lambda_ij

99 """

100

101 # 0.

102 n_occ = 0

103 nbands = len(f_n)

104 while n_occ < nbands and f_n[n_occ] > 1e-10:

105 n_occ += 1

106

107 if occupied_only is True:

108 nbs = n_occ

109 else:

110 nbs = c_nm.shape[0]

111 n_set = c_nm.shape[1]

112

113 timer.start('Construct Gradient Matrix')

114 hc_mn = np.dot(h_mm.conj(), c_nm[:nbs].T)

115 h_mm = np.dot(c_nm[:nbs].conj(), hc_mn)

116 # odd part

117 b_mn = np.zeros(shape=(n_set, nbs), dtype=self.dtype)

118 e_total_sic = np.array([])

119 for n in range(n_occ):

120 F_MM, sic_energy_n =\

121 self.get_orbital_potential_matrix(f_n, c_nm, kpt,

122 wfs, wfs.setups,

123 n, timer)

124

125 b_mn[:, n] = np.dot(c_nm[n], F_MM.conj()).T

126 e_total_sic = np.append(e_total_sic, sic_energy_n, axis=0)

127 l_odd = np.dot(c_nm[:nbs].conj(), b_mn)

128

129 f = f_n[:nbs]

130 grad = f * (h_mm[:nbs, :nbs] + l_odd) - \

131 f[:, np.newaxis] * (h_mm[:nbs, :nbs] + l_odd.T.conj())

132

133 if matrix_exp in ['pade-approx', 'egdecomp2']:

134 grad = np.ascontiguousarray(grad)

135 elif matrix_exp == 'egdecomp':

136 timer.start('Use Eigendecomposition')

137 grad = np.dot(evec.T.conj(), np.dot(grad, evec))

138 grad = grad * d_matrix(evals)

139 grad = np.dot(evec, np.dot(grad, evec.T.conj()))

140 for i in range(grad.shape[0]):

141 grad[i][i] *= 0.5

142 timer.stop('Use Eigendecomposition')

143 else:

144 raise NotImplementedError('Check the keyword '

145 'for matrix_exp. \n'

146 'Must be '

147 '\'pade-approx\' or '

148 '\'egdecomp\'')

149 if self.dtype == float:

150 grad = grad.real

151 if repr_name in ['sparse', 'u-invar']:

152 grad = grad[ind_up]

153

154 timer.stop('Construct Gradient Matrix')

155

156 u = kpt.s * self.n_kps + kpt.q

157 self.e_sic_by_orbitals[u] = \

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

159

160 timer.start('Residual')

161 hc_mn += b_mn

162 h_mm += l_odd

163

164 rhs2 = kpt.S_MM.conj() @ c_nm[:n_occ].T @ h_mm[:n_occ, :n_occ]

165 hc_mn = hc_mn[:, :n_occ] - rhs2[:, :n_occ]

166

167 if constraints:

168 # Zero out the components of the residual that are constrained,

169 # so that the constrained degrees of freedom are always considered

170 # converged

171 for con in constraints:

172 con1 = con[0]

173 hc_mn[:, con1] = 0.0

174

175 norm = []

176 for i in range(n_occ):

177 norm.append(np.dot(hc_mn[:, i].conj(),

178 hc_mn[:, i]).real * kpt.f_n[i])

179

180 error = sum(norm)

181 del rhs2, hc_mn, norm

182 timer.stop('Residual')

183

184 if self.counter == 0:

185 h_mm = 0.5 * (h_mm + h_mm.T.conj())

186 kpt.eps_n[:nbs] = np.linalg.eigh(h_mm)[0]

187 self.counter += 1

188

189 if constraints:

190 constrain_grad(grad, constraints, ind_up)

191

192 return 2.0 * grad, error

193

194 def get_orbital_potential_matrix(self, f_n, C_nM, kpt,

195 wfs, setup, m, timer):

196 """

197 :param f_n:

198 :param C_nM:

199 :param kpt:

200 :param wfs:

201 :param setup:

202 :return: F_i = <m|v_i + PAW_i|n > for occupied

203 F_i = 0 for unoccupied,

204 SI energies

205

206 To calculate this, we need to calculate

207 orbital-depended potential and PAW corrections to it.

208

209 Fot this, we need firstly to get orbitals densities.

210

211 """

212 kpoint = self.n_kps * kpt.s + kpt.q

213 n_set = C_nM.shape[1]

214 F_MM = np.zeros(shape=(n_set, n_set),

215 dtype=self.dtype)

216 # get orbital-density

217 timer.start('Construct Density, Charge, and DM')

218 nt_G, Q_aL, D_ap = \

219 self.get_density(f_n,

220 C_nM, kpt,

221 wfs, setup, m)

222 timer.stop('Construct Density, Charge, and DM')

223

224 # calculate sic energy,

225 # sic pseudo-potential and Hartree

226 timer.start('Get Pseudo Potential')

227 e_sic_m, vt_mG, vHt_g = \

228 self.get_pseudo_pot(nt_G, Q_aL, m, kpoint, timer)

229 timer.stop('Get Pseudo Potential')

230

231 # calculate PAW corrections

232 timer.start('PAW')

233 e_sic_paw_m, dH_ap = \

234 self.get_paw_corrections(D_ap, vHt_g, timer)

235 timer.stop('PAW')

236

237 # total sic:

238 e_sic_m += e_sic_paw_m

239

240 timer.start('ODD Potential Matrices')

241 Vt_MM = np.zeros_like(F_MM)

242 self.bfs.calculate_potential_matrix(vt_mG, Vt_MM, kpt.q)

243 # make matrix hermitian

244 ind_l = np.tril_indices(Vt_MM.shape[0], -1)

245 Vt_MM[(ind_l[1], ind_l[0])] = Vt_MM[ind_l].conj()

246 timer.stop('ODD Potential Matrices')

247

248 # PAW:

249 timer.start('Potential matrix - PAW')

250 for a, dH_p in dH_ap.items():

251 P_Mj = wfs.P_aqMi[a][kpt.q]

252 dH_ij = unpack_hermitian(dH_p)

253

254 if self.dtype == complex:

255 F_MM += P_Mj @ dH_ij @ P_Mj.T.conj()

256 else:

257 F_MM += P_Mj @ dH_ij @ P_Mj.T

258

259 if self.dtype == complex:

260 F_MM += Vt_MM.astype(complex)

261 else:

262 F_MM += Vt_MM

263 if self.sic_coarse_grid:

264 self.cgd.comm.sum(F_MM)

265 else:

266 self.finegd.comm.sum(F_MM)

267 timer.stop('Potential matrix - PAW')

268

269 return F_MM, e_sic_m * f_n[m]

270

271 def get_density(self, f_n, C_nM, kpt,

272 wfs, setup, m):

273

274 # construct orbital density matrix

275 if f_n[m] > 1.0 + 1.0e-4:

276 occup_factor = f_n[m] / (3.0 - wfs.nspins)

277 else:

278 occup_factor = f_n[m]

279 rho_MM = occup_factor * np.outer(C_nM[m].conj(), C_nM[m])

280

281 nt_G = self.cgd.zeros()

282 self.bfs.construct_density(rho_MM, nt_G, kpt.q)

283

284 # calculate atomic density matrix and

285 # compensation charges

286 D_ap = {}

287 Q_aL = {}

288

289 for a in wfs.P_aqMi.keys():

290 P_Mi = wfs.P_aqMi[a][kpt.q]

291 D_ii = np.zeros((wfs.P_aqMi[a].shape[2],

292 wfs.P_aqMi[a].shape[2]),

293 dtype=self.dtype)

294 rhoP_Mi = rho_MM @ P_Mi

295 D_ii = P_Mi.T.conj() @ rhoP_Mi

296 if self.dtype == complex:

297 D_ap[a] = D_p = pack_density(D_ii.real)

298 else:

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

300

301 Q_aL[a] = np.dot(D_p, setup[a].Delta_pL)

302

303 return nt_G, Q_aL, D_ap

304

305 def get_pseudo_pot(self, nt, Q_aL, i, kpoint, timer):

306

307 if self.sic_coarse_grid is False:

308 # change to fine grid

309 vt_sg = self.finegd.zeros(2)

310 vHt_g = self.finegd.zeros()

311 nt_sg = self.finegd.zeros(2)

312 else:

313 vt_sg = self.cgd.zeros(2)

314 vHt_g = self.cgd.zeros()

315 nt_sg = self.cgd.zeros(2)

316

317 if self.sic_coarse_grid is False:

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

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

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

321 else:

322 nt_sg[0] = nt

323

324 timer.start('ODD XC 3D grid')

325 if self.sic_coarse_grid is False:

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

327 else:

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

329 timer.stop('ODD XC 3D grid')

330 vt_sg[0] *= -self.beta_x

331

332 # Hartree

333 if self.sic_coarse_grid is False:

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

335 else:

336 self.ghat_cg.add(nt_sg[0], Q_aL)

337

338 timer.start('ODD Poisson')

339 if self.store_potentials:

340 if self.sic_coarse_grid:

341 vHt_g = self.old_pot[kpoint][i]

342 else:

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

344 vHt_g)

345 self.poiss.solve(vHt_g, nt_sg[0],

346 zero_initial_phi=self.store_potentials,

347 timer=timer)

348 if self.store_potentials:

349 if self.sic_coarse_grid:

350 self.old_pot[kpoint][i] = vHt_g.copy()

351 else:

352 self.restrictor.apply(vHt_g, self.old_pot[kpoint][i])

353

354 timer.stop('ODD Poisson')

355

356 timer.start('ODD Hartree integrate')

357 if self.sic_coarse_grid is False:

358 ec = 0.5 * self.finegd.integrate(nt_sg[0] * vHt_g)

359 else:

360 ec = 0.5 * self.cgd.integrate(nt_sg[0] * vHt_g)

361

362 timer.stop('ODD Hartree integrate')

363 vt_sg[0] -= vHt_g * self.beta_c

364 if self.sic_coarse_grid is False:

365 vt_G = self.cgd.zeros()

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

367 else:

368 vt_G = vt_sg[0]

369

370 return np.array([-ec * self.beta_c,

371 -e_xc * self.beta_x]), vt_G, vHt_g

372

373 def get_paw_corrections(self, D_ap, vHt_g, timer):

374

375 # XC-PAW

376 timer.start('xc-PAW')

377 dH_ap = {}

378 exc = 0.0

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

380 setup = self.setups[a]

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

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

383 exc += self.xc.calculate_paw_correction(setup, D_sp,

384 dH_sp,

385 addcoredensity=False,

386 a=a)

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

388 timer.stop('xc-PAW')

389

390 # Hartree-PAW

391 timer.start('Hartree-PAW')

392 ec = 0.0

393 timer.start('ghat-PAW')

394 if self.sic_coarse_grid is False:

395 W_aL = self.ghat.dict()

396 self.ghat.integrate(vHt_g, W_aL)

397 else:

398 W_aL = self.ghat_cg.dict()

399 self.ghat_cg.integrate(vHt_g, W_aL)

400 timer.stop('ghat-PAW')

401

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

403 setup = self.setups[a]

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

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

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

407 W_aL[a])) * self.beta_c

408 timer.stop('Hartree-PAW')

409

410 timer.start('Wait for sum')

411 if self.sic_coarse_grid is False:

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

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

414 else:

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

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

417 timer.stop('Wait for sum')

418

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

420

421 def get_odd_corrections_to_forces(self, wfs, dens):

422

423 timer = wfs.timer

424

425 timer.start('ODD corrections')

426

427 natoms = len(wfs.setups)

428 F_av = np.zeros((natoms, 3))

429 Ftheta_av = np.zeros_like(F_av)

430 Frho_av = np.zeros_like(F_av)

431 Fatom_av = np.zeros_like(F_av)

432 Fpot_av = np.zeros_like(F_av)

433 Fhart_av = np.zeros_like(F_av)

434

435 ksl = wfs.ksl

436 nao = ksl.nao

437 mynao = ksl.mynao

438 dtype = wfs.dtype

439 manytci = wfs.manytci

440 gd = wfs.gd

441

442 Mstart = ksl.Mstart

443 Mstop = ksl.Mstop

444 n_kps = wfs.kd.nibzkpts

445

446 timer.start('TCI derivative')

447 dThetadR_qvMM, dTdR_qvMM = manytci.O_qMM_T_qMM(

448 gd.comm, Mstart, Mstop, False, derivative=True)

449

450 dPdR_aqvMi = manytci.P_aqMi(

451 self.bfs.my_atom_indices, derivative=True)

452 gd.comm.sum(dThetadR_qvMM)

453 gd.comm.sum(dTdR_qvMM)

454 timer.stop('TCI derivative')

455

456 my_atom_indices = self.bfs.my_atom_indices

457 atom_indices = self.bfs.atom_indices

458

459 def _slices(indices):

460 for a in indices:

461 M1 = self.bfs.M_a[a] - Mstart

462 M2 = M1 + self.setups[a].nao

463 if M2 > 0:

464 yield a, max(0, M1), M2

465

466 def slices():

467 return _slices(atom_indices)

468

469 def my_slices():

470 return _slices(my_atom_indices)

471

472 #

473 # ----- -----

474 # \ -1 \ *

475 # E = ) S H rho = ) c eps f c

476 # mu nu / mu x x z z nu / n mu n n n nu

477 # ----- -----

478 # x z n

479 #

480 # We use the transpose of that matrix. The first form is used

481 # if rho is given, otherwise the coefficients are used.

482

483 for kpt in wfs.kpt_u:

484 u = kpt.s * n_kps + kpt.q

485 f_n = kpt.f_n

486 n_occ = 0

487 for f in f_n:

488 if f > 1.0e-10:

489 n_occ += 1

490

491 for m in range(n_occ):

492

493 # calculate orbital-density matrix

494 rho_xMM = kpt.f_n[m] * np.outer(

495 kpt.C_nM[m].conj(), kpt.C_nM[m]) / (3.0 - wfs.nspins)

496

497 F_MM = self.get_orbital_potential_matrix(

498 f_n, kpt.C_nM, kpt, wfs, self.setups, m, timer)[0]

499

500 sfrhoT_MM = np.linalg.solve(

501 wfs.S_qMM[kpt.q], F_MM @ rho_xMM).T.copy()

502

503 del F_MM

504

505 # Density matrix contribution due to basis overlap

506 #

507 # ----- d Theta

508 # a \ mu nu

509 # F += -2 Re ) ------------ E

510 # / d R nu mu

511 # ----- mu nu

512 # mu in a; nu

513 #

514

515 dThetadRE_vMM = (dThetadR_qvMM[kpt.q] *

516 sfrhoT_MM[np.newaxis]).real

517 for a, M1, M2 in my_slices():

518 Ftheta_av[a, :] += \

519 -2.0 * dThetadRE_vMM[:, M1:M2].sum(-1).sum(-1)

520 del dThetadRE_vMM

521

522 # Density matrix contribution from PAW correction

523 #

524 # ----- -----

525 # a \ a \ b

526 # F += 2 Re ) Z E - 2 Re ) Z E

527 # / mu nu nu mu / mu nu nu mu

528 # ----- -----

529 # mu nu b; mu in a; nu

530 #

531 # with

532 # b*

533 # ----- dP

534 # b \ i mu b b

535 # Z = ) -------- dS P

536 # mu nu / dR ij j nu

537 # ----- b mu

538 # ij

539 #

540 work_MM = np.zeros((mynao, nao), dtype)

541 ZE_MM = None

542 for b in my_atom_indices:

543 setup = self.setups[b]

544 dO_ii = np.asarray(setup.dO_ii, dtype)

545 dOP_iM = dO_ii @ wfs.P_aqMi[b][kpt.q].T.conj()

546

547 for v in range(3):

548 work_MM = \

549 dPdR_aqvMi[b][kpt.q][v][Mstart:Mstop] @ dOP_iM

550 ZE_MM = (work_MM * sfrhoT_MM).real

551 for a, M1, M2 in slices():

552 dE = 2 * ZE_MM[M1:M2].sum()

553 Frho_av[a, v] -= dE # the "b; mu in a; nu" term

554 Frho_av[b, v] += dE # the "mu nu" term

555 del work_MM, ZE_MM

556

557 # Potential contribution

558 #

559 # ----- / d Phi (r)

560 # a \ | mu ~

561 # F += -2 Re ) | ---------- v (r) Phi (r) dr rho

562 # / | d R nu nu mu

563 # ----- / a

564 # mu in a; nu

565 #

566

567 nt_G, Q_aL, D_ap = \

568 self.get_density(f_n, kpt.C_nM, kpt, wfs, self.setups, m)

569

570 e_sic_m, vt_mG, vHt_g = \

571 self.get_pseudo_pot(nt_G, Q_aL, m, u, timer)

572 e_sic_paw_m, dH_ap = \

573 self.get_paw_corrections(D_ap, vHt_g, timer)

574

575 Fpot_av += \

576 self.bfs.calculate_force_contribution(

577 vt_mG, rho_xMM, kpt.q)

578

579 # Atomic density contribution

580 # ----- -----

581 # a \ a \ b

582 # F += -2 Re ) A rho + 2 Re ) A rho

583 # / mu nu nu mu / mununumu

584 # ----- -----

585 # mu nu b; mu in a; nu

586 #

587 # b*

588 # ----- d P

589 # b \ i mu b b

590 # A = ) ------- dH P

591 # mu nu / d R ij j nu

592 # ----- b mu

593 # ij

594 #

595 for b in my_atom_indices:

596 H_ii = np.asarray(unpack_hermitian(dH_ap[b]),

597 dtype)

598 HP_iM = H_ii @ wfs.P_aqMi[b][kpt.q].T.conj()

599 for v in range(3):

600 dPdR_Mi = \

601 dPdR_aqvMi[b][kpt.q][v][Mstart:Mstop]

602 ArhoT_MM = \

603 ((dPdR_Mi @ HP_iM) * rho_xMM.T).real

604 for a, M1, M2 in slices():

605 dE = 2 * ArhoT_MM[M1:M2].sum()

606 Fatom_av[a, v] += dE # the "b; mu in a; nu" term

607 Fatom_av[b, v] -= dE # the "mu nu" term

608

609 # contribution from hartree

610 if self.sic_coarse_grid is False:

611 ghat_aLv = dens.ghat.dict(derivative=True)

612

613 dens.ghat.derivative(vHt_g, ghat_aLv)

614 for a, dF_Lv in ghat_aLv.items():

615 Fhart_av[a] -= self.beta_c * np.dot(Q_aL[a], dF_Lv)

616 else:

617 ghat_aLv = self.ghat_cg.dict(derivative=True)

618

619 self.ghat_cg.derivative(vHt_g, ghat_aLv)

620 for a, dF_Lv in ghat_aLv.items():

621 Fhart_av[a] -= self.beta_c * np.dot(Q_aL[a], dF_Lv)

622

623 F_av += Fpot_av + Ftheta_av + Frho_av + Fatom_av + Fhart_av

624

625 wfs.gd.comm.sum(F_av, 0)

626

627 timer.start('Wait for sum')

628 ksl.orbital_comm.sum(F_av)

629 if wfs.bd.comm.rank == 0:

630 wfs.kd.comm.sum(F_av, 0)

631

632 wfs.world.broadcast(F_av, 0)

633 timer.stop('Wait for sum')

634

635 timer.stop('ODD corrections')

636

637 return F_av * (3.0 - wfs.nspins)

638

639 def get_lagrange_matrices(self, h_mm, c_nm, f_n, kpt,

640 wfs, occupied_only=False,

641 update_eigenvalues=False):

642 n_occ = 0

643 nbands = len(f_n)

644 while n_occ < nbands and f_n[n_occ] > 1e-10:

645 n_occ += 1

646

647 if occupied_only is True:

648 nbs = n_occ

649 else:

650 nbs = c_nm.shape[0]

651 n_set = c_nm.shape[1]

652

653 hc_mn = np.dot(h_mm.conj(), c_nm[:nbs].T)

654 h_mm = np.dot(c_nm[:nbs].conj(), hc_mn)

655 # odd part

656 b_mn = np.zeros(shape=(n_set, nbs), dtype=self.dtype)

657 e_total_sic = np.array([])

658 for n in range(n_occ):

659 F_MM, sic_energy_n =\

660 self.get_orbital_potential_matrix(f_n, c_nm, kpt,

661 wfs, wfs.setups,

662 n, wfs.timer)

663

664 b_mn[:, n] = np.dot(c_nm[n], F_MM.conj()).T

665 e_total_sic = np.append(e_total_sic, sic_energy_n, axis=0)

666 l_odd = np.dot(c_nm[:nbs].conj(), b_mn)

667

668 k = wfs.eigensolver.kpointval(kpt)

669

670 fullham = h_mm + 0.5 * (l_odd + l_odd.T.conj())

671 fullham[:n_occ, n_occ:] = 0.0

672 fullham[n_occ:, :n_occ] = 0.0

673

674 self.lagr_diag_s[k] = np.diagonal(fullham).real

675

676 if update_eigenvalues:

677 eigval, eigvec = np.linalg.eigh(fullham[0:n_occ, 0:n_occ])

678 kpt.eps_n[0:n_occ] = eigval

679 eigval, eigvec = np.linalg.eigh(fullham[n_occ:nbs, n_occ:nbs])

680 kpt.eps_n[n_occ:nbs] = eigval

681

682 return h_mm, l_odd

683

684

685def constrain_grad(grad, constraints, ind_up):

686 """

687 Zero out the components of the gradient that are constrained, so that no

688 optimization step is taken along the constrained degrees of freedom (It

689 would be better not to evaluate these components of the gradient to begin

690 with.).

691 """

692

693 for con in constraints:

694 num = -1

695 for ind1, ind2 in zip(ind_up[0], ind_up[1]):

696 num += 1

697 if con == [ind1, ind2]:

698 grad[num] = 0.0

699 return grad