1import numbers 2import sys 3from functools import partial 4from time import ctime 5 6import gpaw 7import gpaw.mpi as mpi 8import numpy as np 9from ase.units import Ha 10from ase.utils.timing import Timer, timer 11from gpaw.blacs import BlacsDescriptor, BlacsGrid, Redistributor 12from gpaw.bztools import convex_hull_volume 13from gpaw.kpt_descriptor import KPointDescriptor 14from gpaw.response.integrators import PointIntegrator, TetrahedronIntegrator 15from gpaw.response.pair import PairDensity, PWSymmetryAnalyzer 16from gpaw.utilities import devnull 17from gpaw.utilities.blas import gemm 18from gpaw.utilities.memory import maxrss 19from gpaw.wavefunctions.pw import PWDescriptor 20 21 22class ArrayDescriptor: 23 """Describes a single dimensional array.""" 24 25 def __init__(self, data_x): 26 self.data_x = np.array(np.sort(data_x)) 27 self._data_len = len(data_x) 28 29 def __len__(self): 30 return self._data_len 31 32 def get_data(self): 33 return self.data_x 34 35 def get_closest_index(self, scalars_w): 36 """Get closest index. 37 38 Get closest index approximating scalars from below.""" 39 diff_xw = self.data_x[:, np.newaxis] - scalars_w[np.newaxis] 40 return np.argmin(diff_xw, axis=0) 41 42 def get_index_range(self, lim1_m, lim2_m): 43 """Get index range. """ 44 45 i0_m = np.zeros(len(lim1_m), int) 46 i1_m = np.zeros(len(lim2_m), int) 47 48 for m, (lim1, lim2) in enumerate(zip(lim1_m, lim2_m)): 49 i_x = np.logical_and(lim1 <= self.data_x, 50 lim2 >= self.data_x) 51 if i_x.any(): 52 inds = np.argwhere(i_x) 53 i0_m[m] = inds.min() 54 i1_m[m] = inds.max() + 1 55 56 return i0_m, i1_m 57 58 59class FrequencyDescriptor(ArrayDescriptor): 60 61 def __init__(self, domega0, omega2, omegamax): 62 beta = (2**0.5 - 1) * domega0 / omega2 63 wmax = int(omegamax / (domega0 + beta * omegamax)) 64 w = np.arange(wmax + 2) # + 2 is for buffer 65 omega_w = w * domega0 / (1 - beta * w) 66 67 ArrayDescriptor.__init__(self, omega_w) 68 69 self.domega0 = domega0 70 self.omega2 = omega2 71 self.omegamax = omegamax 72 self.omegamin = 0 73 74 self.beta = beta 75 self.wmax = wmax 76 self.omega_w = omega_w 77 self.wmax = wmax 78 self.nw = len(omega_w) 79 80 def get_closest_index(self, o_m): 81 beta = self.beta 82 w_m = (o_m / (self.domega0 + beta * o_m)).astype(int) 83 if isinstance(w_m, np.ndarray): 84 w_m[w_m >= self.wmax] = self.wmax - 1 85 elif isinstance(w_m, numbers.Integral): 86 if w_m >= self.wmax: 87 w_m = self.wmax - 1 88 else: 89 raise TypeError 90 return w_m 91 92 def get_index_range(self, omega1_m, omega2_m): 93 omega1_m = omega1_m.copy() 94 omega2_m = omega2_m.copy() 95 omega1_m[omega1_m < 0] = 0 96 omega2_m[omega2_m < 0] = 0 97 w1_m = self.get_closest_index(omega1_m) 98 w2_m = self.get_closest_index(omega2_m) 99 o1_m = self.omega_w[w1_m] 100 o2_m = self.omega_w[w2_m] 101 w1_m[o1_m < omega1_m] += 1 102 w2_m[o2_m < omega2_m] += 1 103 return w1_m, w2_m 104 105 106def frequency_grid(domega0, omega2, omegamax): 107 beta = (2**0.5 - 1) * domega0 / omega2 108 wmax = int(omegamax / (domega0 + beta * omegamax)) + 2 109 w = np.arange(wmax) 110 omega_w = w * domega0 / (1 - beta * w) 111 return omega_w 112 113 114class Chi0: 115 """Class for calculating non-interacting response functions.""" 116 117 def __init__(self, calc, response='density', 118 frequencies=None, domega0=0.1, omega2=10.0, omegamax=None, 119 ecut=50, gammacentered=False, hilbert=True, nbands=None, 120 timeordered=False, eta=0.2, ftol=1e-6, threshold=1, 121 real_space_derivatives=False, intraband=True, 122 world=mpi.world, txt='-', timer=None, 123 nblocks=1, gate_voltage=None, 124 disable_point_group=False, disable_time_reversal=False, 125 disable_non_symmorphic=True, 126 integrationmode=None, 127 pbc=None, rate=0.0, eshift=0.0, 128 paw_correction='brute-force'): 129 """Construct Chi0 object. 130 131 Parameters 132 ---------- 133 calc : str 134 The groundstate calculation file that the linear response 135 calculation is based on. 136 response : str 137 Type of response function. Currently collinear, scalar options 138 'density', '+-' and '-+' are implemented. 139 frequencies : ndarray or None 140 Array of frequencies to evaluate the response function at. If None, 141 frequencies are determined using the frequency_grid function in 142 gpaw.response.chi0. 143 domega0, omega2, omegamax : float 144 Input parameters for frequency_grid. 145 ecut : float 146 Energy cutoff. 147 gammacentered : bool 148 Center the grid of plane waves around the gamma point or q-vector 149 hilbert : bool 150 Switch for hilbert transform. If True, the full density response 151 is determined from a hilbert transform of its spectral function. 152 This is typically much faster, but does not work for imaginary 153 frequencies. 154 nbands : int 155 Maximum band index to include. 156 timeordered : bool 157 Switch for calculating the time ordered density response function. 158 In this case the hilbert transform cannot be used. 159 eta : float 160 Artificial broadening of spectra. 161 ftol : float 162 Threshold determining whether a band is completely filled 163 (f > 1 - ftol) or completely empty (f < ftol). 164 threshold : float 165 Numerical threshold for the optical limit k dot p perturbation 166 theory expansion (used in gpaw/response/pair.py). 167 real_space_derivatives : bool 168 Switch for calculating nabla matrix elements (in the optical limit) 169 using a real space finite difference approximation. 170 intraband : bool 171 Switch for including the intraband contribution to the density 172 response function. 173 world : MPI comm instance 174 MPI communicator. 175 txt : str 176 Output file. 177 timer : gpaw.utilities.timing.timer instance 178 nblocks : int 179 Divide the response function into nblocks. Useful when the response 180 function is large. 181 gate_voltage : float 182 Shift the fermi level by gate_voltage [Hartree]. 183 disable_point_group : bool 184 Do not use the point group symmetry operators. 185 disable_time_reversal : bool 186 Do not use time reversal symmetry. 187 disable_non_symmorphic : bool 188 Do no use non symmorphic symmetry operators. 189 integrationmode : str 190 Integrator for the kpoint integration. 191 If == 'tetrahedron integration' then the kpoint integral is 192 performed using the linear tetrahedron method. 193 pbc : list 194 Periodic directions of the system. Defaults to [True, True, True]. 195 eshift : float 196 Shift unoccupied bands 197 rate : float,str 198 Phenomenological scattering rate to use in optical limit Drude term 199 (in eV). If rate='eta', then use input artificial broadening eta as 200 rate. Note, for consistency with the formalism the rate is 201 implemented as omegap^2 / (omega + 1j * rate)^2 which differ from 202 some literature by a factor of 2. 203 204 205 Attributes 206 ---------- 207 pair : gpaw.response.pair.PairDensity instance 208 Class for calculating matrix elements of pairs of wavefunctions. 209 210 """ 211 212 self.response = response 213 214 self.timer = timer or Timer() 215 216 self.pair = PairDensity(calc, ecut, self.response, 217 ftol, threshold, 218 real_space_derivatives, world, txt, 219 self.timer, 220 nblocks=nblocks, 221 gate_voltage=gate_voltage, 222 paw_correction=paw_correction) 223 224 self.disable_point_group = disable_point_group 225 self.disable_time_reversal = disable_time_reversal 226 self.disable_non_symmorphic = disable_non_symmorphic 227 self.integrationmode = integrationmode 228 self.eshift = eshift / Ha 229 230 calc = self.pair.calc 231 self.calc = calc 232 233 if world.rank != 0: 234 txt = devnull 235 elif txt == '-': 236 txt = sys.stdout 237 elif isinstance(txt, str): 238 txt = open(txt, 'w') 239 self.fd = txt 240 241 self.vol = abs(np.linalg.det(calc.wfs.gd.cell_cv)) 242 243 self.world = world 244 245 if nblocks == 1: 246 self.blockcomm = self.world.new_communicator([world.rank]) 247 self.kncomm = world 248 else: 249 assert world.size % nblocks == 0, world.size 250 rank1 = world.rank // nblocks * nblocks 251 rank2 = rank1 + nblocks 252 self.blockcomm = self.world.new_communicator(range(rank1, rank2)) 253 ranks = range(world.rank % nblocks, world.size, nblocks) 254 self.kncomm = self.world.new_communicator(ranks) 255 256 self.nblocks = nblocks 257 258 if ecut is not None: 259 ecut /= Ha 260 261 self.ecut = ecut 262 self.gammacentered = gammacentered 263 264 self.eta = eta / Ha 265 if rate == 'eta': 266 self.rate = self.eta 267 else: 268 self.rate = rate / Ha 269 self.domega0 = domega0 / Ha 270 self.omega2 = omega2 / Ha 271 self.omegamax = None if omegamax is None else omegamax / Ha 272 self.nbands = nbands or self.calc.wfs.bd.nbands 273 self.include_intraband = intraband 274 275 omax = self.find_maximum_frequency() 276 277 if frequencies is None: 278 if self.omegamax is None: 279 self.omegamax = omax 280 print('Using nonlinear frequency grid from 0 to %.3f eV' % 281 (self.omegamax * Ha), file=self.fd) 282 self.wd = FrequencyDescriptor(self.domega0, self.omega2, 283 self.omegamax) 284 else: 285 self.wd = ArrayDescriptor(np.asarray(frequencies) / Ha) 286 assert not hilbert 287 288 self.omega_w = self.wd.get_data() 289 self.hilbert = hilbert 290 self.timeordered = bool(timeordered) 291 292 if self.eta == 0.0: 293 assert not hilbert 294 assert not timeordered 295 assert not self.omega_w.real.any() 296 297 self.nocc1 = self.pair.nocc1 # number of completely filled bands 298 self.nocc2 = self.pair.nocc2 # number of non-empty bands 299 300 self.Q_aGii = None 301 302 if pbc is not None: 303 self.pbc = np.array(pbc) 304 else: 305 self.pbc = np.array([True, True, True]) 306 307 if self.pbc is not None and (~self.pbc).any(): 308 assert np.sum((~self.pbc).astype(int)) == 1, \ 309 print('Only one non-periodic direction supported atm.') 310 print('Nonperiodic BC\'s: ', (~self.pbc), 311 file=self.fd) 312 313 if integrationmode is not None: 314 print('Using integration method: ' + self.integrationmode, 315 file=self.fd) 316 else: 317 print('Using integration method: PointIntegrator', file=self.fd) 318 319 def find_maximum_frequency(self): 320 """Determine the maximum electron-hole pair transition energy.""" 321 self.epsmin = 10000.0 322 self.epsmax = -10000.0 323 for kpt in self.calc.wfs.kpt_u: 324 self.epsmin = min(self.epsmin, kpt.eps_n[0]) 325 self.epsmax = max(self.epsmax, kpt.eps_n[self.nbands - 1]) 326 327 print('Minimum eigenvalue: %10.3f eV' % (self.epsmin * Ha), 328 file=self.fd) 329 print('Maximum eigenvalue: %10.3f eV' % (self.epsmax * Ha), 330 file=self.fd) 331 332 return self.epsmax - self.epsmin 333 334 def calculate(self, q_c, spin='all', A_x=None): 335 """Calculate response function. 336 337 Parameters 338 ---------- 339 q_c : list or ndarray 340 Momentum vector. 341 spin : str or int 342 If 'all' then include all spins. 343 If 0 or 1, only include this specific spin. 344 (not used in transverse response functions) 345 A_x : ndarray 346 Output array. If None, the output array is created. 347 348 Returns 349 ------- 350 pd : Planewave descriptor 351 Planewave descriptor for q_c. 352 chi0_wGG : ndarray 353 The response function. 354 chi0_wxvG : ndarray or None 355 (Only in optical limit) Wings of the density response function. 356 chi0_wvv : ndarray or None 357 (Only in optical limit) Head of the density response function. 358 359 """ 360 wfs = self.calc.wfs 361 362 if self.response == 'density': 363 if spin == 'all': 364 spins = range(wfs.nspins) 365 else: 366 assert spin in range(wfs.nspins) 367 spins = [spin] 368 else: 369 if self.response == '+-': 370 spins = [0] 371 elif self.response == '-+': 372 spins = [1] 373 else: 374 raise ValueError('Invalid response %s' % self.response) 375 376 q_c = np.asarray(q_c, dtype=float) 377 optical_limit = np.allclose(q_c, 0.0) and self.response == 'density' 378 379 pd = self.get_PWDescriptor(q_c, self.gammacentered) 380 381 self.print_chi(pd) 382 383 if gpaw.dry_run: 384 print(' Dry run exit', file=self.fd) 385 raise SystemExit 386 387 nG = pd.ngmax + 2 * optical_limit 388 nw = len(self.omega_w) 389 mynG = (nG + self.blockcomm.size - 1) // self.blockcomm.size 390 self.Ga = min(self.blockcomm.rank * mynG, nG) 391 self.Gb = min(self.Ga + mynG, nG) 392 # if self.blockcomm.rank == 0: 393 # assert self.Gb - self.Ga >= 3 394 # assert mynG * (self.blockcomm.size - 1) < nG 395 if A_x is not None: 396 nx = nw * (self.Gb - self.Ga) * nG 397 chi0_wGG = A_x[:nx].reshape((nw, self.Gb - self.Ga, nG)) 398 chi0_wGG[:] = 0.0 399 else: 400 chi0_wGG = np.zeros((nw, self.Gb - self.Ga, nG), complex) 401 402 if optical_limit: 403 chi0_wxvG = np.zeros((len(self.omega_w), 2, 3, nG), complex) 404 chi0_wvv = np.zeros((len(self.omega_w), 3, 3), complex) 405 self.plasmafreq_vv = np.zeros((3, 3), complex) 406 else: 407 chi0_wxvG = None 408 chi0_wvv = None 409 self.plasmafreq_vv = None 410 411 if self.response in ['+-', '-+']: 412 # Do all bands 413 m1 = 0 414 else: 415 # Do all empty bands: 416 m1 = self.nocc1 417 m2 = self.nbands 418 419 pd, chi0_wGG, chi0_wxvG, chi0_wvv = self._calculate(pd, 420 chi0_wGG, 421 chi0_wxvG, 422 chi0_wvv, 423 m1, m2, spins) 424 425 return pd, chi0_wGG, chi0_wxvG, chi0_wvv 426 427 @timer('Calculate CHI_0') 428 def _calculate(self, pd, chi0_wGG, chi0_wxvG, chi0_wvv, m1, m2, spins, 429 extend_head=True): 430 """In-place calculation of the response function. 431 432 Parameters 433 ---------- 434 q_c : list or ndarray 435 Momentum vector.. 436 pd : Planewave descriptor 437 Planewave descriptor for q_c. 438 chi0_wGG : ndarray 439 The response function. 440 chi0_wxvG : ndarray or None 441 Wings of the density response function. 442 chi0_wvv : ndarray or None 443 Head of the density response function. 444 m1 : int 445 Lower band cutoff for band summation 446 m2 : int 447 Upper band cutoff for band summation 448 spins : str or list(ints) 449 If 'all' then include all spins. 450 If [0] or [1], only include this specific spin. 451 extend_head : bool 452 If True: Extend the wings and head of chi in the optical limit to 453 take into account the non-analytic nature of chi. Effectively 454 means that chi has dimension (nw, nG + 2, nG + 2) in the optical 455 limit. This simplifies the code and should only be switched off 456 for parts of the code that do not support this feature i.e., GW 457 RPA total energy and RALDA. 458 """ 459 460 # Parse spins 461 wfs = self.calc.wfs 462 if spins == 'all': 463 spins = range(wfs.nspins) 464 else: 465 for spin in spins: 466 assert spin in range(wfs.nspins) 467 468 # Are we calculating the optical limit. 469 optical_limit = np.allclose(pd.kd.bzk_kc[0], 0.0) and \ 470 self.response == 'density' 471 472 # Use wings in optical limit, if head cannot be extended 473 if optical_limit and not extend_head: 474 wings = True 475 else: 476 wings = False 477 478 # Reset PAW correction in case momentum has change 479 self.Q_aGii = self.pair.initialize_paw_corrections(pd) 480 A_wxx = chi0_wGG # Change notation 481 482 # Initialize integrator. The integrator class is a general class 483 # for brillouin zone integration that can integrate user defined 484 # functions over user defined domains and sum over bands. 485 if self.integrationmode is None or \ 486 self.integrationmode == 'point integration': 487 integrator = PointIntegrator(self.pair.calc.wfs.gd.cell_cv, 488 response=self.response, 489 comm=self.world, 490 timer=self.timer, 491 txt=self.fd, 492 eshift=self.eshift, 493 nblocks=self.nblocks) 494 intnoblock = PointIntegrator(self.pair.calc.wfs.gd.cell_cv, 495 response=self.response, 496 comm=self.world, 497 timer=self.timer, 498 eshift=self.eshift, 499 txt=self.fd) 500 elif self.integrationmode == 'tetrahedron integration': 501 integrator = TetrahedronIntegrator(self.pair.calc.wfs.gd.cell_cv, 502 response=self.response, 503 comm=self.world, 504 timer=self.timer, 505 eshift=self.eshift, 506 txt=self.fd, 507 nblocks=self.nblocks) 508 intnoblock = TetrahedronIntegrator(self.pair.calc.wfs.gd.cell_cv, 509 response=self.response, 510 comm=self.world, 511 timer=self.timer, 512 eshift=self.eshift, 513 txt=self.fd) 514 else: 515 print('Integration mode ' + self.integrationmode + 516 ' not implemented.', file=self.fd) 517 raise NotImplementedError 518 519 # The integration domain is determined by the following function 520 # that reduces the integration domain to the irreducible zone 521 # of the little group of q. 522 bzk_kv, PWSA = self.get_kpoints(pd, 523 integrationmode=self.integrationmode) 524 domain = (bzk_kv, spins) 525 526 if self.integrationmode == 'tetrahedron integration': 527 # If there are non-periodic directions it is possible that the 528 # integration domain is not compatible with the symmetry operations 529 # which essentially means that too large domains will be 530 # integrated. We normalize by vol(BZ) / vol(domain) to make 531 # sure that to fix this. 532 domainvol = convex_hull_volume(bzk_kv) * PWSA.how_many_symmetries() 533 bzvol = (2 * np.pi)**3 / self.vol 534 factor = bzvol / domainvol 535 else: 536 factor = 1 537 538 prefactor = (2 * factor * PWSA.how_many_symmetries() / 539 (wfs.nspins * (2 * np.pi)**3)) # Remember prefactor 540 541 if self.integrationmode is None: 542 if self.calc.wfs.kd.refine_info is not None: 543 nbzkpts = self.calc.wfs.kd.refine_info.mhnbzkpts 544 else: 545 nbzkpts = self.calc.wfs.kd.nbzkpts 546 prefactor *= len(bzk_kv) / nbzkpts 547 548 A_wxx /= prefactor 549 if wings: 550 chi0_wxvG /= prefactor 551 chi0_wvv /= prefactor 552 553 # The functions that are integrated are defined in the bottom 554 # of this file and take a number of constant keyword arguments 555 # which the integrator class accepts through the use of the 556 # kwargs keyword. 557 kd = self.calc.wfs.kd 558 mat_kwargs = {'kd': kd, 'pd': pd, 559 'symmetry': PWSA, 560 'integrationmode': self.integrationmode} 561 eig_kwargs = {'kd': kd, 'pd': pd} 562 563 if not extend_head: 564 mat_kwargs['extend_head'] = False 565 566 # Determine what "kind" of integral to make. 567 extraargs = {} # Initialize extra arguments to integration method. 568 if self.eta == 0: 569 # If eta is 0 then we must be working with imaginary frequencies. 570 # In this case chi is hermitian and it is therefore possible to 571 # reduce the computational costs by a only computing half of the 572 # response function. 573 kind = 'hermitian response function' 574 elif self.hilbert: 575 # The spectral function integrator assumes that the form of the 576 # integrand is a function (a matrix element) multiplied by 577 # a delta function and should return a function of at user defined 578 # x's (frequencies). Thus the integrand is tuple of two functions 579 # and takes an additional argument (x). 580 kind = 'spectral function' 581 else: 582 # Otherwise, we can make no simplifying assumptions of the 583 # form of the response function and we simply perform a brute 584 # force calculation of the response function. 585 kind = 'response function' 586 extraargs['eta'] = self.eta 587 extraargs['timeordered'] = self.timeordered 588 589 # Integrate response function 590 print('Integrating response function.', file=self.fd) 591 # Define band summation 592 if self.response == 'density': 593 bandsum = {'n1': 0, 'n2': self.nocc2, 'm1': m1, 'm2': m2} 594 mat_kwargs.update(bandsum) 595 eig_kwargs.update(bandsum) 596 else: 597 bandsum = {'n1': 0, 'n2': self.nbands, 'm1': m1, 'm2': m2} 598 mat_kwargs.update(bandsum) 599 eig_kwargs.update(bandsum) 600 601 integrator.integrate(kind=kind, # Kind of integral 602 domain=domain, # Integration domain 603 integrand=(self.get_matrix_element, # Integrand 604 self.get_eigenvalues), # Integrand 605 x=self.wd, # Frequency Descriptor 606 kwargs=(mat_kwargs, eig_kwargs), 607 # Arguments for integrand functions 608 out_wxx=A_wxx, # Output array 609 **extraargs) 610 # extraargs: Extra arguments to integration method 611 if wings: 612 mat_kwargs['extend_head'] = True 613 mat_kwargs['block'] = False 614 if self.eta == 0: 615 extraargs['eta'] = self.eta 616 # This is horrible but we need to update the wings manually 617 # in order to make them work with ralda, RPA and GW. This entire 618 # section can be deleted in the future if the ralda and RPA code is 619 # made compatible with the head and wing extension that other parts 620 # of the code is using. 621 chi0_wxvx = np.zeros(np.array(chi0_wxvG.shape) + 622 [0, 0, 0, 2], 623 complex) # Notice the wxv"x" for head extend 624 intnoblock.integrate(kind=kind + ' wings', # kind'o int. 625 domain=domain, # Integration domain 626 integrand=(self.get_matrix_element, # Intgrnd 627 self.get_eigenvalues), # Integrand 628 x=self.wd, # Frequency Descriptor 629 kwargs=(mat_kwargs, eig_kwargs), 630 # Arguments for integrand functions 631 out_wxx=chi0_wxvx, # Output array 632 **extraargs) 633 634 if self.hilbert: 635 # The integrator only returns the spectral function and a Hilbert 636 # transform is performed to return the real part of the density 637 # response function. 638 with self.timer('Hilbert transform'): 639 omega_w = self.wd.get_data() # Get frequencies 640 # Make Hilbert transform 641 ht = HilbertTransform(np.array(omega_w), self.eta, 642 timeordered=self.timeordered) 643 ht(A_wxx) 644 if wings: 645 ht(chi0_wxvx) 646 647 # In the optical limit additional work must be performed 648 # for the intraband response. 649 # Only compute the intraband response if there are partially 650 # unoccupied bands and only if the user has not disabled its 651 # calculation using the include_intraband keyword. 652 if optical_limit and self.nocc1 != self.nocc2: 653 # The intraband response is essentially just the calculation 654 # of the free space Drude plasma frequency. The calculation is 655 # similarly to the interband transitions documented above. 656 mat_kwargs = {'kd': kd, 'symmetry': PWSA, 657 'n1': self.nocc1, 'n2': self.nocc2, 658 'pd': pd} # Integrand arguments 659 eig_kwargs = {'kd': kd, 660 'n1': self.nocc1, 'n2': self.nocc2, 661 'pd': pd} # Integrand arguments 662 domain = (bzk_kv, spins) # Integration domain 663 fermi_level = self.pair.fermi_level # Fermi level 664 665 # Not so elegant solution but it works 666 plasmafreq_wvv = np.zeros((1, 3, 3), complex) # Output array 667 print('Integrating intraband density response.', file=self.fd) 668 669 # Depending on which integration method is used we 670 # have to pass different arguments 671 extraargs = {} 672 if self.integrationmode is None: 673 # Calculate intraband transitions at finite fermi smearing 674 extraargs['intraband'] = True # Calculate intraband 675 elif self.integrationmode == 'tetrahedron integration': 676 # Calculate intraband transitions at T=0 677 extraargs['x'] = ArrayDescriptor([-fermi_level]) 678 679 intnoblock.integrate(kind='spectral function', # Kind of integral 680 domain=domain, # Integration domain 681 # Integrands 682 integrand=(self.get_intraband_response, 683 self.get_intraband_eigenvalue), 684 # Integrand arguments 685 kwargs=(mat_kwargs, eig_kwargs), 686 out_wxx=plasmafreq_wvv, # Output array 687 **extraargs) # Extra args for int. method 688 689 # Again, not so pretty but that's how it is 690 plasmafreq_vv = plasmafreq_wvv[0].copy() 691 if self.include_intraband: 692 if extend_head: 693 va = min(self.Ga, 3) 694 vb = min(self.Gb, 3) 695 A_wxx[:, :vb - va, :3] += (plasmafreq_vv[va:vb] / 696 (self.omega_w[:, np.newaxis, 697 np.newaxis] + 698 1e-10 + self.rate * 1j)**2) 699 elif self.blockcomm.rank == 0: 700 A_wxx[:, 0, 0] += (plasmafreq_vv[2, 2] / 701 (self.omega_w + 1e-10 + 702 self.rate * 1j)**2) 703 704 # Save the plasmafrequency 705 try: 706 self.plasmafreq_vv += 4 * np.pi * plasmafreq_vv * prefactor 707 except AttributeError: 708 self.plasmafreq_vv = 4 * np.pi * plasmafreq_vv * prefactor 709 710 PWSA.symmetrize_wvv(self.plasmafreq_vv[np.newaxis]) 711 print('Plasma frequency:', file=self.fd) 712 print((self.plasmafreq_vv**0.5 * Ha).round(2), 713 file=self.fd) 714 715 # The response function is integrated only over the IBZ. The 716 # chi calculated above must therefore be extended to include the 717 # response from the full BZ. This extension can be performed as a 718 # simple post processing of the response function that makes 719 # sure that the response function fulfills the symmetries of the little 720 # group of q. Due to the specific details of the implementation the chi 721 # calculated above is normalized by the number of symmetries (as seen 722 # below) and then symmetrized. 723 A_wxx *= prefactor 724 725 tmpA_wxx = self.redistribute(A_wxx) 726 if extend_head: 727 PWSA.symmetrize_wxx(tmpA_wxx, 728 optical_limit=optical_limit) 729 else: 730 PWSA.symmetrize_wGG(tmpA_wxx) 731 if wings: 732 chi0_wxvG += chi0_wxvx[..., 2:] 733 chi0_wvv += chi0_wxvx[:, 0, :3, :3] 734 PWSA.symmetrize_wxvG(chi0_wxvG) 735 PWSA.symmetrize_wvv(chi0_wvv) 736 self.redistribute(tmpA_wxx, A_wxx) 737 738 # If point summation was used then the normalization of the 739 # response function is not right and we have to make up for this 740 # fact. 741 742 if wings: 743 chi0_wxvG *= prefactor 744 chi0_wvv *= prefactor 745 746 # In the optical limit, we have extended the wings and the head to 747 # account for their nonanalytic behaviour which means that the size of 748 # the chi0_wGG matrix is nw * (nG + 2)**2. Below we extract these 749 # parameters. 750 751 if optical_limit and extend_head: 752 # The wings are extracted 753 chi0_wxvG[:, 1, :, self.Ga:self.Gb] = np.transpose(A_wxx[..., 0:3], 754 (0, 2, 1)) 755 va = min(self.Ga, 3) 756 vb = min(self.Gb, 3) 757 # print(self.world.rank, va, vb, chi0_wxvG[:, 0, va:vb].shape, 758 # A_wxx[:, va:vb].shape, A_wxx.shape) 759 chi0_wxvG[:, 0, va:vb] = A_wxx[:, :vb - va] 760 761 # Add contributions on different ranks 762 self.blockcomm.sum(chi0_wxvG) 763 chi0_wvv[:] = chi0_wxvG[:, 0, :3, :3] 764 chi0_wxvG = chi0_wxvG[..., 2:] # Jesus, this is complicated 765 766 # The head is extracted 767 # if self.blockcomm.rank == 0: 768 # chi0_wvv[:] = A_wxx[:, :3, :3] 769 # self.blockcomm.broadcast(chi0_wvv, 0) 770 771 # It is easiest to redistribute over freqs to pick body 772 tmpA_wxx = self.redistribute(A_wxx) 773 chi0_wGG = tmpA_wxx[:, 2:, 2:] 774 chi0_wGG = self.redistribute(chi0_wGG) 775 776 elif optical_limit: 777 # Since chi_wGG is nonanalytic in the head 778 # and wings we have to take care that 779 # these are handled correctly. Note that 780 # it is important that the wings are overwritten first. 781 chi0_wGG[:, :, 0] = chi0_wxvG[:, 1, 2, self.Ga:self.Gb] 782 783 if self.blockcomm.rank == 0: 784 chi0_wGG[:, 0, :] = chi0_wxvG[:, 0, 2, :] 785 chi0_wGG[:, 0, 0] = chi0_wvv[:, 2, 2] 786 else: 787 chi0_wGG = A_wxx 788 789 return pd, chi0_wGG, chi0_wxvG, chi0_wvv 790 791 def get_PWDescriptor(self, q_c, gammacentered=False): 792 """Get the planewave descriptor of q_c.""" 793 qd = KPointDescriptor([q_c]) 794 pd = PWDescriptor(self.ecut, self.calc.wfs.gd, 795 complex, qd, gammacentered=gammacentered) 796 return pd 797 798 @timer('Get kpoints') 799 def get_kpoints(self, pd, integrationmode=None): 800 """Get the integration domain.""" 801 # Use symmetries 802 PWSA = PWSymmetryAnalyzer 803 PWSA = PWSA(self.calc.wfs.kd, pd, 804 timer=self.timer, txt=self.fd, 805 disable_point_group=self.disable_point_group, 806 disable_time_reversal=self.disable_time_reversal, 807 disable_non_symmorphic=self.disable_non_symmorphic) 808 809 if integrationmode is None: 810 K_gK = PWSA.group_kpoints() 811 bzk_kc = np.array([self.calc.wfs.kd.bzk_kc[K_K[0]] for 812 K_K in K_gK]) 813 elif integrationmode == 'tetrahedron integration': 814 bzk_kc = PWSA.get_reduced_kd(pbc_c=self.pbc).bzk_kc 815 if (~self.pbc).any(): 816 bzk_kc = np.append(bzk_kc, 817 bzk_kc + (~self.pbc).astype(int), 818 axis=0) 819 820 bzk_kv = np.dot(bzk_kc, pd.gd.icell_cv) * 2 * np.pi 821 822 return bzk_kv, PWSA 823 824 @timer('Get matrix element') 825 def get_matrix_element(self, k_v, s, n1=None, n2=None, 826 m1=None, m2=None, 827 pd=None, kd=None, 828 symmetry=None, integrationmode=None, 829 extend_head=True, block=True): 830 """A function that returns pair-densities. 831 832 A pair density is defined as:: 833 834 <snk| e^(-i (q + G) r) |s'mk+q>, 835 836 where s and s' are spins, n and m are band indices, k is 837 the kpoint and q is the momentum transfer. For dielectric 838 response s'=s, for the transverse magnetic response 839 s' is flipped with respect to s. 840 841 Parameters 842 ---------- 843 k_v : ndarray 844 Kpoint coordinate in cartesian coordinates. 845 s : int 846 Spin index. 847 n1 : int 848 Lower occupied band index. 849 n2 : int 850 Upper occupied band index. 851 m1 : int 852 Lower unoccupied band index. 853 m2 : int 854 Upper unoccupied band index. 855 pd : PlanewaveDescriptor instance 856 kd : KpointDescriptor instance 857 Calculator kpoint descriptor. 858 symmetry: gpaw.response.pair.PWSymmetryAnalyzer instance 859 PWSA object for handling symmetries of the kpoints. 860 integrationmode : str 861 The integration mode employed. 862 extend_head: Bool 863 Extend the head to include non-analytic behaviour 864 865 Return 866 ------ 867 n_nmG : ndarray 868 Pair densities. 869 """ 870 k_c = np.dot(pd.gd.cell_cv, k_v) / (2 * np.pi) 871 872 q_c = pd.kd.bzk_kc[0] 873 874 optical_limit = np.allclose(q_c, 0.0) and self.response == 'density' 875 876 extrapolate_q = False 877 if self.calc.wfs.kd.refine_info is not None: 878 K1 = self.pair.find_kpoint(k_c) 879 label = kd.refine_info.label_k[K1] 880 if label == 'zero': 881 return None 882 elif (kd.refine_info.almostoptical and label == 'mh'): 883 if not hasattr(self, 'pd0'): 884 self.pd0 = self.get_PWDescriptor([0, ] * 3) 885 pd = self.pd0 886 extrapolate_q = True 887 888 nG = pd.ngmax 889 weight = np.sqrt(symmetry.get_kpoint_weight(k_c) / 890 symmetry.how_many_symmetries()) 891 if self.Q_aGii is None: 892 self.Q_aGii = self.pair.initialize_paw_corrections(pd) 893 894 kptpair = self.pair.get_kpoint_pair(pd, s, k_c, n1, n2, 895 m1, m2, block=block) 896 897 m_m = np.arange(m1, m2) 898 n_n = np.arange(n1, n2) 899 900 n_nmG = self.pair.get_pair_density(pd, kptpair, n_n, m_m, 901 Q_aGii=self.Q_aGii, block=block) 902 903 if integrationmode is None: 904 n_nmG *= weight 905 906 df_nm = kptpair.get_occupation_differences(n_n, m_m) 907 if not self.response == 'density': 908 df_nm = np.abs(df_nm) 909 df_nm[df_nm <= 1e-20] = 0.0 910 n_nmG *= df_nm[..., np.newaxis]**0.5 911 912 if extrapolate_q: 913 q_v = np.dot(q_c, pd.gd.icell_cv) * (2 * np.pi) 914 nq_nm = np.dot(n_nmG[:, :, :3], q_v) 915 n_nmG = n_nmG[:, :, 2:] 916 n_nmG[:, :, 0] = nq_nm 917 918 if not extend_head and optical_limit: 919 n_nmG = np.copy(n_nmG[:, :, 2:]) 920 optical_limit = False 921 922 if extend_head and optical_limit: 923 return n_nmG.reshape(-1, nG + 2 * optical_limit) 924 else: 925 return n_nmG.reshape(-1, nG) 926 927 @timer('Get eigenvalues') 928 def get_eigenvalues(self, k_v, s, n1=None, n2=None, 929 m1=None, m2=None, 930 kd=None, pd=None, wfs=None, 931 filter=False): 932 """A function that can return the eigenvalues. 933 934 A simple function describing the integrand of 935 the response function which gives an output that 936 is compatible with the gpaw k-point integration 937 routines.""" 938 if wfs is None: 939 wfs = self.calc.wfs 940 941 kd = wfs.kd 942 k_c = np.dot(pd.gd.cell_cv, k_v) / (2 * np.pi) 943 q_c = pd.kd.bzk_kc[0] 944 K1 = self.pair.find_kpoint(k_c) 945 K2 = self.pair.find_kpoint(k_c + q_c) 946 947 ik1 = kd.bz2ibz_k[K1] 948 ik2 = kd.bz2ibz_k[K2] 949 kpt1 = wfs.kpt_qs[ik1][s] 950 assert wfs.kd.comm.size == 1 951 if self.response in ['+-', '-+']: 952 s2 = 1 - s 953 else: 954 s2 = s 955 kpt2 = wfs.kpt_qs[ik2][s2] 956 deps_nm = np.subtract(kpt1.eps_n[n1:n2][:, np.newaxis], 957 kpt2.eps_n[m1:m2]) 958 959 if filter: 960 fermi_level = self.pair.fermi_level 961 deps_nm[kpt1.eps_n[n1:n2] > fermi_level, :] = np.nan 962 deps_nm[:, kpt2.eps_n[m1:m2] < fermi_level] = np.nan 963 964 return deps_nm.reshape(-1) 965 966 def get_intraband_response(self, k_v, s, n1=None, n2=None, 967 kd=None, symmetry=None, pd=None, 968 integrationmode=None): 969 k_c = np.dot(pd.gd.cell_cv, k_v) / (2 * np.pi) 970 kpt1 = self.pair.get_k_point(s, k_c, n1, n2) 971 n_n = range(n1, n2) 972 973 vel_nv = self.pair.intraband_pair_density(kpt1, n_n) 974 975 if self.integrationmode is None: 976 f_n = kpt1.f_n 977 assert self.calc.wfs.occupations.name in ['fermi-dirac', 978 'zero-width'] 979 width = getattr(self.calc.wfs.occupations, '_width', 0.0) / Ha 980 if width > 1e-15: 981 dfde_n = - 1. / width * (f_n - f_n**2.0) 982 else: 983 dfde_n = np.zeros_like(f_n) 984 vel_nv *= np.sqrt(-dfde_n[:, np.newaxis]) 985 weight = np.sqrt(symmetry.get_kpoint_weight(k_c) / 986 symmetry.how_many_symmetries()) 987 vel_nv *= weight 988 989 return vel_nv 990 991 @timer('Intraband eigenvalue') 992 def get_intraband_eigenvalue(self, k_v, s, 993 n1=None, n2=None, kd=None, pd=None): 994 """A function that can return the eigenvalues. 995 996 A simple function describing the integrand of 997 the response function which gives an output that 998 is compatible with the gpaw k-point integration 999 routines.""" 1000 wfs = self.calc.wfs 1001 kd = wfs.kd 1002 k_c = np.dot(pd.gd.cell_cv, k_v) / (2 * np.pi) 1003 K1 = self.pair.find_kpoint(k_c) 1004 ik = kd.bz2ibz_k[K1] 1005 kpt1 = wfs.kpt_qs[ik][s] 1006 assert wfs.kd.comm.size == 1 1007 1008 return kpt1.eps_n[n1:n2] 1009 1010 @timer('redist') 1011 def redistribute(self, in_wGG, out_x=None): 1012 """Redistribute array. 1013 1014 Switch between two kinds of parallel distributions: 1015 1016 1) parallel over G-vectors (second dimension of in_wGG) 1017 2) parallel over frequency (first dimension of in_wGG) 1018 1019 Returns new array using the memory in the 1-d array out_x. 1020 """ 1021 1022 comm = self.blockcomm 1023 1024 if comm.size == 1: 1025 return in_wGG 1026 1027 nw = len(self.omega_w) 1028 nG = in_wGG.shape[2] 1029 mynw = (nw + comm.size - 1) // comm.size 1030 mynG = (nG + comm.size - 1) // comm.size 1031 1032 bg1 = BlacsGrid(comm, comm.size, 1) 1033 bg2 = BlacsGrid(comm, 1, comm.size) 1034 md1 = BlacsDescriptor(bg1, nw, nG**2, mynw, nG**2) 1035 md2 = BlacsDescriptor(bg2, nw, nG**2, nw, mynG * nG) 1036 1037 if len(in_wGG) == nw: 1038 mdin = md2 1039 mdout = md1 1040 else: 1041 mdin = md1 1042 mdout = md2 1043 1044 r = Redistributor(comm, mdin, mdout) 1045 1046 outshape = (mdout.shape[0], mdout.shape[1] // nG, nG) 1047 if out_x is None: 1048 out_wGG = np.empty(outshape, complex) 1049 else: 1050 out_wGG = out_x[:np.product(outshape)].reshape(outshape) 1051 1052 r.redistribute(in_wGG.reshape(mdin.shape), 1053 out_wGG.reshape(mdout.shape)) 1054 1055 return out_wGG 1056 1057 @timer('dist freq') 1058 def distribute_frequencies(self, chi0_wGG): 1059 """Distribute frequencies to all cores.""" 1060 1061 world = self.world 1062 comm = self.blockcomm 1063 1064 if world.size == 1: 1065 return chi0_wGG 1066 1067 nw = len(self.omega_w) 1068 nG = chi0_wGG.shape[2] 1069 mynw = (nw + world.size - 1) // world.size 1070 mynG = (nG + comm.size - 1) // comm.size 1071 1072 wa = min(world.rank * mynw, nw) 1073 wb = min(wa + mynw, nw) 1074 1075 if self.blockcomm.size == 1: 1076 return chi0_wGG[wa:wb].copy() 1077 1078 if self.kncomm.rank == 0: 1079 bg1 = BlacsGrid(comm, 1, comm.size) 1080 in_wGG = chi0_wGG.reshape((nw, -1)) 1081 else: 1082 bg1 = BlacsGrid(None, 1, 1) 1083 # bg1 = DryRunBlacsGrid(mpi.serial_comm, 1, 1) 1084 in_wGG = np.zeros((0, 0), complex) 1085 md1 = BlacsDescriptor(bg1, nw, nG**2, nw, mynG * nG) 1086 1087 bg2 = BlacsGrid(world, world.size, 1) 1088 md2 = BlacsDescriptor(bg2, nw, nG**2, mynw, nG**2) 1089 1090 r = Redistributor(world, md1, md2) 1091 shape = (wb - wa, nG, nG) 1092 out_wGG = np.empty(shape, complex) 1093 r.redistribute(in_wGG, out_wGG.reshape((wb - wa, nG**2))) 1094 1095 return out_wGG 1096 1097 def print_chi(self, pd): 1098 calc = self.calc 1099 gd = calc.wfs.gd 1100 1101 if gpaw.dry_run: 1102 from gpaw.mpi import SerialCommunicator 1103 size = gpaw.dry_run 1104 world = SerialCommunicator() 1105 world.size = size 1106 else: 1107 world = self.world 1108 1109 q_c = pd.kd.bzk_kc[0] 1110 nw = len(self.omega_w) 1111 ecut = self.ecut * Ha 1112 ns = calc.wfs.nspins 1113 nbands = self.nbands 1114 nk = calc.wfs.kd.nbzkpts 1115 nik = calc.wfs.kd.nibzkpts 1116 ngmax = pd.ngmax 1117 eta = self.eta * Ha 1118 wsize = world.size 1119 knsize = self.kncomm.size 1120 nocc = self.nocc1 1121 npocc = self.nocc2 1122 ngridpoints = gd.N_c[0] * gd.N_c[1] * gd.N_c[2] 1123 nstat = (ns * npocc + world.size - 1) // world.size 1124 occsize = nstat * ngridpoints * 16. / 1024**2 1125 bsize = self.blockcomm.size 1126 chisize = nw * pd.ngmax**2 * 16. / 1024**2 / bsize 1127 1128 p = partial(print, file=self.fd) 1129 1130 p('%s' % ctime()) 1131 p('Called response.chi0.calculate with') 1132 p(' q_c: [%f, %f, %f]' % (q_c[0], q_c[1], q_c[2])) 1133 p(' Number of frequency points: %d' % nw) 1134 p(' Planewave cutoff: %f' % ecut) 1135 p(' Number of spins: %d' % ns) 1136 p(' Number of bands: %d' % nbands) 1137 p(' Number of kpoints: %d' % nk) 1138 p(' Number of irredicible kpoints: %d' % nik) 1139 p(' Number of planewaves: %d' % ngmax) 1140 p(' Broadening (eta): %f' % eta) 1141 p(' world.size: %d' % wsize) 1142 p(' kncomm.size: %d' % knsize) 1143 p(' blockcomm.size: %d' % bsize) 1144 p(' Number of completely occupied states: %d' % nocc) 1145 p(' Number of partially occupied states: %d' % npocc) 1146 p() 1147 p(' Memory estimate of potentially large arrays:') 1148 p(' chi0_wGG: %f M / cpu' % chisize) 1149 p(' Occupied states: %f M / cpu' % occsize) 1150 p(' Memory usage before allocation: %f M / cpu' % (maxrss() / 1151 1024**2)) 1152 p() 1153 1154 1155class HilbertTransform: 1156 1157 def __init__(self, omega_w, eta, timeordered=False, gw=False, 1158 blocksize=500): 1159 """Analytic Hilbert transformation using linear interpolation. 1160 1161 Hilbert transform:: 1162 1163 oo 1164 / 1 1 1165 |dw' (-------------- - --------------) S(w'). 1166 / w - w' + i eta w + w' + i eta 1167 0 1168 1169 With timeordered=True, you get:: 1170 1171 oo 1172 / 1 1 1173 |dw' (-------------- - --------------) S(w'). 1174 / w - w' - i eta w + w' + i eta 1175 0 1176 1177 With gw=True, you get:: 1178 1179 oo 1180 / 1 1 1181 |dw' (-------------- + --------------) S(w'). 1182 / w - w' + i eta w + w' + i eta 1183 0 1184 1185 """ 1186 1187 self.blocksize = blocksize 1188 1189 if timeordered: 1190 self.H_ww = self.H(omega_w, -eta) + self.H(omega_w, -eta, -1) 1191 elif gw: 1192 self.H_ww = self.H(omega_w, eta) - self.H(omega_w, -eta, -1) 1193 else: 1194 self.H_ww = self.H(omega_w, eta) + self.H(omega_w, -eta, -1) 1195 1196 def H(self, o_w, eta, sign=1): 1197 """Calculate transformation matrix. 1198 1199 With s=sign (+1 or -1):: 1200 1201 oo 1202 / dw' 1203 X (w, eta) = | ---------------- S(w'). 1204 s / s w - w' + i eta 1205 0 1206 1207 Returns H_ij so that X_i = np.dot(H_ij, S_j), where:: 1208 1209 X_i = X (omega_w[i]) and S_j = S(omega_w[j]) 1210 s 1211 """ 1212 1213 nw = len(o_w) 1214 H_ij = np.zeros((nw, nw), complex) 1215 do_j = o_w[1:] - o_w[:-1] 1216 for i, o in enumerate(o_w): 1217 d_j = o_w - o * sign 1218 y_j = 1j * np.arctan(d_j / eta) + 0.5 * np.log(d_j**2 + eta**2) 1219 y_j = (y_j[1:] - y_j[:-1]) / do_j 1220 H_ij[i, :-1] = 1 - (d_j[1:] - 1j * eta) * y_j 1221 H_ij[i, 1:] -= 1 - (d_j[:-1] - 1j * eta) * y_j 1222 return H_ij 1223 1224 def __call__(self, S_wx): 1225 """Inplace transform""" 1226 B_wx = S_wx.reshape((len(S_wx), -1)) 1227 nw, nx = B_wx.shape 1228 tmp_wx = np.zeros((nw, min(nx, self.blocksize)), complex) 1229 for x in range(0, nx, self.blocksize): 1230 b_wx = B_wx[:, x:x + self.blocksize] 1231 c_wx = tmp_wx[:, :b_wx.shape[1]] 1232 gemm(1.0, b_wx, self.H_ww, 0.0, c_wx) 1233 b_wx[:] = c_wx 1234 1235 1236if __name__ == '__main__': 1237 do = 0.025 1238 eta = 0.1 1239 omega_w = frequency_grid(do, 10.0, 3) 1240 print(len(omega_w)) 1241 X_w = omega_w * 0j 1242 Xt_w = omega_w * 0j 1243 Xh_w = omega_w * 0j 1244 for o in -np.linspace(2.5, 2.9, 10): 1245 X_w += (1 / (omega_w + o + 1j * eta) - 1246 1 / (omega_w - o + 1j * eta)) / o**2 1247 Xt_w += (1 / (omega_w + o - 1j * eta) - 1248 1 / (omega_w - o + 1j * eta)) / o**2 1249 w = int(-o / do / (1 + 3 * -o / 10)) 1250 o1, o2 = omega_w[w:w + 2] 1251 assert o1 - 1e-12 <= -o <= o2 + 1e-12, (o1, -o, o2) 1252 p = 1 / (o2 - o1)**2 / o**2 1253 Xh_w[w] += p * (o2 - -o) 1254 Xh_w[w + 1] += p * (-o - o1) 1255 1256 ht = HilbertTransform(omega_w, eta, 1) 1257 ht(Xh_w) 1258 1259 import matplotlib.pyplot as plt 1260 plt.plot(omega_w, X_w.imag, label='ImX') 1261 plt.plot(omega_w, X_w.real, label='ReX') 1262 plt.plot(omega_w, Xt_w.imag, label='ImXt') 1263 plt.plot(omega_w, Xt_w.real, label='ReXt') 1264 plt.plot(omega_w, Xh_w.imag, label='ImXh') 1265 plt.plot(omega_w, Xh_w.real, label='ReXh') 1266 plt.legend() 1267 plt.show() 1268