1# Copyright (C) 2018 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 os 36import numpy as np 37import yaml 38try: 39 from yaml import CLoader as Loader 40except ImportError: 41 from yaml import Loader 42 43from phonopy.structure.atoms import PhonopyAtoms 44 45 46def read_cell_yaml(filename, cell_type='unitcell'): 47 ph_yaml = PhonopyYaml() 48 ph_yaml.read(filename) 49 if ph_yaml.unitcell and cell_type == 'unitcell': 50 return ph_yaml.unitcell 51 elif ph_yaml.primitive and cell_type == 'primitive': 52 return ph_yaml.primitive 53 elif ph_yaml.supercell and cell_type == 'supercell': 54 return ph_yaml.supercell 55 else: 56 return None 57 58 59class PhonopyYaml(object): 60 """PhonopyYaml is a container of phonopy setting 61 62 This contains the writer (__str__) and reader (read) of phonopy.yaml type 63 file. 64 65 Methods 66 ------- 67 __str__ 68 Return string of phonopy.yaml. 69 get_yaml_lines 70 Return a list of string lines of phonopy.yaml. 71 read 72 Read specific properties written in phonopy.yaml. 73 set_phonon_info 74 Copy specific properties in Phonopy instance to self. 75 76 Attributes 77 ---------- 78 configuration : dict 79 Phonopy setting tags or options (e.g., {"DIM": "2 2 2", ...}) 80 calculator : str 81 Force calculator. 82 physical_units : dict 83 Set of physical units used in this phonon calculation. 84 unitcell : PhonopyAtoms 85 Unit cell. 86 primitive : PhonopyAtoms 87 Primitive cell. The instance of Primitive class is necessary has to 88 be created from the instance of Supercell class with 89 np.dot(np.linalg.inv(supercell_matrix), primitive_matrix). 90 supercell : PhonopyAtoms 91 Supercell. The instance of Supercell class is necessary has to be 92 created from unitcell with supercel_matrix. 93 dataset 94 supercell_matrix 95 primitive_matrix 96 nac_params 97 force_constants 98 symmetry 99 s2p_map 100 u2p_map 101 frequency_unit_conversion_factor 102 version 103 yaml_filename 104 settings 105 command_name 106 default_filenames 107 108 """ 109 110 command_name = "phonopy" 111 default_filenames = ("phonopy_params.yaml", 112 "phonopy_disp.yaml", 113 "phonopy.yaml") 114 default_settings = {'force_sets': True, 115 'displacements': True, 116 'force_constants': False, 117 'born_effective_charge': True, 118 'dielectric_constant': True} 119 120 def __init__(self, 121 configuration=None, 122 calculator=None, 123 physical_units=None, 124 settings=None): 125 self.configuration = configuration 126 self.calculator = calculator 127 self.physical_units = physical_units 128 129 self.unitcell = None 130 self.primitive = None 131 self.supercell = None 132 self.dataset = None 133 self.supercell_matrix = None 134 self.primitive_matrix = None 135 self.nac_params = None 136 self.force_constants = None 137 138 self.symmetry = None # symmetry of supercell 139 self.s2p_map = None 140 self.u2p_map = None 141 self.frequency_unit_conversion_factor = None 142 self.version = None 143 self.yaml_filename = None 144 145 self.settings = self.default_settings.copy() 146 if type(settings) is dict: 147 self.settings.update(settings) 148 149 self._yaml = None 150 151 def __str__(self): 152 return "\n".join(self.get_yaml_lines()) 153 154 def read(self, filename): 155 self.yaml_filename = filename 156 self._load(filename) 157 158 @property 159 def yaml_data(self): 160 return self._yaml 161 162 @yaml_data.setter 163 def yaml_data(self, yaml_data): 164 self._yaml = yaml_data 165 166 def parse(self): 167 self._parse_transformation_matrices() 168 self._parse_all_cells() 169 self._parse_force_constants() 170 self._parse_dataset() 171 self._parse_nac_params() 172 self._parse_calculator() 173 174 def set_cell(self, cell): 175 self.unitcell = cell 176 177 def set_phonon_info(self, phonopy): 178 self.unitcell = phonopy.unitcell 179 self.primitive = phonopy.primitive 180 self.supercell = phonopy.supercell 181 self.version = phonopy.version 182 self.supercell_matrix = phonopy.supercell_matrix 183 self.symmetry = phonopy.symmetry 184 self.primitive_matrix = phonopy.primitive_matrix 185 s2p_map = self.primitive.s2p_map 186 u2s_map = self.supercell.u2s_map 187 u2u_map = self.supercell.u2u_map 188 s2u_map = self.supercell.s2u_map 189 self.u2p_map = [u2u_map[i] for i in (s2u_map[s2p_map])[u2s_map]] 190 self.nac_params = phonopy.nac_params 191 self.frequency_unit_conversion_factor = phonopy.unit_conversion_factor 192 self.calculator = phonopy.calculator 193 self.force_constants = phonopy.force_constants 194 self.dataset = phonopy.dataset 195 196 def get_yaml_lines(self): 197 lines = self._header_yaml_lines() 198 lines += self._physical_units_yaml_lines() 199 lines += self._symmetry_yaml_lines() 200 lines += self._cell_info_yaml_lines() 201 lines += self._nac_yaml_lines() 202 lines += self._dataset_yaml_lines() 203 lines += self._force_constants_yaml_lines() 204 return lines 205 206 def _header_yaml_lines(self): 207 lines = [] 208 lines.append("%s:" % self.command_name) 209 if self.version is None: 210 from phonopy.version import __version__ 211 lines.append(" version: %s" % __version__) 212 else: 213 lines.append(" version: %s" % self.version) 214 if self.calculator: 215 lines.append(" calculator: %s" % self.calculator) 216 if self.frequency_unit_conversion_factor: 217 lines.append(" frequency_unit_conversion_factor: %f" % 218 self.frequency_unit_conversion_factor) 219 if self.symmetry: 220 lines.append(" symmetry_tolerance: %.5e" % 221 self.symmetry.get_symmetry_tolerance()) 222 if self.nac_params: 223 lines.append(" nac_unit_conversion_factor: %f" 224 % self.nac_params['factor']) 225 if self.configuration is not None: 226 lines.append(" configuration:") 227 for key in self.configuration: 228 val = self.configuration[key] 229 if type(val) is str: 230 val = val.replace('\\', '\\\\') 231 lines.append(" %s: \"%s\"" % (key, val)) 232 lines.append("") 233 return lines 234 235 def _physical_units_yaml_lines(self): 236 lines = [] 237 lines.append("physical_unit:") 238 lines.append(" atomic_mass: \"AMU\"") 239 units = self.physical_units 240 if units is not None: 241 if units['length_unit'] is not None: 242 lines.append(" length: \"%s\"" % units['length_unit']) 243 if (self.command_name == "phonopy" and 244 units['force_constants_unit'] is not None): 245 lines.append(" force_constants: \"%s\"" % 246 units['force_constants_unit']) 247 lines.append("") 248 return lines 249 250 def _symmetry_yaml_lines(self): 251 lines = [] 252 if self.symmetry is not None and self.symmetry.dataset is not None: 253 lines.append("space_group:") 254 lines.append(" type: \"%s\"" % 255 self.symmetry.get_dataset()['international']) 256 lines.append(" number: %d" % 257 self.symmetry.get_dataset()['number']) 258 hall_symbol = self.symmetry.get_dataset()['hall'] 259 if "\"" in hall_symbol: 260 hall_symbol = hall_symbol.replace("\"", "\\\"") 261 lines.append(" Hall_symbol: \"%s\"" % hall_symbol) 262 lines.append("") 263 return lines 264 265 def _cell_info_yaml_lines(self): 266 lines = self._primitive_matrix_yaml_lines( 267 self.primitive_matrix, "primitive_matrix") 268 lines += self._supercell_matrix_yaml_lines( 269 self.supercell_matrix, "supercell_matrix") 270 lines += self._primitive_yaml_lines(self.primitive, "primitive_cell") 271 lines += self._unitcell_yaml_lines() 272 lines += self._supercell_yaml_lines() 273 return lines 274 275 def _primitive_matrix_yaml_lines(self, matrix, name): 276 lines = [] 277 if matrix is not None: 278 lines.append("%s:" % name) 279 for v in matrix: 280 lines.append("- [ %18.15f, %18.15f, %18.15f ]" % tuple(v)) 281 lines.append("") 282 return lines 283 284 def _supercell_matrix_yaml_lines(self, matrix, name): 285 lines = [] 286 if matrix is not None: 287 lines.append("%s:" % name) 288 for v in matrix: 289 lines.append("- [ %3d, %3d, %3d ]" % tuple(v)) 290 lines.append("") 291 return lines 292 293 def _primitive_yaml_lines(self, primitive, name): 294 lines = [] 295 if primitive is not None: 296 lines += self._cell_yaml_lines( 297 self.primitive, name, None) 298 lines.append(" reciprocal_lattice: # without 2pi") 299 rec_lat = np.linalg.inv(primitive.cell) 300 for v, a in zip(rec_lat.T, ('a*', 'b*', 'c*')): 301 lines.append(" - [ %21.15f, %21.15f, %21.15f ] # %s" % 302 (v[0], v[1], v[2], a)) 303 lines.append("") 304 return lines 305 306 def _unitcell_yaml_lines(self): 307 lines = [] 308 if self.unitcell is not None: 309 lines += self._cell_yaml_lines( 310 self.unitcell, "unit_cell", self.u2p_map) 311 lines.append("") 312 return lines 313 314 def _supercell_yaml_lines(self): 315 lines = [] 316 if self.supercell is not None: 317 s2p_map = getattr(self.primitive, 's2p_map', None) 318 lines += self._cell_yaml_lines( 319 self.supercell, "supercell", s2p_map) 320 lines.append("") 321 return lines 322 323 def _cell_yaml_lines(self, cell, name, map_to_primitive): 324 lines = [] 325 lines.append("%s:" % name) 326 count = 0 327 for line in cell.get_yaml_lines(): 328 lines.append(" " + line) 329 if map_to_primitive is not None and "mass" in line: 330 lines.append(" reduced_to: %d" % 331 (map_to_primitive[count] + 1)) 332 count += 1 333 return lines 334 335 def _nac_yaml_lines(self): 336 if self.primitive is None: 337 return [] 338 else: 339 return self._nac_yaml_lines_given_symbols(self.primitive.symbols) 340 341 def _nac_yaml_lines_given_symbols(self, symbols): 342 lines = [] 343 if self.nac_params is not None: 344 if self.settings['born_effective_charge']: 345 lines.append("born_effective_charge:") 346 for i, z in enumerate(self.nac_params['born']): 347 text = "- # %d" % (i + 1) 348 if symbols: 349 text += " (%s)" % symbols[i] 350 lines.append(text) 351 for v in z: 352 lines.append(" - [ %18.15f, %18.15f, %18.15f ]" % 353 tuple(v)) 354 lines.append("") 355 356 if self.settings['dielectric_constant']: 357 lines.append("dielectric_constant:") 358 for v in self.nac_params['dielectric']: 359 lines.append( 360 " - [ %18.15f, %18.15f, %18.15f ]" % tuple(v)) 361 lines.append("") 362 return lines 363 364 def _dataset_yaml_lines(self): 365 lines = [] 366 if self.settings['force_sets'] or self.settings['displacements']: 367 disp_yaml_lines = self._displacements_yaml_lines( 368 with_forces=self.settings['force_sets']) 369 lines += disp_yaml_lines 370 return lines 371 372 def _displacements_yaml_lines(self, with_forces=False): 373 return self._displacements_yaml_lines_2types( 374 self.dataset, with_forces=with_forces) 375 376 def _displacements_yaml_lines_2types(self, dataset, with_forces=False): 377 if dataset is not None: 378 if 'first_atoms' in dataset: 379 return self._displacements_yaml_lines_type1( 380 dataset, with_forces=with_forces) 381 elif 'displacements' in dataset: 382 return self._displacements_yaml_lines_type2( 383 dataset, with_forces=with_forces) 384 return [] 385 386 def _displacements_yaml_lines_type1(self, dataset, with_forces=False): 387 lines = ["displacements:", ] 388 for i, d in enumerate(dataset['first_atoms']): 389 lines.append("- atom: %4d" % (d['number'] + 1)) 390 lines.append(" displacement:") 391 lines.append(" [ %20.16f,%20.16f,%20.16f ]" 392 % tuple(d['displacement'])) 393 if with_forces and 'forces' in d: 394 lines.append(" forces:") 395 for f in d['forces']: 396 lines.append(" - [ %20.16f,%20.16f,%20.16f ]" % tuple(f)) 397 lines.append("") 398 return lines 399 400 def _displacements_yaml_lines_type2(self, dataset, with_forces=False): 401 if 'random_seed' in dataset: 402 lines = ["random_seed: %d" % dataset['random_seed'], 403 "displacements:"] 404 else: 405 lines = ["displacements:", ] 406 for i, dset in enumerate(dataset['displacements']): 407 lines.append("- # %4d" % (i + 1)) 408 for j, d in enumerate(dset): 409 lines.append(" - displacement: # %d" % (j + 1)) 410 lines.append(" [ %20.16f,%20.16f,%20.16f ]" % tuple(d)) 411 if with_forces and 'forces' in dataset: 412 f = dataset['forces'][i][j] 413 lines.append(" force:") 414 lines.append(" [ %20.16f,%20.16f,%20.16f ]" 415 % tuple(f)) 416 lines.append("") 417 return lines 418 419 def _force_constants_yaml_lines(self): 420 lines = [] 421 if (self.settings['force_constants'] and 422 self.force_constants is not None): 423 shape = self.force_constants.shape[:2] 424 lines = ["force_constants:", ] 425 if shape[0] == shape[1]: 426 lines.append(" format: \"full\"") 427 else: 428 lines.append(" format: \"compact\"") 429 lines.append(" shape: [ %d, %d ]" % shape) 430 lines.append(" elements:") 431 for (i, j) in list(np.ndindex(shape)): 432 lines.append(" - # (%d, %d)" % (i + 1, j + 1)) 433 for v in self.force_constants[i, j]: 434 lines.append(" - [ %21.15f, %21.15f, %21.15f ]" 435 % tuple(v)) 436 return lines 437 438 def _load(self, filename): 439 _, ext = os.path.splitext(filename) 440 if ext == '.xz' or ext == '.lzma': 441 try: 442 import lzma 443 except ImportError: 444 raise("Reading a lzma compressed file is not supported " 445 "by this python version.") 446 with lzma.open(filename) as f: 447 self._yaml = yaml.load(f, Loader=Loader) 448 elif ext == '.gz': 449 import gzip 450 with gzip.open(filename) as f: 451 self._yaml = yaml.load(f, Loader=Loader) 452 else: 453 with open(filename, 'r') as f: 454 self._yaml = yaml.load(f, Loader=Loader) 455 456 if type(self._yaml) is str: 457 msg = "Could not open %s's yaml file." % self.command_name 458 raise TypeError(msg) 459 460 self.parse() 461 462 def _parse_transformation_matrices(self): 463 if 'supercell_matrix' in self._yaml: 464 self.supercell_matrix = np.array(self._yaml['supercell_matrix'], 465 dtype='intc', order='C') 466 if 'primitive_matrix' in self._yaml: 467 self.primitive_matrix = np.array(self._yaml['primitive_matrix'], 468 dtype='double', order='C') 469 470 def _parse_all_cells(self): 471 if 'unit_cell' in self._yaml: 472 self.unitcell = self._parse_cell(self._yaml['unit_cell']) 473 if 'primitive_cell' in self._yaml: 474 self.primitive = self._parse_cell(self._yaml['primitive_cell']) 475 if 'supercell' in self._yaml: 476 self.supercell = self._parse_cell(self._yaml['supercell']) 477 if self.unitcell is None: 478 if ('lattice' in self._yaml and 479 ('points' in self._yaml or 'atoms' in self._yaml)): 480 self.unitcell = self._parse_cell(self._yaml) 481 482 def _parse_cell(self, cell_yaml): 483 lattice = None 484 if 'lattice' in cell_yaml: 485 lattice = cell_yaml['lattice'] 486 points = [] 487 symbols = [] 488 masses = [] 489 if 'points' in cell_yaml: 490 for x in cell_yaml['points']: 491 if 'coordinates' in x: 492 points.append(x['coordinates']) 493 if 'symbol' in x: 494 symbols.append(x['symbol']) 495 if 'mass' in x: 496 masses.append(x['mass']) 497 # For version < 1.10.9 498 elif 'atoms' in cell_yaml: 499 for x in cell_yaml['atoms']: 500 if 'coordinates' not in x and 'position' in x: 501 points.append(x['position']) 502 if 'symbol' in x: 503 symbols.append(x['symbol']) 504 if 'mass' in x: 505 masses.append(x['mass']) 506 return self._get_cell(lattice, points, symbols, masses=masses) 507 508 def _get_cell(self, lattice, points, symbols, masses=None): 509 if lattice: 510 _lattice = lattice 511 else: 512 _lattice = None 513 if points: 514 _points = points 515 else: 516 _points = None 517 if symbols: 518 _symbols = symbols 519 else: 520 _symbols = None 521 if masses: 522 _masses = masses 523 else: 524 _masses = None 525 526 if _lattice and _points and _symbols: 527 return PhonopyAtoms(symbols=_symbols, 528 cell=_lattice, 529 masses=_masses, 530 scaled_positions=_points) 531 else: 532 return None 533 534 def _parse_force_constants(self): 535 if 'force_constants' in self._yaml: 536 shape = tuple(self._yaml['force_constants']['shape']) + (3, 3) 537 fc = np.reshape(self._yaml['force_constants']['elements'], shape) 538 self.force_constants = np.array(fc, dtype='double', order='C') 539 540 def _parse_dataset(self): 541 self.dataset = self._get_dataset(self.supercell) 542 543 def _get_dataset(self, supercell): 544 dataset = None 545 if 'displacements' in self._yaml: 546 if supercell is not None: 547 natom = len(supercell) 548 else: 549 natom = None 550 disp = self._yaml['displacements'][0] 551 if type(disp) is dict: # type1 552 dataset = self._parse_force_sets_type1(natom=natom) 553 elif type(disp) is list: # type2 554 if 'displacement' in disp[0]: 555 dataset = self._parse_force_sets_type2() 556 return dataset 557 558 def _parse_force_sets_type1(self, natom=None): 559 with_forces = False 560 if 'forces' in self._yaml['displacements'][0]: 561 with_forces = True 562 dataset = {'natom': len(self._yaml['displacements'][0]['forces'])} 563 elif natom is not None: 564 dataset = {'natom': natom} 565 elif 'natom' in self._yaml: 566 dataset = {'natom': self._yaml['natom']} 567 else: 568 raise RuntimeError( 569 "Number of atoms in supercell could not be found.") 570 571 first_atoms = [] 572 for d in self._yaml['displacements']: 573 data = { 574 'number': d['atom'] - 1, 575 'displacement': np.array(d['displacement'], dtype='double')} 576 if with_forces: 577 data['forces'] = np.array(d['forces'], 578 dtype='double', order='C') 579 first_atoms.append(data) 580 dataset['first_atoms'] = first_atoms 581 582 return dataset 583 584 def _parse_force_sets_type2(self): 585 nsets = len(self._yaml['displacements']) 586 natom = len(self._yaml['displacements'][0]) 587 if 'force' in self._yaml['displacements'][0][0]: 588 with_forces = True 589 forces = np.zeros((nsets, natom, 3), dtype='double', order='C') 590 else: 591 with_forces = False 592 displacements = np.zeros((nsets, natom, 3), dtype='double', order='C') 593 for i, dfset in enumerate(self._yaml['displacements']): 594 for j, df in enumerate(dfset): 595 if with_forces: 596 forces[i, j] = df['force'] 597 displacements[i, j] = df['displacement'] 598 if with_forces: 599 return {'forces': forces, 'displacements': displacements} 600 else: 601 return {'displacements': displacements} 602 603 def _parse_nac_params(self): 604 nac_params = {} 605 if 'born_effective_charge' in self._yaml: 606 nac_params['born'] = np.array(self._yaml['born_effective_charge'], 607 dtype='double', order='C') 608 if 'dielectric_constant' in self._yaml: 609 nac_params['dielectric'] = np.array( 610 self._yaml['dielectric_constant'], dtype='double', order='C') 611 if (self.command_name in self._yaml and 612 'nac_unit_conversion_factor' in self._yaml[self.command_name]): 613 nac_params['factor'] = self._yaml[self.command_name][ 614 'nac_unit_conversion_factor'] 615 if 'born' in nac_params and 'dielectric' in nac_params: 616 self.nac_params = nac_params 617 618 def _parse_calculator(self): 619 if (self.command_name in self._yaml and 620 'calculator' in self._yaml[self.command_name]): 621 self.calculator = self._yaml[self.command_name]['calculator'] 622