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
37from phonopy.interface.calculator import read_crystal_structure
38from phonopy.structure.cells import get_primitive_matrix
39from phonopy.file_IO import (
40    parse_BORN, read_force_constants_hdf5, parse_FORCE_SETS,
41    parse_FORCE_CONSTANTS)
42from phonopy.structure.atoms import PhonopyAtoms
43from phonopy.interface.calculator import get_force_constant_conversion_factor
44
45
46def get_cell_settings(phonopy_yaml=None,
47                      supercell_matrix=None,
48                      primitive_matrix=None,
49                      unitcell=None,
50                      supercell=None,
51                      unitcell_filename=None,
52                      supercell_filename=None,
53                      calculator=None,
54                      symprec=1e-5,
55                      log_level=0):
56    optional_structure_info = None
57    if (primitive_matrix is None or
58        (type(primitive_matrix) is str and primitive_matrix == "auto")):
59        pmat = 'auto'
60    else:
61        pmat = primitive_matrix
62
63    if unitcell_filename is not None:
64        cell, optional_structure_info = _read_crystal_structure(
65            filename=unitcell_filename, interface_mode=calculator)
66        smat = supercell_matrix
67        if log_level:
68            print("Unit cell structure was read from \"%s\"."
69                  % optional_structure_info[0])
70    elif supercell_filename is not None:
71        cell, optional_structure_info = read_crystal_structure(
72            filename=supercell_filename, interface_mode=calculator)
73        smat = np.eye(3, dtype='intc', order='C')
74        if log_level:
75            print("Supercell structure was read from \"%s\"."
76                  % optional_structure_info[0])
77    elif unitcell is not None:
78        cell = PhonopyAtoms(atoms=unitcell)
79        smat = supercell_matrix
80    elif supercell is not None:
81        cell = PhonopyAtoms(atoms=supercell)
82        smat = np.eye(3, dtype='intc', order='C')
83    else:
84        raise RuntimeError("Cell has to be specified.")
85
86    if optional_structure_info is not None and cell is None:
87        filename = optional_structure_info[0]
88        msg = "'%s' could not be found." % filename
89        raise FileNotFoundError(msg)
90
91    pmat = get_primitive_matrix(pmat, symprec=symprec)
92
93    return cell, smat, pmat
94
95
96def get_nac_params(primitive=None,
97                   nac_params=None,
98                   born_filename=None,
99                   is_nac=True,
100                   nac_factor=None,
101                   log_level=0):
102    if born_filename is not None:
103        _nac_params = parse_BORN(primitive, filename=born_filename)
104        if log_level:
105            print("NAC parameters were read from \"%s\"." % born_filename)
106    elif nac_params is not None:  # nac_params input or phonopy_yaml.nac_params
107        _nac_params = nac_params
108    elif is_nac and os.path.isfile("BORN"):
109        _nac_params = parse_BORN(primitive, filename="BORN")
110        if log_level:
111            print("NAC params were read from \"BORN\".")
112    else:
113        _nac_params = None
114
115    if _nac_params is not None:
116        if 'factor' not in _nac_params or _nac_params['factor'] is None:
117            _nac_params['factor'] = nac_factor
118
119    return _nac_params
120
121
122def read_force_constants_from_hdf5(filename='force_constants.hdf5',
123                                   p2s_map=None,
124                                   calculator=None):
125    """Convert force constants physical unit
126
127    Each calculator interface has own default force constants physical unit.
128    This method reads 'physical_unit' in force constants hdf5 file and
129    if this is different from the one for 'calculator', the force constants
130    are converted to have the physical unit of the calculator.
131
132    Note
133    ----
134    This method is also used from phonopy script.
135
136    """
137
138    fc, fc_unit = read_force_constants_hdf5(filename=filename,
139                                            p2s_map=p2s_map,
140                                            return_physical_unit=True)
141    if fc_unit is None:
142        return fc
143    else:
144        factor = get_force_constant_conversion_factor(fc_unit, calculator)
145        return fc * factor
146
147
148def set_dataset_and_force_constants(
149        phonon,
150        dataset,
151        fc,  # From phonopy_yaml
152        force_constants_filename=None,
153        force_sets_filename=None,
154        fc_calculator=None,
155        fc_calculator_options=None,
156        produce_fc=True,
157        symmetrize_fc=True,
158        is_compact_fc=True,
159        log_level=0):
160    natom = len(phonon.supercell)
161
162    # dataset and fc are those obtained from phonopy_yaml unless None.
163    if dataset is not None:
164        phonon.dataset = dataset
165    if fc is not None:
166        phonon.force_constants = fc
167
168    _fc = None
169    _dataset = None
170    if force_constants_filename is not None:
171        _fc = _read_force_constants_file(phonon, force_constants_filename)
172        _force_constants_filename = force_constants_filename
173    elif force_sets_filename is not None:
174        _dataset = parse_FORCE_SETS(natom=natom, filename=force_sets_filename)
175        _force_sets_filename = force_sets_filename
176    elif phonon.forces is None and phonon.force_constants is None:
177        # unless provided these from phonopy_yaml.
178        if os.path.isfile("FORCE_CONSTANTS"):
179            _fc = _read_force_constants_file(phonon, "FORCE_CONSTANTS")
180            _force_constants_filename = "FORCE_CONSTANTS"
181        elif os.path.isfile("force_constants.hdf5"):
182            _fc = _read_force_constants_file(phonon, "force_constants.hdf5")
183            _force_constants_filename = "force_constants.hdf5"
184        elif os.path.isfile("FORCE_SETS"):
185            _dataset = parse_FORCE_SETS(natom=natom)
186            _force_sets_filename = "FORCE_SETS"
187
188    if _fc is not None:
189        phonon.force_constants = _fc
190        if log_level:
191            print("Force constants were read from \"%s\"."
192                  % _force_constants_filename)
193
194    if phonon.force_constants is None:
195        # Overwrite dataset
196        if _dataset is not None:
197            phonon.dataset = _dataset
198            if log_level:
199                print("Force sets were read from \"%s\"."
200                      % _force_sets_filename)
201        if produce_fc:
202            _produce_force_constants(phonon,
203                                     fc_calculator,
204                                     fc_calculator_options,
205                                     symmetrize_fc,
206                                     is_compact_fc,
207                                     log_level)
208
209
210def _read_force_constants_file(phonon, force_constants_filename):
211    dot_split = force_constants_filename.split('.')
212    p2s_map = phonon.primitive.p2s_map
213    if len(dot_split) > 1 and dot_split[-1] == 'hdf5':
214        _fc = read_force_constants_from_hdf5(
215            filename=force_constants_filename,
216            p2s_map=p2s_map,
217            calculator=phonon.calculator)
218    else:
219        _fc = parse_FORCE_CONSTANTS(filename=force_constants_filename,
220                                    p2s_map=p2s_map)
221    return _fc
222
223
224def _produce_force_constants(phonon,
225                             fc_calculator,
226                             fc_calculator_options,
227                             symmetrize_fc,
228                             is_compact_fc,
229                             log_level):
230    phonon.produce_force_constants(
231        calculate_full_force_constants=(not is_compact_fc),
232        fc_calculator=fc_calculator,
233        fc_calculator_options=fc_calculator_options)
234    if symmetrize_fc:
235        phonon.symmetrize_force_constants(show_drift=(log_level > 0))
236        if log_level:
237            print("Force constants were symmetrized.")
238
239
240def _read_crystal_structure(filename=None, interface_mode=None):
241    try:
242        return read_crystal_structure(filename=filename,
243                                      interface_mode=interface_mode)
244    except FileNotFoundError:
245        raise
246    except:
247        print("============================ phonopy.load ============================")
248        print("  Reading crystal structure file failed in phonopy.load.")
249        print("  Maybe phonopy.load(..., calculator='<calculator name>') expected?")
250        print("============================ phonopy.load ============================")
251        raise
252