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 os 36import numpy as np 37import h5py 38 39from phonopy.file_IO import (write_force_constants_to_hdf5, 40 check_force_constants_indices, 41 get_cell_from_disp_yaml) 42from phonopy.cui.load_helper import read_force_constants_from_hdf5 43 44 45def write_cell_yaml(w, supercell): 46 w.write("lattice:\n") 47 for axis in supercell.get_cell(): 48 w.write("- [ %20.15f,%20.15f,%20.15f ]\n" % tuple(axis)) 49 symbols = supercell.get_chemical_symbols() 50 positions = supercell.get_scaled_positions() 51 w.write("atoms:\n") 52 for i, (s, v) in enumerate(zip(symbols, positions)): 53 w.write("- symbol: %-2s # %d\n" % (s, i+1)) 54 w.write(" position: [ %18.14f,%18.14f,%18.14f ]\n" % tuple(v)) 55 56 57def write_disp_fc3_yaml(dataset, supercell, filename='disp_fc3.yaml'): 58 w = open(filename, 'w') 59 w.write("natom: %d\n" % dataset['natom']) 60 61 num_first = len(dataset['first_atoms']) 62 w.write("num_first_displacements: %d\n" % num_first) 63 if 'cutoff_distance' in dataset: 64 w.write("cutoff_distance: %f\n" % dataset['cutoff_distance']) 65 66 num_second = 0 67 num_disp_files = 0 68 for d1 in dataset['first_atoms']: 69 num_disp_files += 1 70 num_second += len(d1['second_atoms']) 71 for d2 in d1['second_atoms']: 72 if 'included' in d2: 73 if d2['included']: 74 num_disp_files += 1 75 else: 76 num_disp_files += 1 77 78 w.write("num_second_displacements: %d\n" % num_second) 79 w.write("num_displacements_created: %d\n" % num_disp_files) 80 81 if 'duplicates' in dataset: 82 w.write("duplicates:\n") 83 for (i, j) in dataset['duplicates']: 84 w.write("- [ %d, %d ]\n" % (i, j)) 85 86 w.write("first_atoms:\n") 87 count1 = 0 88 count2 = len(dataset['first_atoms']) 89 for disp1 in dataset['first_atoms']: 90 count1 += 1 91 disp_cart1 = disp1['displacement'] 92 w.write("- number: %5d\n" % (disp1['number'] + 1)) 93 w.write(" displacement:\n") 94 w.write(" [%20.16f,%20.16f,%20.16f ] # %05d\n" % 95 (disp_cart1[0], disp_cart1[1], disp_cart1[2], count1)) 96 w.write(" displacement_id: %d\n" % count1) 97 w.write(" second_atoms:\n") 98 99 included = None 100 atom2_list = np.array([disp2['number'] 101 for disp2 in disp1['second_atoms']], dtype=int) 102 _, indices = np.unique(atom2_list, return_index=True) 103 for atom2 in atom2_list[indices]: 104 disp2_list = [] 105 for disp2 in disp1['second_atoms']: 106 if disp2['number'] == atom2: 107 disp2_list.append(disp2) 108 109 disp2 = disp2_list[0] 110 atom2 = disp2['number'] 111 if 'included' in disp2: 112 included = disp2['included'] 113 pair_distance = disp2['pair_distance'] 114 w.write(" - number: %5d\n" % (atom2 + 1)) 115 w.write(" distance: %f\n" % pair_distance) 116 if included is not None: 117 if included: 118 w.write(" included: %s\n" % "true") 119 else: 120 w.write(" included: %s\n" % "false") 121 w.write(" displacements:\n") 122 123 for disp2 in disp2_list: 124 count2 += 1 125 126 # Assert all disp2s belonging to same atom2 appear straight. 127 assert disp2['id'] == count2 128 129 disp_cart2 = disp2['displacement'] 130 w.write(" - [%20.16f,%20.16f,%20.16f ] # %05d\n" % 131 (disp_cart2[0], disp_cart2[1], disp_cart2[2], 132 count2)) 133 134 ids = ["%d" % disp2['id'] for disp2 in disp2_list] 135 w.write(" displacement_ids: [ %s ]\n" % ', '.join(ids)) 136 137 write_cell_yaml(w, supercell) 138 139 w.close() 140 141 return num_first + num_second, num_disp_files 142 143 144def write_disp_fc2_yaml(dataset, supercell, filename='disp_fc2.yaml'): 145 w = open(filename, 'w') 146 w.write("natom: %d\n" % dataset['natom']) 147 148 num_first = len(dataset['first_atoms']) 149 w.write("num_first_displacements: %d\n" % num_first) 150 w.write("first_atoms:\n") 151 for i, disp1 in enumerate(dataset['first_atoms']): 152 disp_cart1 = disp1['displacement'] 153 w.write("- number: %5d\n" % (disp1['number'] + 1)) 154 w.write(" displacement:\n") 155 w.write(" [%20.16f,%20.16f,%20.16f ] # %05d\n" % 156 (disp_cart1[0], disp_cart1[1], disp_cart1[2], i + 1)) 157 158 if supercell is not None: 159 write_cell_yaml(w, supercell) 160 161 w.close() 162 163 return num_first 164 165 166def write_FORCES_FC2(disp_dataset, 167 forces_fc2=None, 168 fp=None, 169 filename="FORCES_FC2"): 170 if fp is None: 171 w = open(filename, 'w') 172 else: 173 w = fp 174 175 for i, disp1 in enumerate(disp_dataset['first_atoms']): 176 w.write("# File: %-5d\n" % (i + 1)) 177 w.write("# %-5d " % (disp1['number'] + 1)) 178 w.write("%20.16f %20.16f %20.16f\n" % tuple(disp1['displacement'])) 179 if 'forces' in disp1 and forces_fc2 is None: 180 force_set = disp1['forces'] 181 else: 182 force_set = forces_fc2[i] 183 for forces in force_set: 184 w.write("%15.10f %15.10f %15.10f\n" % tuple(forces)) 185 186 187def write_FORCES_FC3(disp_dataset, 188 forces_fc3=None, 189 fp=None, 190 filename="FORCES_FC3"): 191 if fp is None: 192 w = open(filename, 'w') 193 else: 194 w = fp 195 196 natom = disp_dataset['natom'] 197 num_disp1 = len(disp_dataset['first_atoms']) 198 count = num_disp1 199 file_count = num_disp1 200 201 write_FORCES_FC2(disp_dataset, forces_fc2=forces_fc3, fp=w) 202 203 for i, disp1 in enumerate(disp_dataset['first_atoms']): 204 atom1 = disp1['number'] 205 for disp2 in disp1['second_atoms']: 206 atom2 = disp2['number'] 207 w.write("# File: %-5d\n" % (count + 1)) 208 w.write("# %-5d " % (atom1 + 1)) 209 w.write("%20.16f %20.16f %20.16f\n" % tuple(disp1['displacement'])) 210 w.write("# %-5d " % (atom2 + 1)) 211 w.write("%20.16f %20.16f %20.16f\n" % tuple(disp2['displacement'])) 212 213 # For supercell calculation reduction 214 included = True 215 if 'included' in disp2: 216 included = disp2['included'] 217 if included: 218 if 'forces' in disp2 and forces_fc3 is None: 219 force_set = disp2['forces'] 220 else: 221 force_set = forces_fc3[file_count] 222 for force in force_set: 223 w.write("%15.10f %15.10f %15.10f\n" % tuple(force)) 224 file_count += 1 225 else: 226 # for forces in forces_fc3[i]: 227 # w.write("%15.10f %15.10f %15.10f\n" % (tuple(forces))) 228 for j in range(natom): 229 w.write("%15.10f %15.10f %15.10f\n" % (0, 0, 0)) 230 count += 1 231 232 233def write_fc3_dat(force_constants_third, filename='fc3.dat'): 234 w = open(filename, 'w') 235 for i in range(force_constants_third.shape[0]): 236 for j in range(force_constants_third.shape[1]): 237 for k in range(force_constants_third.shape[2]): 238 tensor3 = force_constants_third[i, j, k] 239 w.write(" %d - %d - %d (%f)\n" % (i + 1, j + 1, k + 1, 240 np.abs(tensor3).sum())) 241 for tensor2 in tensor3: 242 for vec in tensor2: 243 w.write("%20.14f %20.14f %20.14f\n" % tuple(vec)) 244 w.write("\n") 245 246 247def write_fc3_to_hdf5(fc3, 248 filename='fc3.hdf5', 249 p2s_map=None, 250 compression="gzip"): 251 """Write third-order force constants in hdf5 format. 252 253 Parameters 254 ---------- 255 force_constants : ndarray 256 Force constants 257 shape=(n_satom, n_satom, n_satom, 3, 3, 3) or 258 (n_patom, n_satom, n_satom,3,3,3), dtype=double 259 filename : str 260 Filename to be used. 261 p2s_map : ndarray, optional 262 Primitive atom indices in supercell index system 263 shape=(n_patom,), dtype=intc 264 compression : str or int, optional 265 h5py's lossless compression filters (e.g., "gzip", "lzf"). None gives 266 no compression. See the detail at docstring of 267 h5py.Group.create_dataset. Default is "gzip". 268 269 """ 270 271 with h5py.File(filename, 'w') as w: 272 w.create_dataset('fc3', data=fc3, compression=compression) 273 if p2s_map is not None: 274 w.create_dataset('p2s_map', data=p2s_map) 275 276 277def read_fc3_from_hdf5(filename='fc3.hdf5', p2s_map=None): 278 with h5py.File(filename, 'r') as f: 279 fc3 = f['fc3'][:] 280 if 'p2s_map' in f: 281 p2s_map_in_file = f['p2s_map'][:] 282 check_force_constants_indices(fc3.shape[:2], 283 p2s_map_in_file, 284 p2s_map, 285 filename) 286 if fc3.dtype == np.double and fc3.flags.c_contiguous: 287 return fc3 288 else: 289 msg = ("%s has to be read by h5py as numpy ndarray of " 290 "dtype='double' and c_contiguous." % filename) 291 raise TypeError(msg) 292 return None 293 294 295def write_fc2_dat(force_constants, filename='fc2.dat'): 296 w = open(filename, 'w') 297 for i, fcs in enumerate(force_constants): 298 for j, fcb in enumerate(fcs): 299 w.write(" %d - %d\n" % (i+1, j+1)) 300 for vec in fcb: 301 w.write("%20.14f %20.14f %20.14f\n" % tuple(vec)) 302 w.write("\n") 303 304 305def write_fc2_to_hdf5(force_constants, 306 filename='fc2.hdf5', 307 p2s_map=None, 308 physical_unit=None, 309 compression="gzip"): 310 write_force_constants_to_hdf5(force_constants, 311 filename=filename, 312 p2s_map=p2s_map, 313 physical_unit=physical_unit, 314 compression=compression) 315 316 317def read_fc2_from_hdf5(filename='fc2.hdf5', 318 p2s_map=None): 319 return read_force_constants_from_hdf5(filename=filename, 320 p2s_map=p2s_map, 321 calculator='vasp') 322 323 324def write_triplets(triplets, 325 weights, 326 mesh, 327 grid_address, 328 grid_point=None, 329 filename=None): 330 triplets_filename = "triplets" 331 suffix = "-m%d%d%d" % tuple(mesh) 332 if grid_point is not None: 333 suffix += ("-g%d" % grid_point) 334 if filename is not None: 335 suffix += "." + filename 336 suffix += ".dat" 337 triplets_filename += suffix 338 w = open(triplets_filename, 'w') 339 for weight, g3 in zip(weights, triplets): 340 w.write("%4d " % weight) 341 for q3 in grid_address[g3]: 342 w.write("%4d %4d %4d " % tuple(q3)) 343 w.write("\n") 344 w.close() 345 346 347def write_grid_address(grid_address, mesh, filename=None): 348 grid_address_filename = "grid_address" 349 suffix = "-m%d%d%d" % tuple(mesh) 350 if filename is not None: 351 suffix += "." + filename 352 suffix += ".dat" 353 grid_address_filename += suffix 354 355 w = open(grid_address_filename, 'w') 356 w.write("# Grid addresses for %dx%dx%d mesh\n" % tuple(mesh)) 357 w.write("#%9s %8s %8s %8s %8s %8s %8s\n" % 358 ("index", "a", "b", "c", 359 ("a%%%d" % mesh[0]), ("b%%%d" % mesh[1]), ("c%%%d" % mesh[2]))) 360 for i, bz_q in enumerate(grid_address): 361 if i == np.prod(mesh): 362 w.write("#" + "-" * 78 + "\n") 363 q = bz_q % mesh 364 w.write("%10d %8d %8d %8d " % (i, bz_q[0], bz_q[1], bz_q[2])) 365 w.write("%8d %8d %8d\n" % tuple(q)) 366 367 return grid_address_filename 368 369 370def write_grid_address_to_hdf5(grid_address, 371 mesh, 372 grid_mapping_table, 373 compression="gzip", 374 filename=None): 375 suffix = _get_filename_suffix(mesh, filename=filename) 376 full_filename = "grid_address" + suffix + ".hdf5" 377 with h5py.File(full_filename, 'w') as w: 378 w.create_dataset('mesh', data=mesh) 379 w.create_dataset('grid_address', data=grid_address, 380 compression=compression) 381 w.create_dataset('grid_mapping_table', data=grid_mapping_table, 382 compression=compression) 383 return full_filename 384 return None 385 386 387def write_imag_self_energy_at_grid_point(gp, 388 band_indices, 389 mesh, 390 frequencies, 391 gammas, 392 sigma=None, 393 temperature=None, 394 scattering_event_class=None, 395 filename=None, 396 is_mesh_symmetry=True): 397 gammas_filename = "gammas" 398 gammas_filename += "-m%d%d%d-g%d-" % (mesh[0], mesh[1], mesh[2], gp) 399 if sigma is not None: 400 gammas_filename += ("s%f" % sigma).rstrip('0').rstrip(r'\.') + "-" 401 402 if temperature is not None: 403 gammas_filename += ("t%f" % temperature).rstrip('0').rstrip(r'\.') + "-" 404 405 for i in band_indices: 406 gammas_filename += "b%d" % (i + 1) 407 408 if scattering_event_class is not None: 409 gammas_filename += "-c%d" % scattering_event_class 410 411 if filename is not None: 412 gammas_filename += ".%s" % filename 413 elif not is_mesh_symmetry: 414 gammas_filename += ".nosym" 415 gammas_filename += ".dat" 416 417 w = open(gammas_filename, 'w') 418 for freq, g in zip(frequencies, gammas): 419 w.write("%15.7f %20.15e\n" % (freq, g)) 420 w.close() 421 422 return gammas_filename 423 424 425def write_joint_dos(gp, 426 mesh, 427 frequencies, 428 jdos, 429 sigma=None, 430 temperatures=None, 431 filename=None, 432 is_mesh_symmetry=True): 433 if temperatures is None: 434 return _write_joint_dos_at_t(gp, 435 mesh, 436 frequencies, 437 jdos, 438 sigma=sigma, 439 temperature=None, 440 filename=filename, 441 is_mesh_symmetry=is_mesh_symmetry) 442 else: 443 for jdos_at_t, t in zip(jdos, temperatures): 444 return _write_joint_dos_at_t(gp, 445 mesh, 446 frequencies, 447 jdos_at_t, 448 sigma=sigma, 449 temperature=t, 450 filename=filename, 451 is_mesh_symmetry=is_mesh_symmetry) 452 453 454def _write_joint_dos_at_t(grid_point, 455 mesh, 456 frequencies, 457 jdos, 458 sigma=None, 459 temperature=None, 460 filename=None, 461 is_mesh_symmetry=True): 462 suffix = _get_filename_suffix(mesh, 463 grid_point=grid_point, 464 sigma=sigma, 465 filename=filename) 466 jdos_filename = "jdos%s" % suffix 467 if temperature is not None: 468 jdos_filename += ("-t%f" % temperature).rstrip('0').rstrip(r'\.') 469 if not is_mesh_symmetry: 470 jdos_filename += ".nosym" 471 if filename is not None: 472 jdos_filename += ".%s" % filename 473 jdos_filename += ".dat" 474 475 with open(jdos_filename, 'w') as w: 476 for omega, vals in zip(frequencies, jdos): 477 w.write("%15.7f" % omega) 478 w.write((" %20.15e" * len(vals)) % tuple(vals)) 479 w.write("\n") 480 return jdos_filename 481 482 483def write_real_self_energy_at_grid_point(gp, 484 band_indices, 485 frequency_points, 486 deltas, 487 mesh, 488 epsilon, 489 temperature, 490 filename=None, 491 is_mesh_symmetry=True): 492 deltas_filename = "deltas" 493 deltas_filename += _get_filename_suffix(mesh, grid_point=gp) 494 if epsilon > 1e-5: 495 deltas_filename += "-e" + _del_zeros(epsilon) 496 else: 497 deltas_filename += "-e%.3e" % epsilon 498 if temperature is not None: 499 deltas_filename += "-t" + _del_zeros(temperature) + "-" 500 for i in band_indices: 501 deltas_filename += "b%d" % (i + 1) 502 if filename is not None: 503 deltas_filename += ".%s" % filename 504 elif not is_mesh_symmetry: 505 deltas_filename += ".nosym" 506 deltas_filename += ".dat" 507 508 with open(deltas_filename, 'w') as w: 509 for freq, v in zip(frequency_points, deltas): 510 w.write("%15.7f %20.15e\n" % (freq, v)) 511 512 return deltas_filename 513 514 515def write_real_self_energy_to_hdf5(grid_point, 516 band_indices, 517 temperatures, 518 deltas, 519 mesh, 520 epsilon, 521 frequency_points=None, 522 frequencies=None, 523 filename=None): 524 """Wirte real part of self energy (currently only bubble) in hdf5 525 526 deltas : ndarray 527 Real part of self energy. 528 529 With frequency_points: 530 shape=(temperatures, band_indices, frequency_points), 531 dtype='double', order='C' 532 otherwise: 533 shape=(temperatures, band_indices), dtype='double', order='C' 534 535 """ 536 full_filename = "deltas" 537 suffix = _get_filename_suffix(mesh, grid_point=grid_point) 538 _band_indices = np.array(band_indices, dtype='intc') 539 540 full_filename += suffix 541 if epsilon > 1e-5: 542 full_filename += "-e" + _del_zeros(epsilon) 543 else: 544 full_filename += "-e%.3e" % epsilon 545 full_filename += ".hdf5" 546 547 with h5py.File(full_filename, 'w') as w: 548 w.create_dataset('grid_point', data=grid_point) 549 w.create_dataset('mesh', data=mesh) 550 w.create_dataset('band_index', data=_band_indices) 551 w.create_dataset('delta', data=deltas) 552 w.create_dataset('temperature', data=temperatures) 553 w.create_dataset('epsilon', data=epsilon) 554 if frequency_points is not None: 555 w.create_dataset('frequency_points', data=frequency_points) 556 if frequencies is not None: 557 w.create_dataset('frequency', data=frequencies) 558 559 return full_filename 560 561 562def write_spectral_function_at_grid_point(gp, 563 band_indices, 564 frequency_points, 565 spectral_functions, 566 mesh, 567 temperature, 568 sigma=None, 569 filename=None, 570 is_mesh_symmetry=True): 571 spectral_filename = "spectral" 572 spectral_filename += _get_filename_suffix(mesh, grid_point=gp, sigma=sigma) 573 if temperature is not None: 574 spectral_filename += "-t" + _del_zeros(temperature) + "-" 575 for i in band_indices: 576 spectral_filename += "b%d" % (i + 1) 577 if filename is not None: 578 spectral_filename += ".%s" % filename 579 elif not is_mesh_symmetry: 580 spectral_filename += ".nosym" 581 spectral_filename += ".dat" 582 583 with open(spectral_filename, 'w') as w: 584 for freq, v in zip(frequency_points, spectral_functions): 585 w.write("%15.7f %20.15e\n" % (freq, v)) 586 587 return spectral_filename 588 589 590def write_spectral_function_to_hdf5(grid_point, 591 band_indices, 592 temperatures, 593 spectral_functions, 594 shifts, 595 half_linewidths, 596 mesh, 597 sigma=None, 598 frequency_points=None, 599 frequencies=None, 600 filename=None): 601 """Wirte spectral functions (currently only bubble) in hdf5 602 603 spectral_functions : ndarray 604 Spectral functions. 605 shape=(temperature, band_index, frequency_points), 606 dtype='double', order='C' 607 608 """ 609 full_filename = "spectral" 610 suffix = _get_filename_suffix(mesh, grid_point=grid_point, sigma=sigma) 611 _band_indices = np.array(band_indices, dtype='intc') 612 613 full_filename += suffix 614 full_filename += ".hdf5" 615 616 with h5py.File(full_filename, 'w') as w: 617 w.create_dataset('grid_point', data=grid_point) 618 w.create_dataset('mesh', data=mesh) 619 w.create_dataset('band_index', data=_band_indices) 620 w.create_dataset('spectral_function', data=spectral_functions) 621 w.create_dataset('shift', data=shifts) 622 w.create_dataset('half_linewidth', data=half_linewidths) 623 w.create_dataset('temperature', data=temperatures) 624 if frequency_points is not None: 625 w.create_dataset('frequency_point', data=frequency_points) 626 if frequencies is not None: 627 w.create_dataset('frequency', data=frequencies) 628 629 return full_filename 630 631 632def write_collision_to_hdf5(temperature, 633 mesh, 634 gamma=None, 635 gamma_isotope=None, 636 collision_matrix=None, 637 grid_point=None, 638 band_index=None, 639 sigma=None, 640 sigma_cutoff=None, 641 filename=None): 642 if band_index is None: 643 band_indices = None 644 else: 645 band_indices = [band_index] 646 suffix = _get_filename_suffix(mesh, 647 grid_point=grid_point, 648 band_indices=band_indices, 649 sigma=sigma, 650 sigma_cutoff=sigma_cutoff, 651 filename=filename) 652 full_filename = "collision" + suffix + ".hdf5" 653 with h5py.File(full_filename, 'w') as w: 654 w.create_dataset('temperature', data=temperature) 655 if gamma is not None: 656 w.create_dataset('gamma', data=gamma) 657 if gamma_isotope is not None: 658 w.create_dataset('gamma_isotope', data=gamma_isotope) 659 if collision_matrix is not None: 660 w.create_dataset('collision_matrix', data=collision_matrix) 661 if grid_point is not None: 662 w.create_dataset('grid_point', data=grid_point) 663 if band_index is not None: 664 w.create_dataset('band_index', data=(band_index + 1)) 665 if sigma is not None: 666 w.create_dataset('sigma', data=sigma) 667 if sigma_cutoff is not None: 668 w.create_dataset('sigma_cutoff_width', data=sigma_cutoff) 669 670 text = "Collisions " 671 if grid_point is not None: 672 text += "at grid adress %d " % grid_point 673 if sigma is not None: 674 if grid_point is not None: 675 text += "and " 676 else: 677 text += "at " 678 text += "sigma %s " % _del_zeros(sigma) 679 text += "were written into " 680 if sigma is not None: 681 text += "\n" 682 text += "\"%s\"." % ("collision" + suffix + ".hdf5") 683 print(text) 684 685 return full_filename 686 687 688def write_full_collision_matrix(collision_matrix, filename='fcm.hdf5'): 689 with h5py.File(filename, 'w') as w: 690 w.create_dataset('collision_matrix', data=collision_matrix) 691 692 693def write_unitary_matrix_to_hdf5(temperature, 694 mesh, 695 unitary_matrix=None, 696 sigma=None, 697 sigma_cutoff=None, 698 solver=None, 699 filename=None, 700 verbose=False): 701 """Write eigenvectors of collision matrices at temperatures. 702 703 Depending on the choice of the solver, eigenvectors are sotred in 704 either column-wise or row-wise. 705 706 """ 707 708 suffix = _get_filename_suffix(mesh, 709 sigma=sigma, 710 sigma_cutoff=sigma_cutoff, 711 filename=filename) 712 hdf5_filename = "unitary" + suffix + ".hdf5" 713 with h5py.File(hdf5_filename, 'w') as w: 714 w.create_dataset('temperature', data=temperature) 715 if unitary_matrix is not None: 716 w.create_dataset('unitary_matrix', data=unitary_matrix) 717 if solver is not None: 718 w.create_dataset('solver', data=solver) 719 720 if verbose: 721 if len(temperature) > 1: 722 text = "Unitary matrices " 723 else: 724 text = "Unitary matrix " 725 if sigma is not None: 726 text += "at sigma %s " % _del_zeros(sigma) 727 if sigma_cutoff is not None: 728 text += "(%4.2f SD) " % sigma_cutoff 729 if len(temperature) > 1: 730 text += "were written into " 731 else: 732 text += "was written into " 733 if sigma is not None: 734 text += "\n" 735 text += "\"%s\"." % hdf5_filename 736 print(text) 737 738 739def write_collision_eigenvalues_to_hdf5(temperatures, 740 mesh, 741 collision_eigenvalues, 742 sigma=None, 743 sigma_cutoff=None, 744 filename=None, 745 verbose=True): 746 suffix = _get_filename_suffix(mesh, 747 sigma=sigma, 748 sigma_cutoff=sigma_cutoff, 749 filename=filename) 750 with h5py.File("coleigs" + suffix + ".hdf5", 'w') as w: 751 w.create_dataset('temperature', data=temperatures) 752 w.create_dataset('collision_eigenvalues', data=collision_eigenvalues) 753 w.close() 754 755 if verbose: 756 text = "Eigenvalues of collision matrix " 757 if sigma is not None: 758 text += "with sigma %s\n" % sigma 759 text += "were written into " 760 text += "\"%s\"" % ("coleigs" + suffix + ".hdf5") 761 print(text) 762 763 764def write_kappa_to_hdf5(temperature, 765 mesh, 766 frequency=None, 767 group_velocity=None, 768 gv_by_gv=None, 769 mean_free_path=None, 770 heat_capacity=None, 771 kappa=None, 772 mode_kappa=None, 773 kappa_RTA=None, # RTA calculated in LBTE 774 mode_kappa_RTA=None, # RTA calculated in LBTE 775 f_vector=None, 776 gamma=None, 777 gamma_isotope=None, 778 gamma_N=None, 779 gamma_U=None, 780 averaged_pp_interaction=None, 781 qpoint=None, 782 weight=None, 783 mesh_divisors=None, 784 grid_point=None, 785 band_index=None, 786 sigma=None, 787 sigma_cutoff=None, 788 kappa_unit_conversion=None, 789 compression="gzip", 790 filename=None, 791 verbose=True): 792 if band_index is None: 793 band_indices = None 794 else: 795 band_indices = [band_index] 796 suffix = _get_filename_suffix(mesh, 797 mesh_divisors=mesh_divisors, 798 grid_point=grid_point, 799 band_indices=band_indices, 800 sigma=sigma, 801 sigma_cutoff=sigma_cutoff, 802 filename=filename) 803 full_filename = "kappa" + suffix + ".hdf5" 804 with h5py.File(full_filename, 'w') as w: 805 w.create_dataset('temperature', data=temperature) 806 w.create_dataset('mesh', data=mesh) 807 if frequency is not None: 808 w.create_dataset('frequency', data=frequency, 809 compression=compression) 810 if group_velocity is not None: 811 w.create_dataset('group_velocity', data=group_velocity, 812 compression=compression) 813 if gv_by_gv is not None: 814 w.create_dataset('gv_by_gv', data=gv_by_gv) 815 # if mean_free_path is not None: 816 # w.create_dataset('mean_free_path', data=mean_free_path, 817 # compression=compression) 818 if heat_capacity is not None: 819 w.create_dataset('heat_capacity', data=heat_capacity, 820 compression=compression) 821 if kappa is not None: 822 w.create_dataset('kappa', data=kappa) 823 if mode_kappa is not None: 824 w.create_dataset('mode_kappa', data=mode_kappa, 825 compression=compression) 826 if kappa_RTA is not None: 827 w.create_dataset('kappa_RTA', data=kappa_RTA) 828 if mode_kappa_RTA is not None: 829 w.create_dataset('mode_kappa_RTA', data=mode_kappa_RTA, 830 compression=compression) 831 if f_vector is not None: 832 w.create_dataset('f_vector', data=f_vector, 833 compression=compression) 834 if gamma is not None: 835 w.create_dataset('gamma', data=gamma, 836 compression=compression) 837 if gamma_isotope is not None: 838 w.create_dataset('gamma_isotope', data=gamma_isotope, 839 compression=compression) 840 if gamma_N is not None: 841 w.create_dataset('gamma_N', data=gamma_N, 842 compression=compression) 843 if gamma_U is not None: 844 w.create_dataset('gamma_U', data=gamma_U, 845 compression=compression) 846 if averaged_pp_interaction is not None: 847 w.create_dataset('ave_pp', data=averaged_pp_interaction, 848 compression=compression) 849 if qpoint is not None: 850 w.create_dataset('qpoint', data=qpoint, 851 compression=compression) 852 if weight is not None: 853 w.create_dataset('weight', data=weight, 854 compression=compression) 855 if grid_point is not None: 856 w.create_dataset('grid_point', data=grid_point) 857 if band_index is not None: 858 w.create_dataset('band_index', data=(band_index + 1)) 859 if sigma is not None: 860 w.create_dataset('sigma', data=sigma) 861 if sigma_cutoff is not None: 862 w.create_dataset('sigma_cutoff_width', data=sigma_cutoff) 863 if kappa_unit_conversion is not None: 864 w.create_dataset('kappa_unit_conversion', 865 data=kappa_unit_conversion) 866 867 if verbose: 868 text = "" 869 if kappa is not None: 870 text += "Thermal conductivity and related properties " 871 else: 872 text += "Thermal conductivity related properties " 873 if grid_point is not None: 874 text += "at gp-%d " % grid_point 875 if band_index is not None: 876 text += "and band_index-%d\n" % (band_index + 1) 877 if sigma is not None: 878 if grid_point is not None: 879 text += "and " 880 else: 881 text += "at " 882 text += "sigma %s" % sigma 883 if sigma_cutoff is None: 884 text += "\n" 885 else: 886 text += "(%4.2f SD)\n" % sigma_cutoff 887 text += "were written into " 888 else: 889 text += "were written into " 890 if band_index is None: 891 text += "\n" 892 text += "\"%s\"." % full_filename 893 print(text) 894 895 return full_filename 896 897 898def read_gamma_from_hdf5(mesh, 899 mesh_divisors=None, 900 grid_point=None, 901 band_index=None, 902 sigma=None, 903 sigma_cutoff=None, 904 filename=None, 905 verbose=True): 906 if band_index is None: 907 band_indices = None 908 else: 909 band_indices = [band_index] 910 suffix = _get_filename_suffix(mesh, 911 mesh_divisors=mesh_divisors, 912 grid_point=grid_point, 913 band_indices=band_indices, 914 sigma=sigma, 915 sigma_cutoff=sigma_cutoff, 916 filename=filename) 917 full_filename = "kappa" + suffix + ".hdf5" 918 if not os.path.exists(full_filename): 919 if verbose: 920 print("%s not found." % full_filename) 921 return None 922 923 read_data = {} 924 925 with h5py.File(full_filename, 'r') as f: 926 read_data['gamma'] = f['gamma'][:] 927 for key in ('gamma_isotope', 928 'ave_pp', 929 'gamma_N', 930 'gamma_U'): 931 if key in f.keys(): 932 if len(f[key].shape) > 0: 933 read_data[key] = f[key][:] 934 else: 935 read_data[key] = f[key][()] 936 if verbose: 937 print("Read data from %s." % full_filename) 938 939 return read_data 940 941 942def read_collision_from_hdf5(mesh, 943 indices=None, 944 grid_point=None, 945 band_index=None, 946 sigma=None, 947 sigma_cutoff=None, 948 filename=None, 949 verbose=True): 950 if band_index is None: 951 band_indices = None 952 else: 953 band_indices = [band_index] 954 suffix = _get_filename_suffix(mesh, 955 grid_point=grid_point, 956 band_indices=band_indices, 957 sigma=sigma, 958 sigma_cutoff=sigma_cutoff, 959 filename=filename) 960 full_filename = "collision" + suffix + ".hdf5" 961 if not os.path.exists(full_filename): 962 if verbose: 963 print("%s not found." % full_filename) 964 return None 965 966 with h5py.File(full_filename, 'r') as f: 967 if indices == 'all': 968 colmat_shape = (1,) + f['collision_matrix'].shape 969 collision_matrix = np.zeros(colmat_shape, 970 dtype='double', order='C') 971 gamma = np.array(f['gamma'][:], dtype='double', order='C') 972 collision_matrix[0] = f['collision_matrix'][:] 973 temperatures = np.array(f['temperature'][:], dtype='double') 974 else: 975 colmat_shape = (1, len(indices)) + f['collision_matrix'].shape[1:] 976 collision_matrix = np.zeros(colmat_shape, dtype='double') 977 gamma = np.array(f['gamma'][indices], dtype='double', order='C') 978 collision_matrix[0] = f['collision_matrix'][indices] 979 temperatures = np.array(f['temperature'][indices], dtype='double') 980 981 if verbose: 982 text = "Collisions " 983 if band_index is None: 984 if grid_point is not None: 985 text += "at grid point %d " % grid_point 986 else: 987 if grid_point is not None: 988 text += ("at (grid point %d, band index %d) " % 989 (grid_point, band_index)) 990 if sigma is not None: 991 if grid_point is not None: 992 text += "and " 993 else: 994 text += "at " 995 text += "sigma %s" % _del_zeros(sigma) 996 if sigma_cutoff is not None: 997 text += "(%4.2f SD)" % sigma_cutoff 998 if band_index is None and grid_point is not None: 999 text += " were read from " 1000 text += "\n" 1001 else: 1002 text += "\n" 1003 text += "were read from " 1004 text += "\"%s\"." % full_filename 1005 print(text) 1006 1007 return collision_matrix, gamma, temperatures 1008 1009 return None 1010 1011 1012def write_pp_to_hdf5(mesh, 1013 pp=None, 1014 g_zero=None, 1015 grid_point=None, 1016 triplet=None, 1017 weight=None, 1018 triplet_map=None, 1019 triplet_all=None, 1020 sigma=None, 1021 sigma_cutoff=None, 1022 filename=None, 1023 verbose=True, 1024 check_consistency=False, 1025 compression="gzip"): 1026 suffix = _get_filename_suffix(mesh, 1027 grid_point=grid_point, 1028 sigma=sigma, 1029 sigma_cutoff=sigma_cutoff, 1030 filename=filename) 1031 full_filename = "pp" + suffix + ".hdf5" 1032 1033 with h5py.File(full_filename, 'w') as w: 1034 if pp is not None: 1035 if g_zero is None: 1036 w.create_dataset('pp', data=pp, 1037 compression=compression) 1038 if triplet is not None: 1039 w.create_dataset('triplet', data=triplet, 1040 compression=compression) 1041 if weight is not None: 1042 w.create_dataset('weight', data=weight, 1043 compression=compression) 1044 if triplet_map is not None: 1045 w.create_dataset('triplet_map', data=triplet_map, 1046 compression=compression) 1047 if triplet_all is not None: 1048 w.create_dataset('triplet_all', data=triplet_all, 1049 compression=compression) 1050 else: 1051 x = g_zero.ravel() 1052 nonzero_pp = np.array(pp.ravel()[x == 0], dtype='double') 1053 bytelen = len(x) // 8 1054 remlen = len(x) % 8 1055 y = x[:bytelen * 8].reshape(-1, 8) 1056 z = np.packbits(y) 1057 if remlen != 0: 1058 z_rem = np.packbits(x[bytelen * 8:]) 1059 1060 # No compression for pp because already almost random. 1061 w.create_dataset('nonzero_pp', data=nonzero_pp, 1062 compression=None) 1063 w.create_dataset('pp_shape', data=pp.shape, 1064 compression=compression) 1065 w.create_dataset('g_zero_bits', data=z, 1066 compression=compression) 1067 if remlen != 0: 1068 w.create_dataset('g_zero_bits_reminder', data=z_rem) 1069 1070 # This is only for the test and coupled with read_pp_from_hdf5. 1071 if check_consistency: 1072 w.create_dataset('pp', data=pp, 1073 compression=compression) 1074 w.create_dataset('g_zero', data=g_zero, 1075 compression=compression) 1076 1077 if verbose: 1078 text = "" 1079 text += "Ph-ph interaction strength " 1080 if grid_point is not None: 1081 text += "at gp-%d " % grid_point 1082 if sigma is not None: 1083 if grid_point is not None: 1084 text += "and " 1085 else: 1086 text += "at " 1087 text += "sigma %s" % sigma 1088 if sigma_cutoff is None: 1089 text += "\n" 1090 else: 1091 text += "(%4.2f SD)\n" % sigma_cutoff 1092 text += "were written into " 1093 else: 1094 text += "were written into " 1095 text += "\n" 1096 text += "\"%s\"." % full_filename 1097 print(text) 1098 1099 return full_filename 1100 1101 1102def read_pp_from_hdf5(mesh, 1103 grid_point=None, 1104 sigma=None, 1105 sigma_cutoff=None, 1106 filename=None, 1107 verbose=True, 1108 check_consistency=False): 1109 suffix = _get_filename_suffix(mesh, 1110 grid_point=grid_point, 1111 sigma=sigma, 1112 sigma_cutoff=sigma_cutoff, 1113 filename=filename) 1114 full_filename = "pp" + suffix + ".hdf5" 1115 if not os.path.exists(full_filename): 1116 if verbose: 1117 print("%s not found." % full_filename) 1118 return None 1119 1120 with h5py.File(full_filename, 'r') as f: 1121 if 'nonzero_pp' in f: 1122 nonzero_pp = f['nonzero_pp'][:] 1123 pp_shape = f['pp_shape'][:] 1124 z = f['g_zero_bits'][:] 1125 bytelen = np.prod(pp_shape) // 8 1126 remlen = 0 1127 if 'g_zero_bits_reminder' in f: 1128 z_rem = f['g_zero_bits_reminder'][:] 1129 remlen = np.prod(pp_shape) - bytelen * 8 1130 1131 bits = np.unpackbits(z) 1132 if not bits.flags['C_CONTIGUOUS']: 1133 bits = np.array(bits, dtype='uint8') 1134 1135 g_zero = np.zeros(pp_shape, dtype='byte', order='C') 1136 b = g_zero.ravel() 1137 b[:(bytelen * 8)] = bits 1138 if remlen != 0: 1139 b[-remlen:] = np.unpackbits(z_rem)[:remlen] 1140 1141 pp = np.zeros(pp_shape, dtype='double', order='C') 1142 pp_ravel = pp.ravel() 1143 pp_ravel[g_zero.ravel() == 0] = nonzero_pp 1144 1145 # check_consistency==True in write_pp_to_hdf5 required. 1146 if check_consistency and g_zero is not None: 1147 if verbose: 1148 print("Checking consistency of ph-ph interanction " 1149 "strength.") 1150 assert (g_zero == f['g_zero'][:]).all() 1151 assert np.allclose(pp, f['pp'][:]) 1152 else: 1153 pp = np.zeros(f['pp'].shape, dtype='double', order='C') 1154 pp[:] = f['pp'][:] 1155 g_zero = None 1156 1157 if verbose: 1158 print("Ph-ph interaction strength was read from \"%s\"." % 1159 full_filename) 1160 1161 return pp, g_zero 1162 1163 return None 1164 1165 1166def write_gamma_detail_to_hdf5(temperature, 1167 mesh, 1168 gamma_detail=None, 1169 grid_point=None, 1170 triplet=None, 1171 weight=None, 1172 triplet_map=None, 1173 triplet_all=None, 1174 frequency_points=None, 1175 band_index=None, 1176 sigma=None, 1177 sigma_cutoff=None, 1178 compression="gzip", 1179 filename=None, 1180 verbose=True): 1181 if band_index is None: 1182 band_indices = None 1183 else: 1184 band_indices = [band_index] 1185 suffix = _get_filename_suffix(mesh, 1186 grid_point=grid_point, 1187 band_indices=band_indices, 1188 sigma=sigma, 1189 sigma_cutoff=sigma_cutoff, 1190 filename=filename) 1191 full_filename = "gamma_detail" + suffix + ".hdf5" 1192 1193 with h5py.File(full_filename, 'w') as w: 1194 w.create_dataset('temperature', data=temperature) 1195 w.create_dataset('mesh', data=mesh) 1196 if gamma_detail is not None: 1197 w.create_dataset('gamma_detail', data=gamma_detail, 1198 compression=compression) 1199 if triplet is not None: 1200 w.create_dataset('triplet', data=triplet, 1201 compression=compression) 1202 if weight is not None: 1203 w.create_dataset('weight', data=weight, 1204 compression=compression) 1205 if triplet_map is not None: 1206 w.create_dataset('triplet_map', data=triplet_map, 1207 compression=compression) 1208 if triplet_all is not None: 1209 w.create_dataset('triplet_all', data=triplet_all, 1210 compression=compression) 1211 if grid_point is not None: 1212 w.create_dataset('grid_point', data=grid_point) 1213 if band_index is not None: 1214 w.create_dataset('band_index', data=(band_index + 1)) 1215 if sigma is not None: 1216 w.create_dataset('sigma', data=sigma) 1217 if sigma_cutoff is not None: 1218 w.create_dataset('sigma_cutoff_width', data=sigma_cutoff) 1219 if frequency_points is not None: 1220 w.create_dataset('frequency_point', data=frequency_points) 1221 1222 if verbose: 1223 text = "" 1224 text += "Phonon triplets contributions to Gamma " 1225 if grid_point is not None: 1226 text += "at gp-%d " % grid_point 1227 if band_index is not None: 1228 text += "and band_index-%d\n" % (band_index + 1) 1229 if sigma is not None: 1230 if grid_point is not None: 1231 text += "and " 1232 else: 1233 text += "at " 1234 text += "sigma %s" % sigma 1235 if sigma_cutoff is None: 1236 text += "\n" 1237 else: 1238 text += "(%4.2f SD)\n" % sigma_cutoff 1239 text += "were written into " 1240 else: 1241 text += "were written into " 1242 if band_index is None: 1243 text += "\n" 1244 text += "\"%s\"." % full_filename 1245 print(text) 1246 1247 return full_filename 1248 1249 return None 1250 1251 1252def write_phonon_to_hdf5(frequency, 1253 eigenvector, 1254 grid_address, 1255 mesh, 1256 compression="gzip", 1257 filename=None): 1258 suffix = _get_filename_suffix(mesh, filename=filename) 1259 full_filename = "phonon" + suffix + ".hdf5" 1260 1261 with h5py.File(full_filename, 'w') as w: 1262 w.create_dataset('mesh', data=mesh) 1263 w.create_dataset('grid_address', data=grid_address, 1264 compression=compression) 1265 w.create_dataset('frequency', data=frequency, 1266 compression=compression) 1267 w.create_dataset('eigenvector', data=eigenvector, 1268 compression=compression) 1269 return full_filename 1270 1271 return None 1272 1273 1274def read_phonon_from_hdf5(mesh, 1275 filename=None, 1276 verbose=True): 1277 suffix = _get_filename_suffix(mesh, filename=filename) 1278 full_filename = "phonon" + suffix + ".hdf5" 1279 if not os.path.exists(full_filename): 1280 if verbose: 1281 print("%s not found." % full_filename) 1282 return None 1283 1284 with h5py.File(full_filename, 'r') as f: 1285 frequencies = np.array(f['frequency'][:], dtype='double', order='C') 1286 itemsize = frequencies.itemsize 1287 eigenvectors = np.array(f['eigenvector'][:], 1288 dtype=("c%d" % (itemsize * 2)), order='C') 1289 mesh_in_file = np.array(f['mesh'][:], dtype='intc') 1290 grid_address = np.array(f['grid_address'][:], dtype='intc', order='C') 1291 1292 assert (mesh_in_file == mesh).all(), "Mesh numbers are inconsistent." 1293 1294 if verbose: 1295 print("Phonons are read from \"%s\"." % full_filename) 1296 1297 return frequencies, eigenvectors, grid_address 1298 1299 return None 1300 1301 1302def write_ir_grid_points(mesh, 1303 mesh_divs, 1304 grid_points, 1305 coarse_grid_weights, 1306 grid_address, 1307 primitive_lattice): 1308 w = open("ir_grid_points.yaml", 'w') 1309 w.write("mesh: [ %d, %d, %d ]\n" % tuple(mesh)) 1310 if mesh_divs is not None: 1311 w.write("mesh_divisors: [ %d, %d, %d ]\n" % tuple(mesh_divs)) 1312 w.write("reciprocal_lattice:\n") 1313 for vec, axis in zip(primitive_lattice.T, ('a*', 'b*', 'c*')): 1314 w.write("- [ %12.8f, %12.8f, %12.8f ] # %2s\n" 1315 % (tuple(vec) + (axis,))) 1316 w.write("num_reduced_ir_grid_points: %d\n" % len(grid_points)) 1317 w.write("ir_grid_points: # [address, weight]\n") 1318 1319 for g, weight in zip(grid_points, coarse_grid_weights): 1320 w.write("- grid_point: %d\n" % g) 1321 w.write(" weight: %d\n" % weight) 1322 w.write(" grid_address: [ %12d, %12d, %12d ]\n" % 1323 tuple(grid_address[g])) 1324 w.write(" q-point: [ %12.7f, %12.7f, %12.7f ]\n" % 1325 tuple(grid_address[g].astype('double') / mesh)) 1326 1327 1328def parse_disp_fc2_yaml(filename="disp_fc2.yaml", return_cell=False): 1329 dataset = _parse_yaml(filename) 1330 natom = dataset['natom'] 1331 new_dataset = {} 1332 new_dataset['natom'] = natom 1333 new_first_atoms = [] 1334 for first_atoms in dataset['first_atoms']: 1335 first_atoms['number'] -= 1 1336 atom1 = first_atoms['number'] 1337 disp1 = first_atoms['displacement'] 1338 new_first_atoms.append({'number': atom1, 'displacement': disp1}) 1339 new_dataset['first_atoms'] = new_first_atoms 1340 1341 if return_cell: 1342 cell = get_cell_from_disp_yaml(dataset) 1343 return new_dataset, cell 1344 else: 1345 return new_dataset 1346 1347 1348def parse_disp_fc3_yaml(filename="disp_fc3.yaml", return_cell=False): 1349 dataset = _parse_yaml(filename) 1350 natom = dataset['natom'] 1351 new_dataset = {} 1352 new_dataset['natom'] = natom 1353 if 'cutoff_distance' in dataset: 1354 new_dataset['cutoff_distance'] = dataset['cutoff_distance'] 1355 new_first_atoms = [] 1356 for first_atoms in dataset['first_atoms']: 1357 atom1 = first_atoms['number'] - 1 1358 disp1 = first_atoms['displacement'] 1359 new_second_atoms = [] 1360 for second_atom in first_atoms['second_atoms']: 1361 disp2_dataset = {'number': second_atom['number'] - 1} 1362 if 'included' in second_atom: 1363 disp2_dataset.update({'included': second_atom['included']}) 1364 if 'distance' in second_atom: 1365 disp2_dataset.update( 1366 {'pair_distance': second_atom['distance']}) 1367 for disp2 in second_atom['displacements']: 1368 disp2_dataset.update({'displacement': disp2}) 1369 new_second_atoms.append(disp2_dataset.copy()) 1370 new_first_atoms.append({'number': atom1, 1371 'displacement': disp1, 1372 'second_atoms': new_second_atoms}) 1373 new_dataset['first_atoms'] = new_first_atoms 1374 1375 if return_cell: 1376 cell = get_cell_from_disp_yaml(dataset) 1377 return new_dataset, cell 1378 else: 1379 return new_dataset 1380 1381 1382def parse_FORCES_FC2(disp_dataset, 1383 filename="FORCES_FC2", 1384 unit_conversion_factor=None): 1385 num_atom = disp_dataset['natom'] 1386 num_disp = len(disp_dataset['first_atoms']) 1387 forces_fc2 = [] 1388 with open(filename, 'r') as f2: 1389 for i in range(num_disp): 1390 forces = _parse_force_lines(f2, num_atom) 1391 if forces is None: 1392 return [] 1393 else: 1394 forces_fc2.append(forces) 1395 1396 for i, disp1 in enumerate(disp_dataset['first_atoms']): 1397 if unit_conversion_factor is not None: 1398 disp1['forces'] = forces_fc2[i] * unit_conversion_factor 1399 else: 1400 disp1['forces'] = forces_fc2[i] 1401 1402 1403def parse_FORCES_FC3(disp_dataset, 1404 filename="FORCES_FC3", 1405 use_loadtxt=False, 1406 unit_conversion_factor=None): 1407 """Parse type1 FORCES_FC3 and store forces in disp_dataset""" 1408 1409 num_atom = disp_dataset['natom'] 1410 num_disp = len(disp_dataset['first_atoms']) 1411 for disp1 in disp_dataset['first_atoms']: 1412 num_disp += len(disp1['second_atoms']) 1413 1414 if use_loadtxt: 1415 forces_fc3 = np.loadtxt(filename).reshape((num_disp, -1, 3)) 1416 else: 1417 forces_fc3 = np.zeros((num_disp, num_atom, 3), 1418 dtype='double', order='C') 1419 with open(filename, 'r') as f3: 1420 for i in range(num_disp): 1421 forces = _parse_force_lines(f3, num_atom) 1422 if forces is None: 1423 raise RuntimeError("Failed to parse %s." % filename) 1424 else: 1425 forces_fc3[i] = forces 1426 1427 if unit_conversion_factor is not None: 1428 forces_fc3 *= unit_conversion_factor 1429 1430 i = 0 1431 for disp1 in disp_dataset['first_atoms']: 1432 disp1['forces'] = forces_fc3[i] 1433 i += 1 1434 1435 for disp1 in disp_dataset['first_atoms']: 1436 for disp2 in disp1['second_atoms']: 1437 disp2['forces'] = forces_fc3[i] 1438 i += 1 1439 1440 1441def parse_QPOINTS3(filename='QPOINTS3'): 1442 f = open(filename) 1443 num = int(f.readline().strip()) 1444 count = 0 1445 qpoints3 = [] 1446 for line in f: 1447 line_array = [float(x) for x in line.strip().split()] 1448 1449 if len(line_array) < 9: 1450 raise RuntimeError("Failed to parse %s." % filename) 1451 else: 1452 qpoints3.append(line_array[0:9]) 1453 1454 count += 1 1455 if count == num: 1456 break 1457 1458 return np.array(qpoints3) 1459 1460 1461def parse_fc3(num_atom, filename='fc3.dat'): 1462 f = open(filename) 1463 fc3 = np.zeros((num_atom, num_atom, num_atom, 3, 3, 3), dtype=float) 1464 for i in range(num_atom): 1465 for j in range(num_atom): 1466 for k in range(num_atom): 1467 f.readline() 1468 for l in range(3): 1469 fc3[i, j, k, l] = [ 1470 [float(x) for x in f.readline().split()], 1471 [float(x) for x in f.readline().split()], 1472 [float(x) for x in f.readline().split()]] 1473 f.readline() 1474 return fc3 1475 1476 1477def parse_fc2(num_atom, filename='fc2.dat'): 1478 f = open(filename) 1479 fc2 = np.zeros((num_atom, num_atom, 3, 3), dtype=float) 1480 for i in range(num_atom): 1481 for j in range(num_atom): 1482 f.readline() 1483 fc2[i, j] = [[float(x) for x in f.readline().split()], 1484 [float(x) for x in f.readline().split()], 1485 [float(x) for x in f.readline().split()]] 1486 f.readline() 1487 1488 return fc2 1489 1490 1491def parse_triplets(filename): 1492 f = open(filename) 1493 triplets = [] 1494 weights = [] 1495 for line in f: 1496 if line.strip()[0] == "#": 1497 continue 1498 1499 line_array = [int(x) for x in line.split()] 1500 triplets.append(line_array[:3]) 1501 weights.append(line_array[3]) 1502 1503 return np.array(triplets), np.array(weights) 1504 1505 1506def parse_grid_address(filename): 1507 f = open(filename, 'r') 1508 grid_address = [] 1509 for line in f: 1510 if line.strip()[0] == "#": 1511 continue 1512 1513 line_array = [int(x) for x in line.split()] 1514 grid_address.append(line_array[1:4]) 1515 1516 return np.array(grid_address) 1517 1518 1519def get_filename_suffix(mesh, 1520 mesh_divisors=None, 1521 grid_point=None, 1522 band_indices=None, 1523 sigma=None, 1524 sigma_cutoff=None, 1525 temperature=None, 1526 filename=None): 1527 return _get_filename_suffix(mesh, 1528 mesh_divisors=mesh_divisors, 1529 grid_point=grid_point, 1530 band_indices=band_indices, 1531 sigma=sigma, 1532 sigma_cutoff=sigma_cutoff, 1533 temperature=temperature, 1534 filename=filename) 1535 1536 1537def _get_filename_suffix(mesh, 1538 mesh_divisors=None, 1539 grid_point=None, 1540 band_indices=None, 1541 sigma=None, 1542 sigma_cutoff=None, 1543 temperature=None, 1544 filename=None): 1545 suffix = "-m%d%d%d" % tuple(mesh) 1546 if mesh_divisors is not None: 1547 if (np.array(mesh_divisors, dtype=int) != 1).any(): 1548 suffix += "-d%d%d%d" % tuple(mesh_divisors) 1549 if grid_point is not None: 1550 suffix += ("-g%d" % grid_point) 1551 if band_indices is not None: 1552 suffix += "-" 1553 for bi in band_indices: 1554 suffix += "b%d" % (bi + 1) 1555 if sigma is not None: 1556 suffix += "-s" + _del_zeros(sigma) 1557 if sigma_cutoff is not None: 1558 sigma_cutoff_str = _del_zeros(sigma_cutoff) 1559 suffix += "-sd" + sigma_cutoff_str 1560 if temperature is not None: 1561 suffix += "-t" + _del_zeros(temperature) 1562 if filename is not None: 1563 suffix += "." + filename 1564 1565 return suffix 1566 1567 1568def _del_zeros(val): 1569 return ("%f" % val).rstrip('0').rstrip(r'\.') 1570 1571 1572def _parse_yaml(file_yaml): 1573 import yaml 1574 try: 1575 from yaml import CLoader as Loader 1576 from yaml import CDumper as Dumper 1577 except ImportError: 1578 from yaml import Loader, Dumper 1579 1580 with open(file_yaml) as f: 1581 string = f.read() 1582 data = yaml.load(string, Loader=Loader) 1583 return data 1584 1585 1586def _parse_force_lines(forcefile, num_atom): 1587 forces = [] 1588 for line in forcefile: 1589 if line.strip() == '': 1590 continue 1591 if line.strip()[0] == '#': 1592 continue 1593 forces.append([float(x) for x in line.strip().split()]) 1594 if len(forces) == num_atom: 1595 break 1596 1597 if not len(forces) == num_atom: 1598 return None 1599 else: 1600 return np.array(forces) 1601 1602 1603def _parse_force_constants_lines(fcthird_file, num_atom): 1604 fc2 = [] 1605 for line in fcthird_file: 1606 if line.strip() == '': 1607 continue 1608 if line.strip()[0] == '#': 1609 continue 1610 fc2.append([float(x) for x in line.strip().split()]) 1611 if len(fc2) == num_atom ** 2 * 3: 1612 break 1613 1614 if not len(fc2) == num_atom ** 2 * 3: 1615 return None 1616 else: 1617 return np.array(fc2).reshape(num_atom, num_atom, 3, 3) 1618 1619 1620def get_length_of_first_line(f): 1621 for line in f: 1622 if line.strip() == '': 1623 continue 1624 elif line.strip()[0] == '#': 1625 continue 1626 else: 1627 f.seek(0) 1628 return len(line.split()) 1629 1630 raise RuntimeError("File doesn't contain relevant infomration.") 1631