Coverage for gpaw/response/g0w0.py: 95%
641 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1from __future__ import annotations
2import pickle
3import warnings
4from math import pi, isclose
5from pathlib import Path
6from collections.abc import Iterable
8import numpy as np
10from ase.parallel import paropen
11from ase.units import Ha
13from gpaw import GPAW, debug
14import gpaw.mpi as mpi
15from gpaw.hybrids.eigenvalues import non_self_consistent_eigenvalues
16from gpaw.pw.descriptor import (count_reciprocal_vectors, PWMapping)
17from gpaw.utilities.progressbar import ProgressBar
19from gpaw.response import ResponseContext, ResponseGroundStateAdapter
20from gpaw.response.chi0 import Chi0Calculator, get_frequency_descriptor
21from gpaw.response.pair import phase_shifted_fft_indices
22from gpaw.response.qpd import SingleQPWDescriptor
23from gpaw.response.pw_parallelization import Blocks1D
24from gpaw.response.screened_interaction import (initialize_w_calculator,
25 GammaIntegrationMode)
26from gpaw.response.coulomb_kernels import CoulombKernel
27from gpaw.response import timer
28from gpaw.response.mpa_sampling import mpa_frequency_sampling
29from gpaw.mpi import broadcast_exception
31from ase.utils.filecache import MultiFileJSONCache as FileCache
32from contextlib import ExitStack
33from ase.parallel import broadcast
36def compare_inputs(inp1, inp2, rel_tol=1e-14, abs_tol=1e-14):
37 """
38 Compares nested structures of dictionarys, lists, etc. and
39 makes sure the nested structure is the same, and also that all
40 floating points match within the given tolerances.
42 :params inp1: Structure 1 to compare.
43 :params inp2: Structure 2 to compare.
44 :params rel_tol: Maximum difference for being considered "close",
45 relative to the magnitude of the input values as defined by math.isclose().
46 :params abs_tol: Maximum difference for being considered "close",
47 regardless of the magnitude of the input values as defined by
48 math.isclose().
50 :returns: bool indicating if structures don't match (False) or do match
51 (True)
52 """
53 if isinstance(inp1, dict):
54 if inp1.keys() != inp2.keys():
55 return False
56 for key in inp1.keys() & inp2.keys():
57 val1 = inp1[key]
58 val2 = inp2[key]
59 if not compare_inputs(val1, val2,
60 rel_tol=rel_tol, abs_tol=abs_tol):
61 return False
62 elif isinstance(inp1, float):
63 if not isclose(inp1, inp2, rel_tol=rel_tol, abs_tol=abs_tol):
64 return False
65 elif not isinstance(inp1, str) and isinstance(inp1, Iterable):
66 if len(inp1) != len(inp2):
67 return False
68 for val1, val2 in zip(inp1, inp2):
69 if not compare_inputs(val1, val2,
70 rel_tol=rel_tol, abs_tol=abs_tol):
71 return False
72 else:
73 if inp1 != inp2:
74 return False
76 return True
79class Sigma:
80 def __init__(self, iq, q_c, fxc, esknshape, nw, **inputs):
81 """Inputs are used for cache invalidation, and are stored for each
82 file.
83 """
84 self.iq = iq
85 self.q_c = q_c
86 self.fxc = fxc
87 # We might as well allocate both from same array
88 # in order to add and communicate to them faster.
89 self._buf = np.zeros((2, *esknshape))
90 # self-energies and derivatives:
91 self.sigma_eskn, self.dsigma_eskn = self._buf
93 eskwnshape = (*esknshape[:3], nw, esknshape[3])
94 self.sigma_eskwn = np.zeros(eskwnshape, dtype=complex)
95 self.inputs = inputs
97 def sum(self, comm):
98 comm.sum(self._buf)
99 comm.sum(self.sigma_eskwn)
101 def __iadd__(self, other):
102 self.validate_inputs(other.inputs)
103 self._buf += other._buf
104 self.sigma_eskwn += other.sigma_eskwn
105 return self
107 def validate_inputs(self, inputs):
108 equals = compare_inputs(inputs, self.inputs, rel_tol=1e-12,
109 abs_tol=1e-12)
110 if not equals:
111 raise RuntimeError('There exists a cache with mismatching input '
112 f'parameters: {inputs} != {self.inputs}.')
114 @classmethod
115 def fromdict(cls, dct):
116 instance = cls(dct['iq'], dct['q_c'], dct['fxc'],
117 dct['sigma_eskn'].shape, dct['sigma_eskwn'].shape[3],
118 **dct['inputs'])
119 instance.sigma_eskn[:] = dct['sigma_eskn']
120 instance.dsigma_eskn[:] = dct['dsigma_eskn']
121 instance.sigma_eskwn[:] = dct['sigma_eskwn']
122 return instance
124 def todict(self):
125 return {'iq': self.iq,
126 'q_c': self.q_c,
127 'fxc': self.fxc,
128 'sigma_eskn': self.sigma_eskn,
129 'sigma_eskwn': self.sigma_eskwn,
130 'dsigma_eskn': self.dsigma_eskn,
131 'inputs': self.inputs}
134class G0W0Outputs:
135 def __init__(self, context, shape, ecut_e, sigma_eskn, dsigma_eskn,
136 sigma_eskwn, eps_skn, vxc_skn, exx_skn, f_skn):
137 self.extrapolate(context, shape, ecut_e, sigma_eskn, dsigma_eskn)
138 self.Z_skn = 1 / (1 - self.dsigma_skn)
140 # G0W0 single-step.
141 # If we want GW0 again, we need to grab the expressions
142 # from e.g. e73917fca5b9dc06c899f00b26a7c46e7d6fa749
143 # or earlier and use qp correctly.
144 self.qp_skn = eps_skn + self.Z_skn * (
145 -vxc_skn + exx_skn + self.sigma_skn)
147 self.sigma_eskn = sigma_eskn
148 self.dsigma_eskn = dsigma_eskn
150 self.eps_skn = eps_skn
151 self.vxc_skn = vxc_skn
152 self.exx_skn = exx_skn
153 self.f_skn = f_skn
154 self.sigma_eskwn = sigma_eskwn
156 def extrapolate(self, context, shape, ecut_e, sigma_eskn, dsigma_eskn):
157 if len(ecut_e) == 1:
158 self.sigma_skn = sigma_eskn[0]
159 self.dsigma_skn = dsigma_eskn[0]
160 self.sigr2_skn = None
161 self.dsigr2_skn = None
162 return
164 from scipy.stats import linregress
166 # Do linear fit of selfenergy vs. inverse of number of plane waves
167 # to extrapolate to infinite number of plane waves
169 context.print('', flush=False)
170 context.print('Extrapolating selfenergy to infinite energy cutoff:',
171 flush=False)
172 context.print(' Performing linear fit to %d points' % len(ecut_e))
173 self.sigr2_skn = np.zeros(shape)
174 self.dsigr2_skn = np.zeros(shape)
175 self.sigma_skn = np.zeros(shape)
176 self.dsigma_skn = np.zeros(shape)
177 invN_i = ecut_e**(-3. / 2)
178 for m in range(np.prod(shape)):
179 s, k, n = np.unravel_index(m, shape)
181 slope, intercept, r_value, p_value, std_err = \
182 linregress(invN_i, sigma_eskn[:, s, k, n])
184 self.sigr2_skn[s, k, n] = r_value**2
185 self.sigma_skn[s, k, n] = intercept
187 slope, intercept, r_value, p_value, std_err = \
188 linregress(invN_i, dsigma_eskn[:, s, k, n])
190 self.dsigr2_skn[s, k, n] = r_value**2
191 self.dsigma_skn[s, k, n] = intercept
193 if np.any(self.sigr2_skn < 0.9) or np.any(self.dsigr2_skn < 0.9):
194 context.print(' Warning: Bad quality of linear fit for some ('
195 'n,k). ', flush=False)
196 context.print(' Higher cutoff might be necessary.',
197 flush=False)
199 context.print(' Minimum R^2 = %1.4f. (R^2 Should be close to 1)' %
200 min(np.min(self.sigr2_skn), np.min(self.dsigr2_skn)))
202 def get_results_eV(self):
203 results = {
204 'f': self.f_skn,
205 'eps': self.eps_skn * Ha,
206 'vxc': self.vxc_skn * Ha,
207 'exx': self.exx_skn * Ha,
208 'sigma': self.sigma_skn * Ha,
209 'dsigma': self.dsigma_skn,
210 'Z': self.Z_skn,
211 'qp': self.qp_skn * Ha}
213 results.update(
214 sigma_eskn=self.sigma_eskn * Ha,
215 dsigma_eskn=self.dsigma_eskn,
216 sigma_eskwn=self.sigma_eskwn * Ha)
218 if self.sigr2_skn is not None:
219 assert self.dsigr2_skn is not None
220 results['sigr2_skn'] = self.sigr2_skn
221 results['dsigr2_skn'] = self.dsigr2_skn
223 return results
226class QSymmetryOp:
227 def __init__(self, symno, U_cc, sign):
228 self.symno = symno
229 self.U_cc = U_cc
230 self.sign = sign
232 def apply(self, q_c):
233 return self.sign * (self.U_cc @ q_c)
235 def check_q_Q_symmetry(self, Q_c, q_c):
236 d_c = self.apply(q_c) - Q_c
237 assert np.allclose(d_c.round(), d_c)
239 def get_M_vv(self, cell_cv):
240 # We'll be inverting these cells a lot.
241 # Should have an object with the cell and its inverse which does this.
242 return cell_cv.T @ self.U_cc.T @ np.linalg.inv(cell_cv).T
244 @classmethod
245 def get_symops(cls, qd, iq, q_c):
246 # Loop over all k-points in the BZ and find those that are
247 # related to the current IBZ k-point by symmetry
248 Q1 = qd.ibz2bz_k[iq]
249 done = set()
250 for Q2 in qd.bz2bz_ks[Q1]:
251 if Q2 >= 0 and Q2 not in done:
252 time_reversal = qd.time_reversal_k[Q2]
253 symno = qd.sym_k[Q2]
254 Q_c = qd.bzk_kc[Q2]
256 symop = cls(
257 symno=symno,
258 U_cc=qd.symmetry.op_scc[symno],
259 sign=1 - 2 * time_reversal)
261 symop.check_q_Q_symmetry(Q_c, q_c)
262 # Q_c, symop = QSymmetryOp.from_qd(qd, Q2, q_c)
263 yield Q_c, symop
264 done.add(Q2)
266 @classmethod
267 def get_symop_from_kpair(cls, kd, qd, kpt1, kpt2):
268 # from k-point pair kpt1, kpt2 get Q_c = kpt2-kpt1, corrsponding IBZ
269 # k-point q_c, indexes iQ, iq and symmetry transformation relating
270 # Q_c to q_c
271 Q_c = kd.bzk_kc[kpt2.K] - kd.bzk_kc[kpt1.K]
272 iQ = qd.where_is_q(Q_c, qd.bzk_kc)
273 iq = qd.bz2ibz_k[iQ]
274 q_c = qd.ibzk_kc[iq]
276 # Find symmetry that transforms Q_c into q_c
277 sym = qd.sym_k[iQ]
278 U_cc = qd.symmetry.op_scc[sym]
279 time_reversal = qd.time_reversal_k[iQ]
280 sign = 1 - 2 * time_reversal
281 symop = cls(sym, U_cc, sign)
282 symop.check_q_Q_symmetry(Q_c, q_c)
284 return symop, iq
286 def apply_symop_q(self, qpd, pawcorr, kpt1, kpt2):
287 # returns necessary quantities to get symmetry transformed
288 # density matrix
289 Q_G = phase_shifted_fft_indices(kpt1.k_c, kpt2.k_c, qpd,
290 coordinate_transformation=self.apply)
292 qG_Gv = qpd.get_reciprocal_vectors(add_q=True)
293 M_vv = self.get_M_vv(qpd.gd.cell_cv)
294 mypawcorr = pawcorr.remap_by_symop(self, qG_Gv, M_vv)
296 return mypawcorr, Q_G
299def get_nmG(kpt1, kpt2, mypawcorr, n, qpd, I_G, pair_calc, timer=None):
300 if timer:
301 timer.start('utcc and pawcorr multiply')
302 ut1cc_R = kpt1.ut_nR[n].conj()
303 C1_aGi = mypawcorr.multiply(kpt1.P_ani, band=n)
304 if timer:
305 timer.stop('utcc and pawcorr multiply')
306 n_mG = pair_calc.calculate_pair_density(
307 ut1cc_R, C1_aGi, kpt2, qpd, I_G)
308 return n_mG
311gw_logo = """\
312 ___ _ _ _
313 | || | | |
314 | | || | | |
315 |__ ||_____|
316 |___|
317"""
320def get_max_nblocks(world, calc, ecut):
321 nblocks = world.size
322 if not isinstance(calc, (str, Path)):
323 raise Exception('Using a calulator is not implemented at '
324 'the moment, load from file!')
325 # nblocks_calc = calc
326 else:
327 nblocks_calc = GPAW(calc)
328 ngmax = []
329 for q_c in nblocks_calc.wfs.kd.bzk_kc:
330 qpd = SingleQPWDescriptor.from_q(q_c, np.min(ecut) / Ha,
331 nblocks_calc.wfs.gd)
332 ngmax.append(qpd.ngmax)
333 nG = np.min(ngmax)
335 while nblocks > nG**0.5 + 1 or world.size % nblocks != 0:
336 nblocks -= 1
338 mynG = (nG + nblocks - 1) // nblocks
339 assert mynG * (nblocks - 1) < nG
340 return nblocks
343def get_frequencies(frequencies: dict | None,
344 domega0: float | None, omega2: float | None):
345 if domega0 is not None or omega2 is not None:
346 assert frequencies is None
347 frequencies = {'type': 'nonlinear',
348 'domega0': 0.025 if domega0 is None else domega0,
349 'omega2': 10.0 if omega2 is None else omega2}
350 warnings.warn(f'Please use frequencies={frequencies}')
351 elif frequencies is None:
352 frequencies = {'type': 'nonlinear',
353 'domega0': 0.025,
354 'omega2': 10.0}
355 else:
356 assert frequencies['type'] == 'nonlinear'
357 return frequencies
360def choose_ecut_things(ecut, ecut_extrapolation):
361 if ecut_extrapolation is True:
362 pct = 0.8
363 necuts = 3
364 ecut_e = ecut * (1 + (1. / pct - 1) * np.arange(necuts)[::-1] /
365 (necuts - 1))**(-2 / 3)
366 elif isinstance(ecut_extrapolation, (list, np.ndarray)):
367 ecut_e = np.array(np.sort(ecut_extrapolation))
368 if not np.allclose(ecut, ecut_e[-1]):
369 raise ValueError('ecut parameter must be the largest value'
370 'of ecut_extrapolation, when it is a list.')
371 ecut = ecut_e[-1]
372 else:
373 ecut_e = np.array([ecut])
374 return ecut, ecut_e
377def select_kpts(kpts, kd):
378 """Function to process input parameters that take a list of k-points given
379 in different format and returns a list of indices of the corresponding
380 k-points in the IBZ."""
382 if kpts is None:
383 # Do all k-points in the IBZ:
384 return np.arange(kd.nibzkpts)
386 if np.asarray(kpts).ndim == 1:
387 return kpts
389 # Find k-points:
390 bzk_Kc = kd.bzk_kc
391 indices = []
392 for k_c in kpts:
393 d_Kc = bzk_Kc - k_c
394 d_Kc -= d_Kc.round()
395 K = abs(d_Kc).sum(1).argmin()
396 if not np.allclose(d_Kc[K], 0):
397 raise ValueError('Could not find k-point: {k_c}'
398 .format(k_c=k_c))
399 k = kd.bz2ibz_k[K]
400 indices.append(k)
401 return indices
404class PairDistribution:
405 def __init__(self, kptpair_factory, blockcomm, mysKn1n2):
406 self.get_k_point = kptpair_factory.get_k_point
407 self.kd = kptpair_factory.gs.kd
408 self.blockcomm = blockcomm
409 self.mysKn1n2 = mysKn1n2
410 self.mykpts = [self.get_k_point(s, K, n1, n2)
411 for s, K, n1, n2 in self.mysKn1n2]
413 def kpt_pairs_by_q(self, q_c, m1, m2):
414 mykpts = self.mykpts
415 for u, kpt1 in enumerate(mykpts):
416 progress = u / len(mykpts)
417 K2 = self.kd.find_k_plus_q(q_c, [kpt1.K])[0]
418 kpt2 = self.get_k_point(kpt1.s, K2, m1, m2,
419 blockcomm=self.blockcomm)
421 yield progress, kpt1, kpt2
424def distribute_k_points_and_bands(chi0_body_calc, band1, band2, kpts=None):
425 """Distribute spins, k-points and bands.
427 The attribute self.mysKn1n2 will be set to a list of (s, K, n1, n2)
428 tuples that this process handles.
429 """
430 gs = chi0_body_calc.gs
431 blockcomm = chi0_body_calc.blockcomm
432 kncomm = chi0_body_calc.kncomm
434 if kpts is None:
435 kpts = np.arange(gs.kd.nbzkpts)
437 # nbands is the number of bands for each spin/k-point combination.
438 nbands = band2 - band1
439 size = kncomm.size
440 rank = kncomm.rank
441 ns = gs.nspins
442 nk = len(kpts)
443 n = (ns * nk * nbands + size - 1) // size
444 i1 = min(rank * n, ns * nk * nbands)
445 i2 = min(i1 + n, ns * nk * nbands)
447 mysKn1n2 = []
448 i = 0
449 for s in range(ns):
450 for K in kpts:
451 n1 = min(max(0, i1 - i), nbands)
452 n2 = min(max(0, i2 - i), nbands)
453 if n1 != n2:
454 mysKn1n2.append((s, K, n1 + band1, n2 + band1))
455 i += nbands
457 p = chi0_body_calc.context.print
458 p('BZ k-points:', gs.kd, flush=False)
459 p('Distributing spins, k-points and bands (%d x %d x %d)' %
460 (ns, nk, nbands), 'over %d process%s' %
461 (kncomm.size, ['es', ''][kncomm.size == 1]),
462 flush=False)
463 p('Number of blocks:', blockcomm.size)
465 return PairDistribution(
466 chi0_body_calc.kptpair_factory, blockcomm, mysKn1n2)
469class G0W0Calculator:
470 def __init__(self, filename='gw', *,
471 wd,
472 chi0calc,
473 wcalc,
474 kpts, bands, nbands=None,
475 fxc_modes,
476 eta,
477 ecut_e,
478 frequencies=None,
479 exx_vxc_calculator,
480 qcache,
481 ppa=False,
482 mpa=None,
483 evaluate_sigma=None):
484 """G0W0 calculator, initialized through G0W0 object.
486 The G0W0 calculator is used to calculate the quasi
487 particle energies through the G0W0 approximation for a number
488 of states.
490 Parameters
491 ----------
492 filename: str
493 Base filename of output files.
494 wcalc: WCalculator object
495 Defines the calculator for computing the screened interaction
496 kpts: list
497 List of indices of the IBZ k-points to calculate the quasi particle
498 energies for.
499 bands:
500 Range of band indices, like (n1, n2), to calculate the quasi
501 particle energies for. Bands n where n1<=n<n2 will be
502 calculated. Note that the second band index is not included.
503 frequencies:
504 Input parameters for frequency_grid.
505 Can be array of frequencies to evaluate the response function at
506 or dictionary of parameters for build-in nonlinear grid
507 (see :ref:`frequency grid`).
508 ecut_e: array(float)
509 Plane wave cut-off energies in eV. Defined with choose_ecut_things
510 nbands: int
511 Number of bands to use in the calculation. If None, the number will
512 be determined from :ecut: to yield a number close to the number of
513 plane waves used.
514 do_GW_too: bool
515 When carrying out a calculation including vertex corrections, it
516 is possible to get the standard GW results at the same time
517 (almost for free).
518 ppa: bool
519 Use Godby-Needs plasmon-pole approximation for screened interaction
520 and self-energy (reformulated as mpa with npoles = 1)
521 mpa: dict
522 Use multipole approximation for screened interaction
523 and self-energy [PRB 104, 115157 (2021)]
524 This method uses a sampling along one or two lines in the complex
525 frequency plane.
527 MPA parameters
528 ----------
529 npoles: Number of poles (positive integer generally lower than 15)
530 parallel_lines: How many (1-2) parallel lines to the real frequency
531 axis the sampling has.
532 wrange: Real interval defining the range of energy along the real
533 frequency axis.
534 alpha: exponent of the power distribution of points along the real
535 frequency axis [PRB 107, 155130 (2023)]
536 varpi: Distance of the second line to the real axis.
537 eta0: Imaginary part of the first point of the first line.
538 eta_rest: Imaginary part of the rest of the points of the first
539 line.
540 evaluate_sigma: array(float)
541 List of frequencies (in eV), where to evaluate the frequency
542 dependent self energy for each k-point and band involved in the
543 sigma-evaluation. This will be done in addition to evaluating the
544 normal self-energy quasiparticle matrix elements in G0W0
545 approximation.
546 """
547 self.chi0calc = chi0calc
548 self.wcalc = wcalc
549 self.context = self.wcalc.context
550 self.ppa = ppa
551 self.mpa = mpa
552 if evaluate_sigma is None:
553 evaluate_sigma = np.array([])
554 self.evaluate_sigma = evaluate_sigma
555 self.qcache = qcache
557 # Note: self.wd should be our only representation of the frequencies.
558 # We should therefore get rid of self.frequencies.
559 # It is currently only used by the restart code,
560 # so should be easy to remove after some further adaptation.
561 self.wd = wd
562 self.frequencies = frequencies
564 self.ecut_e = ecut_e / Ha
566 self.context.print(gw_logo)
568 if self.chi0calc.gs.metallic:
569 self.context.print('WARNING: \n'
570 'The current GW implementation cannot'
571 ' handle intraband screening. \n'
572 'This results in poor k-point'
573 ' convergence for metals')
575 self.fxc_modes = fxc_modes
577 if self.fxc_modes[0] != 'GW':
578 assert self.wcalc.xckernel.xc != 'RPA'
580 if len(self.fxc_modes) == 2:
581 # With multiple fxc_modes, we previously could do only
582 # GW plus one other fxc_mode. Now we can have any set
583 # of modes, but whether things are consistent or not may
584 # depend on how wcalc is configured.
585 assert 'GW' in self.fxc_modes
586 assert self.wcalc.xckernel.xc != 'RPA'
588 self.filename = filename
589 self.eta = eta / Ha
591 self.kpts = kpts
592 self.bands = bands
594 b1, b2 = self.bands
595 self.shape = (self.wcalc.gs.nspins, len(self.kpts), b2 - b1)
597 self.nbands = nbands
599 if self.wcalc.gs.nspins != 1:
600 for fxc_mode in self.fxc_modes:
601 if fxc_mode != 'GW':
602 raise RuntimeError('Including a xc kernel does not '
603 'currently work for spin-polarized '
604 f'systems. Invalid fxc_mode {fxc_mode}.'
605 )
607 self.pair_distribution = distribute_k_points_and_bands(
608 self.chi0calc.chi0_body_calc, b1, b2,
609 self.chi0calc.gs.kd.ibz2bz_k[self.kpts])
611 self.print_parameters(kpts, b1, b2)
613 self.exx_vxc_calculator = exx_vxc_calculator
615 p = self.context.print
616 if self.ppa:
617 p('Using Godby-Needs plasmon-pole approximation:')
618 p(' Fitting energy: i*E0, E0 = '
619 f'{self.wd.omega_w[1].imag:.3f} Hartree')
620 elif self.mpa:
621 omega_w = self.chi0calc.wd.omega_w
622 p('Using multipole approximation:')
623 p(f' Number of poles: {len(omega_w) // 2}')
624 p(f' Energy range: Re(E[-1]) = {omega_w[-1].real:.3f} Hartree')
625 p(' Imaginary range: Im(E[-1]) = '
626 f'{self.wd.omega_w[-1].imag:.3f} Hartree')
627 p(' Imaginary shift: Im(E[1]) = '
628 f'{self.wd.omega_w[1].imag:.3f} Hartree')
629 p(' Imaginary Origin shift: Im(E[0])'
630 f'= {self.wd.omega_w[0].imag:.3f} Hartree')
631 else:
632 self.context.print('Using full-frequency real axis integration')
634 def print_parameters(self, kpts, b1, b2):
635 isl = ['',
636 'Quasi particle states:']
637 if kpts is None:
638 isl.append('All k-points in IBZ')
639 else:
640 kptstxt = ', '.join([f'{k:d}' for k in self.kpts])
641 isl.append(f'k-points (IBZ indices): [{kptstxt}]')
642 isl.extend([f'Band range: ({b1:d}, {b2:d})',
643 '',
644 'Computational parameters:'])
645 if len(self.ecut_e) == 1:
646 isl.append(
647 'Plane wave cut-off: '
648 f'{self.chi0calc.chi0_body_calc.ecut * Ha:g} eV')
649 else:
650 assert len(self.ecut_e) > 1
651 isl.append('Extrapolating to infinite plane wave cut-off using '
652 'points at:')
653 for ec in self.ecut_e:
654 isl.append(f' {ec * Ha:.3f} eV')
655 isl.extend([f'Number of bands: {self.nbands:d}',
656 f'Coulomb cutoff: {self.wcalc.coulomb.truncation}',
657 f'Broadening: {self.eta * Ha:g} eV',
658 '',
659 f'fxc modes: {", ".join(sorted(self.fxc_modes))}',
660 f'Kernel: {self.wcalc.xckernel.xc}'])
661 self.context.print('\n'.join(isl))
663 def get_eps_and_occs(self):
664 eps_skn = np.empty(self.shape) # KS-eigenvalues
665 f_skn = np.empty(self.shape) # occupation numbers
667 nspins = self.wcalc.gs.nspins
668 b1, b2 = self.bands
669 for i, k in enumerate(self.kpts):
670 for s in range(nspins):
671 u = s + k * nspins
672 kpt = self.wcalc.gs.kpt_u[u]
673 eps_skn[s, i] = kpt.eps_n[b1:b2]
674 f_skn[s, i] = kpt.f_n[b1:b2] / kpt.weight
676 return eps_skn, f_skn
678 @timer('G0W0')
679 def calculate(self, qpoints=None):
680 """Starts the G0W0 calculation.
682 qpoints: list[int]
683 Set of q-points to calculate.
685 Returns a dict with the results with the following key/value pairs:
687 =========== =============================================
688 key value
689 =========== =============================================
690 ``f`` Occupation numbers
691 ``eps`` Kohn-Sham eigenvalues in eV
692 ``vxc`` Exchange-correlation
693 contributions in eV
694 ``exx`` Exact exchange contributions in eV
695 ``sigma`` Self-energy contributions in eV
696 ``dsigma`` Self-energy derivatives
697 ``sigma_e`` Self-energy contributions in eV
698 used for ecut extrapolation
699 ``Z`` Renormalization factors
700 ``qp`` Quasi particle (QP) energies in eV
701 ``iqp`` GW0/GW: QP energies for each iteration in eV
702 =========== =============================================
704 All the values are ``ndarray``'s of shape
705 (spins, IBZ k-points, bands)."""
707 qpoints = set(qpoints) if qpoints else None
709 if qpoints is None:
710 self.context.print('Summing all q:')
711 else:
712 qpt_str = ' '.join(map(str, qpoints))
713 self.context.print(f'Calculating following q-points: {qpt_str}')
714 self.calculate_q_points(qpoints=qpoints)
715 if qpoints is not None:
716 return f'A partial result of q-points: {qpt_str}'
717 sigmas = self.read_sigmas()
718 self.all_results = self.postprocess(sigmas)
719 # Note: self.results is a pointer pointing to one of the results,
720 # for historical reasons.
722 self.savepckl()
723 return self.results
725 def postprocess(self, sigmas):
726 all_results = {}
727 for fxc_mode, sigma in sigmas.items():
728 all_results[fxc_mode] = self.postprocess_single(fxc_mode, sigma)
730 self.print_results(all_results)
731 return all_results
733 def read_sigmas(self):
734 if self.context.comm.rank == 0:
735 sigmas = self._read_sigmas()
736 else:
737 sigmas = None
739 return broadcast(sigmas, comm=self.context.comm)
741 def _read_sigmas(self):
742 assert self.context.comm.rank == 0
744 # Integrate over all q-points, and accumulate the quasiparticle shifts
745 for iq, q_c in enumerate(self.wcalc.qd.ibzk_kc):
746 key = str(iq)
748 sigmas_contrib = self.get_sigmas_dict(key)
750 if iq == 0:
751 sigmas = sigmas_contrib
752 else:
753 for fxc_mode in self.fxc_modes:
754 sigmas[fxc_mode] += sigmas_contrib[fxc_mode]
756 return sigmas
758 def get_sigmas_dict(self, key):
759 assert self.context.comm.rank == 0
760 return {fxc_mode: Sigma.fromdict(sigma)
761 for fxc_mode, sigma in self.qcache[key].items()}
763 def postprocess_single(self, fxc_name, sigma):
764 output = self.calculate_g0w0_outputs(sigma)
765 return output.get_results_eV()
767 def savepckl(self):
768 """Save outputs to pckl files and return paths to those files."""
769 # Note: this is always called, but the paths aren't returned
770 # to the caller. Calling it again then overwrites the files.
771 #
772 # TODO:
773 # * Replace with JSON
774 # * Save to different files or same file?
775 # * Move this functionality to g0w0 result object
776 paths = {}
777 for fxc_mode in self.fxc_modes:
778 path = Path(f'{self.filename}_results_{fxc_mode}.pckl')
779 with paropen(path, 'wb', comm=self.context.comm) as fd:
780 pickle.dump(self.all_results[fxc_mode], fd, 2)
781 paths[fxc_mode] = path
783 # Do not return paths to caller before we know they all exist:
784 self.context.comm.barrier()
785 return paths
787 @property
788 def nqpts(self):
789 """Returns the number of q-points in the system."""
790 return len(self.wcalc.qd.ibzk_kc)
792 @timer('evaluate sigma')
793 def calculate_q(self, ie, k, kpt1, kpt2, qpd, Wdict,
794 *, symop, sigmas, blocks1d, pawcorr):
795 """Calculates the contribution to the self-energy and its derivative
796 for a given set of k-points, kpt1 and kpt2."""
797 mypawcorr, I_G = symop.apply_symop_q(qpd, pawcorr, kpt1, kpt2)
798 if debug:
799 N_c = qpd.gd.N_c
800 i_cG = symop.apply(np.unravel_index(qpd.Q_qG[0], N_c))
801 bzk_kc = self.wcalc.gs.kd.bzk_kc
802 Q_c = bzk_kc[kpt2.K] - bzk_kc[kpt1.K]
803 shift0_c = Q_c - symop.apply(qpd.q_c)
804 self.check(ie, i_cG, shift0_c, N_c, Q_c, mypawcorr)
806 for n in range(kpt1.n2 - kpt1.n1):
807 eps1 = kpt1.eps_n[n]
808 self.context.timer.start('get_nmG')
809 n_mG = get_nmG(kpt1, kpt2, mypawcorr,
810 n, qpd, I_G, self.chi0calc.pair_calc)
811 self.context.timer.stop('get_nmG')
813 if symop.sign == 1:
814 n_mG = n_mG.conj()
816 f_m = kpt2.f_n
817 deps_m = eps1 - kpt2.eps_n
819 nn = kpt1.n1 + n - self.bands[0]
821 assert set(Wdict) == set(sigmas)
823 for fxc_mode in self.fxc_modes:
824 sigma = sigmas[fxc_mode]
825 Wmodel = Wdict[fxc_mode]
827 # m is band index of all (both unoccupied and occupied) wave
828 # functions in G
829 for m, (deps, f, n_G) in enumerate(zip(deps_m, f_m, n_mG)):
830 # 2 * f - 1 will be used to select the branch of Hilbert
831 # transform, see get_HW of screened_interaction.py
832 # at FullFrequencyHWModel class.
834 nc_G = n_G.conj()
835 myn_G = n_G[blocks1d.myslice]
837 if self.evaluate_sigma is not None:
838 for w, omega in enumerate(self.evaluate_sigma):
839 S_GG, _ = Wmodel.get_HW(deps - eps1 + omega, f)
840 if S_GG is None:
841 continue
842 # print(myn_G.shape, S_GG.shape, nc_G.shape)
843 sigma.sigma_eskwn[ie, kpt1.s, k, w, nn] += \
844 myn_G @ S_GG @ nc_G
846 self.context.timer.start('Wmodel.get_HW')
847 S_GG, dSdw_GG = Wmodel.get_HW(deps, f)
848 self.context.timer.stop('Wmodel.get_HW')
849 if S_GG is None:
850 continue
852 # ie: ecut index for extrapolation
853 # kpt1.s: spin index of *
854 # k: k-point index of *
855 # nn: band index of *
856 # * wave function, where the sigma expectation value is
857 # evaluated
858 slot = ie, kpt1.s, k, nn
859 self.context.timer.start('n_G @ S_GG @ n_G')
860 sigma.sigma_eskn[slot] += (myn_G @ S_GG @ nc_G).real
861 sigma.dsigma_eskn[slot] += (myn_G @ dSdw_GG @ nc_G).real
862 self.context.timer.stop('n_G @ S_GG @ n_G')
864 def check(self, ie, i_cG, shift0_c, N_c, Q_c, pawcorr):
865 # Can we delete this check? XXX
866 assert np.allclose(shift0_c.round(), shift0_c)
867 shift0_c = shift0_c.round().astype(int)
868 I0_G = np.ravel_multi_index(i_cG - shift0_c[:, None], N_c, 'wrap')
869 qpd = SingleQPWDescriptor.from_q(Q_c, self.ecut_e[ie],
870 self.wcalc.gs.gd)
871 G_I = np.empty(N_c.prod(), int)
872 G_I[:] = -1
873 I1_G = qpd.Q_qG[0]
874 G_I[I1_G] = np.arange(len(I0_G))
875 G_G = G_I[I0_G]
876 # This indexing magic should definitely be moved to a method.
877 # What on earth is it really?
879 assert len(I0_G) == len(I1_G)
880 assert (G_G >= 0).all()
881 pairden_paw_corr = self.wcalc.gs.pair_density_paw_corrections
882 pawcorr_wcalc1 = pairden_paw_corr(qpd)
883 assert pawcorr.almost_equal(pawcorr_wcalc1, G_G)
885 def calculate_q_points(self, qpoints):
886 """Main loop over irreducible Brillouin zone points.
887 Handles restarts of individual qpoints using FileCache from ASE,
888 and subsequently calls calculate_q."""
890 pb = ProgressBar(self.context.fd)
892 self.context.timer.start('W')
893 self.context.print('\nCalculating screened Coulomb potential')
894 self.context.print(self.wcalc.coulomb.description())
896 chi0calc = self.chi0calc
897 self.context.print(self.wd)
899 # Find maximum size of chi-0 matrices:
900 nGmax = max(count_reciprocal_vectors(chi0calc.chi0_body_calc.ecut,
901 self.wcalc.gs.gd, q_c)
902 for q_c in self.wcalc.qd.ibzk_kc)
903 nw = len(self.wd)
905 size = self.chi0calc.chi0_body_calc.integrator.blockcomm.size
907 mynGmax = (nGmax + size - 1) // size
908 mynw = (nw + size - 1) // size
910 # some memory sizes...
911 if self.context.comm.rank == 0:
912 siz = (nw * mynGmax * nGmax +
913 max(mynw * nGmax, nw * mynGmax) * nGmax) * 16
914 sizA = (nw * nGmax * nGmax + nw * nGmax * nGmax) * 16
915 self.context.print(
916 ' memory estimate for chi0: local=%.2f MB, global=%.2f MB'
917 % (siz / 1024**2, sizA / 1024**2))
919 if self.context.comm.rank == 0 and qpoints is None:
920 self.context.print('Removing empty qpoint cache files...')
921 self.qcache.strip_empties()
923 self.context.comm.barrier()
925 # Need to pause the timer in between iterations
926 self.context.timer.stop('W')
928 with broadcast_exception(self.context.comm):
929 if self.context.comm.rank == 0:
930 for key, sigmas in self.qcache.items():
931 if qpoints and int(key) not in qpoints:
932 continue
933 sigmas = {fxc_mode: Sigma.fromdict(sigma)
934 for fxc_mode, sigma in sigmas.items()}
935 for fxc_mode, sigma in sigmas.items():
936 sigma.validate_inputs(self.get_validation_inputs())
938 for iq, q_c in enumerate(self.wcalc.qd.ibzk_kc):
939 # If a list of q-points is specified,
940 # skip the q-points not in the list
941 if qpoints and (iq not in qpoints):
942 continue
943 with ExitStack() as stack:
944 if self.context.comm.rank == 0:
945 qhandle = stack.enter_context(self.qcache.lock(str(iq)))
946 skip = qhandle is None
947 else:
948 skip = False
950 skip = broadcast(skip, comm=self.context.comm)
952 if skip:
953 continue
955 result = self.calculate_q_point(iq, q_c, pb, chi0calc)
957 if self.context.comm.rank == 0:
958 qhandle.save(result)
959 pb.finish()
961 def calculate_q_point(self, iq, q_c, pb, chi0calc):
962 # Reset calculation
963 sigmashape = (len(self.ecut_e), *self.shape)
964 sigmas = {fxc_mode: Sigma(iq, q_c, fxc_mode, sigmashape,
965 len(self.evaluate_sigma),
966 **self.get_validation_inputs())
967 for fxc_mode in self.fxc_modes}
969 chi0 = chi0calc.create_chi0(q_c)
971 m1 = chi0calc.gs.nocc1
972 for ie, ecut in enumerate(self.ecut_e):
973 self.context.timer.start('W')
975 # First time calculation
976 if ecut == chi0.qpd.ecut:
977 # Nothing to cut away:
978 m2 = self.nbands
979 else:
980 m2 = int(self.wcalc.gs.volume * ecut**1.5
981 * 2**0.5 / 3 / pi**2)
982 if m2 > self.nbands:
983 raise ValueError(f'Trying to extrapolate ecut to'
984 f'larger number of bands ({m2})'
985 f' than there are bands '
986 f'({self.nbands}).')
987 qpdi, Wdict, blocks1d, pawcorr = self.calculate_w(
988 chi0calc, q_c, chi0,
989 m1, m2, ecut, iq)
990 m1 = m2
992 self.context.timer.stop('W')
994 for nQ, (bzq_c, symop) in enumerate(QSymmetryOp.get_symops(
995 self.wcalc.qd, iq, q_c)):
997 for (progress, kpt1, kpt2)\
998 in self.pair_distribution.kpt_pairs_by_q(bzq_c, 0, m2):
999 pb.update((nQ + progress) / self.wcalc.qd.mynk)
1001 k1 = self.wcalc.gs.kd.bz2ibz_k[kpt1.K]
1002 i = self.kpts.index(k1)
1003 self.calculate_q(ie, i, kpt1, kpt2, qpdi, Wdict,
1004 symop=symop,
1005 sigmas=sigmas,
1006 blocks1d=blocks1d,
1007 pawcorr=pawcorr)
1009 for sigma in sigmas.values():
1010 sigma.sum(self.context.comm)
1012 return sigmas
1014 def get_validation_inputs(self):
1015 return {'kpts': self.kpts,
1016 'bands': list(self.bands),
1017 'nbands': self.nbands,
1018 'ecut_e': list(self.ecut_e),
1019 'frequencies': self.frequencies,
1020 'fxc_modes': self.fxc_modes,
1021 'integrate_gamma': repr(self.wcalc.integrate_gamma)}
1023 @timer('calculate_w')
1024 def calculate_w(self, chi0calc, q_c, chi0,
1025 m1, m2, ecut,
1026 iq):
1027 """Calculates the screened potential for a specified q-point."""
1029 chi0calc.chi0_body_calc.print_info(chi0.qpd)
1030 chi0calc.update_chi0(chi0, m1=m1, m2=m2,
1031 spins=range(self.wcalc.gs.nspins))
1033 Wdict = {}
1035 for fxc_mode in self.fxc_modes:
1036 rqpd = chi0.qpd.copy_with(ecut=ecut) # reduced qpd
1037 rchi0 = chi0.copy_with_reduced_pd(rqpd)
1038 Wdict[fxc_mode] = self.wcalc.get_HW_model(rchi0,
1039 fxc_mode=fxc_mode)
1040 if (chi0calc.chi0_body_calc.pawcorr is not None and
1041 rqpd.ecut < chi0.qpd.ecut):
1042 pw_map = PWMapping(rqpd, chi0.qpd)
1044 """This is extremely bad behaviour! G0W0Calculator
1045 should not change properties on the
1046 Chi0BodyCalculator! Change in the future! XXX"""
1047 chi0calc.chi0_body_calc.pawcorr = \
1048 chi0calc.chi0_body_calc.pawcorr.reduce_ecut(pw_map.G2_G1)
1050 # Create a blocks1d for the reduced plane-wave description
1051 blocks1d = Blocks1D(chi0.body.blockdist.blockcomm, rqpd.ngmax)
1053 return rqpd, Wdict, blocks1d, chi0calc.chi0_body_calc.pawcorr
1055 @timer('calculate_vxc_and_exx')
1056 def calculate_vxc_and_exx(self):
1057 return self.exx_vxc_calculator.calculate(
1058 n1=self.bands[0], n2=self.bands[1],
1059 kpt_indices=self.kpts)
1061 def print_results(self, results):
1062 description = ['f: Occupation numbers',
1063 'eps: KS-eigenvalues [eV]',
1064 'vxc: KS vxc [eV]',
1065 'exx: Exact exchange [eV]',
1066 'sigma: Self-energies [eV]',
1067 'dsigma: Self-energy derivatives',
1068 'Z: Renormalization factors',
1069 'qp: QP-energies [eV]']
1071 self.context.print('\nResults:')
1072 for line in description:
1073 self.context.print(line)
1075 b1, b2 = self.bands
1076 names = [line.split(':', 1)[0] for line in description]
1077 ibzk_kc = self.wcalc.gs.kd.ibzk_kc
1078 for s in range(self.wcalc.gs.nspins):
1079 for i, ik in enumerate(self.kpts):
1080 self.context.print(
1081 '\nk-point ' + '{} ({}): ({:.3f}, {:.3f}, '
1082 '{:.3f})'.format(i, ik, *ibzk_kc[ik]) +
1083 ' ' + self.fxc_modes[0])
1084 self.context.print('band' + ''.join(f'{name:>8}'
1085 for name in names))
1087 def actually_print_results(resultset):
1088 for n in range(b2 - b1):
1089 self.context.print(
1090 f'{n + b1:4}' +
1091 ''.join('{:8.3f}'.format(
1092 resultset[name][s, i, n]) for name in names))
1094 for fxc_mode in results:
1095 self.context.print(fxc_mode.rjust(69))
1096 actually_print_results(results[fxc_mode])
1098 self.context.write_timer()
1100 def calculate_g0w0_outputs(self, sigma):
1101 eps_skn, f_skn = self.get_eps_and_occs()
1102 vxc_skn, exx_skn = self.calculate_vxc_and_exx()
1103 kwargs = dict(
1104 context=self.context,
1105 shape=self.shape,
1106 ecut_e=self.ecut_e,
1107 eps_skn=eps_skn,
1108 vxc_skn=vxc_skn,
1109 exx_skn=exx_skn,
1110 f_skn=f_skn)
1112 return G0W0Outputs(sigma_eskn=sigma.sigma_eskn,
1113 dsigma_eskn=sigma.dsigma_eskn,
1114 sigma_eskwn=sigma.sigma_eskwn,
1115 **kwargs)
1118def choose_bands(bands, relbands, nvalence, nocc):
1119 if bands is not None and relbands is not None:
1120 raise ValueError('Use bands or relbands!')
1122 if relbands is not None:
1123 bands = [nvalence // 2 + b for b in relbands]
1125 if bands is None:
1126 bands = [0, nocc]
1128 return bands
1131class G0W0(G0W0Calculator):
1132 def __init__(self, calc, filename='gw',
1133 ecut=150.0,
1134 ecut_extrapolation=False,
1135 xc='RPA',
1136 ppa=False,
1137 mpa=None,
1138 E0=Ha,
1139 eta=0.1,
1140 nbands=None,
1141 bands=None,
1142 relbands=None,
1143 frequencies=None,
1144 domega0=None, # deprecated
1145 omega2=None, # deprecated
1146 nblocks=1,
1147 nblocksmax=False,
1148 kpts=None,
1149 world=mpi.world,
1150 timer=None,
1151 fxc_mode='GW',
1152 fxc_modes=None,
1153 truncation=None,
1154 integrate_gamma='sphere',
1155 q0_correction=False,
1156 do_GW_too=False,
1157 output_prefix=None,
1158 **kwargs):
1159 """G0W0 calculator wrapper.
1161 The G0W0 calculator is used to calculate the quasi
1162 particle energies through the G0W0 approximation for a number
1163 of states.
1165 Parameters
1166 ----------
1167 calc:
1168 Filename of saved calculator object.
1169 filename: str
1170 Base filename (a prefix) of output files.
1171 kpts: list
1172 List of indices of the IBZ k-points to calculate the quasi particle
1173 energies for.
1174 bands:
1175 Range of band indices, like (n1, n2), to calculate the quasi
1176 particle energies for. Bands n where n1<=n<n2 will be
1177 calculated. Note that the second band index is not included.
1178 relbands:
1179 Same as *bands* except that the numbers are relative to the
1180 number of occupied bands.
1181 E.g. (-1, 1) will use HOMO+LUMO.
1182 frequencies:
1183 Input parameters for the nonlinear frequency descriptor.
1184 ecut: float
1185 Plane wave cut-off energy in eV.
1186 ecut_extrapolation: bool or list
1187 If set to True an automatic extrapolation of the selfenergy to
1188 infinite cutoff will be performed based on three points
1189 for the cutoff energy.
1190 If an array is given, the extrapolation will be performed based on
1191 the cutoff energies given in the array.
1192 nbands: int
1193 Number of bands to use in the calculation. If None, the number will
1194 be determined from :ecut: to yield a number close to the number of
1195 plane waves used.
1196 ppa: bool
1197 Sets whether the Godby-Needs plasmon-pole approximation for the
1198 dielectric function should be used.
1199 mpa: dict
1200 Sets whether the multipole approximation for the response
1201 function should be used.
1202 xc: str
1203 Kernel to use when including vertex corrections.
1204 fxc_mode: str
1205 Where to include the vertex corrections; polarizability and/or
1206 self-energy. 'GWP': Polarizability only, 'GWS': Self-energy only,
1207 'GWG': Both.
1208 do_GW_too: bool
1209 When carrying out a calculation including vertex corrections, it
1210 is possible to get the standard GW results at the same time
1211 (almost for free).
1212 truncation: str
1213 Coulomb truncation scheme. Can be either 2D, 1D, or 0D.
1214 integrate_gamma: str or dict
1215 Method to integrate the Coulomb interaction.
1217 The default is 'sphere'. If 'reduced' key is not given,
1218 it defaults to False.
1220 {'type': 'sphere'} or 'sphere':
1221 Analytical integration of q=0, G=0 1/q^2 integrand in a sphere
1222 matching the volume of a single q-point.
1223 Used to be integrate_gamma=0.
1225 {'type': 'reciprocal'} or 'reciprocal':
1226 Numerical integration of q=0, G=0 1/q^2 integral in a volume
1227 resembling the reciprocal cell (parallelpiped).
1228 Used to be integrate_gamma=1.
1230 {'type': 'reciprocal', 'reduced':True} or 'reciprocal2D':
1231 Numerical integration of q=0, G=0 1/q^2 integral in a area
1232 resembling the reciprocal 2D cell (parallelogram) to be used
1233 to be usedwith 2D systems.
1234 Used to be integrate_gamma=2.
1236 {'type': '1BZ'} or '1BZ':
1237 Numerical integration of q=0, G=0 1/q^2 integral in a volume
1238 resembling the Wigner-Seitz cell of the reciprocal lattice
1239 (voronoi). More accurate than 'reciprocal'.
1241 A. Guandalini, P. D’Amico, A. Ferretti and D. Varsano:
1242 npj Computational Materials volume 9, Article number: 44 (2023)
1244 {'type': '1BZ', 'reduced': True} or '1BZ2D':
1245 Same as above, but everything is done in 2D (for 2D systems).
1247 {'type': 'WS'} or 'WS':
1248 The most accurate method to use for bulk systems.
1249 Instead of numerically integrating only q=0, G=0, all (q,G)-
1250 pairs participate to the truncation, which is done in real
1251 space utilizing the Wigner-Seitz truncation in the
1252 Born-von-Karmann supercell of the system.
1254 Numerical integration of q=0, G=0 1/q^2 integral in a volume
1255 resembling the Wigner-Seitz cell of the reciprocal lattice
1256 (Voronoi). More accurate than 'reciprocal'.
1258 R. Sundararaman and T. A. Arias: Phys. Rev. B 87, 165122 (2013)
1259 E0: float
1260 Energy (in eV) used for fitting in the plasmon-pole approximation.
1261 q0_correction: bool
1262 Analytic correction to the q=0 contribution applicable to 2D
1263 systems.
1264 nblocks: int
1265 Number of blocks chi0 should be distributed in so each core
1266 does not have to store the entire matrix. This is to reduce
1267 memory requirement. nblocks must be less than or equal to the
1268 number of processors.
1269 nblocksmax: bool
1270 Cuts chi0 into as many blocks as possible to reduce memory
1271 requirements as much as possible.
1272 output_prefix: None | str
1273 Where to direct the txt output. If set to None (default),
1274 will be deduced from filename (the default output prefix).
1275 This is to allow multiple processes to work on same cache
1276 (given by filename-prefix), while writing to different out
1277 files.
1278 """
1279 if fxc_mode:
1280 assert fxc_modes is None
1281 if fxc_modes:
1282 assert fxc_mode is None
1284 frequencies = get_frequencies(frequencies, domega0, omega2)
1286 integrate_gamma = GammaIntegrationMode(integrate_gamma)
1288 # We pass a serial communicator because the parallel handling
1289 # is somewhat wonky, we'd rather do that ourselves:
1290 try:
1291 qcache = FileCache(f'qcache_{filename}',
1292 comm=mpi.SerialCommunicator())
1293 except TypeError as err:
1294 raise RuntimeError(
1295 'File cache requires ASE master '
1296 'from September 20 2022 or newer. '
1297 'You may need to pull newest ASE.') from err
1299 mode = 'a' if qcache.filecount() > 1 else 'w'
1301 # (calc can not actually be a calculator at all.)
1302 gpwfile = Path(calc)
1304 output_prefix = filename or output_prefix
1305 context = ResponseContext(txt=output_prefix + '.txt',
1306 comm=world, timer=timer)
1307 gs = ResponseGroundStateAdapter.from_gpw_file(gpwfile)
1309 # Check if nblocks is compatible, adjust if not
1310 if nblocksmax:
1311 nblocks = get_max_nblocks(context.comm, gpwfile, ecut)
1313 kpts = list(select_kpts(kpts, gs.kd))
1315 ecut, ecut_e = choose_ecut_things(ecut, ecut_extrapolation)
1317 if nbands is None:
1318 nbands = int(gs.volume * (ecut / Ha)**1.5 * 2**0.5 / 3 / pi**2)
1319 else:
1320 if ecut_extrapolation:
1321 raise RuntimeError(
1322 'nbands cannot be supplied with ecut-extrapolation.')
1324 if ppa:
1325 # ppa reformulated as mpa with one pole
1326 mpa = {'npoles': 1, 'wrange': [0, 0], 'varpi': E0,
1327 'eta0': 1e-6, 'eta_rest': Ha, 'alpha': 1}
1329 if mpa:
1331 frequencies = mpa_frequency_sampling(**mpa)
1333 parameters = {'eta': 1e-6,
1334 'hilbert': False,
1335 'timeordered': False}
1337 else:
1338 # use nonlinear frequency grid
1339 frequencies = get_frequencies(frequencies, domega0, omega2)
1341 parameters = {'eta': eta,
1342 'hilbert': True,
1343 'timeordered': True}
1344 wd = get_frequency_descriptor(frequencies, gs=gs, nbands=nbands)
1346 wcontext = context.with_txt(output_prefix + '.w.txt', mode=mode)
1348 chi0calc = Chi0Calculator(
1349 gs, wcontext, nblocks=nblocks,
1350 wd=wd,
1351 nbands=nbands,
1352 ecut=ecut,
1353 intraband=False,
1354 **parameters)
1356 bands = choose_bands(bands, relbands, gs.nvalence, chi0calc.gs.nocc2)
1358 coulomb = CoulombKernel.from_gs(gs, truncation=truncation)
1359 # XXX eta needs to be converted to Hartree here,
1360 # XXX and it is also converted to Hartree at superclass constructor
1361 # XXX called below. This needs to be cleaned up.
1362 wcalc = initialize_w_calculator(chi0calc, wcontext,
1363 mpa=mpa,
1364 xc=xc,
1365 E0=E0, eta=eta / Ha, coulomb=coulomb,
1366 integrate_gamma=integrate_gamma,
1367 q0_correction=q0_correction)
1369 if fxc_mode:
1370 fxc_modes = [fxc_mode]
1372 if do_GW_too:
1373 fxc_modes.append('GW')
1375 exx_vxc_calculator = EXXVXCCalculator(
1376 gpwfile,
1377 snapshotfile_prefix=filename)
1379 super().__init__(filename=filename,
1380 wd=wd,
1381 chi0calc=chi0calc,
1382 wcalc=wcalc,
1383 ecut_e=ecut_e,
1384 eta=eta,
1385 fxc_modes=fxc_modes,
1386 nbands=nbands,
1387 bands=bands,
1388 frequencies=frequencies,
1389 kpts=kpts,
1390 exx_vxc_calculator=exx_vxc_calculator,
1391 qcache=qcache,
1392 ppa=ppa,
1393 mpa=mpa,
1394 **kwargs)
1396 @property
1397 def results_GW(self):
1398 # Compatibility with old "do_GW_too" behaviour
1399 if 'GW' in self.fxc_modes and self.fxc_modes[0] != 'GW':
1400 return self.all_results['GW']
1402 @property
1403 def results(self):
1404 return self.all_results[self.fxc_modes[0]]
1407class EXXVXCCalculator:
1408 """EXX and Kohn-Sham XC contribution."""
1410 def __init__(self, gpwfile, snapshotfile_prefix):
1411 self._gpwfile = gpwfile
1412 self._snapshotfile_prefix = snapshotfile_prefix
1414 def calculate(self, n1, n2, kpt_indices):
1415 _, vxc_skn, exx_skn = non_self_consistent_eigenvalues(
1416 self._gpwfile,
1417 'EXX',
1418 n1, n2,
1419 kpt_indices=kpt_indices,
1420 snapshot=f'{self._snapshotfile_prefix}-vxc-exx.json',
1421 )
1422 return vxc_skn / Ha, exx_skn / Ha