1# Copyright (C) 2011 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phonopy. 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 36try: 37 from StringIO import StringIO 38except ImportError: 39 from io import StringIO 40import numpy as np 41 42import yaml 43try: 44 from yaml import CLoader as Loader 45except ImportError: 46 from yaml import Loader 47 48from phonopy.structure.atoms import PhonopyAtoms 49from phonopy.cui.settings import fracval 50from phonopy.structure.dataset import get_displacements_and_forces 51from phonopy.structure.symmetry import Symmetry, elaborate_borns_and_epsilon 52from phonopy.harmonic.force_constants import similarity_transformation 53 54 55# 56# FORCE_SETS 57# 58def write_FORCE_SETS(dataset, filename='FORCE_SETS'): 59 lines = get_FORCE_SETS_lines(dataset) 60 with open(filename, 'w') as w: 61 w.write("\n".join(lines)) 62 w.write("\n") 63 64 65def get_FORCE_SETS_lines(dataset, forces=None): 66 """Generate FORCE_SETS string 67 68 See the format of dataset in the docstring of 69 Phonopy.set_displacement_dataset. Optionally, sets of forces of supercells 70 can be given. In this case, these forces are unnecessary to be stored 71 in the dataset. 72 73 """ 74 75 if 'first_atoms' in dataset: 76 return _get_FORCE_SETS_lines_type1(dataset, forces=forces) 77 elif 'displacements' in dataset: 78 if forces is not None: 79 dataset['forces'] = forces 80 return _get_FORCE_SETS_lines_type2(dataset) 81 82 83def _get_FORCE_SETS_lines_type1(dataset, forces=None): 84 num_atom = dataset['natom'] 85 displacements = dataset['first_atoms'] 86 if forces is None: 87 _forces = [x['forces'] for x in dataset['first_atoms']] 88 else: 89 _forces = forces 90 91 lines = [] 92 lines.append("%-5d" % num_atom) 93 lines.append("%-5d" % len(displacements)) 94 for count, disp in enumerate(displacements): 95 lines.append("") 96 lines.append("%-5d" % (disp['number'] + 1)) 97 lines.append("%20.16f %20.16f %20.16f" % tuple(disp['displacement'])) 98 for f in _forces[count]: 99 lines.append("%15.10f %15.10f %15.10f" % tuple(f)) 100 101 return lines 102 103 104def _get_FORCE_SETS_lines_type2(dataset): 105 lines = [] 106 for displacements, forces in zip(dataset['displacements'], 107 dataset['forces']): 108 for d, f in zip(displacements, forces): 109 lines.append(("%15.8f" * 6) % (tuple(d) + tuple(f))) 110 111 return lines 112 113 114def parse_FORCE_SETS(natom=None, 115 is_translational_invariance=False, 116 filename="FORCE_SETS", 117 to_type2=False): 118 """ 119 120 to_type2 : bool 121 dataset of type2 is returned when True. 122 123 """ 124 125 with open(filename, 'r') as f: 126 return _get_dataset( 127 f, 128 natom=natom, 129 is_translational_invariance=is_translational_invariance, 130 to_type2=to_type2) 131 132 133def parse_FORCE_SETS_from_strings(strings, 134 natom=None, 135 is_translational_invariance=False, 136 to_type2=False): 137 return _get_dataset( 138 StringIO(strings), 139 natom=natom, 140 is_translational_invariance=is_translational_invariance, 141 to_type2=to_type2) 142 143 144def _get_dataset(f, 145 natom=None, 146 is_translational_invariance=False, 147 to_type2=False): 148 first_line_ary = _get_line_ignore_blank(f).split() 149 f.seek(0) 150 if len(first_line_ary) == 1: 151 if natom is None or int(first_line_ary[0]) == natom: 152 dataset = _get_dataset_type1(f, is_translational_invariance) 153 else: 154 msg = "Number of forces is not consistent with supercell setting." 155 raise RuntimeError(msg) 156 157 if to_type2: 158 disps, forces = get_displacements_and_forces(dataset) 159 return {'displacements': disps, 'forces': forces} 160 else: 161 return dataset 162 163 elif len(first_line_ary) == 6: 164 return _get_dataset_type2(f, natom) 165 166 167def _get_dataset_type1(f, is_translational_invariance): 168 set_of_forces = [] 169 num_atom = int(_get_line_ignore_blank(f)) 170 num_displacements = int(_get_line_ignore_blank(f)) 171 172 for i in range(num_displacements): 173 line = _get_line_ignore_blank(f) 174 atom_number = int(line) 175 line = _get_line_ignore_blank(f).split() 176 displacement = np.array([float(x) for x in line]) 177 forces_tmp = [] 178 for j in range(num_atom): 179 line = _get_line_ignore_blank(f).split() 180 forces_tmp.append(np.array([float(x) for x in line])) 181 forces_tmp = np.array(forces_tmp, dtype='double') 182 183 if is_translational_invariance: 184 forces_tmp -= np.sum(forces_tmp, axis=0) / len(forces_tmp) 185 186 forces = {'number': atom_number - 1, 187 'displacement': displacement, 188 'forces': forces_tmp} 189 set_of_forces.append(forces) 190 191 dataset = {'natom': num_atom, 192 'first_atoms': set_of_forces} 193 194 return dataset 195 196 197def get_dataset_type2(f, natom): 198 return _get_dataset_type2(f, natom) 199 200 201def _get_dataset_type2(f, natom): 202 data = np.loadtxt(f, dtype='double') 203 if data.shape[1] != 6 or (natom and data.shape[0] % natom != 0): 204 msg = "Data shape of forces and displacements is incorrect." 205 raise RuntimeError(msg) 206 if natom: 207 data = data.reshape(-1, natom, 6) 208 displacements = data[:, :, :3] 209 forces = data[:, :, 3:] 210 else: 211 displacements = data[:, :3] 212 forces = data[:, 3:] 213 dataset = {'displacements': 214 np.array(displacements, dtype='double', order='C'), 215 'forces': np.array(forces, dtype='double', order='C')} 216 return dataset 217 218 219def _get_line_ignore_blank(f): 220 line = f.readline().strip() 221 if line == '': 222 line = _get_line_ignore_blank(f) 223 return line 224 225 226def collect_forces(f, num_atom, hook, force_pos, word=None): 227 for line in f: 228 if hook in line: 229 break 230 231 forces = [] 232 for line in f: 233 if line.strip() == '': 234 continue 235 if word is not None: 236 if word not in line: 237 continue 238 239 elems = line.split() 240 if len(elems) > force_pos[2]: 241 try: 242 forces.append([float(elems[i]) for i in force_pos]) 243 except ValueError: 244 forces = [] 245 break 246 else: 247 return False 248 249 if len(forces) == num_atom: 250 break 251 252 return forces 253 254 255def iter_collect_forces(filename, 256 num_atom, 257 hook, 258 force_pos, 259 word=None, 260 max_iter=1000): 261 with open(filename) as f: 262 forces = [] 263 prev_forces = [] 264 265 for i in range(max_iter): 266 forces = collect_forces(f, num_atom, hook, force_pos, word=word) 267 if not forces: 268 forces = prev_forces[:] 269 break 270 else: 271 prev_forces = forces[:] 272 273 if i == max_iter - 1: 274 sys.stderr.write("Reached to max number of iterations (%d).\n" % 275 max_iter) 276 277 return forces 278 279 280# 281# FORCE_CONSTANTS, force_constants.hdf5 282# 283def write_FORCE_CONSTANTS(force_constants, 284 filename='FORCE_CONSTANTS', 285 p2s_map=None): 286 """Write force constants in text file format. 287 288 Parameters 289 ---------- 290 force_constants: ndarray 291 Force constants 292 shape=(n_satom,n_satom,3,3) or (n_patom,n_satom,3,3) 293 dtype=double 294 filename: str 295 Filename to be saved 296 p2s_map: ndarray 297 Primitive atom indices in supercell index system 298 dtype=intc 299 300 """ 301 302 lines = get_FORCE_CONSTANTS_lines(force_constants, p2s_map=p2s_map) 303 with open(filename, 'w') as w: 304 w.write("\n".join(lines)) 305 306 307def get_FORCE_CONSTANTS_lines(force_constants, p2s_map=None): 308 if p2s_map is not None and len(p2s_map) == force_constants.shape[0]: 309 indices = p2s_map 310 else: 311 indices = np.arange(force_constants.shape[0], dtype='intc') 312 313 lines = [] 314 fc_shape = force_constants.shape 315 lines.append("%4d %4d" % fc_shape[:2]) 316 for i, s_i in enumerate(indices): 317 for j in range(fc_shape[1]): 318 lines.append("%d %d" % (s_i + 1, j + 1)) 319 for vec in force_constants[i][j]: 320 lines.append(("%22.15f" * 3) % tuple(vec)) 321 322 return lines 323 324 325def write_force_constants_to_hdf5(force_constants, 326 filename='force_constants.hdf5', 327 p2s_map=None, 328 physical_unit=None, 329 compression=None): 330 """Write force constants in hdf5 format. 331 332 Parameters 333 ---------- 334 force_constants: ndarray 335 Force constants 336 shape=(n_satom,n_satom,3,3) or (n_patom,n_satom,3,3) 337 dtype=double 338 filename: str 339 Filename to be saved 340 p2s_map: ndarray 341 Primitive atom indices in supercell index system 342 shape=(n_patom,) 343 dtype=intc 344 physical_unit : str, optional 345 Physical unit used for force contants. Default is None. 346 compression : str or int, optional 347 h5py's lossless compression filters (e.g., "gzip", "lzf"). 348 See the detail at docstring of h5py.Group.create_dataset. Default is 349 None. 350 351 """ 352 353 try: 354 import h5py 355 except ImportError: 356 raise ModuleNotFoundError("You need to install python-h5py.") 357 358 with h5py.File(filename, 'w') as w: 359 w.create_dataset('force_constants', data=force_constants, 360 compression=compression) 361 if p2s_map is not None: 362 w.create_dataset('p2s_map', data=p2s_map) 363 if physical_unit is not None: 364 dset = w.create_dataset('physical_unit', (1,), 365 dtype='S%d' % len(physical_unit)) 366 dset[0] = np.string_(physical_unit) 367 368 369def parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS", 370 p2s_map=None): 371 with open(filename) as fcfile: 372 idx1 = [] 373 374 line = fcfile.readline() 375 idx = [int(x) for x in line.split()] 376 if len(idx) == 1: 377 idx = [idx[0], idx[0]] 378 force_constants = np.zeros((idx[0], idx[1], 3, 3), dtype='double') 379 for i in range(idx[0]): 380 for j in range(idx[1]): 381 s_i = int(fcfile.readline().split()[0]) - 1 382 if s_i not in idx1: 383 idx1.append(s_i) 384 tensor = [] 385 for k in range(3): 386 tensor.append([float(x) 387 for x in fcfile.readline().split()]) 388 force_constants[i, j] = tensor 389 390 check_force_constants_indices(idx, idx1, p2s_map, filename) 391 392 return force_constants 393 394 395def read_force_constants_hdf5(filename="force_constants.hdf5", 396 p2s_map=None, 397 return_physical_unit=False): 398 try: 399 import h5py 400 except ImportError: 401 raise ModuleNotFoundError("You need to install python-h5py.") 402 403 with h5py.File(filename, 'r') as f: 404 if 'fc2' in f: 405 key = 'fc2' 406 elif 'force_constants' in f: 407 key = 'force_constants' 408 else: 409 raise RuntimeError("%s doesn't contain necessary information" % 410 filename) 411 412 fc = f[key][:] 413 if 'p2s_map' in f: 414 p2s_map_in_file = f['p2s_map'][:] 415 check_force_constants_indices(fc.shape[:2], 416 p2s_map_in_file, 417 p2s_map, 418 filename) 419 420 if return_physical_unit: 421 if 'physical_unit' in f: 422 physical_unit = f['physical_unit'][0].decode('utf-8') 423 else: 424 physical_unit = None 425 return fc, physical_unit 426 else: 427 return fc 428 429 430def check_force_constants_indices(shape, indices, p2s_map, filename): 431 if shape[0] != shape[1] and p2s_map is not None: 432 if len(p2s_map) != len(indices) or (p2s_map != indices).any(): 433 text = ("%s file is inconsistent with the calculation setting. " 434 "PRIMITIVE_AXIS may not be set correctly.") % filename 435 raise RuntimeError(text) 436 437 438def parse_disp_yaml(filename="disp.yaml", return_cell=False): 439 """Read disp.yaml or phonopy_disp.yaml 440 441 This method was originally made for parsing disp.yaml. Later this 442 started to work for phonopy_disp.yaml, too. But now this method is not 443 allowed to read phonopy_disp.yaml because of existance of PhonopyYaml 444 class. 445 446 """ 447 448 with open(filename) as f: 449 new_dataset = {} 450 dataset = yaml.load(f, Loader=Loader) 451 if 'phonopy' in dataset and 'calculator' in dataset['phonopy']: 452 new_dataset['calculator'] = dataset['phonopy']['calculator'] 453 if 'natom' in dataset: 454 natom = dataset['natom'] 455 elif 'supercell' and 'points' in dataset['supercell']: 456 natom = len(dataset['supercell']['points']) 457 else: 458 raise RuntimeError("%s doesn't contain necessary information.") 459 new_dataset['natom'] = natom 460 new_first_atoms = [] 461 462 try: 463 displacements = dataset['displacements'] 464 except KeyError: 465 raise 466 467 if type(displacements[0]) is dict: 468 for first_atoms in displacements: 469 first_atoms['atom'] -= 1 470 atom1 = first_atoms['atom'] 471 disp1 = first_atoms['displacement'] 472 new_first_atoms.append({'number': atom1, 473 'displacement': disp1}) 474 new_dataset['first_atoms'] = new_first_atoms 475 476 if return_cell: 477 cell = get_cell_from_disp_yaml(dataset) 478 return new_dataset, cell 479 else: 480 return new_dataset 481 482 483def write_disp_yaml_from_dataset(dataset, supercell, filename='disp.yaml'): 484 displacements = [(d['number'],) + tuple(d['displacement']) 485 for d in dataset['first_atoms']] 486 write_disp_yaml(displacements, supercell, filename=filename) 487 488 489def write_disp_yaml(displacements, supercell, filename='disp.yaml'): 490 lines = [] 491 lines.append("natom: %4d" % supercell.get_number_of_atoms()) 492 lines += get_disp_yaml_lines(displacements, supercell) 493 lines.append(str(supercell)) 494 495 with open(filename, 'w') as w: 496 w.write("\n".join(lines)) 497 498 499def get_disp_yaml_lines(displacements, supercell): 500 lines = [] 501 lines.append("displacements:") 502 for i, disp in enumerate(displacements): 503 lines.append("- atom: %4d" % (disp[0] + 1)) 504 lines.append(" displacement:") 505 lines.append(" [ %20.16f,%20.16f,%20.16f ]" % tuple(disp[1:4])) 506 return lines 507 508 509# 510# DISP (old phonopy displacement format) 511# 512def parse_DISP(filename='DISP'): 513 with open(filename) as disp: 514 displacements = [] 515 for line in disp: 516 if line.strip() != '': 517 a = line.split() 518 displacements.append( 519 [int(a[0])-1, float(a[1]), float(a[2]), float(a[3])]) 520 return displacements 521 522 523# 524# Parse supercell in disp.yaml 525# 526def get_cell_from_disp_yaml(dataset): 527 if 'lattice' in dataset: 528 lattice = dataset['lattice'] 529 if 'points' in dataset: 530 data_key = 'points' 531 pos_key = 'coordinates' 532 elif 'atoms' in dataset: 533 data_key = 'atoms' 534 pos_key = 'position' 535 else: 536 data_key = None 537 pos_key = None 538 539 try: 540 positions = [x[pos_key] for x in dataset[data_key]] 541 except KeyError: 542 msg = ("\"disp.yaml\" format is too old. " 543 "Please re-create it as \"phonopy_disp.yaml\" to contain " 544 "supercell crystal structure information.") 545 raise RuntimeError(msg) 546 symbols = [x['symbol'] for x in dataset[data_key]] 547 cell = PhonopyAtoms(cell=lattice, 548 scaled_positions=positions, 549 symbols=symbols, 550 pbc=True) 551 return cell 552 else: 553 return get_cell_from_disp_yaml(dataset['supercell']) 554 555 556# 557# QPOINTS 558# 559def parse_QPOINTS(filename="QPOINTS"): 560 with open(filename, 'r') as f: 561 num_qpoints = int(f.readline().strip()) 562 qpoints = [] 563 for i in range(num_qpoints): 564 qpoints.append([fracval(x) for x in f.readline().strip().split()]) 565 return np.array(qpoints) 566 567 568# 569# BORN 570# 571def write_BORN(primitive, borns, epsilon, filename="BORN"): 572 lines = get_BORN_lines(primitive, borns, epsilon) 573 with open(filename, 'w') as w: 574 w.write('\n'.join(lines)) 575 576 577def get_BORN_lines(unitcell, borns, epsilon, 578 factor=None, 579 primitive_matrix=None, 580 supercell_matrix=None, 581 symprec=1e-5): 582 borns, epsilon, atom_indices = elaborate_borns_and_epsilon( 583 unitcell, borns, epsilon, symmetrize_tensors=True, 584 primitive_matrix=primitive_matrix, 585 supercell_matrix=supercell_matrix, 586 symprec=symprec) 587 588 text = "# epsilon and Z* of atoms " 589 text += ' '.join(["%d" % n for n in atom_indices + 1]) 590 lines = [text, ] 591 lines.append(("%13.8f " * 9) % tuple(epsilon.flatten())) 592 for z in borns: 593 lines.append(("%13.8f " * 9) % tuple(z.flatten())) 594 return lines 595 596 597def parse_BORN(primitive, symprec=1e-5, is_symmetry=True, filename="BORN"): 598 with open(filename, 'r') as f: 599 return _parse_BORN_from_file_object(f, primitive, symprec, is_symmetry) 600 601 602def parse_BORN_from_strings(strings, primitive, 603 symprec=1e-5, is_symmetry=True): 604 f = StringIO(strings) 605 return _parse_BORN_from_file_object(f, primitive, symprec, is_symmetry) 606 607 608def _parse_BORN_from_file_object(f, primitive, symprec, is_symmetry): 609 symmetry = Symmetry(primitive, symprec=symprec, is_symmetry=is_symmetry) 610 return get_born_parameters(f, primitive, symmetry) 611 612 613def get_born_parameters(f, primitive, prim_symmetry): 614 line_arr = f.readline().split() 615 if len(line_arr) < 1: 616 print("BORN file format of line 1 is incorrect") 617 return None 618 619 factor = None 620 G_cutoff = None 621 Lambda = None 622 623 if len(line_arr) > 0: 624 try: 625 factor = float(line_arr[0]) 626 except (ValueError, TypeError): 627 factor = None 628 if len(line_arr) > 1: 629 try: 630 G_cutoff = float(line_arr[1]) 631 except (ValueError, TypeError): 632 G_cutoff = None 633 if len(line_arr) > 2: 634 try: 635 Lambda = float(line_arr[2]) 636 except (ValueError, TypeError): 637 Lambda = None 638 639 # Read dielectric constant 640 line = f.readline().split() 641 if not len(line) == 9: 642 print("BORN file format of line 2 is incorrect") 643 return None 644 dielectric = np.reshape([float(x) for x in line], (3, 3)) 645 646 # Read Born effective charge 647 independent_atoms = prim_symmetry.get_independent_atoms() 648 borns = np.zeros((primitive.get_number_of_atoms(), 3, 3), 649 dtype='double', order='C') 650 651 for i in independent_atoms: 652 line = f.readline().split() 653 if len(line) == 0: 654 print("Number of lines for Born effect charge is not enough.") 655 return None 656 if not len(line) == 9: 657 print("BORN file format of line %d is incorrect" % (i + 3)) 658 return None 659 borns[i] = np.reshape([float(x) for x in line], (3, 3)) 660 661 # Check that the number of atoms in the BORN file was correct 662 line = f.readline().split() 663 if len(line) > 0: 664 print("Too many atoms in the BORN file (it should only contain " 665 "symmetry-independent atoms)") 666 return None 667 668 _expand_borns(borns, primitive, prim_symmetry) 669 non_anal = {'born': borns, 670 'factor': factor, 671 'dielectric': dielectric} 672 if G_cutoff is not None: 673 non_anal['G_cutoff'] = G_cutoff 674 if Lambda is not None: 675 non_anal['Lambda'] = Lambda 676 677 return non_anal 678 679 680def _expand_borns(borns, primitive, prim_symmetry): 681 # Expand Born effective charges to all atoms in the primitive cell 682 rotations = prim_symmetry.get_symmetry_operations()['rotations'] 683 map_operations = prim_symmetry.get_map_operations() 684 map_atoms = prim_symmetry.get_map_atoms() 685 686 for i in range(primitive.get_number_of_atoms()): 687 # R_cart = L R L^-1 688 rot_cartesian = similarity_transformation( 689 primitive.get_cell().transpose(), rotations[map_operations[i]]) 690 # R_cart^T B R_cart^-T (inverse rotation is required to transform) 691 borns[i] = similarity_transformation(rot_cartesian.transpose(), 692 borns[map_atoms[i]]) 693 694 695# 696# phonopy.yaml 697# 698def is_file_phonopy_yaml(filename, keyword='phonopy'): 699 with open(filename, 'r') as f: 700 try: 701 data = yaml.load(f, Loader=Loader) 702 if data is None: 703 return False 704 if keyword in data: 705 return True 706 else: 707 return False 708 except yaml.YAMLError: 709 return False 710 711 712# 713# e-v.dat, thermal_properties.yaml 714# 715def read_thermal_properties_yaml(filenames): 716 thermal_properties = [] 717 num_modes = [] 718 num_integrated_modes = [] 719 for filename in filenames: 720 with open(filename) as f: 721 tp_yaml = yaml.load(f, Loader=Loader) 722 thermal_properties.append(tp_yaml['thermal_properties']) 723 if 'num_modes' in tp_yaml and 'num_integrated_modes' in tp_yaml: 724 num_modes.append(tp_yaml['num_modes']) 725 num_integrated_modes.append(tp_yaml['num_integrated_modes']) 726 727 temperatures = [v['temperature'] for v in thermal_properties[0]] 728 temp = [] 729 cv = [] 730 entropy = [] 731 fe_phonon = [] 732 for i, tp in enumerate(thermal_properties): 733 temp.append([v['temperature'] for v in tp]) 734 if not np.allclose(temperatures, temp): 735 msg = ['', ] 736 msg.append("Check your input files") 737 msg.append("Disagreement of temperature range or step") 738 for t, fname in zip(temp, filenames): 739 msg.append("%s: Range [ %d, %d ], Step %f" % 740 (fname, int(t[0]), int(t[-1]), t[1] - t[0])) 741 msg.append('') 742 msg.append("Stop phonopy-qha") 743 raise RuntimeError(msg) 744 cv.append([v['heat_capacity'] for v in tp]) 745 entropy.append([v['entropy'] for v in tp]) 746 fe_phonon.append([v['free_energy'] for v in tp]) 747 748 # shape=(temperatures, volumes) 749 cv = np.array(cv).T 750 entropy = np.array(entropy).T 751 fe_phonon = np.array(fe_phonon).T 752 753 return (temperatures, cv, entropy, fe_phonon, num_modes, 754 num_integrated_modes) 755 756 757def read_v_e(filename): 758 data = _parse_QHA_data(filename) 759 if data.shape[1] != 2: 760 msg = ("File format of %s is incorrect for reading e-v data." % 761 filename) 762 raise RuntimeError(msg) 763 volumes, electronic_energies = data.T 764 return volumes, electronic_energies 765 766 767def read_efe(filename): 768 data = _parse_QHA_data(filename) 769 temperatures = data[:, 0] 770 free_energies = data[:, 1:] 771 return temperatures, free_energies 772 773 774def _parse_QHA_data(filename): 775 data = [] 776 with open(filename) as f: 777 for line in f: 778 if line.strip() == '' or line.strip()[0] == '#': 779 continue 780 if '#' in line: 781 data.append([float(x) for x in line.split('#')[0].split()]) 782 else: 783 data.append([float(x) for x in line.split()]) 784 return np.array(data) 785