Coverage for gpaw/xc/hybrid.py: 92%
404 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1# Copyright (C) 2003 CAMP
2# Please see the accompanying LICENSE file for further information.
4"""This module provides all the classes and functions associated with the
5evaluation of exact exchange.
6"""
8from math import exp, ceil
9import numpy as np
11from gpaw.atom.configurations import core_states
12from gpaw.gaunt import gaunt
13from gpaw.lfc import LFC
14from gpaw.poisson import PoissonSolver
15from gpaw.helmholtz import HelmholtzSolver
16from gpaw.transformers import Transformer
17from gpaw.utilities import (hartree, pack_density, pack_hermitian,
18 packed_index, unpack_density, unpack_hermitian)
19from gpaw.utilities.blas import mmm
20from gpaw.utilities.tools import symmetrize
21from gpaw.xc import XC
22from gpaw.xc.functional import XCFunctional
23from gpaw.xc.kernel import XCNull
24from gpaw.setup import CachedYukawaInteractions
27class HybridXCBase(XCFunctional):
28 orbital_dependent = True
30 def __init__(self, name, stencil=2, hybrid=None, xc=None, omega=None):
31 """Mix standard functionals with exact exchange.
33 name: str
34 Name of hybrid functional.
35 hybrid: float
36 Fraction of exact exchange.
37 xc: str or XCFunctional object
38 Standard DFT functional with scaled down exchange.
39 """
41 rsf_functionals = { # Parameters can also be taken from libxc
42 'CAMY-BLYP': { # Akinaga, Ten-no CPL 462 (2008) 348-351
43 'alpha': 0.2,
44 'beta': 0.8,
45 'omega': 0.44,
46 'cam': True,
47 'rsf': 'Yukawa',
48 'xc': 'HYB_GGA_XC_CAMY_BLYP'
49 },
50 'CAMY-B3LYP': { # Seth, Ziegler JCTC 8 (2012) 901-907
51 'alpha': 0.19,
52 'beta': 0.46,
53 'omega': 0.34,
54 'cam': True,
55 'rsf': 'Yukawa',
56 'xc': 'HYB_GGA_XC_CAMY_B3LYP'
57 },
58 'LCY-BLYP': { # Seth, Ziegler JCTC 8 (2012) 901-907
59 'alpha': 0.0,
60 'beta': 1.0,
61 'omega': 0.75,
62 'cam': False,
63 'rsf': 'Yukawa',
64 'xc': 'HYB_GGA_XC_LCY_BLYP'
65 },
66 'LCY-PBE': { # Seth, Ziegler JCTC 8 (2012) 901-907
67 'alpha': 0.0,
68 'beta': 1.0,
69 'omega': 0.75,
70 'cam': False,
71 'rsf': 'Yukawa',
72 'xc': 'HYB_GGA_XC_LCY_PBE'
73 }
74 }
75 self.omega = None
76 self.cam_alpha = None
77 self.cam_beta = None
78 self.is_cam = False
79 self.rsf = None
81 def _xc(name):
82 return {'name': name, 'stencil': stencil}
84 if name == 'EXX':
85 hybrid = 1.0
86 xc = XC(XCNull())
87 elif name == 'PBE0':
88 hybrid = 0.25
89 xc = XC(_xc('HYB_GGA_XC_PBEH'))
90 elif name == 'B3LYP':
91 hybrid = 0.2
92 xc = XC(_xc('HYB_GGA_XC_B3LYP'))
93 elif name in rsf_functionals:
94 rsf_functional = rsf_functionals[name]
95 self.cam_alpha = rsf_functional['alpha']
96 self.cam_beta = rsf_functional['beta']
97 self.omega = rsf_functional['omega']
98 self.is_cam = rsf_functional['cam']
99 self.rsf = rsf_functional['rsf']
100 xc = XC(rsf_functional['xc'])
101 hybrid = self.cam_alpha + self.cam_beta
103 if isinstance(xc, (str, dict)):
104 xc = XC(xc)
106 self.hybrid = float(hybrid)
107 self.xc = xc
108 if omega is not None:
109 omega = float(omega)
110 if self.omega is not None and self.omega != omega:
111 self.xc.kernel.set_omega(omega)
112 # Needed to tune omega for RSF
113 self.omega = omega
114 XCFunctional.__init__(self, name, xc.type)
116 def todict(self):
117 return {'name': self.name,
118 'hybrid': self.hybrid,
119 'excitation': self.excitation,
120 'excited': self.excited,
121 'xc': self.xc.todict(),
122 'omega': self.omega}
124 def tostring(self):
125 """Return string suitable to generate xc from string."""
126 xc_dict = self.todict()
127 for test_key in ['name', 'xc', 'kernel', 'type']:
128 if test_key in xc_dict:
129 del xc_dict[test_key]
130 return self.name + ':' + ':'.join([(k + '=' + repr(v))
131 for k, v in xc_dict.items()])
133 def get_setup_name(self):
134 return 'PBE'
137class HybridXC(HybridXCBase):
138 def __init__(self, name, hybrid=None, xc=None,
139 finegrid=False, unocc=False, omega=None,
140 excitation=None, excited=0, stencil=2):
141 """Mix standard functionals with exact exchange.
143 finegrid: boolean
144 Use fine grid for energy functional evaluations ?
145 unocc: boolean
146 Apply vxx also to unoccupied states ?
147 omega: float
148 RSF mixing parameter
149 excitation: string:
150 Apply operator for improved virtual orbitals
151 to unocc states? Possible modes:
152 singlet: excitations to singlets
153 triplet: excitations to triplets
154 average: average between singlets and tripletts
155 see f.e. https://doi.org/10.1021/acs.jctc.8b00238
156 excited: number
157 Band to excite from - counted from HOMO downwards
159 """
160 self.finegrid = finegrid
161 self.unocc = unocc
162 self.excitation = excitation
163 self.excited = excited
164 HybridXCBase.__init__(self, name, hybrid=hybrid, xc=xc, omega=omega,
165 stencil=stencil)
167 # Note: self.omega may not be identical to omega!
168 self.yukawa_interactions = CachedYukawaInteractions(self.omega)
170 def calculate_paw_correction(self, setup, D_sp, dEdD_sp=None,
171 addcoredensity=True, a=None):
172 return self.xc.calculate_paw_correction(setup, D_sp, dEdD_sp,
173 addcoredensity, a)
175 def initialize(self, density, hamiltonian, wfs):
176 assert wfs.kd.gamma
177 self.xc.initialize(density, hamiltonian, wfs)
178 self.kpt_comm = wfs.kd.comm
179 self.nspins = wfs.nspins
180 self.setups = wfs.setups
181 self.density = density
182 self.kpt_u = wfs.kpt_u
183 self.exx_s = np.zeros(self.nspins)
184 self.ekin_s = np.zeros(self.nspins)
185 self.nocc_s = np.empty(self.nspins, int)
187 self.gd = density.gd
188 self.redistributor = density.redistributor
190 use_charge_center = hamiltonian.poisson.use_charge_center
191 # XXX How do we construct a copy of the Poisson solver of the
192 # Hamiltonian? We don't know what class it is, etc., but gd
193 # may differ.
194 # XXX One might consider using a charged centered compensation
195 # charge for the PoissonSolver in the case of EXX as standard
196 self.poissonsolver = PoissonSolver(
197 'fd', eps=1e-12, use_charge_center=use_charge_center)
198 # self.poissonsolver = hamiltonian.poisson
200 if self.finegrid:
201 self.finegd = self.gd.refine()
202 # XXX Taking restrictor from Hamiltonian will not work in PW mode,
203 # will it? I think this supports only real-space mode.
204 # self.restrictor = hamiltonian.restrictor
205 self.restrictor = Transformer(self.finegd, self.gd, 3)
206 self.interpolator = Transformer(self.gd, self.finegd, 3)
207 else:
208 self.finegd = self.gd
210 self.ghat = LFC(self.finegd,
211 [setup.ghat_l for setup in density.setups],
212 integral=np.sqrt(4 * np.pi), forces=True)
213 self.poissonsolver.set_grid_descriptor(self.finegd)
214 if self.rsf == 'Yukawa':
215 omega2 = self.omega**2
216 self.screened_poissonsolver = HelmholtzSolver(
217 k2=-omega2, eps=1e-12, nn=3,
218 use_charge_center=use_charge_center)
219 self.screened_poissonsolver.set_grid_descriptor(self.finegd)
221 def set_positions(self, spos_ac):
222 self.ghat.set_positions(spos_ac)
224 def calculate(self, gd, n_sg, v_sg=None, e_g=None):
225 # Normal XC contribution:
226 exc = self.xc.calculate(gd, n_sg, v_sg, e_g)
227 # Note that the quantities passed are on the
228 # density/Hamiltonian grids!
229 # They may be distributed differently from own quantities.
230 self.ekin = self.kpt_comm.sum_scalar(self.ekin_s.sum())
231 return exc + self.kpt_comm.sum_scalar(self.exx_s.sum())
233 def calculate_exx(self):
234 for kpt in self.kpt_u:
235 self.apply_orbital_dependent_hamiltonian(kpt, kpt.psit_nG)
237 def apply_orbital_dependent_hamiltonian(self, kpt, psit_nG,
238 Htpsit_nG=None, dH_asp=None):
239 if kpt.f_n is None:
240 return
242 deg = 2 // self.nspins # Spin degeneracy
243 hybrid = self.hybrid
244 P_ani = kpt.P_ani
245 setups = self.setups
246 is_cam = self.is_cam
248 vt_g = self.finegd.empty()
249 if self.gd is not self.finegd:
250 vt_G = self.gd.empty()
251 if self.rsf == 'Yukawa':
252 y_vt_g = self.finegd.empty()
253 # if self.gd is not self.finegd:
254 # y_vt_G = self.gd.empty()
256 nocc = int(ceil(kpt.f_n.sum())) // (3 - self.nspins)
257 if self.excitation is not None:
258 ex_band = nocc - self.excited - 1
259 if self.excitation == 'singlet':
260 ex_weight = -1
261 elif self.excitation == 'triplet':
262 ex_weight = +1
263 else:
264 ex_weight = 0
266 if self.unocc or self.excitation is not None:
267 nbands = len(kpt.f_n)
268 else:
269 nbands = nocc
270 self.nocc_s[kpt.s] = nocc
272 if Htpsit_nG is not None:
273 kpt.vt_nG = self.gd.empty(nbands)
274 kpt.vxx_ani = {}
275 kpt.vxx_anii = {}
276 for a, P_ni in P_ani.items():
277 I = P_ni.shape[1]
278 kpt.vxx_ani[a] = np.zeros((nbands, I))
279 kpt.vxx_anii[a] = np.zeros((nbands, I, I))
281 exx = 0.0
282 ekin = 0.0
284 # XXXX nbands can be different numbers on different cpus!
285 # That means some will execute the loop and others not.
286 # And deadlocks with augment-grids.
288 # Determine pseudo-exchange
289 for n1 in range(nbands):
290 psit1_G = psit_nG[n1]
291 f1 = kpt.f_n[n1] / deg
292 for n2 in range(n1, nbands):
293 psit2_G = psit_nG[n2]
294 f2 = kpt.f_n[n2] / deg
295 if n1 != n2 and f1 == 0 and f1 == f2:
296 continue # Don't work on double unocc. bands
297 # Double count factor:
298 dc = (1 + (n1 != n2)) * deg
299 nt_G, rhot_g = self.calculate_pair_density(n1, n2, psit_nG,
300 P_ani)
301 vt_g[:] = 0.0
302 # XXXXX This will go wrong because we are solving the
303 # Poisson equation on the distribution of gd, not finegd
304 # Or maybe it's fixed now
306 self.poissonsolver.solve(vt_g, -rhot_g,
307 charge=-float(n1 == n2),
308 zero_initial_phi=True)
309 vt_g *= hybrid
310 if self.rsf == 'Yukawa':
311 y_vt_g[:] = 0.0
312 self.screened_poissonsolver.solve(
313 y_vt_g, -rhot_g, charge=-float(n1 == n2),
314 zero_initial_phi=True)
315 if is_cam: # Cam like correction
316 y_vt_g *= self.cam_beta
317 else:
318 y_vt_g *= hybrid
319 vt_g -= y_vt_g
320 if self.gd is self.finegd:
321 vt_G = vt_g
322 else:
323 self.restrictor.apply(vt_g, vt_G)
325 # Integrate the potential on fine and coarse grids
326 int_fine = self.finegd.integrate(vt_g * rhot_g)
327 int_coarse = self.gd.integrate(vt_G * nt_G)
328 if self.gd.comm.rank == 0: # only add to energy on master CPU
329 exx += 0.5 * dc * f1 * f2 * int_fine
330 ekin -= dc * f1 * f2 * int_coarse
331 if Htpsit_nG is not None:
332 Htpsit_nG[n1] += f2 * vt_G * psit2_G
333 if n1 == n2:
334 kpt.vt_nG[n1] = f1 * vt_G
335 if self.excitation is not None and n1 == ex_band:
336 Htpsit_nG[nocc:] += f1 * vt_G * psit_nG[nocc:]
337 else:
338 if self.excitation is None or n1 != ex_band \
339 or n2 < nocc:
340 Htpsit_nG[n2] += f1 * vt_G * psit1_G
341 else:
342 Htpsit_nG[n2] += f1 * ex_weight * vt_G * psit1_G
344 # Update the vxx_uni and vxx_unii vectors of the nuclei,
345 # used to determine the atomic hamiltonian, and the
346 # residuals
347 v_aL = self.ghat.dict()
348 self.ghat.integrate(vt_g, v_aL)
349 for a, v_L in v_aL.items():
350 v_ii = unpack_hermitian(
351 np.dot(setups[a].Delta_pL, v_L))
352 v_ni = kpt.vxx_ani[a]
353 v_nii = kpt.vxx_anii[a]
354 P_ni = P_ani[a]
355 v_ni[n1] += f2 * np.dot(v_ii, P_ni[n2])
356 if n1 != n2:
357 if self.excitation is None or n1 != ex_band or \
358 n2 < nocc:
359 v_ni[n2] += f1 * np.dot(v_ii, P_ni[n1])
360 else:
361 v_ni[n2] += f1 * ex_weight * \
362 np.dot(v_ii, P_ni[n1])
363 else:
364 # XXX Check this:
365 v_nii[n1] = f1 * v_ii
366 if self.excitation is not None and n1 == ex_band:
367 for nuoc in range(nocc, nbands):
368 v_ni[nuoc] += f1 * \
369 np.dot(v_ii, P_ni[nuoc])
371 def calculate_vv(ni, D_ii, M_pp, weight, addme=False):
372 """Calculate the local corrections depending on Mpp."""
373 dexx = 0
374 dekin = 0
375 if not addme:
376 addsign = -2.0
377 else:
378 addsign = 2.0
379 for i1 in range(ni):
380 for i2 in range(ni):
381 A = 0.0
382 for i3 in range(ni):
383 p13 = packed_index(i1, i3, ni)
384 for i4 in range(ni):
385 p24 = packed_index(i2, i4, ni)
386 A += M_pp[p13, p24] * D_ii[i3, i4]
387 p12 = packed_index(i1, i2, ni)
388 if Htpsit_nG is not None:
389 dH_p[p12] += addsign * weight / \
390 deg * A / ((i1 != i2) + 1)
391 dekin += 2 * weight / deg * D_ii[i1, i2] * A
392 dexx -= weight / deg * D_ii[i1, i2] * A
393 return (dexx, dekin)
395 # Apply the atomic corrections to the energy and the Hamiltonian
396 # matrix
397 for a, P_ni in P_ani.items():
398 setup = setups[a]
400 if Htpsit_nG is not None:
401 # Add non-trivial corrections the Hamiltonian matrix
402 h_nn = symmetrize(np.inner(P_ni[:nbands],
403 kpt.vxx_ani[a][:nbands]))
404 ekin -= np.dot(kpt.f_n[:nbands], h_nn.diagonal())
406 dH_p = dH_asp[a][kpt.s]
408 # Get atomic density and Hamiltonian matrices
409 D_p = self.density.D_asp[a][kpt.s]
410 D_ii = unpack_density(D_p)
411 ni = len(D_ii)
413 # Add atomic corrections to the valence-valence exchange energy
414 # --
415 # > D C D
416 # -- ii iiii ii
417 (dexx, dekin) = calculate_vv(ni, D_ii, setup.M_pp, hybrid)
418 ekin += dekin
419 exx += dexx
421 if self.rsf is not None:
422 Mg_pp = self.yukawa_interactions.get_Mg_pp(setup)
423 if is_cam:
424 (dexx, dekin) = calculate_vv(
425 ni, D_ii, Mg_pp, self.cam_beta, addme=True)
426 else:
427 (dexx, dekin) = calculate_vv(
428 ni, D_ii, Mg_pp, hybrid, addme=True)
429 ekin -= dekin
430 exx -= dexx
431 # Add valence-core exchange energy
432 # --
433 # > X D
434 # -- ii ii
435 if setup.X_p is not None:
436 exx -= hybrid * np.dot(D_p, setup.X_p)
437 if Htpsit_nG is not None:
438 dH_p -= hybrid * setup.X_p
439 ekin += hybrid * np.dot(D_p, setup.X_p)
441 if self.rsf == 'Yukawa' and setup.X_pg is not None:
442 if is_cam:
443 thybrid = self.cam_beta # 0th order
444 else:
445 thybrid = hybrid
446 exx += thybrid * np.dot(D_p, setup.X_pg)
447 if Htpsit_nG is not None:
448 dH_p += thybrid * setup.X_pg
449 ekin -= thybrid * np.dot(D_p, setup.X_pg)
450 elif self.rsf == 'Yukawa' and setup.X_pg is None:
451 thybrid = exp(-3.62e-2 * self.omega) # educated guess
452 if is_cam:
453 thybrid *= self.cam_beta
454 else:
455 thybrid *= hybrid
456 exx += thybrid * np.dot(D_p, setup.X_p)
457 if Htpsit_nG is not None:
458 dH_p += thybrid * setup.X_p
459 ekin -= thybrid * np.dot(D_p, setup.X_p)
460 # Add core-core exchange energy
461 if kpt.s == 0:
462 if self.rsf is None or is_cam:
463 if is_cam:
464 exx += self.cam_alpha * setup.ExxC
465 else:
466 exx += hybrid * setup.ExxC
468 self.exx_s[kpt.s] = self.gd.comm.sum_scalar(exx)
469 self.ekin_s[kpt.s] = self.gd.comm.sum_scalar(ekin)
471 def correct_hamiltonian_matrix(self, kpt, H_nn):
472 if not hasattr(kpt, 'vxx_ani'):
473 return
475 # if self.gd.comm.rank > 0:
476 # H_nn[:] = 0.0
478 nocc = self.nocc_s[kpt.s]
479 nbands = len(kpt.vt_nG)
480 for a, P_ni in kpt.P_ani.items():
481 H_nn[:nbands, :nbands] += symmetrize(np.inner(P_ni[:nbands],
482 kpt.vxx_ani[a]))
483 # self.gd.comm.sum(H_nn)
485 if not self.unocc or self.excitation is not None:
486 H_nn[:nocc, nocc:] = 0.0
487 H_nn[nocc:, :nocc] = 0.0
489 def calculate_pair_density(self, n1, n2, psit_nG, P_ani):
490 Q_aL = {}
491 for a, P_ni in P_ani.items():
492 P1_i = P_ni[n1]
493 P2_i = P_ni[n2]
494 D_ii = np.outer(P1_i, P2_i.conj()).real
495 D_p = pack_density(D_ii)
496 Q_aL[a] = np.dot(D_p, self.setups[a].Delta_pL)
498 nt_G = psit_nG[n1] * psit_nG[n2]
500 if self.finegd is self.gd:
501 nt_g = nt_G
502 else:
503 nt_g = self.finegd.empty()
504 self.interpolator.apply(nt_G, nt_g)
506 rhot_g = nt_g.copy()
507 self.ghat.add(rhot_g, Q_aL)
509 return nt_G, rhot_g
511 def add_correction(self, kpt, psit_xG, Htpsit_xG, P_axi, c_axi, n_x,
512 calculate_change=False):
513 if kpt.f_n is None:
514 return
516 if self.unocc or self.excitation is not None:
517 nocc = len(kpt.vt_nG)
518 else:
519 nocc = self.nocc_s[kpt.s]
521 if calculate_change:
522 for x, n in enumerate(n_x):
523 if n < nocc:
524 Htpsit_xG[x] += kpt.vt_nG[n] * psit_xG[x]
525 for a, P_xi in P_axi.items():
526 c_axi[a][x] += np.dot(kpt.vxx_anii[a][n], P_xi[x])
527 else:
528 for a, c_xi in c_axi.items():
529 c_xi[:nocc] += kpt.vxx_ani[a][:nocc]
531 def rotate(self, kpt, U_nn):
532 if kpt.f_n is None:
533 return
535 U_nn = U_nn.T.copy()
536 nocc = self.nocc_s[kpt.s]
537 if len(kpt.vt_nG) == nocc:
538 U_nn = U_nn[:nocc, :nocc]
539 # gemm(1.0, kpt.vt_nG.copy(), U_nn, 0.0, kpt.vt_nG)
540 n, G1, G2, G3 = kpt.vt_nG.shape
541 vt_nG = kpt.vt_nG.reshape((n, G1 * G2 * G3))
542 mmm(1.0, U_nn, 'N', vt_nG.copy(), 'N', 0.0, vt_nG)
543 for v_ni in kpt.vxx_ani.values():
544 # gemm(1.0, v_ni.copy(), U_nn, 0.0, v_ni)
545 mmm(1.0, U_nn, 'N', v_ni.copy(), 'N', 0.0, v_ni)
546 for v_nii in kpt.vxx_anii.values():
547 # gemm(1.0, v_nii.copy(), U_nn, 0.0, v_nii)
548 n, i, i = v_nii.shape
549 v_nii = v_nii.reshape((n, i**2))
550 mmm(1.0, U_nn, 'N', v_nii.copy(), 'N', 0.0, v_nii)
553def atomic_exact_exchange(atom, type='all'):
554 """Returns the exact exchange energy of the atom defined by the
555 instantiated AllElectron object 'atom'
556 """
557 G_LLL = gaunt(lmax=max(atom.l_j)) # Make gaunt coeff. list
558 Nj = len(atom.n_j) # The total number of orbitals
560 # determine relevant states for chosen type of exchange contribution
561 if type == 'all':
562 nstates = mstates = range(Nj)
563 else:
564 Njcore = core_states(atom.symbol) # The number of core orbitals
565 if type == 'val-val':
566 nstates = mstates = range(Njcore, Nj)
567 elif type == 'core-core':
568 nstates = mstates = range(Njcore)
569 elif type == 'val-core':
570 nstates = range(Njcore, Nj)
571 mstates = range(Njcore)
572 else:
573 raise RuntimeError('Unknown type of exchange: ', type)
575 # Arrays for storing the potential (times radius)
576 vr = np.zeros(atom.N)
577 vrl = np.zeros(atom.N)
579 # do actual calculation of exchange contribution
580 Exx = 0.0
581 for j1 in nstates:
582 # angular momentum of first state
583 l1 = atom.l_j[j1]
585 for j2 in mstates:
586 # angular momentum of second state
587 l2 = atom.l_j[j2]
589 # joint occupation number
590 f12 = 0.5 * (atom.f_j[j1] / (2. * l1 + 1) *
591 atom.f_j[j2] / (2. * l2 + 1))
593 # electron density times radius times length element
594 nrdr = atom.u_j[j1] * atom.u_j[j2] * atom.dr
595 nrdr[1:] /= atom.r[1:]
597 # potential times radius
598 vr[:] = 0.0
600 # L summation
601 for l in range(l1 + l2 + 1):
602 # get potential for current l-value
603 hartree(l, nrdr, atom.r, vrl)
605 # take all m1 m2 and m values of Gaunt matrix of the form
606 # G(L1,L2,L) where L = {l,m}
607 G2 = G_LLL[l1**2:(l1 + 1)**2,
608 l2**2:(l2 + 1)**2,
609 l**2:(l + 1)**2]**2
611 # add to total potential
612 vr += vrl * np.sum(G2)
614 # add to total exchange the contribution from current two states
615 Exx += -.5 * f12 * np.dot(vr, nrdr)
617 # double energy if mixed contribution
618 if type == 'val-core':
619 Exx *= 2.
621 # return exchange energy
622 return Exx
625def constructX(gen, gamma=0):
626 """Construct the X_p^a matrix for the given atom.
628 The X_p^a matrix describes the valence-core interactions of the
629 partial waves.
630 """
631 # initialize attributes
632 uv_j = gen.vu_j # soft valence states * r:
633 lv_j = gen.vl_j # their repective l quantum numbers
634 Nvi = 0
635 for l in lv_j:
636 Nvi += 2 * l + 1 # total number of valence states (including m)
638 # number of core and valence orbitals (j only, i.e. not m-number)
639 Njcore = gen.njcore
640 Njval = len(lv_j)
642 # core states * r:
643 uc_j = gen.u_j[:Njcore]
644 r, dr, N = gen.r, gen.dr, gen.N
645 r2 = r**2
647 # potential times radius
648 vr = np.zeros(N)
650 # initialize X_ii matrix
651 X_ii = np.zeros((Nvi, Nvi))
653 # make gaunt coeff. list
654 lmax = max(gen.l_j[:Njcore] + gen.vl_j)
655 G_LLL = gaunt(lmax=lmax)
657 # sum over core states
658 for jc in range(Njcore):
659 lc = gen.l_j[jc]
661 # sum over first valence state index
662 i1 = 0
663 for jv1 in range(Njval):
664 lv1 = lv_j[jv1]
666 # electron density 1 times radius times length element
667 n1c = uv_j[jv1] * uc_j[jc] * dr
668 n1c[1:] /= r[1:]
670 # sum over second valence state index
671 i2 = 0
672 for jv2 in range(Njval):
673 lv2 = lv_j[jv2]
675 # electron density 2
676 n2c = uv_j[jv2] * uc_j[jc]
677 n2c[1:] /= r2[1:]
679 # sum expansion in angular momenta
680 for l in range(min(lv1, lv2) + lc + 1):
681 # Int density * potential * r^2 * dr:
682 if gamma == 0:
683 vr = gen.rgd.poisson(n2c, l)
684 else:
685 vr = gen.rgd.yukawa(n2c, l, gamma)
686 nv = np.dot(n1c, vr)
688 # expansion coefficients
689 A_mm = X_ii[i1:i1 + 2 * lv1 + 1, i2:i2 + 2 * lv2 + 1]
690 for mc in range(2 * lc + 1):
691 for m in range(2 * l + 1):
692 G1c = G_LLL[lv1**2:(lv1 + 1)**2,
693 lc**2 + mc, l**2 + m]
694 G2c = G_LLL[lv2**2:(lv2 + 1)**2,
695 lc**2 + mc, l**2 + m]
696 A_mm += nv * np.outer(G1c, G2c)
697 i2 += 2 * lv2 + 1
698 i1 += 2 * lv1 + 1
700 # pack X_ii matrix
701 X_p = pack_hermitian(X_ii)
702 return X_p
705def H_coulomb_val_core(paw, u=0):
706 """Short description here.
708 ::
710 core * *
711 // -- i(r) k(r') k(r) j (r')
712 H = || drdr' > ----------------------
713 ij // -- |r - r'|
714 k
715 """
716 H_nn = np.zeros((paw.wfs.bd.nbands, paw.wfs.bd.nbands),
717 dtype=paw.wfs.dtype)
718 for a, P_ni in paw.wfs.kpt_u[u].P_ani.items():
719 X_ii = unpack_hermitian(paw.wfs.setups[a].X_p)
720 H_nn += np.dot(P_ni.conj(), np.dot(X_ii, P_ni.T))
721 paw.wfs.gd.comm.sum(H_nn)
722 from ase.units import Hartree
723 return H_nn * Hartree