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 numpy as np 36import spglib 37from phonopy.structure.symmetry import Symmetry 38from phonopy.structure.tetrahedron_method import TetrahedronMethod 39from phonopy.structure.grid_points import extract_ir_grid_points 40from phono3py.phonon.func import gaussian 41 42 43def get_triplets_at_q(grid_point, 44 mesh, 45 point_group, # real space point group of space group 46 reciprocal_lattice, # column vectors 47 is_time_reversal=True, 48 swappable=True, 49 stores_triplets_map=False): 50 """Parameters 51 ---------- 52 grid_point : int 53 A grid point 54 mesh : array_like 55 Mesh numbers 56 shape=(3,), dtype='intc' 57 point_group : array_like 58 Rotation matrices in real space. Note that those in reciprocal space 59 mean these matrices transposed (local terminology). 60 shape=(n_rot, 3, 3), dtype='intc', order='C' 61 reciprocal_lattice : array_like 62 Reciprocal primitive basis vectors given as column vectors 63 shape=(3, 3), dtype='double', order='C' 64 is_time_reversal : bool, optional 65 Inversion symemtry is added if it doesn't exist. Default is True. 66 swappable : bool, optional 67 q1 and q2 among (q0, q1, q2) can be swapped. Deafult is True. 68 69 Returns 70 ------- 71 triplets_at_q : ndarray 72 Symmetry reduced number of triplets are stored as grid point 73 integer numbers. 74 shape=(n_triplets, 3), dtype='uintp' 75 weights : ndarray 76 Weights of triplets in Brillouin zone 77 shape=(n_triplets,), dtype='intc' 78 bz_grid_address : ndarray 79 Integer grid address of the points in Brillouin zone including 80 surface. The first prod(mesh) numbers of points are 81 independent. But the rest of points are 82 translational-symmetrically equivalent to some other points. 83 shape=(n_grid_points, 3), dtype='intc', order='C' 84 bz_map : ndarray 85 Grid point mapping table containing BZ surface. See more 86 detail in spglib docstring. 87 shape=(prod(mesh*2),), dtype='uintp' 88 map_tripelts : ndarray or None 89 Returns when stores_triplets_map=True, otherwise None is 90 returned. Mapping table of all triplets to symmetrically 91 independent tripelts. More precisely, this gives a list of 92 index mapping from all q-points to independent q' of 93 q+q'+q''=G. Considering q' is enough because q is fixed and 94 q''=G-q-q' where G is automatically determined to choose 95 smallest |G|. 96 shape=(prod(mesh),), dtype='uintp' 97 map_q : ndarray or None 98 Returns when stores_triplets_map=True, otherwise None is 99 returned. Irreducible q-points stabilized by q-point of 100 specified grid_point. 101 shape=(prod(mesh),), dtype='uintp' 102 103 """ 104 105 map_triplets, map_q, grid_address = _get_triplets_reciprocal_mesh_at_q( 106 grid_point, 107 mesh, 108 point_group, 109 is_time_reversal=is_time_reversal, 110 swappable=swappable) 111 bz_grid_address, bz_map = spglib.relocate_BZ_grid_address( 112 grid_address, 113 mesh, 114 reciprocal_lattice, 115 is_dense=True) 116 triplets_at_q, weights = _get_BZ_triplets_at_q( 117 grid_point, 118 bz_grid_address, 119 bz_map, 120 map_triplets, 121 mesh) 122 123 assert np.prod(mesh) == weights.sum(), \ 124 "Num grid points %d, sum of weight %d" % ( 125 np.prod(mesh), weights.sum()) 126 127 # These maps are required for collision matrix calculation. 128 if not stores_triplets_map: 129 map_triplets = None 130 map_q = None 131 132 return triplets_at_q, weights, bz_grid_address, bz_map, map_triplets, map_q 133 134 135def get_all_triplets(grid_point, 136 bz_grid_address, 137 bz_map, 138 mesh): 139 triplets_at_q, _ = _get_BZ_triplets_at_q( 140 grid_point, 141 bz_grid_address, 142 bz_map, 143 np.arange(np.prod(mesh), dtype=bz_map.dtype), 144 mesh) 145 146 return triplets_at_q 147 148 149def get_nosym_triplets_at_q(grid_point, 150 mesh, 151 reciprocal_lattice, 152 stores_triplets_map=False): 153 grid_address = get_grid_address(mesh) 154 bz_grid_address, bz_map = spglib.relocate_BZ_grid_address( 155 grid_address, 156 mesh, 157 reciprocal_lattice, 158 is_dense=True) 159 map_triplets = np.arange(len(grid_address), dtype=bz_map.dtype) 160 triplets_at_q, weights = _get_BZ_triplets_at_q( 161 grid_point, 162 bz_grid_address, 163 bz_map, 164 map_triplets, 165 mesh) 166 167 if not stores_triplets_map: 168 map_triplets = None 169 map_q = None 170 else: 171 map_q = map_triplets.copy() 172 173 return triplets_at_q, weights, bz_grid_address, bz_map, map_triplets, map_q 174 175 176def get_grid_address(mesh): 177 grid_mapping_table, grid_address = spglib.get_stabilized_reciprocal_mesh( 178 mesh, 179 [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]], 180 is_time_reversal=False, 181 is_dense=True) 182 183 return grid_address 184 185 186def get_bz_grid_address(mesh, reciprocal_lattice, with_boundary=False): 187 grid_address = get_grid_address(mesh) 188 bz_grid_address, bz_map = spglib.relocate_BZ_grid_address( 189 grid_address, 190 mesh, 191 reciprocal_lattice, 192 is_dense=True) 193 if with_boundary: 194 return bz_grid_address, bz_map 195 else: 196 return bz_grid_address[:np.prod(mesh)] 197 198 199def get_grid_point_from_address_py(address, mesh): 200 # X runs first in XYZ 201 # (*In spglib, Z first is possible with MACRO setting.) 202 m = mesh 203 return (address[0] % m[0] + 204 (address[1] % m[1]) * m[0] + 205 (address[2] % m[2]) * m[0] * m[1]) 206 207 208def get_grid_point_from_address(address, mesh): 209 """Grid point number is given by grid address. 210 211 Parameters 212 ---------- 213 address : array_like 214 Grid address. 215 shape=(3,), dtype='intc' 216 mesh : array_like 217 Mesh numbers. 218 shape=(3,), dtype='intc' 219 220 Returns 221 ------- 222 int 223 Grid point number. 224 225 """ 226 227 return spglib.get_grid_point_from_address(address, mesh) 228 229 230def get_ir_grid_points(mesh, rotations, mesh_shifts=None): 231 if mesh_shifts is None: 232 mesh_shifts = [False, False, False] 233 grid_mapping_table, grid_address = spglib.get_stabilized_reciprocal_mesh( 234 mesh, 235 rotations, 236 is_shift=np.where(mesh_shifts, 1, 0), 237 is_dense=True) 238 (ir_grid_points, 239 ir_grid_weights) = extract_ir_grid_points(grid_mapping_table) 240 241 return ir_grid_points, ir_grid_weights, grid_address, grid_mapping_table 242 243 244def get_grid_points_by_rotations(grid_point, 245 reciprocal_rotations, 246 mesh, 247 mesh_shifts=None): 248 if mesh_shifts is None: 249 mesh_shifts = [False, False, False] 250 return spglib.get_grid_points_by_rotations( 251 grid_point, 252 reciprocal_rotations, 253 mesh, 254 is_shift=np.where(mesh_shifts, 1, 0), 255 is_dense=True) 256 257 258def reduce_grid_points(mesh_divisors, 259 grid_address, 260 dense_grid_points, 261 dense_grid_weights=None, 262 coarse_mesh_shifts=None): 263 divisors = np.array(mesh_divisors, dtype='intc') 264 if (divisors == 1).all(): 265 coarse_grid_points = np.array(dense_grid_points, dtype='uintp') 266 if dense_grid_weights is not None: 267 coarse_grid_weights = np.array(dense_grid_weights, dtype='intc') 268 else: 269 if coarse_mesh_shifts is None: 270 shift = [0, 0, 0] 271 else: 272 shift = np.where(coarse_mesh_shifts, divisors // 2, [0, 0, 0]) 273 modulo = grid_address[dense_grid_points] % divisors 274 condition = (modulo == shift).all(axis=1) 275 coarse_grid_points = np.extract(condition, dense_grid_points) 276 if dense_grid_weights is not None: 277 coarse_grid_weights = np.extract(condition, dense_grid_weights) 278 279 if dense_grid_weights is None: 280 return coarse_grid_points 281 else: 282 return coarse_grid_points, coarse_grid_weights 283 284 285def from_coarse_to_dense_grid_points(dense_mesh, 286 mesh_divisors, 287 coarse_grid_points, 288 coarse_grid_address, 289 coarse_mesh_shifts=None): 290 if coarse_mesh_shifts is None: 291 coarse_mesh_shifts = [False, False, False] 292 shifts = np.where(coarse_mesh_shifts, 1, 0) 293 dense_grid_points = [] 294 for cga in coarse_grid_address[coarse_grid_points]: 295 dense_address = cga * mesh_divisors + shifts * (mesh_divisors // 2) 296 dense_grid_points.append(get_grid_point_from_address(dense_address, 297 dense_mesh)) 298 return np.array(dense_grid_points, dtype='uintp') 299 300 301def get_coarse_ir_grid_points(primitive, 302 mesh, 303 mesh_divisors, 304 coarse_mesh_shifts, 305 is_kappa_star=True, 306 symprec=1e-5): 307 mesh = np.array(mesh, dtype='intc') 308 309 symmetry = Symmetry(primitive, symprec) 310 point_group = symmetry.get_pointgroup_operations() 311 312 if mesh_divisors is None: 313 (ir_grid_points, 314 ir_grid_weights, 315 grid_address, 316 grid_mapping_table) = get_ir_grid_points(mesh, point_group) 317 else: 318 mesh_divs = np.array(mesh_divisors, dtype='intc') 319 coarse_mesh = mesh // mesh_divs 320 if coarse_mesh_shifts is None: 321 coarse_mesh_shifts = [False, False, False] 322 323 if not is_kappa_star: 324 coarse_grid_address = get_grid_address(coarse_mesh) 325 coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='uintp') 326 else: 327 (coarse_ir_grid_points, 328 coarse_ir_grid_weights, 329 coarse_grid_address, 330 coarse_grid_mapping_table) = get_ir_grid_points( 331 coarse_mesh, 332 point_group, 333 mesh_shifts=coarse_mesh_shifts) 334 ir_grid_points = from_coarse_to_dense_grid_points( 335 mesh, 336 mesh_divs, 337 coarse_grid_points, 338 coarse_grid_address, 339 coarse_mesh_shifts=coarse_mesh_shifts) 340 grid_address = get_grid_address(mesh) 341 ir_grid_weights = ir_grid_weights 342 343 reciprocal_lattice = np.linalg.inv(primitive.get_cell()) 344 bz_grid_address, bz_map = spglib.relocate_BZ_grid_address( 345 grid_address, 346 mesh, 347 reciprocal_lattice, 348 is_dense=True) 349 350 return (ir_grid_points, 351 ir_grid_weights, 352 bz_grid_address, 353 grid_mapping_table) 354 355 356def get_number_of_triplets(primitive, 357 mesh, 358 grid_point, 359 swappable=True, 360 symprec=1e-5): 361 mesh = np.array(mesh, dtype='intc') 362 symmetry = Symmetry(primitive, symprec) 363 point_group = symmetry.get_pointgroup_operations() 364 reciprocal_lattice = np.linalg.inv(primitive.get_cell()) 365 triplets_at_q, _, _, _, _, _ = get_triplets_at_q( 366 grid_point, 367 mesh, 368 point_group, 369 reciprocal_lattice, 370 swappable=swappable) 371 372 return len(triplets_at_q) 373 374 375def get_triplets_integration_weights(interaction, 376 frequency_points, 377 sigma, 378 sigma_cutoff=None, 379 is_collision_matrix=False, 380 neighboring_phonons=False, 381 lang='C'): 382 triplets = interaction.get_triplets_at_q()[0] 383 frequencies = interaction.get_phonons()[0] 384 num_band = frequencies.shape[1] 385 g_zero = None 386 387 if is_collision_matrix: 388 g = np.empty( 389 (3, len(triplets), len(frequency_points), num_band, num_band), 390 dtype='double', order='C') 391 else: 392 g = np.empty( 393 (2, len(triplets), len(frequency_points), num_band, num_band), 394 dtype='double', order='C') 395 g[:] = 0 396 397 if sigma: 398 if lang == 'C': 399 import phono3py._phono3py as phono3c 400 g_zero = np.zeros(g.shape[1:], dtype='byte', order='C') 401 if sigma_cutoff is None: 402 cutoff = -1 403 else: 404 cutoff = float(sigma_cutoff) 405 # cutoff < 0 disables g_zero feature. 406 phono3c.triplets_integration_weights_with_sigma( 407 g, 408 g_zero, 409 frequency_points, 410 triplets, 411 frequencies, 412 sigma, 413 cutoff) 414 else: 415 for i, tp in enumerate(triplets): 416 f1s = frequencies[tp[1]] 417 f2s = frequencies[tp[2]] 418 for j, k in list(np.ndindex((num_band, num_band))): 419 f1 = f1s[j] 420 f2 = f2s[k] 421 g0 = gaussian(frequency_points - f1 - f2, sigma) 422 g[0, i, :, j, k] = g0 423 g1 = gaussian(frequency_points + f1 - f2, sigma) 424 g2 = gaussian(frequency_points - f1 + f2, sigma) 425 g[1, i, :, j, k] = g1 - g2 426 if len(g) == 3: 427 g[2, i, :, j, k] = g0 + g1 + g2 428 else: 429 if lang == 'C': 430 g_zero = np.zeros(g.shape[1:], dtype='byte', order='C') 431 _set_triplets_integration_weights_c( 432 g, 433 g_zero, 434 interaction, 435 frequency_points, 436 neighboring_phonons=neighboring_phonons) 437 else: 438 _set_triplets_integration_weights_py( 439 g, interaction, frequency_points) 440 441 return g, g_zero 442 443 444def get_tetrahedra_vertices(relative_address, 445 mesh, 446 triplets_at_q, 447 bz_grid_address, 448 bz_map): 449 bzmesh = mesh * 2 450 grid_order = [1, mesh[0], mesh[0] * mesh[1]] 451 bz_grid_order = [1, bzmesh[0], bzmesh[0] * bzmesh[1]] 452 num_triplets = len(triplets_at_q) 453 vertices = np.zeros((num_triplets, 2, 24, 4), dtype='uintp') 454 for i, tp in enumerate(triplets_at_q): 455 for j, adrs_shift in enumerate( 456 (relative_address, -relative_address)): 457 adrs = bz_grid_address[tp[j + 1]] + adrs_shift 458 bz_gp = np.dot(adrs % bzmesh, bz_grid_order) 459 gp = np.dot(adrs % mesh, grid_order) 460 vgp = bz_map[bz_gp] 461 vertices[i, j] = vgp + (vgp == -1) * (gp + 1) 462 return vertices 463 464 465def _get_triplets_reciprocal_mesh_at_q(fixed_grid_number, 466 mesh, 467 rotations, 468 is_time_reversal=True, 469 swappable=True): 470 """Search symmetry reduced triplets fixing one q-point 471 472 Triplets of (q0, q1, q2) are searched. 473 474 Parameters 475 ---------- 476 fixed_grid_number : int 477 Grid point of q0 478 mesh : array_like 479 Mesh numbers 480 dtype='intc' 481 shape=(3,) 482 rotations : array_like 483 Rotation matrices in real space. Note that those in reciprocal space 484 mean these matrices transposed (local terminology). 485 dtype='intc' 486 shape=(n_rot, 3, 3) 487 is_time_reversal : bool 488 Inversion symemtry is added if it doesn't exist. 489 swappable : bool 490 q1 and q2 can be swapped. By this number of triplets decreases. 491 492 """ 493 494 import phono3py._phono3py as phono3c 495 496 map_triplets = np.zeros(np.prod(mesh), dtype='uintp') 497 map_q = np.zeros(np.prod(mesh), dtype='uintp') 498 grid_address = np.zeros((np.prod(mesh), 3), dtype='intc') 499 500 phono3c.triplets_reciprocal_mesh_at_q( 501 map_triplets, 502 map_q, 503 grid_address, 504 fixed_grid_number, 505 np.array(mesh, dtype='intc'), 506 is_time_reversal * 1, 507 np.array(rotations, dtype='intc', order='C'), 508 swappable * 1) 509 510 return map_triplets, map_q, grid_address 511 512 513def _get_BZ_triplets_at_q(grid_point, 514 bz_grid_address, 515 bz_map, 516 map_triplets, 517 mesh): 518 import phono3py._phono3py as phono3c 519 520 weights = np.zeros(len(map_triplets), dtype='intc') 521 for g in map_triplets: 522 weights[g] += 1 523 ir_weights = np.extract(weights > 0, weights) 524 triplets = np.zeros((len(ir_weights), 3), dtype=bz_map.dtype) 525 # triplets are overwritten. 526 num_ir_ret = phono3c.BZ_triplets_at_q(triplets, 527 grid_point, 528 bz_grid_address, 529 bz_map, 530 map_triplets, 531 np.array(mesh, dtype='intc')) 532 assert num_ir_ret == len(ir_weights) 533 534 return triplets, np.array(ir_weights, dtype='intc') 535 536 537def _set_triplets_integration_weights_c(g, 538 g_zero, 539 interaction, 540 frequency_points, 541 neighboring_phonons=False): 542 import phono3py._phono3py as phono3c 543 544 reciprocal_lattice = np.linalg.inv(interaction.primitive.cell) 545 mesh = interaction.mesh_numbers 546 thm = TetrahedronMethod(reciprocal_lattice, mesh=mesh) 547 grid_address = interaction.grid_address 548 bz_map = interaction.bz_map 549 triplets_at_q = interaction.get_triplets_at_q()[0] 550 551 if neighboring_phonons: 552 unique_vertices = thm.get_unique_tetrahedra_vertices() 553 for i, j in zip((1, 2), (1, -1)): 554 neighboring_grid_points = np.zeros( 555 len(unique_vertices) * len(triplets_at_q), dtype=bz_map.dtype) 556 phono3c.neighboring_grid_points( 557 neighboring_grid_points, 558 np.array(triplets_at_q[:, i], dtype='uintp').ravel(), 559 np.array(j * unique_vertices, dtype='intc', order='C'), 560 mesh, 561 grid_address, 562 bz_map) 563 interaction.run_phonon_solver(np.unique(neighboring_grid_points)) 564 565 frequencies = interaction.get_phonons()[0] 566 phono3c.triplets_integration_weights( 567 g, 568 g_zero, 569 frequency_points, # f0 570 np.array(thm.get_tetrahedra(), dtype='intc', order='C'), 571 mesh, 572 triplets_at_q, 573 frequencies, # f1 574 frequencies, # f2 575 grid_address, 576 bz_map, 577 g.shape[0]) 578 579 580def _set_triplets_integration_weights_py(g, interaction, frequency_points): 581 reciprocal_lattice = np.linalg.inv(interaction.get_primitive().get_cell()) 582 mesh = interaction.get_mesh_numbers() 583 thm = TetrahedronMethod(reciprocal_lattice, mesh=mesh) 584 grid_address = interaction.get_grid_address() 585 bz_map = interaction.get_bz_map() 586 triplets_at_q = interaction.get_triplets_at_q()[0] 587 tetrahedra_vertices = get_tetrahedra_vertices( 588 thm.get_tetrahedra(), 589 mesh, 590 triplets_at_q, 591 grid_address, 592 bz_map) 593 interaction.run_phonon_solver(np.unique(tetrahedra_vertices)) 594 frequencies = interaction.get_phonons()[0] 595 num_band = frequencies.shape[1] 596 for i, vertices in enumerate(tetrahedra_vertices): 597 for j, k in list(np.ndindex((num_band, num_band))): 598 f1_v = frequencies[vertices[0], j] 599 f2_v = frequencies[vertices[1], k] 600 thm.set_tetrahedra_omegas(f1_v + f2_v) 601 thm.run(frequency_points) 602 g0 = thm.get_integration_weight() 603 g[0, i, :, j, k] = g0 604 thm.set_tetrahedra_omegas(-f1_v + f2_v) 605 thm.run(frequency_points) 606 g1 = thm.get_integration_weight() 607 thm.set_tetrahedra_omegas(f1_v - f2_v) 608 thm.run(frequency_points) 609 g2 = thm.get_integration_weight() 610 g[1, i, :, j, k] = g1 - g2 611 if len(g) == 3: 612 g[2, i, :, j, k] = g0 + g1 + g2 613