1# Copyright (C) 2020 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phono3py. 5# 6# Redistribution and use in source and binary forms, with or without 7# modification, are permitted provided that the following conditions 8# are met: 9# 10# * Redistributions of source code must retain the above copyright 11# notice, this list of conditions and the following disclaimer. 12# 13# * Redistributions in binary form must reproduce the above copyright 14# notice, this list of conditions and the following disclaimer in 15# the documentation and/or other materials provided with the 16# distribution. 17# 18# * Neither the name of the phonopy project nor the names of its 19# contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 33# POSSIBILITY OF SUCH DAMAGE. 34 35import sys 36import numpy as np 37from phonopy.units import Hbar, EV, THz 38from phonopy.phonon.degeneracy import degenerate_sets 39from phono3py.phonon3.triplets import get_triplets_integration_weights 40from phono3py.phonon.func import bose_einstein 41from phono3py.file_IO import (write_gamma_detail_to_hdf5, 42 write_imag_self_energy_at_grid_point) 43 44 45def get_imag_self_energy(interaction, 46 grid_points, 47 temperatures, 48 sigmas=None, 49 frequency_points=None, 50 frequency_step=None, 51 num_frequency_points=None, 52 num_points_in_batch=None, 53 scattering_event_class=None, # class 1 or 2 54 write_gamma_detail=False, 55 return_gamma_detail=False, 56 output_filename=None, 57 log_level=0): 58 """Imaginary part of self energy at frequency points 59 60 Band indices to be calculated at are kept in Interaction instance. 61 62 Parameters 63 ---------- 64 interaction : Interaction 65 Ph-ph interaction. 66 grid_points : array_like 67 Grid-point indices where imag-self-energeis are caclculated. 68 dtype=int, shape=(grid_points,) 69 temperatures : array_like 70 Temperatures where imag-self-energies are calculated. 71 dtype=float, shape=(temperatures,) 72 sigmas : array_like, optional 73 A set of sigmas. simgas=[None, ] means to use tetrahedron method, 74 otherwise smearing method with real positive value of sigma. 75 For example, sigmas=[None, 0.01, 0.03] is possible. Default is None, 76 which results in [None, ]. 77 dtype=float, shape=(sigmas,) 78 frequency_points : array_like, optional 79 Frequency sampling points. Default is None. In this case, 80 num_frequency_points or frequency_step is used to generate uniform 81 frequency sampling points. 82 dtype=float, shape=(frequency_points,) 83 frequency_step : float, optional 84 Uniform pitch of frequency sampling points. Default is None. This 85 results in using num_frequency_points. 86 num_frequency_points: int, optional 87 Number of sampling sampling points to be used instead of 88 frequency_step. This number includes end points. Default is None, 89 which gives 201. 90 num_points_in_batch: int, optional 91 Number of sampling points in one batch. This is for the frequency 92 sampling mode and the sampling points are divided into batches. 93 Lager number provides efficient use of multi-cores but more 94 memory demanding. Default is None, which give the number of 10. 95 scattering_event_class : int, optional 96 Specific choice of scattering event class, 1 or 2 that is specified 97 1 or 2, respectively. The result is stored in gammas. Therefore 98 usual gammas are not stored in the variable. Default is None, which 99 doesn't specify scattering_event_class. 100 write_gamma_detail : bool, optional 101 Detailed gammas are written into a file in hdf5. Default is False. 102 return_gamma_detail : bool, optional 103 With True, detailed gammas are returned. Default is False. 104 log_level: int 105 Log level. Default is 0. 106 107 Returns 108 ------- 109 tuple : 110 (frequency_points, gammas) are returned. With return_gamma_detail=True, 111 (frequency_points, gammas, detailed_gammas) are returned. 112 113 """ 114 115 if sigmas is None: 116 _sigmas = [None, ] 117 else: 118 _sigmas = sigmas 119 120 if (interaction.get_phonons()[2] == 0).any(): 121 if log_level: 122 print("Running harmonic phonon calculations...") 123 interaction.run_phonon_solver() 124 125 mesh = interaction.mesh_numbers 126 frequencies = interaction.get_phonons()[0] 127 max_phonon_freq = np.amax(frequencies) 128 129 _frequency_points = get_frequency_points( 130 max_phonon_freq=max_phonon_freq, 131 sigmas=_sigmas, 132 frequency_points=frequency_points, 133 frequency_step=frequency_step, 134 num_frequency_points=num_frequency_points) 135 136 ise = ImagSelfEnergy( 137 interaction, with_detail=(write_gamma_detail or return_gamma_detail)) 138 gamma = np.zeros( 139 (len(grid_points), len(_sigmas), len(temperatures), 140 len(interaction.band_indices), len(_frequency_points)), 141 dtype='double', order='C') 142 detailed_gamma = [] 143 144 for i, gp in enumerate(grid_points): 145 ise.set_grid_point(gp) 146 if log_level: 147 weights = interaction.get_triplets_at_q()[1] 148 print("------------------- Imaginary part of self energy (%d/%d) " 149 "-------------------" % (i + 1, len(grid_points))) 150 print("Grid point: %d" % gp) 151 print("Number of ir-triplets: " 152 "%d / %d" % (len(weights), weights.sum())) 153 154 ise.run_interaction() 155 156 if log_level: 157 adrs = interaction.grid_address[gp] 158 q = adrs.astype('double') / mesh 159 print("q-point: %s" % q) 160 print("Phonon frequency:") 161 text = "[ " 162 for bi, freq in enumerate(frequencies[gp]): 163 if bi % 6 == 0 and bi != 0: 164 text += "\n" 165 text += "%8.4f " % freq 166 text += "]" 167 print(text) 168 sys.stdout.flush() 169 170 if write_gamma_detail or return_gamma_detail: 171 (triplets, weights, 172 map_triplets, _) = interaction.get_triplets_at_q() 173 num_band0 = len(interaction.band_indices) 174 num_band = frequencies.shape[1] 175 detailed_gamma_at_gp = np.zeros( 176 (len(_sigmas), len(temperatures), len(_frequency_points), 177 len(weights), num_band0, num_band, num_band), 178 dtype='double') 179 else: 180 detailed_gamma_at_gp = None 181 182 for j, sigma in enumerate(_sigmas): 183 if log_level: 184 if sigma: 185 print("Sigma: %s" % sigma) 186 else: 187 print("Tetrahedron method is used for BZ integration.") 188 189 ise.set_sigma(sigma) 190 191 # Run one by one at frequency points 192 if detailed_gamma_at_gp is None: 193 detailed_gamma_at_gp_at_j = None 194 else: 195 detailed_gamma_at_gp_at_j = detailed_gamma_at_gp[j] 196 run_ise_at_frequency_points_batch( 197 _frequency_points, 198 ise, 199 temperatures, 200 gamma[i, j], 201 write_gamma_detail=write_gamma_detail, 202 return_gamma_detail=return_gamma_detail, 203 detailed_gamma_at_gp=detailed_gamma_at_gp_at_j, 204 scattering_event_class=scattering_event_class, 205 nelems_in_batch=num_points_in_batch, 206 log_level=log_level) 207 208 if write_gamma_detail: 209 full_filename = write_gamma_detail_to_hdf5( 210 temperatures, 211 mesh, 212 gamma_detail=detailed_gamma_at_gp[j], 213 grid_point=gp, 214 triplet=triplets, 215 weight=weights, 216 triplet_map=map_triplets, 217 sigma=sigma, 218 frequency_points=_frequency_points, 219 filename=output_filename) 220 221 if log_level: 222 print("Contribution of each triplet to imaginary part of " 223 "self energy is written in\n\"%s\"." % full_filename) 224 225 if return_gamma_detail: 226 detailed_gamma.append(detailed_gamma_at_gp) 227 228 if return_gamma_detail: 229 return _frequency_points, gamma, detailed_gamma 230 else: 231 return _frequency_points, gamma 232 233 234def get_frequency_points(max_phonon_freq=None, 235 sigmas=None, 236 frequency_points=None, 237 frequency_step=None, 238 num_frequency_points=None): 239 if frequency_points is None: 240 if sigmas is not None: 241 sigma_vals = [sigma for sigma in sigmas if sigma is not None] 242 else: 243 sigma_vals = [] 244 if sigma_vals: 245 fmax = max_phonon_freq * 2 + np.max(sigma_vals) * 4 246 else: 247 fmax = max_phonon_freq * 2 248 fmax *= 1.005 249 fmin = 0 250 _frequency_points = _sample_frequency_points( 251 fmin, 252 fmax, 253 frequency_step=frequency_step, 254 num_frequency_points=num_frequency_points) 255 else: 256 _frequency_points = np.array(frequency_points, dtype='double') 257 258 return _frequency_points 259 260 261def _sample_frequency_points(f_min, 262 f_max, 263 frequency_step=None, 264 num_frequency_points=None): 265 if num_frequency_points is None: 266 if frequency_step is not None: 267 frequency_points = np.arange( 268 f_min, f_max, frequency_step, dtype='double') 269 else: 270 frequency_points = np.array(np.linspace( 271 f_min, f_max, 201), dtype='double') 272 else: 273 frequency_points = np.array(np.linspace( 274 f_min, f_max, num_frequency_points), dtype='double') 275 276 return frequency_points 277 278 279def write_imag_self_energy(imag_self_energy, 280 mesh, 281 grid_points, 282 band_indices, 283 frequency_points, 284 temperatures, 285 sigmas, 286 scattering_event_class=None, 287 output_filename=None, 288 is_mesh_symmetry=True, 289 log_level=0): 290 for gp, ise_sigmas in zip(grid_points, imag_self_energy): 291 for sigma, ise_temps in zip(sigmas, ise_sigmas): 292 for t, ise in zip(temperatures, ise_temps): 293 for i, bi in enumerate(band_indices): 294 pos = 0 295 for j in range(i): 296 pos += len(band_indices[j]) 297 filename = write_imag_self_energy_at_grid_point( 298 gp, 299 bi, 300 mesh, 301 frequency_points, 302 ise[pos:(pos + len(bi))].sum(axis=0) / len(bi), 303 sigma=sigma, 304 temperature=t, 305 scattering_event_class=scattering_event_class, 306 filename=output_filename, 307 is_mesh_symmetry=is_mesh_symmetry) 308 if log_level: 309 print("Imaginary parts of self-energies were " 310 "written to \"%s\"." % filename) 311 312 313def average_by_degeneracy(imag_self_energy, band_indices, freqs_at_gp): 314 deg_sets = degenerate_sets(freqs_at_gp) 315 imag_se = np.zeros_like(imag_self_energy) 316 for dset in deg_sets: 317 bi_set = [] 318 for i, bi in enumerate(band_indices): 319 if bi in dset: 320 bi_set.append(i) 321 for i in bi_set: 322 if imag_self_energy.ndim == 1: 323 imag_se[i] = (imag_self_energy[bi_set].sum() / 324 len(bi_set)) 325 else: 326 imag_se[:, i] = ( 327 imag_self_energy[:, bi_set].sum(axis=1) / 328 len(bi_set)) 329 return imag_se 330 331 332class ImagSelfEnergy(object): 333 def __init__(self, 334 interaction, 335 frequency_points=None, 336 temperature=None, 337 sigma=None, 338 sigma_cutoff=None, 339 with_detail=False, 340 unit_conversion=None, 341 lang='C'): 342 self._pp = interaction 343 self._sigma = None 344 self.set_sigma(sigma, sigma_cutoff=sigma_cutoff) 345 self._temperature = None 346 self.set_temperature(temperature) 347 self._frequency_points = None 348 self.set_frequency_points(frequency_points) 349 self._grid_point = None 350 351 self._lang = lang 352 self._imag_self_energy = None 353 self._detailed_imag_self_energy = None 354 self._pp_strength = None 355 self._frequencies = None 356 self._triplets_at_q = None 357 self._weights_at_q = None 358 self._with_detail = with_detail 359 self._unit_conversion = None 360 self._cutoff_frequency = interaction.cutoff_frequency 361 362 self._g = None # integration weights 363 self._g_zero = None # Necessary elements of interaction strength 364 self._g_zero_frequency_points = None 365 self._g_zero_zeros = None # always zeros for frequency sampling mode 366 self._mesh = self._pp.mesh_numbers 367 self._is_collision_matrix = False 368 369 # Unit to THz of Gamma 370 if unit_conversion is None: 371 self._unit_conversion = (18 * np.pi / (Hbar * EV) ** 2 372 / (2 * np.pi * THz) ** 2 373 * EV ** 2) 374 else: 375 self._unit_conversion = unit_conversion 376 377 def run(self): 378 if self._pp_strength is None: 379 self.run_interaction() 380 381 num_band0 = self._pp_strength.shape[1] 382 if self._frequency_points is None: 383 self._imag_self_energy = np.zeros(num_band0, dtype='double') 384 if self._with_detail: 385 self._detailed_imag_self_energy = np.empty_like( 386 self._pp_strength) 387 self._detailed_imag_self_energy[:] = 0 388 self._ise_N = np.zeros_like(self._imag_self_energy) 389 self._ise_U = np.zeros_like(self._imag_self_energy) 390 self._run_with_band_indices() 391 else: 392 self._imag_self_energy = np.zeros( 393 (len(self._frequency_points), num_band0), 394 order='C', dtype='double') 395 if self._with_detail: 396 self._detailed_imag_self_energy = np.zeros( 397 (len(self._frequency_points),) + self._pp_strength.shape, 398 order='C', dtype='double') 399 self._ise_N = np.zeros_like(self._imag_self_energy) 400 self._ise_U = np.zeros_like(self._imag_self_energy) 401 self._run_with_frequency_points() 402 403 def run_interaction(self, is_full_pp=True): 404 if is_full_pp or self._frequency_points is not None: 405 self._pp.run(lang=self._lang) 406 else: 407 self._pp.run(lang=self._lang, g_zero=self._g_zero) 408 self._pp_strength = self._pp.interaction_strength 409 410 def set_integration_weights(self, scattering_event_class=None): 411 if self._frequency_points is None: 412 bi = self._pp.band_indices 413 f_points = self._frequencies[self._grid_point][bi] 414 else: 415 f_points = self._frequency_points 416 417 self._g, _g_zero = get_triplets_integration_weights( 418 self._pp, 419 np.array(f_points, dtype='double'), 420 self._sigma, 421 self._sigma_cutoff, 422 is_collision_matrix=self._is_collision_matrix) 423 if self._frequency_points is None: 424 self._g_zero = _g_zero 425 else: 426 # g_zero feature can not be used in frequency sampling mode. 427 # zero values of the following array shape is used in C-routine. 428 # shape = [num_triplets, num_band0, num_band, num_band] 429 shape = list(self._g.shape[1:]) 430 shape[1] = len(self._pp.band_indices) 431 self._g_zero_zeros = np.zeros(shape=shape, dtype='byte', order='C') 432 self._g_zero_frequency_points = _g_zero 433 434 if scattering_event_class == 1 or scattering_event_class == 2: 435 self._g[scattering_event_class - 1] = 0 436 437 def get_imag_self_energy(self): 438 if self._cutoff_frequency is None: 439 return self._imag_self_energy 440 else: 441 return self._average_by_degeneracy(self._imag_self_energy) 442 443 def get_imag_self_energy_N_and_U(self): 444 if self._cutoff_frequency is None: 445 return self._ise_N, self._ise_U 446 else: 447 return (self._average_by_degeneracy(self._ise_N), 448 self._average_by_degeneracy(self._ise_U)) 449 450 def get_detailed_imag_self_energy(self): 451 return self._detailed_imag_self_energy 452 453 def get_integration_weights(self): 454 return self._g, self._g_zero 455 456 def get_unit_conversion_factor(self): 457 return self._unit_conversion 458 459 def set_grid_point(self, grid_point=None, stores_triplets_map=False): 460 if grid_point is None: 461 self._grid_point = None 462 else: 463 self._pp.set_grid_point(grid_point, 464 stores_triplets_map=stores_triplets_map) 465 self._pp_strength = None 466 (self._triplets_at_q, 467 self._weights_at_q) = self._pp.get_triplets_at_q()[:2] 468 self._grid_point = grid_point 469 self._frequencies, self._eigenvectors, _ = self._pp.get_phonons() 470 471 def set_sigma(self, sigma, sigma_cutoff=None): 472 if sigma is None: 473 self._sigma = None 474 else: 475 self._sigma = float(sigma) 476 477 if sigma_cutoff is None: 478 self._sigma_cutoff = None 479 else: 480 self._sigma_cutoff = float(sigma_cutoff) 481 482 self.delete_integration_weights() 483 484 def set_frequency_points(self, frequency_points): 485 if frequency_points is None: 486 self._frequency_points = None 487 else: 488 self._frequency_points = np.array(frequency_points, dtype='double') 489 490 def set_temperature(self, temperature): 491 if temperature is None: 492 self._temperature = None 493 else: 494 self._temperature = float(temperature) 495 496 def set_averaged_pp_interaction(self, ave_pp): 497 num_triplets = len(self._triplets_at_q) 498 num_band = self._pp.get_primitive().get_number_of_atoms() * 3 499 num_grid = np.prod(self._mesh) 500 bi = self._pp.get_band_indices() 501 self._pp_strength = np.zeros( 502 (num_triplets, len(bi), num_band, num_band), dtype='double') 503 504 for i, v_ave in enumerate(ave_pp): 505 self._pp_strength[:, i, :, :] = v_ave / num_grid 506 507 def set_interaction_strength(self, pp_strength): 508 self._pp_strength = pp_strength 509 self._pp.set_interaction_strength(pp_strength, g_zero=self._g_zero) 510 511 def delete_integration_weights(self): 512 self._g = None 513 self._g_zero = None 514 self._pp_strength = None 515 516 def _run_with_band_indices(self): 517 if self._g is not None: 518 if self._lang == 'C': 519 if self._with_detail: 520 # self._detailed_imag_self_energy.shape = 521 # (num_triplets, num_band0, num_band, num_band) 522 # self._imag_self_energy is also set. 523 self._run_c_detailed_with_band_indices_with_g() 524 else: 525 # self._imag_self_energy.shape = (num_band0,) 526 self._run_c_with_band_indices_with_g() 527 else: 528 print("Running into _run_py_with_band_indices_with_g()") 529 print("This routine is super slow and only for the test.") 530 self._run_py_with_band_indices_with_g() 531 else: 532 print("get_triplets_integration_weights must be executed " 533 "before calling this method.") 534 import sys 535 sys.exit(1) 536 537 def _run_with_frequency_points(self): 538 if self._g is not None: 539 if self._lang == 'C': 540 if self._with_detail: 541 self._run_c_detailed_with_frequency_points_with_g() 542 else: 543 self._run_c_with_frequency_points_with_g() 544 else: 545 print("Running into _run_py_with_frequency_points_with_g()") 546 print("This routine is super slow and only for the test.") 547 self._run_py_with_frequency_points_with_g() 548 else: 549 print("get_triplets_integration_weights must be executed " 550 "before calling this method.") 551 import sys 552 sys.exit(1) 553 554 def _run_c_with_band_indices_with_g(self): 555 import phono3py._phono3py as phono3c 556 557 if self._g_zero is None: 558 _g_zero = np.zeros(self._pp_strength.shape, 559 dtype='byte', order='C') 560 else: 561 _g_zero = self._g_zero 562 563 phono3c.imag_self_energy_with_g(self._imag_self_energy, 564 self._pp_strength, 565 self._triplets_at_q, 566 self._weights_at_q, 567 self._frequencies, 568 self._temperature, 569 self._g, 570 _g_zero, 571 self._cutoff_frequency, 572 -1) 573 self._imag_self_energy *= self._unit_conversion 574 575 def _run_c_detailed_with_band_indices_with_g(self): 576 import phono3py._phono3py as phono3c 577 578 if self._g_zero is None: 579 _g_zero = np.zeros(self._pp_strength.shape, 580 dtype='byte', order='C') 581 else: 582 _g_zero = self._g_zero 583 584 phono3c.detailed_imag_self_energy_with_g( 585 self._detailed_imag_self_energy, 586 self._ise_N, # Normal 587 self._ise_U, # Umklapp 588 self._pp_strength, 589 self._triplets_at_q, 590 self._weights_at_q, 591 self._pp.get_grid_address(), 592 self._frequencies, 593 self._temperature, 594 self._g, 595 _g_zero, 596 self._cutoff_frequency) 597 598 self._detailed_imag_self_energy *= self._unit_conversion 599 self._ise_N *= self._unit_conversion 600 self._ise_U *= self._unit_conversion 601 self._imag_self_energy = self._ise_N + self._ise_U 602 603 def _run_c_with_frequency_points_with_g(self): 604 import phono3py._phono3py as phono3c 605 num_band0 = self._pp_strength.shape[1] 606 ise_at_f = np.zeros(num_band0, dtype='double') 607 608 for i in range(len(self._frequency_points)): 609 phono3c.imag_self_energy_with_g(ise_at_f, 610 self._pp_strength, 611 self._triplets_at_q, 612 self._weights_at_q, 613 self._frequencies, 614 self._temperature, 615 self._g, 616 self._g_zero_frequency_points, 617 self._cutoff_frequency, 618 i) 619 self._imag_self_energy[i] = ise_at_f 620 self._imag_self_energy *= self._unit_conversion 621 622 def _run_c_detailed_with_frequency_points_with_g(self): 623 import phono3py._phono3py as phono3c 624 num_band0 = self._pp_strength.shape[1] 625 g_shape = list(self._g.shape) 626 g_shape[2] = num_band0 627 g = np.zeros((2,) + self._pp_strength.shape, order='C', dtype='double') 628 detailed_ise_at_f = np.zeros( 629 self._detailed_imag_self_energy.shape[1:5], 630 order='C', dtype='double') 631 ise_at_f_N = np.zeros(num_band0, dtype='double') 632 ise_at_f_U = np.zeros(num_band0, dtype='double') 633 _g_zero = np.zeros(g_shape, dtype='byte', order='C') 634 635 for i in range(len(self._frequency_points)): 636 for j in range(g.shape[2]): 637 g[:, :, j, :, :] = self._g[:, :, i, :, :] 638 phono3c.detailed_imag_self_energy_with_g( 639 detailed_ise_at_f, 640 ise_at_f_N, 641 ise_at_f_U, 642 self._pp_strength, 643 self._triplets_at_q, 644 self._weights_at_q, 645 self._pp.grid_address, 646 self._frequencies, 647 self._temperature, 648 g, 649 _g_zero, 650 self._cutoff_frequency) 651 self._detailed_imag_self_energy[i] = (detailed_ise_at_f * 652 self._unit_conversion) 653 self._ise_N[i] = ise_at_f_N * self._unit_conversion 654 self._ise_U[i] = ise_at_f_U * self._unit_conversion 655 self._imag_self_energy[i] = self._ise_N[i] + self._ise_U[i] 656 657 def _run_py_with_band_indices_with_g(self): 658 if self._temperature > 0: 659 self._ise_thm_with_band_indices() 660 else: 661 self._ise_thm_with_band_indices_0K() 662 663 def _ise_thm_with_band_indices(self): 664 freqs = self._frequencies[self._triplets_at_q[:, [1, 2]]] 665 freqs = np.where(freqs > self._cutoff_frequency, freqs, 1) 666 n = bose_einstein(freqs, self._temperature) 667 for i, (tp, w, interaction) in enumerate(zip(self._triplets_at_q, 668 self._weights_at_q, 669 self._pp_strength)): 670 for j, k in list(np.ndindex(interaction.shape[1:])): 671 f1 = self._frequencies[tp[1]][j] 672 f2 = self._frequencies[tp[2]][k] 673 if (f1 > self._cutoff_frequency and 674 f2 > self._cutoff_frequency): 675 n2 = n[i, 0, j] 676 n3 = n[i, 1, k] 677 g1 = self._g[0, i, :, j, k] 678 g2_g3 = self._g[1, i, :, j, k] # g2 - g3 679 self._imag_self_energy[:] += ( 680 (n2 + n3 + 1) * g1 + 681 (n2 - n3) * (g2_g3)) * interaction[:, j, k] * w 682 683 self._imag_self_energy *= self._unit_conversion 684 685 def _ise_thm_with_band_indices_0K(self): 686 for i, (w, interaction) in enumerate(zip(self._weights_at_q, 687 self._pp_strength)): 688 for j, k in list(np.ndindex(interaction.shape[1:])): 689 g1 = self._g[0, i, :, j, k] 690 self._imag_self_energy[:] += g1 * interaction[:, j, k] * w 691 692 self._imag_self_energy *= self._unit_conversion 693 694 def _run_py_with_frequency_points_with_g(self): 695 if self._temperature > 0: 696 self._ise_thm_with_frequency_points() 697 else: 698 self._ise_thm_with_frequency_points_0K() 699 700 def _ise_thm_with_frequency_points(self): 701 for i, (tp, w, interaction) in enumerate(zip(self._triplets_at_q, 702 self._weights_at_q, 703 self._pp_strength)): 704 for j, k in list(np.ndindex(interaction.shape[1:])): 705 f1 = self._frequencies[tp[1]][j] 706 f2 = self._frequencies[tp[2]][k] 707 if (f1 > self._cutoff_frequency and 708 f2 > self._cutoff_frequency): 709 n2 = bose_einstein(f1, self._temperature) 710 n3 = bose_einstein(f2, self._temperature) 711 g1 = self._g[0, i, :, j, k] 712 g2_g3 = self._g[1, i, :, j, k] # g2 - g3 713 for l in range(len(interaction)): 714 self._imag_self_energy[:, l] += ( 715 (n2 + n3 + 1) * g1 + 716 (n2 - n3) * (g2_g3)) * interaction[l, j, k] * w 717 718 self._imag_self_energy *= self._unit_conversion 719 720 def _ise_thm_with_frequency_points_0K(self): 721 for i, (w, interaction) in enumerate(zip(self._weights_at_q, 722 self._pp_strength)): 723 for j, k in list(np.ndindex(interaction.shape[1:])): 724 g1 = self._g[0, i, :, j, k] 725 for l in range(len(interaction)): 726 self._imag_self_energy[:, l] += g1 * interaction[l, j, k] * w 727 728 self._imag_self_energy *= self._unit_conversion 729 730 def _average_by_degeneracy(self, imag_self_energy): 731 return average_by_degeneracy(imag_self_energy, 732 self._pp.band_indices, 733 self._frequencies[self._grid_point]) 734 735 736def run_ise_at_frequency_points_batch( 737 _frequency_points, 738 ise, 739 temperatures, 740 gamma, 741 write_gamma_detail=False, 742 return_gamma_detail=False, 743 detailed_gamma_at_gp=None, 744 scattering_event_class=None, 745 nelems_in_batch=50, 746 log_level=0): 747 if nelems_in_batch is None: 748 _nelems_in_batch = 10 749 else: 750 _nelems_in_batch = nelems_in_batch 751 752 batches = _get_batches(len(_frequency_points), _nelems_in_batch) 753 754 if log_level: 755 print("Calculations at %d frequency points are devided into " 756 "%d batches." % (len(_frequency_points), len(batches))) 757 758 for bi, fpts_batch in enumerate(batches): 759 if log_level: 760 print("%d/%d: %s" % (bi + 1, len(batches), fpts_batch + 1)) 761 sys.stdout.flush() 762 763 ise.set_frequency_points(_frequency_points[fpts_batch]) 764 ise.set_integration_weights( 765 scattering_event_class=scattering_event_class) 766 for l, t in enumerate(temperatures): 767 ise.set_temperature(t) 768 ise.run() 769 gamma[l, :, fpts_batch] = ise.get_imag_self_energy() 770 if write_gamma_detail or return_gamma_detail: 771 detailed_gamma_at_gp[l, fpts_batch] = ( 772 ise.get_detailed_imag_self_energy()) 773 774 775def _get_batches(tot_nelems, nelems=10): 776 nbatch = tot_nelems // nelems 777 batches = [np.arange(i * nelems, (i + 1) * nelems) 778 for i in range(nbatch)] 779 if tot_nelems % nelems > 0: 780 batches.append(np.arange(nelems * nbatch, tot_nelems)) 781 return batches 782