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 sys
37import numpy as np
38from phonopy.harmonic.force_constants import (
39    show_drift_force_constants,
40    symmetrize_force_constants,
41    symmetrize_compact_force_constants)
42from phonopy.file_IO import get_dataset_type2, parse_FORCE_SETS
43from phonopy.cui.phonopy_script import print_error, file_exists
44from phonopy.interface.calculator import get_default_physical_units
45from phono3py.phonon3.fc3 import show_drift_fc3
46from phono3py.file_IO import (
47    parse_disp_fc3_yaml, parse_disp_fc2_yaml, parse_FORCES_FC2,
48    parse_FORCES_FC3, read_fc3_from_hdf5, read_fc2_from_hdf5,
49    write_fc3_to_hdf5, write_fc2_to_hdf5, get_length_of_first_line)
50from phono3py.cui.show_log import show_phono3py_force_constants_settings
51from phono3py.phonon3.fc3 import (
52    set_permutation_symmetry_fc3, set_translational_invariance_fc3)
53from phono3py.interface.phono3py_yaml import Phono3pyYaml
54
55
56def create_phono3py_force_constants(phono3py,
57                                    settings,
58                                    ph3py_yaml=None,
59                                    input_filename=None,
60                                    output_filename=None,
61                                    phono3py_yaml_filename=None,
62                                    log_level=1):
63    if settings.fc_calculator is None:
64        symmetrize_fc3r = (settings.is_symmetrize_fc3_r or
65                           settings.fc_symmetry)
66        symmetrize_fc2 = (settings.is_symmetrize_fc2 or
67                          settings.fc_symmetry)
68    else:  # Rely on fc calculator the symmetrization of fc.
69        symmetrize_fc2 = False
70        symmetrize_fc3r = False
71
72    if log_level:
73        show_phono3py_force_constants_settings(settings)
74
75    #######
76    # fc3 #
77    #######
78    if (settings.is_joint_dos or
79        (settings.is_isotope and
80         not (settings.is_bterta or settings.is_lbte)) or
81        settings.read_gamma or
82        settings.read_pp or
83        (not settings.is_bterta and settings.write_phonon) or
84        settings.constant_averaged_pp_interaction is not None):
85        pass
86    else:
87        if settings.read_fc3:
88            _read_phono3py_fc3(phono3py,
89                               symmetrize_fc3r,
90                               input_filename,
91                               log_level)
92        else:  # fc3 from FORCES_FC3 or ph3py_yaml
93            _create_phono3py_fc3(phono3py,
94                                 ph3py_yaml,
95                                 symmetrize_fc3r,
96                                 symmetrize_fc2,
97                                 input_filename,
98                                 output_filename,
99                                 settings.is_compact_fc,
100                                 settings.cutoff_pair_distance,
101                                 settings.fc_calculator,
102                                 settings.fc_calculator_options,
103                                 log_level)
104            if output_filename is None:
105                filename = 'fc3.hdf5'
106            else:
107                filename = 'fc3.' + output_filename + '.hdf5'
108            if log_level:
109                print("Writing fc3 to \"%s\"." % filename)
110            write_fc3_to_hdf5(phono3py.fc3,
111                              filename=filename,
112                              p2s_map=phono3py.primitive.p2s_map,
113                              compression=settings.hdf5_compression)
114
115        cutoff_distance = settings.cutoff_fc3_distance
116        if cutoff_distance is not None and cutoff_distance > 0:
117            if log_level:
118                print("Cutting-off fc3 by zero (cut-off distance: %f)" %
119                      cutoff_distance)
120            phono3py.cutoff_fc3_by_zero(cutoff_distance)
121
122        if log_level:
123            show_drift_fc3(phono3py.fc3, primitive=phono3py.primitive)
124
125    #######
126    # fc2 #
127    #######
128    phonon_primitive = phono3py.phonon_primitive
129    p2s_map = phonon_primitive.p2s_map
130    if settings.read_fc2:
131        _read_phono3py_fc2(phono3py,
132                           symmetrize_fc2,
133                           input_filename,
134                           log_level)
135    else:
136        if phono3py.phonon_supercell_matrix is None:
137            if (settings.fc_calculator == 'alm' and phono3py.fc2 is not None):
138                if log_level:
139                    print("fc2 that was fit simultaneously with fc3 "
140                          "by ALM is used.")
141            else:
142                _create_phono3py_fc2(phono3py,
143                                     ph3py_yaml,
144                                     symmetrize_fc2,
145                                     input_filename,
146                                     settings.is_compact_fc,
147                                     settings.fc_calculator,
148                                     settings.fc_calculator_options,
149                                     log_level)
150        else:
151            _create_phono3py_phonon_fc2(phono3py,
152                                        ph3py_yaml,
153                                        symmetrize_fc2,
154                                        input_filename,
155                                        settings.is_compact_fc,
156                                        settings.fc_calculator,
157                                        settings.fc_calculator_options,
158                                        log_level)
159        if output_filename is None:
160            filename = 'fc2.hdf5'
161        else:
162            filename = 'fc2.' + output_filename + '.hdf5'
163        if log_level:
164            print("Writing fc2 to \"%s\"." % filename)
165        write_fc2_to_hdf5(phono3py.fc2,
166                          filename=filename,
167                          p2s_map=p2s_map,
168                          physical_unit='eV/angstrom^2',
169                          compression=settings.hdf5_compression)
170
171    if log_level:
172        show_drift_force_constants(phono3py.fc2,
173                                   primitive=phonon_primitive,
174                                   name='fc2')
175
176
177def parse_forces(phono3py,
178                 ph3py_yaml=None,
179                 cutoff_pair_distance=None,
180                 force_filename="FORCES_FC3",
181                 disp_filename=None,
182                 fc_type=None,
183                 log_level=0):
184    filename_read_from = None
185
186    if fc_type == 'phonon_fc2':
187        natom = len(phono3py.phonon_supercell)
188    else:
189        natom = len(phono3py.supercell)
190
191    # Get dataset from ph3py_yaml. dataset can be None.
192    dataset = _extract_datast_from_ph3py_yaml(ph3py_yaml, fc_type)
193    if dataset:
194        filename_read_from = ph3py_yaml.yaml_filename
195
196    # Try to read FORCES_FC* if type-2 and return dataset.
197    # None is returned unless type-2.
198    # can emit FileNotFoundError.
199    if (dataset is None or
200        dataset is not None and not forces_in_dataset(dataset)):
201        _dataset = _get_type2_dataset(natom,
202                                      phono3py.calculator,
203                                      filename=force_filename,
204                                      log_level=log_level)
205        # Do not overwrite dataset when _dataset is None.
206        if _dataset:
207            filename_read_from = force_filename
208            dataset = _dataset
209
210    if dataset is None:
211        # Displacement dataset is obtained from disp_filename.
212        # can emit FileNotFoundError.
213        dataset = _read_disp_fc_yaml(disp_filename, fc_type)
214        filename_read_from = disp_filename
215
216    if 'natom' in dataset and dataset['natom'] != natom:
217        msg = ("Number of atoms in supercell is not consistent with "
218               "\"%s\"." % filename_read_from)
219        raise RuntimeError(msg)
220
221    if log_level and filename_read_from is not None:
222        print("Displacement dataset for %s was read from \"%s\"."
223              % (fc_type, filename_read_from))
224
225    if cutoff_pair_distance:
226        if ('cutoff_distance' not in dataset or
227            ('cutoff_distance' in dataset and
228             cutoff_pair_distance < dataset['cutoff_distance'])):
229            dataset['cutoff_distance'] = cutoff_pair_distance
230            if log_level:
231                print("Cutoff-pair-distance: %f" % cutoff_pair_distance)
232
233    # Type-1 FORCES_FC*.
234    # dataset comes either from disp_fc*.yaml or phono3py*.yaml.
235    if not forces_in_dataset(dataset):
236        if fc_type == 'phonon_fc2':
237            parse_FORCES_FC2(dataset, filename=force_filename)
238        else:
239            parse_FORCES_FC3(dataset, filename=force_filename)
240
241        if log_level:
242            print("Sets of supercell forces were read from \"%s\"."
243                  % force_filename)
244            sys.stdout.flush()
245
246    _convert_unit_in_dataset(dataset, phono3py.calculator)
247
248    return dataset
249
250
251def forces_in_dataset(dataset):
252    return ('forces' in dataset or
253            ('first_atoms' in dataset and
254             'forces' in dataset['first_atoms'][0]))
255
256
257def _read_disp_fc_yaml(disp_filename, fc_type):
258    if fc_type == 'phonon_fc2':
259        dataset = parse_disp_fc2_yaml(filename=disp_filename)
260    else:
261        dataset = parse_disp_fc3_yaml(filename=disp_filename)
262
263    return dataset
264
265
266def _read_phono3py_fc3(phono3py,
267                       symmetrize_fc3r,
268                       input_filename,
269                       log_level):
270    if input_filename is None:
271        filename = 'fc3.hdf5'
272    else:
273        filename = 'fc3.' + input_filename + '.hdf5'
274    file_exists(filename, log_level)
275    if log_level:
276        print("Reading fc3 from \"%s\"." % filename)
277
278    p2s_map = phono3py.primitive.p2s_map
279    try:
280        fc3 = read_fc3_from_hdf5(filename=filename, p2s_map=p2s_map)
281    except RuntimeError:
282        import traceback
283        traceback.print_exc()
284        if log_level:
285            print_error()
286        sys.exit(1)
287    num_atom = phono3py.supercell.get_number_of_atoms()
288    if fc3.shape[1] != num_atom:
289        print("Matrix shape of fc3 doesn't agree with supercell size.")
290        if log_level:
291            print_error()
292        sys.exit(1)
293
294    if symmetrize_fc3r:
295        set_translational_invariance_fc3(fc3)
296        set_permutation_symmetry_fc3(fc3)
297
298    phono3py.fc3 = fc3
299
300
301def _read_phono3py_fc2(phono3py,
302                       symmetrize_fc2,
303                       input_filename,
304                       log_level):
305    if input_filename is None:
306        filename = 'fc2.hdf5'
307    else:
308        filename = 'fc2.' + input_filename + '.hdf5'
309    file_exists(filename, log_level)
310    if log_level:
311        print("Reading fc2 from \"%s\"." % filename)
312
313    num_atom = phono3py.phonon_supercell.get_number_of_atoms()
314    p2s_map = phono3py.phonon_primitive.p2s_map
315    try:
316        phonon_fc2 = read_fc2_from_hdf5(filename=filename, p2s_map=p2s_map)
317    except RuntimeError:
318        import traceback
319        traceback.print_exc()
320        if log_level:
321            print_error()
322        sys.exit(1)
323
324    if phonon_fc2.shape[1] != num_atom:
325        print("Matrix shape of fc2 doesn't agree with supercell size.")
326        if log_level:
327            print_error()
328        sys.exit(1)
329
330    if symmetrize_fc2:
331        if phonon_fc2.shape[0] == phonon_fc2.shape[1]:
332            symmetrize_force_constants(phonon_fc2)
333        else:
334            symmetrize_compact_force_constants(phonon_fc2,
335                                               phono3py.phonon_primitive)
336
337    phono3py.fc2 = phonon_fc2
338
339
340def _get_type2_dataset(natom, calculator, filename="FORCES_FC3", log_level=0):
341    if not os.path.isfile(filename):
342        return None
343
344    with open(filename, 'r') as f:
345        len_first_line = get_length_of_first_line(f)
346        if len_first_line == 6:
347            dataset = get_dataset_type2(f, natom)
348            if log_level:
349                print("%d snapshots were found in %s."
350                      % (len(dataset['displacements']), "FORCES_FC3"))
351        else:
352            dataset = None
353    return dataset
354
355
356def _create_phono3py_fc3(phono3py,
357                         ph3py_yaml,
358                         symmetrize_fc3r,
359                         symmetrize_fc2,
360                         input_filename,
361                         output_filename,
362                         is_compact_fc,
363                         cutoff_pair_distance,
364                         fc_calculator,
365                         fc_calculator_options,
366                         log_level):
367    """
368
369    Note
370    ----
371    cutoff_pair_distance is the parameter to determine each displaced
372    supercell is included to the computation of fc3. It is assumed that
373    cutoff_pair_distance is stored in the step to create sets of
374    displacements and the value is stored n the displacement dataset and
375    also as the parameter 'included': True or False for each displacement.
376    The parameter cutoff_pair_distance here can be used in the step to
377    create fc3 by overwriting original cutoff_pair_distance value only
378    when the former value is smaller than the later.
379
380    """
381    if input_filename is None:
382        disp_filename = 'disp_fc3.yaml'
383    else:
384        disp_filename = 'disp_fc3.' + input_filename + '.yaml'
385
386    _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml)
387
388    try:
389        dataset = parse_forces(phono3py,
390                               ph3py_yaml=_ph3py_yaml,
391                               cutoff_pair_distance=cutoff_pair_distance,
392                               force_filename="FORCES_FC3",
393                               disp_filename=disp_filename,
394                               fc_type='fc3',
395                               log_level=log_level)
396    except RuntimeError as e:
397        # from _parse_forces_type1
398        if log_level:
399            print(str(e))
400            print_error()
401        sys.exit(1)
402    except FileNotFoundError as e:
403        # from _get_type2_dataset
404        file_exists(e.filename, log_level)
405
406    phono3py.dataset = dataset
407    phono3py.produce_fc3(symmetrize_fc3r=symmetrize_fc3r,
408                         is_compact_fc=is_compact_fc,
409                         fc_calculator=fc_calculator,
410                         fc_calculator_options=fc_calculator_options)
411
412
413def _create_phono3py_fc2(phono3py,
414                         ph3py_yaml,
415                         symmetrize_fc2,
416                         input_filename,
417                         is_compact_fc,
418                         fc_calculator,
419                         fc_calculator_options,
420                         log_level):
421    if input_filename is None:
422        disp_filename = 'disp_fc3.yaml'
423    else:
424        disp_filename = 'disp_fc3.' + input_filename + '.yaml'
425
426    _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml)
427
428    try:
429        dataset = parse_forces(phono3py,
430                               ph3py_yaml=_ph3py_yaml,
431                               force_filename="FORCES_FC3",
432                               disp_filename=disp_filename,
433                               fc_type='fc2',
434                               log_level=log_level)
435    except RuntimeError as e:
436        if log_level:
437            print(str(e))
438            print_error()
439        sys.exit(1)
440    except FileNotFoundError as e:
441        file_exists(e.filename, log_level)
442
443    phono3py.phonon_dataset = dataset
444    phono3py.produce_fc2(
445        symmetrize_fc2=symmetrize_fc2,
446        is_compact_fc=is_compact_fc,
447        fc_calculator=fc_calculator,
448        fc_calculator_options=fc_calculator_options)
449
450
451def _get_ph3py_yaml(disp_filename, ph3py_yaml):
452    _ph3py_yaml = ph3py_yaml
453    # Try to use phono3py.phonon_dataset when the disp file not found
454    if not os.path.isfile(disp_filename):
455        disp_filename = None
456        if _ph3py_yaml is None and os.path.isfile('phono3py_disp.yaml'):
457            _ph3py_yaml = Phono3pyYaml()
458            _ph3py_yaml.read('phono3py_disp.yaml')
459    return _ph3py_yaml
460
461
462def _create_phono3py_phonon_fc2(phono3py,
463                                ph3py_yaml,
464                                symmetrize_fc2,
465                                input_filename,
466                                is_compact_fc,
467                                fc_calculator,
468                                fc_calculator_options,
469                                log_level):
470    if input_filename is None:
471        disp_filename = 'disp_fc2.yaml'
472    else:
473        disp_filename = 'disp_fc2.' + input_filename + '.yaml'
474
475    _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml)
476
477    try:
478        dataset = parse_forces(phono3py,
479                               ph3py_yaml=_ph3py_yaml,
480                               force_filename="FORCES_FC2",
481                               disp_filename=disp_filename,
482                               fc_type='phonon_fc2',
483                               log_level=log_level)
484    except RuntimeError as e:
485        if log_level:
486            print(str(e))
487            print_error()
488        sys.exit(1)
489    except FileNotFoundError as e:
490        file_exists(e.filename, log_level)
491
492    phono3py.phonon_dataset = dataset
493    phono3py.produce_fc2(
494        symmetrize_fc2=symmetrize_fc2,
495        is_compact_fc=is_compact_fc,
496        fc_calculator=fc_calculator,
497        fc_calculator_options=fc_calculator_options)
498
499
500def _convert_unit_in_dataset(dataset, calculator):
501    physical_units = get_default_physical_units(calculator)
502    force_to_eVperA = physical_units['force_to_eVperA']
503    distance_to_A = physical_units['distance_to_A']
504
505    if 'first_atoms' in dataset:
506        for d1 in dataset['first_atoms']:
507            if distance_to_A is not None:
508                disp = _to_ndarray(d1['displacement'])
509                d1['displacement'] = disp * distance_to_A
510            if force_to_eVperA is not None and 'forces' in d1:
511                forces = _to_ndarray(d1['forces'])
512                d1['forces'] = forces * force_to_eVperA
513            if 'second_atoms' in d1:
514                for d2 in d1['second_atoms']:
515                    if distance_to_A is not None:
516                        disp = _to_ndarray(d2['displacement'])
517                        d2['displacement'] = disp * distance_to_A
518                    if force_to_eVperA is not None and 'forces' in d2:
519                        forces = _to_ndarray(d2['forces'])
520                        d2['forces'] = forces * force_to_eVperA
521    else:
522        if distance_to_A is not None and 'displacements' in dataset:
523            disp = _to_ndarray(dataset['displacements'])
524            dataset['displacements'] = disp * distance_to_A
525        if force_to_eVperA is not None and 'forces' in dataset:
526            forces = _to_ndarray(dataset['forces'])
527            dataset['forces'] = forces * force_to_eVperA
528
529
530def _to_ndarray(array, dtype='double'):
531    if type(array) is not np.ndarray:
532        return np.array(array, dtype=dtype, order='C')
533    else:
534        return array
535
536
537def _extract_datast_from_ph3py_yaml(ph3py_yaml, fc_type):
538    dataset = None
539    if fc_type == 'phonon_fc2':
540        if ph3py_yaml and ph3py_yaml.phonon_dataset is not None:
541            # copy dataset
542            # otherwise unit conversion can be applied multiple times.
543            _ph3py_yaml = Phono3pyYaml()
544            _ph3py_yaml.yaml_data = ph3py_yaml.yaml_data
545            _ph3py_yaml.parse()
546            dataset = _ph3py_yaml.phonon_dataset
547    else:
548        if ph3py_yaml and ph3py_yaml.dataset is not None:
549            # copy dataset
550            # otherwise unit conversion can be applied multiple times.
551            _ph3py_yaml = Phono3pyYaml()
552            _ph3py_yaml.yaml_data = ph3py_yaml.yaml_data
553            _ph3py_yaml.parse()
554            dataset = _ph3py_yaml.dataset
555    return dataset
556