Coverage for gpaw/response/screened_interaction.py: 85%
240 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1import numpy as np
2from math import pi
3from gpaw.response.q0_correction import Q0Correction
4from ase.units import Ha
5from ase.dft.kpoints import monkhorst_pack
6from gpaw.kpt_descriptor import KPointDescriptor
7from gpaw.response.temp import DielectricFunctionCalculator
8from gpaw.response.hilbert import GWHilbertTransforms
9from gpaw.response.mpa_interpolation import RESolver
10from gpaw.cgpaw import evaluate_mpa_poly
13class GammaIntegrationMode:
14 def __init__(self, gamma_integration):
15 if isinstance(gamma_integration, GammaIntegrationMode):
16 self.type = gamma_integration.type
17 self.reduced = gamma_integration.reduced
18 self._N = gamma_integration._N
19 return
21 defaults = {'sphere': {'type': 'sphere'},
22 'reciprocal': {'type': 'reciprocal'},
23 'reciprocal2D': {'type': 'reciprocal', 'reduced': True},
24 '1BZ': {'type': '1BZ'},
25 '1BZ2D': {'type': '1BZ', 'reduced': True},
26 'WS': {'type': 'WS'}}
28 if isinstance(gamma_integration, int):
29 raise TypeError("gamma_integration=INT is no longer supported. "
30 "Please start using the new notations, as is given"
31 " in the documentation in gpaw/response/g0w0.py"
32 " of __init__ of class G0W0.")
34 if isinstance(gamma_integration, str):
35 gamma_integration = defaults[gamma_integration]
37 self.type = gamma_integration['type']
38 self.reduced = gamma_integration.get('reduced', False)
39 self._N = gamma_integration.get('N', 100)
41 if self.type not in {'sphere', 'reciprocal', '1BZ', 'WS'}:
42 raise TypeError('type in gamma_integration should be one of sphere'
43 ', reciprocal, 1BZ, or WS.')
45 if not self.is_numerical:
46 if gamma_integration.get('reduced', False):
47 raise TypeError('reduced key being True is only supported for '
48 'type reciprocal or 1BZ.')
50 def __repr__(self):
51 return f'type: {self.type} reduced: {self.reduced}'
53 @property
54 def is_analytical(self):
55 return self.type == 'sphere'
57 @property
58 def is_numerical(self):
59 return self.type in {'reciprocal', '1BZ'}
61 @property
62 def is_Wigner_Seitz(self):
63 return self.type == 'WS'
65 @property
66 def to_1bz(self):
67 return self.type == '1BZ'
69 @property
70 def N(self):
71 return self._N
74class QPointDescriptor(KPointDescriptor):
76 @staticmethod
77 def from_gs(gs):
78 kd, atoms = gs.kd, gs.atoms
79 # Find q-vectors and weights in the IBZ:
80 assert -1 not in kd.bz2bz_ks
81 offset_c = 0.5 * ((kd.N_c + 1) % 2) / kd.N_c
82 bzq_qc = monkhorst_pack(kd.N_c) + offset_c
83 qd = KPointDescriptor(bzq_qc)
84 qd.set_symmetry(atoms, kd.symmetry)
85 return qd
88def initialize_w_calculator(chi0calc, context, *,
89 coulomb,
90 xc='RPA', # G0W0Kernel arguments
91 mpa=None, E0=Ha, eta=None,
92 integrate_gamma=GammaIntegrationMode('sphere'),
93 q0_correction=False):
94 """Initialize a WCalculator from a Chi0Calculator.
96 Parameters
97 ----------
98 chi0calc : Chi0Calculator
99 xc : str
100 Kernel to use when including vertex corrections.
102 Remaining arguments: See WCalculator
103 """
104 from gpaw.response.g0w0_kernels import G0W0Kernel
106 gs = chi0calc.gs
107 qd = QPointDescriptor.from_gs(gs)
109 xckernel = G0W0Kernel(xc=xc, ecut=chi0calc.chi0_body_calc.ecut,
110 gs=gs, qd=qd,
111 context=context)
113 kwargs = dict()
114 if mpa:
115 wcalc_cls = MPACalculator
116 kwargs['mpa'] = mpa
117 else:
118 wcalc_cls = WCalculator
120 return wcalc_cls(gs, context, qd=qd,
121 coulomb=coulomb, xckernel=xckernel,
122 integrate_gamma=integrate_gamma, eta=eta,
123 q0_correction=q0_correction, **kwargs)
126class WBaseCalculator():
128 def __init__(self, gs, context, *, qd,
129 coulomb, xckernel,
130 integrate_gamma=GammaIntegrationMode('sphere'),
131 eta=None,
132 q0_correction=False):
133 """
134 Base class for W Calculator including basic initializations and Gamma
135 Gamma handling.
137 Parameters
138 ----------
139 gs : ResponseGroundStateAdapter
140 context : ResponseContext
141 qd : QPointDescriptor
142 coulomb : CoulombKernel
143 xckernel : G0W0Kernel
144 integrate_gamma: GammaIntegrationMode
145 q0_correction : bool
146 Analytic correction to the q=0 contribution applicable to 2D
147 systems.
148 """
149 self.gs = gs
150 self.context = context
151 self.qd = qd
152 self.coulomb = coulomb
153 self.xckernel = xckernel
154 self.integrate_gamma = integrate_gamma
155 self.eta = eta
157 if q0_correction:
158 assert self.coulomb.truncation == '2D'
159 self.q0_corrector = Q0Correction(
160 cell_cv=self.gs.gd.cell_cv,
161 bzk_kc=self.gs.kd.bzk_kc,
162 N_c=self.qd.N_c)
164 npts_c = self.q0_corrector.npts_c
165 self.context.print('Applying analytical 2D correction to W:',
166 flush=False)
167 self.context.print(' Evaluating Gamma point contribution to W '
168 + 'on a %dx%dx%d grid' % tuple(npts_c))
169 else:
170 self.q0_corrector = None
172 def get_V0sqrtV0(self, chi0):
173 """
174 Integrated Coulomb kernels.
175 """
176 V0 = None
177 sqrtV0 = None
178 if self.integrate_gamma.is_numerical:
179 reduced = self.integrate_gamma.reduced
180 tofirstbz = self.integrate_gamma.to_1bz
181 N = self.integrate_gamma.N
182 V0, sqrtV0 = self.coulomb.integrated_kernel(qpd=chi0.qpd,
183 reduced=reduced,
184 tofirstbz=tofirstbz,
185 N=N)
186 elif self.integrate_gamma.is_analytical:
187 if chi0.optical_limit:
188 # The volume of reciprocal cell occupied by a single q-point
189 bzvol = (2 * np.pi)**3 / self.gs.volume / self.qd.nbzkpts
190 # Radius of a sphere with a volume of the bzvol above
191 Rq0 = (3 * bzvol / (4 * np.pi))**(1. / 3.)
192 # Analytical integral of Coulomb interaction over the sphere
193 # defined above centered at q=0.
194 # V0 = 1/|bzvol| int_bzvol dq 4 pi / |q^2|
195 V0 = 16 * np.pi**2 * Rq0 / bzvol
196 # Analytical integral of square root of Coulomb interaction
197 # over the same sphere
198 # sqrtV0 = 1/|bzvol| int_bzvol dq sqrt(4 pi / |q^2|)
199 sqrtV0 = (4 * np.pi)**(1.5) * Rq0**2 / bzvol / 2
200 else:
201 raise KeyError('Unknown integrate_gamma option:'
202 f'{self.integrate_gamma}.')
203 return V0, sqrtV0
205 def apply_gamma_correction(self, W_GG, einv_GG, V0, sqrtV0, sqrtV_G):
206 """
207 Replacing q=0, (G,G')= (0,0), (0,:), (:,0) with corresponding
208 matrix elements calculated with an average of the (diverging)
209 Coulomb interaction.
210 XXX: Understand and document exact expressions
211 """
212 W_GG[0, 0] = einv_GG[0, 0] * V0
213 W_GG[0, 1:] = einv_GG[0, 1:] * sqrtV_G[1:] * sqrtV0
214 W_GG[1:, 0] = einv_GG[1:, 0] * sqrtV0 * sqrtV_G[1:]
217class WCalculator(WBaseCalculator):
218 def get_HW_model(self, chi0, fxc_mode, only_correlation=True):
219 assert only_correlation
220 W_wGG = self.calculate_W_WgG(chi0,
221 fxc_mode=fxc_mode,
222 only_correlation=True)
223 # HT used to calculate convulution between time-ordered G and W
224 hilbert_transform = GWHilbertTransforms(chi0.wd.omega_w, self.eta)
225 with self.context.timer('Hilbert'):
226 W_xwGG = hilbert_transform(W_wGG)
228 factor = 1.0 / (self.qd.nbzkpts * 2 * pi * self.gs.volume)
229 return FullFrequencyHWModel(chi0.wd, W_xwGG, factor)
231 def calculate_W_WgG(self, chi0,
232 fxc_mode='GW',
233 only_correlation=False):
234 """Calculate the screened interaction in W_wGG or W_WgG representation.
236 Additional Parameters
237 ----------
238 only_correlation: bool
239 if true calculate Wc otherwise calculate full W
240 out_dist: str
241 specifices output distribution of W array (wGG or WgG)
242 """
243 W_wGG = self.calculate_W_wGG(chi0, fxc_mode,
244 only_correlation=only_correlation)
246 W_WgG = chi0.body.blockdist.distribute_as(W_wGG, chi0.body.nw, 'WgG')
247 return W_WgG
249 def calculate_W_wGG(self, chi0, fxc_mode='GW',
250 only_correlation=False):
251 """In-place calculation of the screened interaction."""
252 dfc = DielectricFunctionCalculator(chi0, self.coulomb,
253 self.xckernel, fxc_mode)
254 self.context.timer.start('Dyson eq.')
256 if self.integrate_gamma.is_Wigner_Seitz:
257 from gpaw.hybrids.wstc import WignerSeitzTruncatedCoulomb
258 wstc = WignerSeitzTruncatedCoulomb(chi0.qpd.gd.cell_cv,
259 dfc.coulomb.N_c)
260 sqrtV_G = wstc.get_potential(chi0.qpd)**0.5
261 else:
262 sqrtV_G = dfc.sqrtV_G
263 V0, sqrtV0 = self.get_V0sqrtV0(chi0)
265 einv_wGG = dfc.get_epsinv_wGG(only_correlation=False)
266 W_wGG = np.empty_like(einv_wGG)
267 for iw, (einv_GG, W_GG) in enumerate(zip(einv_wGG, W_wGG)):
268 # If only_correlation = True function spits out
269 # W^c = sqrt(V)(epsinv - delta_GG')sqrt(V). However, full epsinv
270 # is still needed for q0_corrector.
271 einvt_GG = (einv_GG - dfc.I_GG) if only_correlation else einv_GG
272 W_GG[:] = einvt_GG * (sqrtV_G *
273 sqrtV_G[:, np.newaxis])
274 if self.q0_corrector is not None and chi0.optical_limit:
275 W = dfc.wblocks.a + iw
276 self.q0_corrector.add_q0_correction(chi0.qpd, W_GG,
277 einv_GG,
278 chi0.chi0_WxvG[W],
279 chi0.chi0_Wvv[W],
280 sqrtV_G)
281 elif (self.integrate_gamma.is_analytical and chi0.optical_limit) \
282 or self.integrate_gamma.is_numerical:
283 self.apply_gamma_correction(W_GG, einvt_GG,
284 V0, sqrtV0, dfc.sqrtV_G)
286 self.context.timer.stop('Dyson eq.')
287 return W_wGG
289 def dyson_and_W_new(self, iq, q_c, chi0, ecut, coulomb):
290 # assert not self.do_GW_too
291 assert ecut == chi0.qpd.ecut
292 assert self.fxc_mode == 'GW'
293 assert not np.allclose(q_c, 0)
295 nW = len(self.wd)
296 nG = chi0.qpd.ngmax
298 from gpaw.response.wgg import Grid
300 WGG = (nW, nG, nG)
301 WgG_grid = Grid(
302 comm=self.blockcomm,
303 shape=WGG,
304 cpugrid=(1, self.blockcomm.size, 1))
305 assert chi0.chi0_wGG.shape == WgG_grid.myshape
307 my_gslice = WgG_grid.myslice[1]
309 dielectric_WgG = chi0.chi0_wGG # XXX
310 for iw, chi0_GG in enumerate(chi0.chi0_wGG):
311 sqrtV_G = coulomb.sqrtV(chi0.qpd, q_v=None)
312 e_GG = np.eye(nG) - chi0_GG * sqrtV_G * sqrtV_G[:, np.newaxis]
313 e_gG = e_GG[my_gslice]
315 dielectric_WgG[iw, :, :] = e_gG
317 wgg_grid = Grid(comm=self.blockcomm, shape=WGG)
319 dielectric_wgg = wgg_grid.zeros(dtype=complex)
320 WgG_grid.redistribute(wgg_grid, dielectric_WgG, dielectric_wgg)
322 assert np.allclose(dielectric_wgg, dielectric_WgG)
324 wgg_grid.invert_inplace(dielectric_wgg)
326 wgg_grid.redistribute(WgG_grid, dielectric_wgg, dielectric_WgG)
327 inveps_WgG = dielectric_WgG
329 self.context.timer.start('Dyson eq.')
331 for iw, inveps_gG in enumerate(inveps_WgG):
332 inveps_gG -= np.identity(nG)[my_gslice]
333 thing_GG = sqrtV_G * sqrtV_G[:, np.newaxis]
334 inveps_gG *= thing_GG[my_gslice]
336 W_WgG = inveps_WgG
337 Wp_wGG = W_WgG.copy()
338 Wm_wGG = W_WgG.copy()
339 return chi0.qpd, Wm_wGG, Wp_wGG # not Hilbert transformed yet
342class HWModel:
343 """
344 Hilbert Transformed W Model.
345 """
347 def get_HW(self, omega, occ):
348 """
349 Get Hilbert transformed W at frequency omega.
351 occ: The occupation number for the orbital of the Greens function.
352 """
353 raise NotImplementedError
356class FullFrequencyHWModel(HWModel):
357 def __init__(self, wd, HW_swGG, factor):
358 self.wd = wd
359 self.HW_swGG = HW_swGG
360 self.factor = factor
362 def get_HW(self, omega, occ):
363 # For more information about how fsign and wsign works, see
364 # https://backend.orbit.dtu.dk/ws/portalfiles/portal/93075765/hueser_PhDthesis.pdf
365 # eq. 2.2 endind up to eq. 2.11
366 # Effectively, the symmetry of time ordered W is used,
367 # i.e. W(w) = -W(-w). To allow that data is only stored for w>=0.
368 # Hence, the interpolation happends always to the positive side, but
369 # the information of true w is keps tract using wsign.
370 # In addition, whether the orbital in question at G is occupied or
371 # unoccupied, which then again affects, which Hilbert transform of
372 # W is chosen, is kept track with fsign.
373 fsign = np.sign(2 * occ - 1)
374 o = abs(omega)
375 wsign = np.sign(omega + 1e-15)
376 wd = self.wd
377 # Pick +i*eta or -i*eta:
378 s = (1 + wsign * np.sign(-fsign)).astype(int) // 2
379 w = wd.get_floor_index(o, safe=False)
381 # Interpolation indexes w + 1, therefore - 2 here
382 if w > len(wd) - 2:
383 return None, None
385 o1 = wd.omega_w[w]
386 o2 = wd.omega_w[w + 1]
388 C1_GG = self.HW_swGG[s][w]
389 C2_GG = self.HW_swGG[s][w + 1]
390 p = self.factor * wsign
392 sigma_GG = ((o - o1) * C2_GG + (o2 - o) * C1_GG) / (o2 - o1)
393 dsigma_GG = wsign * (C2_GG - C1_GG) / (o2 - o1)
394 return -1j * p * sigma_GG, -1j * p * dsigma_GG
397class MPAHWModel(HWModel):
398 def __init__(self, W_nGG, omegat_nGG, eta, factor):
399 self.W_nGG = np.ascontiguousarray(W_nGG)
400 self.omegat_nGG = np.ascontiguousarray(omegat_nGG)
401 self.eta = eta
402 self.factor = factor
404 def get_HW(self, omega, occ):
405 x_GG = np.empty(self.omegat_nGG.shape[1:], dtype=complex)
406 dx_GG = np.empty(self.omegat_nGG.shape[1:], dtype=complex)
407 evaluate_mpa_poly(x_GG, dx_GG, omega, occ, self.omegat_nGG, self.W_nGG,
408 self.eta, self.factor)
410 return x_GG.conj(), dx_GG.conj() # Why do we have to do a conjugate
413class MPACalculator(WBaseCalculator):
414 def __init__(self, gs, context, *, eta, mpa, **kwargs):
415 super().__init__(gs, context, **kwargs)
416 self.eta = eta
417 self.mpa = mpa
419 def get_HW_model(self, chi0,
420 fxc_mode='GW'):
421 """Calculate the MPA parametrization of screened interaction.
422 """
424 dfc = DielectricFunctionCalculator(chi0,
425 self.coulomb,
426 self.xckernel,
427 fxc_mode)
429 self.context.timer.start('Dyson eq.')
430 einv_wGG = dfc.get_epsinv_wGG(only_correlation=True)
431 einv_WgG = chi0.body.blockdist.distribute_as(einv_wGG, chi0.nw, 'WgG')
433 solver = RESolver(chi0.wd.omega_w)
434 E_pGG, R_pGG = solver.solve(einv_WgG)
435 E_pGG -= 1j * self.eta # DALV: This is just to match the FF results
437 R_pGG = chi0.body.blockdist.distribute_as(R_pGG, self.mpa['npoles'],
438 'wGG')
439 E_pGG = chi0.body.blockdist.distribute_as(E_pGG,
440 self.mpa['npoles'], 'wGG')
442 if self.integrate_gamma.is_Wigner_Seitz:
443 from gpaw.hybrids.wstc import WignerSeitzTruncatedCoulomb
444 wstc = WignerSeitzTruncatedCoulomb(chi0.qpd.gd.cell_cv,
445 dfc.coulomb.N_c)
446 sqrtV_G = wstc.get_potential(chi0.qpd)**0.5
447 else:
448 sqrtV_G = dfc.sqrtV_G
450 W_pGG = pi * R_pGG * sqrtV_G[np.newaxis, :, np.newaxis] \
451 * sqrtV_G[np.newaxis, np.newaxis, :]
453 assert self.q0_corrector is None
454 if (self.integrate_gamma.is_analytical and chi0.optical_limit)\
455 or self.integrate_gamma.is_numerical:
456 V0, sqrtV0 = self.get_V0sqrtV0(chi0)
457 for W_GG, R_GG in zip(W_pGG, R_pGG):
458 self.apply_gamma_correction(W_GG, pi * R_GG,
459 V0, sqrtV0,
460 dfc.sqrtV_G)
462 W_pGG = np.transpose(W_pGG, axes=(0, 2, 1)) # Why the transpose
463 E_pGG = np.transpose(E_pGG, axes=(0, 2, 1))
465 W_pGG = chi0.body.blockdist.distribute_as(W_pGG, self.mpa['npoles'],
466 'WgG')
467 E_pGG = chi0.body.blockdist.distribute_as(E_pGG,
468 self.mpa['npoles'], 'WgG')
470 self.context.timer.stop('Dyson eq.')
472 factor = 1.0 / (self.qd.nbzkpts * 2 * pi * self.gs.volume)
473 return MPAHWModel(W_pGG, E_pGG, self.eta, factor)