1#!/usr/bin/env python 2 3import sys 4import numpy as np 5from phonopy.phonon.tetrahedron_mesh import TetrahedronMesh 6from phonopy.harmonic.force_constants import similarity_transformation 7from phonopy.interface.calculator import read_crystal_structure 8from phono3py.phonon3.triplets import (get_ir_grid_points, 9 get_grid_points_by_rotations, 10 get_grid_address) 11from phonopy.phonon.dos import NormalDistribution 12 13epsilon = 1.0e-8 14 15 16def show_tensor(kdos, temperatures, sampling_points, args): 17 for i, kdos_t in enumerate(kdos): 18 if not args.gv: 19 print("# %d K" % temperatures[i]) 20 21 for f, k in zip(sampling_points[i], kdos_t): # show kappa_xx 22 if args.average: 23 print(("%13.5f " * 3) % 24 (f, k[0][:3].sum() / 3, k[1][:3].sum() / 3)) 25 elif args.trace: 26 print(("%13.5f " * 3) % (f, k[0][:3].sum(), k[1][:3].sum())) 27 else: 28 print(("%f " * 13) % ((f,) + tuple(k[0]) + tuple(k[1]))) 29 30 print('') 31 print('') 32 33 34def show_scalar(gdos, temperatures, sampling_points, args): 35 if args.pqj or args.gruneisen or args.gv_norm: 36 for f, g in zip(sampling_points, gdos[0]): 37 print("%f %e %e" % (f, g[0], g[1])) 38 else: 39 for i, gdos_t in enumerate(gdos): 40 print("# %d K" % temperatures[i]) 41 for f, g in zip(sampling_points, gdos_t): 42 print("%f %f %f" % (f, g[0], g[1])) 43 print('') 44 print('') 45 46 47def set_T_target(temperatures, 48 mode_prop, 49 T_target, 50 mean_freepath=None): 51 for i, t in enumerate(temperatures): 52 if np.abs(t - args.temperature) < epsilon: 53 temperatures = temperatures[i:i+1] 54 mode_prop = mode_prop[i:i+1, :, :] 55 if mean_freepath is not None: 56 mean_freepath = mean_freepath[i:i+1] 57 return temperatures, mode_prop, mean_freepath 58 else: 59 return temperatures, mode_prop 60 61 62def run_prop_dos(frequencies, 63 mode_prop, 64 primitive, 65 mesh, 66 grid_address, 67 grid_mapping_table, 68 ir_grid_points, 69 num_sampling_points): 70 kappa_dos = KappaDOS(mode_prop, 71 primitive, 72 frequencies, 73 mesh, 74 grid_address, 75 grid_mapping_table, 76 ir_grid_points, 77 num_sampling_points=num_sampling_points) 78 freq_points, kdos = kappa_dos.get_kdos() 79 sampling_points = np.tile(freq_points, (len(kdos), 1)) 80 81 return kdos, sampling_points 82 83 84def get_mfp(g, gv): 85 g = np.where(g > 0, g, -1) 86 gv_norm = np.sqrt((gv ** 2).sum(axis=2)) 87 mean_freepath = np.where(g > 0, gv_norm / (2 * 2 * np.pi * g), 0) 88 return mean_freepath 89 90 91def run_mfp_dos(mean_freepath, 92 mode_prop, 93 primitive, 94 mesh, 95 grid_address, 96 grid_mapping_table, 97 ir_grid_points, 98 num_sampling_points): 99 kdos = [] 100 sampling_points = [] 101 for i, mfp in enumerate(mean_freepath): 102 kappa_dos = KappaDOS(mode_prop[i:i+1, :, :], 103 primitive, 104 mfp, 105 mesh, 106 grid_address, 107 grid_mapping_table, 108 ir_grid_points, 109 num_sampling_points=num_sampling_points) 110 sampling_points_at_T, kdos_at_T = kappa_dos.get_kdos() 111 kdos.append(kdos_at_T[0]) 112 sampling_points.append(sampling_points_at_T) 113 kdos = np.array(kdos) 114 sampling_points = np.array(sampling_points) 115 116 return kdos, sampling_points 117 118 119def fracval(frac): 120 if frac.find('/') == -1: 121 return float(frac) 122 else: 123 x = frac.split('/') 124 return float(x[0]) / float(x[1]) 125 126 127def get_grid_symmetry(primitive, mesh, weights, qpoints): 128 symmetry = Symmetry(primitive) 129 rotations = symmetry.get_pointgroup_operations() 130 (ir_grid_points, 131 weights_for_check, 132 grid_address, 133 grid_mapping_table) = get_ir_grid_points(mesh, rotations) 134 135 try: 136 np.testing.assert_array_equal(weights, weights_for_check) 137 except: 138 print("*******************************") 139 print("** Might forget --pa option? **") 140 print("*******************************") 141 raise 142 143 qpoints_for_check = grid_address[ir_grid_points] / mesh.astype('double') 144 diff_q = qpoints - qpoints_for_check 145 np.testing.assert_almost_equal(diff_q, np.rint(diff_q)) 146 147 return ir_grid_points, grid_address, grid_mapping_table 148 149 150def get_gv_by_gv(gv, 151 symmetry, 152 primitive, 153 mesh, 154 grid_points, 155 grid_address): 156 point_operations = symmetry.get_reciprocal_operations() 157 rec_lat = np.linalg.inv(primitive.get_cell()) 158 rotations_cartesian = np.array( 159 [similarity_transformation(rec_lat, r) 160 for r in point_operations], dtype='double') 161 162 num_band = gv.shape[1] 163 gv_sum2 = np.zeros((gv.shape[0], num_band, 6), dtype='double') 164 for i, gp in enumerate(grid_points): 165 rotation_map = get_grid_points_by_rotations( 166 grid_address[gp], 167 point_operations, 168 mesh) 169 gv_by_gv = np.zeros((num_band, 3, 3), dtype='double') 170 for r in rotations_cartesian: 171 gvs_rot = np.dot(gv[i], r.T) 172 gv_by_gv += [np.outer(r_gv, r_gv) for r_gv in gvs_rot] 173 gv_by_gv /= len(rotation_map) // len(np.unique(rotation_map)) 174 for j, vxv in enumerate( 175 ([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])): 176 gv_sum2[i, :, j] = gv_by_gv[:, vxv[0], vxv[1]] 177 178 return gv_sum2 179 180 181class KappaDOS(object): 182 def __init__(self, 183 mode_kappa, 184 cell, 185 frequencies, 186 mesh, 187 grid_address, 188 grid_mapping_table, 189 ir_grid_points, 190 grid_order=None, 191 num_sampling_points=100): 192 self._mode_kappa = mode_kappa 193 self._tetrahedron_mesh = TetrahedronMesh( 194 cell, 195 frequencies, 196 mesh, 197 grid_address, 198 grid_mapping_table, 199 ir_grid_points) 200 201 min_freq = min(frequencies.ravel()) 202 max_freq = max(frequencies.ravel()) + epsilon 203 self._frequency_points = np.linspace(min_freq, 204 max_freq, 205 num_sampling_points) 206 self._kdos = np.zeros( 207 (len(mode_kappa), len(self._frequency_points), 2, 6), 208 dtype='double') 209 self._run_tetrahedron_method() 210 211 def get_kdos(self): 212 return self._frequency_points, self._kdos 213 214 def _run_tetrahedron_method(self): 215 thm = self._tetrahedron_mesh 216 for j, value in enumerate(('J', 'I')): 217 thm.set(value=value, frequency_points=self._frequency_points) 218 for i, iw in enumerate(thm): 219 # kdos[temp, freq_points, IJ, tensor_elem] 220 # iw[freq_points, band] 221 # mode_kappa[temp, ir_gp, band, tensor_elem] 222 self._kdos[:, :, j] += np.transpose( 223 np.dot(iw, self._mode_kappa[:, i]), axes=(1, 0, 2)) 224 225 226class GammaDOS(object): 227 def __init__(self, 228 gamma, 229 frequencies, 230 ir_grid_weights, 231 num_fpoints=200): 232 self._gamma = gamma 233 self._frequencies = frequencies 234 self._ir_grid_weights = ir_grid_weights 235 self._num_fpoints = num_fpoints 236 self._set_frequency_points() 237 self._gdos = np.zeros( 238 (len(gamma), len(self._frequency_points), 2), dtype='double') 239 240 def get_gdos(self): 241 return self._frequency_points, self._gdos 242 243 def _set_frequency_points(self): 244 min_freq = np.min(self._frequencies) 245 max_freq = np.max(self._frequencies) + epsilon 246 self._frequency_points = np.linspace(min_freq, 247 max_freq, 248 self._num_fpoints) 249 250 251class GammaDOSsmearing(GammaDOS): 252 def __init__(self, 253 gamma, 254 frequencies, 255 ir_grid_weights, 256 sigma=None, 257 num_fpoints=200): 258 GammaDOS.__init__(self, 259 gamma, 260 frequencies, 261 ir_grid_weights, 262 num_fpoints=num_fpoints) 263 if sigma is None: 264 self._sigma = (max(self._frequency_points) - 265 min(self._frequency_points)) / 100 266 else: 267 self._sigma = 0.1 268 self._smearing_function = NormalDistribution(self._sigma) 269 self._run_smearing_method() 270 271 def _run_smearing_method(self): 272 self._dos = [] 273 num_gp = np.sum(self._ir_grid_weights) 274 for i, f in enumerate(self._frequency_points): 275 dos = self._smearing_function.calc(self._frequencies - f) 276 for j, g_t in enumerate(self._gamma): 277 self._gdos[j, i, 1] = np.sum(np.dot(self._ir_grid_weights, 278 dos * g_t)) / num_gp 279 280 281class GammaDOStetrahedron(GammaDOS): 282 def __init__(self, 283 gamma, 284 cell, 285 frequencies, 286 mesh, 287 grid_address, 288 grid_mapping_table, 289 ir_grid_points, 290 ir_grid_weights, 291 num_fpoints=200): 292 GammaDOS.__init__(self, 293 gamma, 294 frequencies, 295 ir_grid_weights, 296 num_fpoints=num_fpoints) 297 self._cell = cell 298 self._mesh = mesh 299 self._grid_address = grid_address 300 self._grid_mapping_table = grid_mapping_table 301 self._ir_grid_points = ir_grid_points 302 303 self._set_tetrahedron_method() 304 self._run_tetrahedron_method() 305 306 def _set_tetrahedron_method(self): 307 self._tetrahedron_mesh = TetrahedronMesh( 308 self._cell, 309 self._frequencies, 310 self._mesh, 311 self._grid_address, 312 self._grid_mapping_table, 313 self._ir_grid_points) 314 315 def _run_tetrahedron_method(self): 316 thm = self._tetrahedron_mesh 317 for j, value in enumerate(('J', 'I')): 318 thm.set(value=value, frequency_points=self._frequency_points) 319 for i, iw in enumerate(thm): 320 # gdos[temp, freq_points, IJ] 321 # iw[freq_points, band] 322 # gamma[temp, ir_gp, band] 323 self._gdos[:, :, j] += np.dot( 324 self._gamma[:, i] * self._ir_grid_weights[i], iw.T) 325 326 327if __name__ == '__main__': 328 """Incremental kappa with respect to frequency and the derivative""" 329 330 import h5py 331 from phonopy.structure.cells import get_primitive 332 from phonopy.structure.symmetry import Symmetry 333 import argparse 334 335 parser = argparse.ArgumentParser(description="Show unit cell volume") 336 parser.add_argument( 337 "--pa", dest="primitive_matrix", default="1 0 0 0 1 0 0 0 1", 338 help="Primitive matrix") 339 parser.add_argument( 340 "--mesh", dest="mesh", default="1 1 1", 341 help="Mesh numbers") 342 parser.add_argument( 343 "-c", "--cell", dest="cell_filename", 344 help="Unit cell filename") 345 parser.add_argument( 346 '--gv', action='store_true', 347 help='Calculate for gv_x_gv (tensor)') 348 parser.add_argument( 349 '--pqj', action='store_true', 350 help='Calculate for Pqj (scalar)') 351 parser.add_argument( 352 '--cv', action='store_true', 353 help='Calculate for Cv (scalar)') 354 parser.add_argument( 355 '--tau', action='store_true', 356 help='Calculate for lifetimes (scalar)') 357 parser.add_argument( 358 '--gamma', action='store_true', 359 help='Calculate for Gamma (scalar)') 360 parser.add_argument( 361 '--gruneisen', action='store_true', 362 help='Calculate for mode-Gruneisen parameters squared (scalar)') 363 parser.add_argument( 364 '--gv-norm', action='store_true', 365 help='Calculate for |g_v| (scalar)') 366 parser.add_argument( 367 '--mfp', action='store_true', 368 help='Mean free path is used instead of frequency') 369 parser.add_argument( 370 '--temperature', type=float, dest='temperature', 371 help='Temperature to output data at') 372 parser.add_argument( 373 '--nsp', '--num-sampling-points', type=int, dest='num_sampling_points', 374 default=100, 375 help="Number of sampling points in frequency or MFP axis") 376 parser.add_argument( 377 '--average', action='store_true', 378 help=("Output the traces of the tensors divided by 3 " 379 "rather than the unique elements")) 380 parser.add_argument( 381 '--trace', action='store_true', 382 help=("Output the traces of the tensors " 383 "rather than the unique elements")) 384 parser.add_argument( 385 '--smearing', action='store_true', 386 help='Use smearing method (only for scalar density)') 387 parser.add_argument( 388 '--qe', '--pwscf', dest="qe_mode", 389 action="store_true", help="Invoke Pwscf mode") 390 parser.add_argument( 391 '--crystal', dest="crystal_mode", 392 action="store_true", help="Invoke CRYSTAL mode") 393 parser.add_argument( 394 '--abinit', dest="abinit_mode", 395 action="store_true", help="Invoke Abinit mode") 396 parser.add_argument( 397 '--turbomole', dest="turbomole_mode", 398 action="store_true", help="Invoke TURBOMOLE mode") 399 parser.add_argument( 400 "--noks", "--no-kappa-stars", 401 dest="no_kappa_stars", action="store_true", 402 help="Deactivate summation of partial kappa at q-stars") 403 parser.add_argument('filenames', nargs='*') 404 args = parser.parse_args() 405 406 interface_mode = None 407 if args.qe_mode: 408 interface_mode = 'qe' 409 elif args.crystal_mode: 410 interface_mode = 'crystal' 411 elif args.abinit_mode: 412 interface_mode = 'abinit' 413 elif args.turbomole_mode: 414 interface_mode = 'turbomole' 415 if len(args.filenames) > 1: 416 cell, _ = read_crystal_structure(args.filenames[0], 417 interface_mode=interface_mode) 418 f = h5py.File(args.filenames[1], 'r') 419 else: 420 cell, _ = read_crystal_structure(args.cell_filename, 421 interface_mode=interface_mode) 422 f = h5py.File(args.filenames[0], 'r') 423 424 primitive_matrix = np.reshape( 425 [fracval(x) for x in args.primitive_matrix.split()], (3, 3)) 426 primitive = get_primitive(cell, primitive_matrix) 427 428 if 'mesh' in f: 429 mesh = np.array(f['mesh'][:], dtype='intc') 430 else: 431 mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc') 432 433 if 'temperature' in f: 434 temperatures = f['temperature'][:] 435 else: 436 temperatures = None 437 weights = f['weight'][:] 438 439 if args.no_kappa_stars or (weights == 1).all(): 440 grid_address = get_grid_address(mesh) 441 ir_grid_points = np.arange(np.prod(mesh), dtype='uintp') 442 grid_mapping_table = np.arange(np.prod(mesh), dtype='uintp') 443 else: 444 (ir_grid_points, 445 grid_address, 446 grid_mapping_table) = get_grid_symmetry(primitive, 447 mesh, 448 weights, 449 f['qpoint'][:]) 450 451 ################ 452 # Set property # 453 ################ 454 if args.gv: 455 if 'gv_by_gv' in f: 456 gv_sum2 = f['gv_by_gv'][:] 457 else: # For backward compatibility. This will be removed someday. 458 gv = f['group_velocity'][:] 459 symmetry = Symmetry(primitive) 460 gv_sum2 = get_gv_by_gv(gv, 461 symmetry, 462 primitive, 463 mesh, 464 ir_grid_points, 465 grid_address) 466 467 # gv x gv is divied by primitive cell volume. 468 unit_conversion = primitive.get_volume() 469 mode_prop = gv_sum2.reshape((1,) + gv_sum2.shape) / unit_conversion 470 else: 471 if 'mode_kappa' in f: 472 mode_prop = f['mode_kappa'][:] 473 else: 474 mode_prop = None 475 476 frequencies = f['frequency'][:] 477 conditions = frequencies > 0 478 if np.logical_not(conditions).sum() > 3: 479 sys.stderr.write("# Imaginary frequencies are found. " 480 "They are set to be zero.\n") 481 frequencies = np.where(conditions, frequencies, 0) 482 483 ####### 484 # Run # 485 ####### 486 if (args.gamma or args.gruneisen or args.pqj or 487 args.cv or args.tau or args.gv_norm): 488 if args.pqj: 489 mode_prop = f['ave_pp'][:].reshape((1,) + f['ave_pp'].shape) 490 elif args.cv: 491 mode_prop = f['heat_capacity'][:] 492 elif args.tau: 493 g = f['gamma'][:] 494 g = np.where(g > 0, g, -1) 495 mode_prop = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0) # tau 496 elif args.gv_norm: 497 mode_prop = np.sqrt( 498 (f['group_velocity'][:, :, :] ** 2).sum(axis=2)) 499 mode_prop = mode_prop.reshape((1,) + mode_prop.shape) 500 elif args.gamma: 501 mode_prop = f['gamma'][:] 502 elif args.gruneisen: 503 mode_prop = f['gruneisen'][:].reshape((1,) + f['gruneisen'].shape) 504 mode_prop **= 2 505 else: 506 raise 507 if (args.temperature is not None and 508 not (args.gv_norm or args.pqj or args.gruneisen)): 509 temperatures, mode_prop = set_T_target(temperatures, 510 mode_prop, 511 args.temperature) 512 if args.smearing: 513 mode_prop_dos = GammaDOSsmearing( 514 mode_prop, 515 frequencies, 516 weights, 517 num_fpoints=args.num_sampling_points) 518 sampling_points, gdos = mode_prop_dos.get_gdos() 519 else: 520 mode_prop_dos = GammaDOStetrahedron( 521 mode_prop, 522 primitive, 523 frequencies, 524 mesh, 525 grid_address, 526 grid_mapping_table, 527 ir_grid_points, 528 weights, 529 num_fpoints=args.num_sampling_points) 530 sampling_points, gdos = mode_prop_dos.get_gdos() 531 for i, gdos_t in enumerate(gdos): 532 total = np.dot(weights, mode_prop[i]).sum() / weights.sum() 533 assert np.isclose(gdos_t[-1][0], total) 534 show_scalar(gdos, temperatures, sampling_points, args) 535 else: 536 if args.mfp: 537 if 'mean_free_path' in f: 538 mfp = f['mean_free_path'][:] 539 mean_freepath = np.sqrt((mfp ** 2).sum(axis=3)) 540 else: 541 mean_freepath = get_mfp(f['gamma'][:], f['group_velocity'][:]) 542 if args.temperature is not None: 543 (temperatures, 544 mode_prop, 545 mean_freepath) = set_T_target(temperatures, 546 mode_prop, 547 args.temperature, 548 mean_freepath=mean_freepath) 549 550 kdos, sampling_points = run_mfp_dos(mean_freepath, 551 mode_prop, 552 primitive, 553 mesh, 554 grid_address, 555 grid_mapping_table, 556 ir_grid_points, 557 args.num_sampling_points) 558 show_tensor(kdos, temperatures, sampling_points, args) 559 else: 560 if args.temperature is not None and not args.gv: 561 temperatures, mode_prop = set_T_target(temperatures, 562 mode_prop, 563 args.temperature) 564 kdos, sampling_points = run_prop_dos(frequencies, 565 mode_prop, 566 primitive, 567 mesh, 568 grid_address, 569 grid_mapping_table, 570 ir_grid_points, 571 args.num_sampling_points) 572 show_tensor(kdos, temperatures, sampling_points, args) 573