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