#!/usr/bin/env python

# Copyright (C) 2011 Atsushi Togo
# All rights reserved.
#
# This file is part of phonopy.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
#   notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
#   notice, this list of conditions and the following disclaimer in
#   the documentation and/or other materials provided with the
#   distribution.
#
# * Neither the name of the phonopy project nor the names of its
#   contributors may be used to endorse or promote products derived
#   from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

import sys
import numpy as np
from phonopy.file_IO import (get_born_parameters, write_FORCE_SETS,
                             parse_FORCE_SETS)
from phonopy.units import VaspToTHz, Bohr, Rydberg, Hartree
from phonopy.cui.collect_cell_info import collect_cell_info
from phonopy.interface.calculator import (
    get_force_sets, get_default_physical_units, get_interface_mode)
from phonopy.cui.phonopy_argparse import show_deprecated_option_warnings
from phonopy.phonon.band_structure import get_band_qpoints
from phono3py.version import __version__
from phono3py.file_IO import (parse_disp_fc2_yaml, parse_disp_fc3_yaml,
                              write_disp_fc2_yaml, read_phonon_from_hdf5,
                              write_phonon_to_hdf5,
                              parse_FORCES_FC2,
                              write_FORCES_FC2, write_FORCES_FC3,
                              get_lenghth_of_first_line)
from phono3py.cui.settings import Phono3pyConfParser
from phono3py import Phono3py, Phono3pyJointDos, Phono3pyIsotope
from phono3py.phonon3.gruneisen import get_gruneisen_parameters
from phono3py.cui.phono3py_argparse import get_parser
from phono3py.cui.show_log import (print_phono3py, print_version, print_end,
                                   print_error, print_error_message,
                                   show_general_settings,
                                   show_phono3py_settings,
                                   show_phono3py_cells, file_exists)
from phono3py.cui.triplets_info import write_grid_points, show_num_triplets
from phono3py.cui.translate_settings import get_phono3py_configurations
from phono3py.cui.create_supercells import create_phono3py_supercells
from phono3py.cui.create_force_constants import create_phono3py_force_constants
from phono3py.cui.phono3py_yaml import Phono3pyYaml


# import logging
# logging.basicConfig()
# logging.getLogger("phono3py.phonon3.fc3").setLevel(level=logging.DEBUG)


def finalize_phono3py(log_level,
                      phonop3y_conf,
                      phono3py,
                      interface_mode,
                      filename="phono3py.yaml"):
    if log_level > 0:
        units = get_default_physical_units(interface_mode)
        ph3py_yaml = Phono3pyYaml(configuration=phono3py_conf.get_configures(),
                                  calculator=interface_mode,
                                  physical_units=units)
        ph3py_yaml.set_phonon_info(phono3py)
        with open(filename, 'w') as w:
            w.write(str(ph3py_yaml))
        print_end()
    sys.exit(0)


def get_run_mode(settings):
    run_mode = None
    if settings.get_is_gruneisen():
        run_mode = "gruneisen"
    elif settings.get_is_joint_dos():
        run_mode = "jdos"
    elif (settings.get_is_isotope() and
          not (settings.get_is_bterta() or settings.get_is_lbte())):
        run_mode = "isotope"
    elif settings.get_is_imag_self_energy():
        run_mode = "imag_self_energy"
    elif settings.get_is_frequency_shift():
        run_mode = "frequency_shift"
    elif settings.get_is_bterta():
        run_mode = "conductivity-RTA"
    elif settings.get_is_lbte():
        run_mode = "conductivity-LBTE"
    elif settings.get_create_displacements():
        run_mode = "displacements"
    elif settings.get_write_phonon():
        run_mode = "phonon"
    return run_mode


# Parse arguments
parser, deprecated = get_parser()
args = parser.parse_args()

# Log level
log_level = 1
if args.verbose:
    log_level = 2
if args.quiet:
    log_level = 0
if args.log_level is not None:
    log_level = args.log_level

# Title
if log_level:
    print_phono3py()
    print_version(__version__)
    print("Python version %d.%d.%d" % sys.version_info[:3])
    import phonopy.structure.spglib as spglib
    print("Spglib version %d.%d.%d" % spglib.get_version())

if deprecated:
    show_deprecated_option_warnings(deprecated)

#
# Phono3py interface mode
#
# Physical units: energy,  distance,  atomic mass, force
# vasp          : eV,      Angstrom,  AMU,         eV/Angstrom
# pwscf         : Ry,      au,        AMU,         Ry/au
# CRYSTAL       : eV,      Angstrom,  AMU,         eV/Angstrom
# abinit        : hartree, au,        AMU,         eV/Angstrom
# TURBOMOLE     : hartree, au,        AMU,         hartree/au
#
if args.qe_mode:
    interface_mode = 'qe'
    try:
        from phonopy.interface.qe import write_supercells_with_displacements
    except ImportError:  # to be compatible to older phonopy
        from phonopy.interface.pwscf import write_supercells_with_displacements
elif args.crystal_mode:
    interface_mode = 'crystal'
    from phonopy.interface.crystal import write_supercells_with_displacements
elif args.abinit_mode:
    interface_mode = 'abinit'
    from phonopy.interface.abinit import write_supercells_with_displacements
elif args.turbomole_mode:
    interface_mode = 'turbomole'
    from phonopy.interface.turbomole import write_supercells_with_displacements
else:
    interface_mode = 'vasp'
    from phonopy.interface.vasp import write_supercells_with_displacements

units = get_default_physical_units(interface_mode)
distance_to_A = units['distance_to_A']
force_to_eVperA = units['force_to_eVperA']


# Input and output filename extension
input_filename = args.input_filename
output_filename = args.output_filename
if args.input_output_filename is not None:
    input_filename = args.input_output_filename
    output_filename = args.input_output_filename

# HDF5 compression filter
try:
    compression = int(args.hdf5_compression)
except ValueError:  # str
    compression = args.hdf5_compression
    if compression.lower() == "none":
        compression = None
except TypeError:  # None
    compression = args.hdf5_compression

#####################
# Create FORCES_FC3 #
#####################
if args.forces_fc3 or args.forces_fc3_file:
    if input_filename is None:
        disp_filename = 'disp_fc3.yaml'
    else:
        disp_filename = 'disp_fc3.' + input_filename + '.yaml'
    file_exists(disp_filename, log_level)
    if log_level:
        print("Displacement dataset is read from %s." % disp_filename)
    disp_dataset = parse_disp_fc3_yaml(filename=disp_filename)
    num_atoms = disp_dataset['natom']
    num_disps = len(disp_dataset['first_atoms'])
    for d1 in disp_dataset['first_atoms']:
        for d2 in d1['second_atoms']:
            if 'included' not in d2 or d2['included']:
                num_disps += 1
    if args.forces_fc3_file:
        file_exists(args.forces_fc3_file, log_level)
        force_filenames = [x.strip() for x in open(args.forces_fc3_file)]
    else:
        force_filenames = args.forces_fc3

    for filename in force_filenames:
        file_exists(filename, log_level)

    if log_level > 0:
        print("Number of displacements: %d" % num_disps)
        print("Number of supercell files: %d" % len(force_filenames))

    force_sets = get_force_sets(interface_mode,
                                num_atoms,
                                num_disps,
                                force_filenames,
                                disp_filename=disp_filename,
                                check_number_of_files=True,
                                verbose=(log_level > 0))

    if args.forces_fcz:
        force_filename = args.forces_fcz
        file_exists(force_filename, log_level)
        force_set_zero = get_force_sets(interface_mode,
                                        num_atoms,
                                        1,
                                        [force_filename, ],
                                        verbose=(log_level > 0))[0]
        for fs in force_sets:
            fs -= force_set_zero

        if log_level > 0:
            print("Forces in \'%s\' were subtracted from supercell forces."
                  % force_filename)

    if force_sets:
        write_FORCES_FC3(disp_dataset, force_sets, filename="FORCES_FC3")
        if log_level:
            print("")
            print("%s has been created." % "FORCES_FC3")
            print_end()
        sys.exit(0)
    else:
        if log_level:
            print("")
            print("%s could not be created." % "FORCES_FC3")
            print_error()
        sys.exit(1)

#####################
# Create FORCES_FC2 #
#####################
if args.forces_fc2:
    if input_filename is None:
        disp_filename = 'disp_fc2.yaml'
    else:
        disp_filename = 'disp_fc2.' + input_filename + '.yaml'
    file_exists(disp_filename, log_level)
    if log_level:
        print("Displacement dataset is read from %s." % disp_filename)
    disp_dataset = parse_disp_fc2_yaml(filename=disp_filename)
    num_atoms = disp_dataset['natom']
    num_disps = len(disp_dataset['first_atoms'])
    force_filenames = args.forces_fc2
    for filename in force_filenames:
        file_exists(filename, log_level)

    if log_level > 0:
        print("Number of displacements: %d" % num_disps)
        print("Number of supercell files: %d" % len(force_filenames))
    force_sets = get_force_sets(interface_mode,
                                num_atoms,
                                num_disps,
                                force_filenames,
                                disp_filename,
                                verbose=(log_level > 0))

    if args.forces_fcz:
        force_filename = args.forces_fcz
        file_exists(force_filename, log_level)
        force_set_zero = get_force_sets(interface_mode,
                                        num_atoms,
                                        1,
                                        [force_filename, ],
                                        verbose=(log_level > 0))[0]
        for fs in force_sets:
            fs -= force_set_zero

        if log_level > 0:
            print("Forces in \'%s\' were subtracted from supercell forces."
                  % force_filename)

    if force_sets:
        write_FORCES_FC2(disp_dataset,
                         forces_fc2=force_sets,
                         filename="FORCES_FC2")
        if log_level:
            print("")
            print("%s has been created." % "FORCES_FC2")
            print_end()
        sys.exit(0)
    else:
        if log_level:
            print("")
            print("%s could not be created." % "FORCES_FC2")
            print_error()
        sys.exit(1)

if args.force_sets_to_forces_fc2_mode:
    filename = 'FORCE_SETS'
    file_exists(filename, log_level)
    disp_dataset = parse_FORCE_SETS(filename=filename)
    write_FORCES_FC2(disp_dataset)
    write_disp_fc2_yaml(disp_dataset, None)

    if log_level:
        print("FORCES_FC2 and disp_fc2.yaml have been created from "
              "FORCE_SETS.")
        print_end()
    sys.exit(0)

#####################################
# Create FORCE_SETS from FORCES_FC* #
#####################################
if args.force_sets_mode:
    if args.phonon_supercell_dimension is not None:
        if input_filename is None:
            disp_filename = 'disp_fc2.yaml'
        else:
            disp_filename = 'disp_fc2.' + input_filename + '.yaml'
        forces_filename = "FORCES_FC2"
    else:
        if input_filename is None:
            disp_filename = 'disp_fc3.yaml'
        else:
            disp_filename = 'disp_fc3.' + input_filename + '.yaml'
        forces_filename = "FORCES_FC3"

    with open(forces_filename, 'r') as f:
        len_first_line = get_lenghth_of_first_line(f)

    if len_first_line == 3:
        file_exists(disp_filename, log_level)
        disp_dataset = parse_disp_fc2_yaml(filename=disp_filename)
        file_exists(forces_filename, log_level)
        parse_FORCES_FC2(disp_dataset, filename=forces_filename)

        if log_level:
            print("Displacement dataset is read from %s." % disp_filename)

        write_FORCE_SETS(disp_dataset)

        if log_level:
            print("FORCE_SETS has been created.")
            print_end()
    else:
        if log_level:
            print("The file format of %s is already readable by phonopy."
                  % forces_filename)
            print_end()

    sys.exit(0)

##################
# Parse settings #
##################
if len(args.conf_file) > 0:
    phono3py_conf = Phono3pyConfParser(filename=args.conf_file[0], args=args)
    settings = phono3py_conf.get_settings()
else:
    phono3py_conf = Phono3pyConfParser(args=args)
    settings = phono3py_conf.get_settings()

##################################################
# Set calculator interface and crystal structure #
##################################################
cell_info = collect_cell_info(
    supercell_matrix=settings.get_supercell_matrix(),
    primitive_matrix=settings.get_primitive_matrix(),
    interface_mode=get_interface_mode(args),
    cell_filename=settings.get_cell_filename(),
    chemical_symbols=settings.get_chemical_symbols(),
    command_name="phono3py",
    symprec=args.symprec)
if type(cell_info) is str:
    print_error_message(cell_info)
    if log_level > 0:
        print_error()
    sys.exit(1)
(unitcell, supercell_matrix, primitive_matrix,
 unitcell_filename, optional_structure_info, interface_mode,
 has_read_phonopy_yaml) = cell_info

if has_read_phonopy_yaml:
    ph3py = Phono3pyYaml()
    ph3py.read(unitcell_filename)
    phonon_supercell_matrix = ph3py.phonon_supercell_matrix
else:
    phonon_supercell_matrix = None

if unitcell is None:
    print_error_message("Crystal structure file of %s could not be found." %
                        unitcell_filename)
    if log_level > 0:
        print_error()
    sys.exit(1)

# Check unit cell
if np.linalg.det(unitcell.get_cell()) < 0.0:
    print_error_message("Determinant of the lattice vector matrix "
                        "has to be positive.")
    if log_level > 0:
        print_end()
    sys.exit(0)

######################
# Translate settings #
######################
conf = get_phono3py_configurations(settings)
if phonon_supercell_matrix is None:
    phonon_supercell_matrix = conf['phonon_supercell_matrix']
masses = conf['masses']
mesh = conf['mesh']
mesh_divs = conf['mesh_divs']
coarse_mesh_shifts = conf['coarse_mesh_shifts']
grid_points = conf['grid_points']
band_indices = conf['band_indices']
sigmas = conf['sigmas']
sigma_cutoff = conf['sigma_cutoff']
temperature_points = conf['temperature_points']
temperatures = conf['temperatures']
frequency_factor_to_THz = conf['frequency_factor_to_THz']
num_frequency_points = conf['num_frequency_points']
frequency_step = conf['frequency_step']
frequency_scale_factor = conf['frequency_scale_factor']
cutoff_frequency = conf['cutoff_frequency']

symprec = args.symprec

# Check supercell size
if supercell_matrix is None:
    print_error_message("Supercell dimension (--dim) has to be specified.")
    if log_level > 0:
        print_end()
    sys.exit(0)

#################################################
# Create supercells with displacements and exit #
#################################################
if settings.get_create_displacements():
    phono3py = create_phono3py_supercells(
        unitcell,
        supercell_matrix,
        phonon_supercell_matrix,
        settings.get_displacement_distance(),
        settings.get_is_plusminus_displacement(),
        settings.get_is_diagonal_displacement(),
        settings.get_cutoff_pair_distance(),
        write_supercells_with_displacements,
        optional_structure_info,
        settings.get_is_symmetry(),
        symprec,
        output_filename=output_filename,
        interface_mode=interface_mode,
        log_level=log_level)

    finalize_phono3py(log_level,
                      phono3py_conf,
                      phono3py,
                      interface_mode,
                      filename="phono3py_disp.yaml")

#################################################
# Change unit of lattice parameters to Angstrom #
#################################################
if distance_to_A is not None:
    lattice = unitcell.get_cell()
    lattice *= distance_to_A
    unitcell.set_cell(lattice)

#####################
# Initiate phono3py #
#####################
phono3py = Phono3py(
    unitcell,
    supercell_matrix,
    primitive_matrix=primitive_matrix,
    phonon_supercell_matrix=phonon_supercell_matrix,
    masses=masses,
    mesh=mesh,
    band_indices=band_indices,
    sigmas=sigmas,
    sigma_cutoff=sigma_cutoff,
    cutoff_frequency=cutoff_frequency,
    frequency_factor_to_THz=frequency_factor_to_THz,
    is_symmetry=settings.get_is_symmetry(),
    is_mesh_symmetry=settings.get_is_mesh_symmetry(),
    symmetrize_fc3q=settings.get_is_symmetrize_fc3_q(),
    symprec=symprec,
    calculator=interface_mode,
    log_level=log_level,
    lapack_zheev_uplo=args.uplo)

supercell = phono3py.supercell
primitive = phono3py.primitive
phonon_supercell = phono3py.phonon_supercell
phonon_primitive = phono3py.phonon_primitive
symmetry = phono3py.symmetry
mesh_numbers = phono3py.mesh_numbers

run_mode = get_run_mode(settings)
if log_level:
    show_general_settings(settings,
                          run_mode,
                          (type(primitive_matrix) is str and
                           primitive_matrix == 'auto'),
                          phono3py.primitive_matrix,
                          supercell_matrix,
                          phonon_supercell_matrix,
                          unitcell_filename,
                          input_filename,
                          output_filename,
                          compression,
                          interface_mode)

if log_level > 1:
    show_phono3py_cells(symmetry,
                        primitive,
                        supercell,
                        phonon_primitive,
                        phonon_supercell,
                        settings)
else:
    print("Spacegroup: %s" %
          phono3py.symmetry.get_international_table())
    print("Use -v option to watch primitive cell, unit cell, "
          "and supercell structures.")

#####################################################
# Write ir-grid points and grid addresses, and exit #
#####################################################
if args.write_grid_points:
    write_grid_points(primitive,
                      mesh_numbers,
                      mesh_divs=mesh_divs,
                      band_indices=band_indices,
                      sigmas=sigmas,
                      temperatures=temperatures,
                      coarse_mesh_shifts=coarse_mesh_shifts,
                      is_kappa_star=settings.get_is_kappa_star(),
                      is_lbte=(settings.get_write_collision() or
                               settings.get_is_lbte()),
                      compression=compression,
                      symprec=symprec)

    if log_level:
        print_end()
    sys.exit(0)

##################################################
# Show reduced number of triplets at grid points #
##################################################
if args.show_num_triplets:
    show_num_triplets(primitive,
                      mesh_numbers,
                      mesh_divs=mesh_divs,
                      band_indices=band_indices,
                      grid_points=grid_points,
                      coarse_mesh_shifts=coarse_mesh_shifts,
                      is_kappa_star=settings.get_is_kappa_star(),
                      symprec=symprec)

    if log_level:
        print_end()
    sys.exit(0)

##################################
# Non-analytical term correction #
##################################
if settings.get_is_nac() or has_read_phonopy_yaml:
    def read_nac_params_from_phonopy_yaml(unitcell_filename, log_level):
        ph3_yaml = Phono3pyYaml()
        ph3_yaml.read(unitcell_filename)
        if ph3_yaml.nac_params is not None:
            if log_level:
                print("NAC parameters are read from \"%s\"."
                      % unitcell_filename)
            return ph3_yaml.nac_params
        else:
            return None

    def read_BORN(phono3py, log_level):
        with open("BORN") as f:
            nac_params = get_born_parameters(
                f, phono3py.phonon_primitive, phono3py.primitive_symmetry)
            if log_level:
                print("NAC parameters are read from \"%s\"." % "BORN")
            return nac_params

    if has_read_phonopy_yaml:
        nac_params = read_nac_params_from_phonopy_yaml(unitcell_filename,
                                                       log_level)
        if nac_params is None and file_exists("BORN", log_level, is_any=True):
            nac_params = read_BORN(phono3py, log_level)
    elif file_exists("BORN", log_level):
        nac_params = read_BORN(phono3py, log_level)
        if not nac_params:
            error_text = "BORN file could not be read correctly."
            print_error_message(error_text)
            if log_level > 0:
                print_error()
            sys.exit(1)
    else:
        nac_params = None  # This would not happen though.

    nac_q_direction = settings.get_nac_q_direction()

    if nac_params is not None:
        if nac_params['factor'] is None:
            nac_params['factor'] = Hartree * Bohr
        if settings.get_nac_method() is not None:
            nac_params['method'] = settings.get_nac_method()

        if log_level > 1:
            print("-" * 27 + " Dielectric constant " + "-" * 28)
            for v in nac_params['dielectric']:
                print("         %12.7f %12.7f %12.7f" % tuple(v))
            print("-" * 26 + " Born effective charges " + "-" * 26)
            symbols = primitive.get_chemical_symbols()
            for i, (z, s) in enumerate(zip(nac_params['born'], symbols)):
                for j, v in enumerate(z):
                    if j == 0:
                        text = "%5d %-2s" % (i + 1, s)
                    else:
                        text = "        "
                    print("%s %12.7f %12.7f %12.7f" % ((text,) + tuple(v)))
else:
    nac_params = None
    nac_q_direction = None

# if settings.get_is_nac():
#     file_exists('BORN', log_level)
#     nac_params = parse_BORN(phonon_primitive)
#     nac_factor = Hartree * Bohr
#     if nac_params['factor'] is not None:
#         if log_level:
#             print("-" * 22 + " Non-analytical term correction " + "-" * 22)
#             print("Default NAC unit conversion factor is %6.4f." % nac_factor)
#             print("But instead the value in BORN file %f is used." %
#                   nac_params['factor'])
#     else:
#         nac_params['factor'] = nac_factor
#     if settings.get_nac_method() is not None:
#         nac_params['method'] = settings.get_nac_method()

#     if log_level > 1:
#         print("-" * 27 + " Dielectric constant " + "-" * 28)
#         for v in nac_params['dielectric']:
#             print("         %12.7f %12.7f %12.7f" % tuple(v))
#         print("-" * 26 + " Born effective charges " + "-" * 26)
#         symbols = phonon_primitive.get_chemical_symbols()
#         for i, (z, s) in enumerate(zip(nac_params['born'], symbols)):
#             for j, v in enumerate(z):
#                 if j == 0:
#                     text = "%5d %-2s" % (i + 1, s)
#                 else:
#                     text = "        "
#                 print("%s %12.7f %12.7f %12.7f" % ((text,) + tuple(v)))

#     nac_q_direction = settings.get_nac_q_direction()
# else:
#     nac_params = None
#     nac_q_direction = None

###################
# Force constants #
###################
create_phono3py_force_constants(phono3py,
                                phonon_supercell_matrix,
                                settings,
                                force_to_eVperA=force_to_eVperA,
                                distance_to_A=distance_to_A,
                                compression=compression,
                                input_filename=input_filename,
                                output_filename=output_filename,
                                log_level=log_level)

##############################
# Phonon Gruneisen parameter #
##############################
if settings.get_is_gruneisen():
    if (mesh_numbers is None and
        settings.get_band_paths() is None and
        settings.get_qpoints() is None):
        print("An option of --mesh, --band, or --qpoints has to be specified.")
        if log_level:
            print_error()
        sys.exit(1)

    if len(phono3py.get_fc2()) != len(phono3py.get_fc3()):
        print("Supercells used for fc2 and fc3 have to be same.")
        if log_level:
            print_error()
        sys.exit(1)

    if settings.get_band_paths() is not None:
        if settings.get_band_points() is None:
            npoints = 51
        else:
            npoints = settings.get_band_points()
        band_paths = settings.get_band_paths()
        bands = get_band_qpoints(band_paths, npoints=npoints)
    else:
        bands = None

    rotations = phono3py.primitive_symmetry.get_pointgroup_operations()
    gruneisen = get_gruneisen_parameters(
        phono3py.get_fc2(),
        phono3py.get_fc3(),
        supercell,
        primitive,
        bands,
        mesh_numbers,
        rotations,
        settings.get_qpoints(),
        nac_params=nac_params,
        nac_q_direction=nac_q_direction,
        ion_clamped=settings.get_ion_clamped(),
        factor=VaspToTHz,
        symprec=symprec,
        output_filename=output_filename,
        log_level=log_level)

    if log_level:
        print_end()
    sys.exit(0)

##################
# Check settings #
##################
if mesh_numbers is None and run_mode is not None:
    print("")
    print("Q-point mesh has to be specified with --mesh option.")
    print("")
    if log_level:
        print_error()
    sys.exit(1)

if run_mode in ("imag_self_energy", "jdos", "isotope", "frequency_shift"):
    if grid_points is None:
        print("")
        print("Grid point(s) has to be specified with --gp or --ga option.")
        print("")
        if log_level:
            print_error()
        sys.exit(1)

#################
# Show settings #
#################
if log_level and run_mode is not None:
    show_phono3py_settings(settings,
                           mesh_numbers,
                           mesh_divs,
                           band_indices,
                           sigmas,
                           sigma_cutoff,
                           temperatures,
                           temperature_points,
                           grid_points,
                           cutoff_frequency,
                           frequency_factor_to_THz,
                           frequency_scale_factor,
                           frequency_step,
                           num_frequency_points,
                           nac_params,
                           log_level)

#############
# Joint DOS #
#############
if run_mode == "jdos":
    joint_dos = Phono3pyJointDos(
        phonon_supercell,
        phonon_primitive,
        mesh_numbers,
        phono3py.get_fc2(),
        nac_params=nac_params,
        nac_q_direction=nac_q_direction,
        sigmas=sigmas,
        cutoff_frequency=cutoff_frequency,
        frequency_step=frequency_step,
        num_frequency_points=num_frequency_points,
        temperatures=temperature_points,
        frequency_factor_to_THz=frequency_factor_to_THz,
        frequency_scale_factor=frequency_scale_factor,
        is_mesh_symmetry=settings.get_is_mesh_symmetry(),
        symprec=symprec,
        output_filename=output_filename,
        log_level=log_level)
    joint_dos.run(grid_points)
    if log_level:
        print_end()
    sys.exit(0)

#############################
# Phonon-isotope scattering #
#############################
if settings.get_is_isotope() and settings.get_mass_variances() is None:
    from phonopy.structure.atoms import isotope_data
    symbols = phonon_primitive.get_chemical_symbols()
    in_database = True
    for s in set(symbols):
        if s not in isotope_data:
            print("%s is not in the list of isotope databese" % s)
            print("(not implemented).")
            print("Use --mass_variances option.")
            in_database = False
    if not in_database:
        if log_level:
            print_end()
        sys.exit(0)

###########################
# Phonon-isotope lifetime #
###########################
if run_mode == "isotope":
    mass_variances = settings.get_mass_variances()
    if band_indices is not None:
        band_indices = np.hstack(band_indices).astype('intc')
    iso = Phono3pyIsotope(
        mesh_numbers,
        phonon_primitive,
        mass_variances=mass_variances,
        band_indices=band_indices,
        sigmas=sigmas,
        frequency_factor_to_THz=frequency_factor_to_THz,
        symprec=symprec,
        cutoff_frequency=settings.get_cutoff_frequency(),
        lapack_zheev_uplo=args.uplo)
    iso.set_dynamical_matrix(phono3py.get_fc2(),
                             phonon_supercell,
                             phonon_primitive,
                             nac_params=nac_params,
                             frequency_scale_factor=frequency_scale_factor)
    iso.run(grid_points)
    if log_level:
        print_end()
    sys.exit(0)

#############################
# Phonon-phonon interaction #
#############################
if run_mode is not None:
    ave_pp = settings.get_constant_averaged_pp_interaction()
    phono3py.set_phph_interaction(
        nac_params=nac_params,
        nac_q_direction=nac_q_direction,
        constant_averaged_interaction=ave_pp,
        frequency_scale_factor=frequency_scale_factor,
        unit_conversion=settings.get_pp_conversion_factor(),
        solve_dynamical_matrices=(not settings.get_read_phonon()))

    if settings.get_write_phonon():
        freqs, eigvecs, grid_address = phono3py.get_phonon_data()
        filename = write_phonon_to_hdf5(freqs,
                                        eigvecs,
                                        grid_address,
                                        mesh_numbers,
                                        compression=compression,
                                        filename=output_filename)
        if filename:
            if log_level:
                print("Phonons are written into \"%s\"." % filename)
                print_end()
            sys.exit(0)
        else:
            print("Writing phonons failed.")
            if log_level:
                print_error()
            sys.exit(1)

    if settings.get_read_phonon():
        phonons = read_phonon_from_hdf5(mesh_numbers,
                                        filename=input_filename,
                                        verbose=(log_level > 0))
        if phonons is None:
            print("Reading phonons failed.")
            if log_level:
                print_error()
            sys.exit(1)

        frequencies = phonons[0]
        eigenvectors = phonons[1]
        grid_address = phonons[2]
        if phono3py.set_phonon_data(frequencies, eigenvectors, grid_address):
            pass
        else:
            if log_level:
                print_error()
            sys.exit(1)

if run_mode == "imag_self_energy":
    phono3py.run_imag_self_energy(
        grid_points,
        frequency_step=frequency_step,
        num_frequency_points=num_frequency_points,
        temperatures=temperature_points,
        scattering_event_class=settings.get_scattering_event_class(),
        write_gamma_detail=settings.get_write_gamma_detail(),
        output_filename=output_filename)
    phono3py.write_imag_self_energy(filename=output_filename)
elif run_mode == "frequency_shift":
    phono3py.get_frequency_shift(
        grid_points,
        temperatures=temperatures,
        output_filename=output_filename)
elif run_mode == "conductivity-RTA" or run_mode == "conductivity-LBTE":
    phono3py.run_thermal_conductivity(
        is_LBTE=settings.get_is_lbte(),
        temperatures=temperatures,
        is_isotope=settings.get_is_isotope(),
        mass_variances=settings.get_mass_variances(),
        grid_points=grid_points,
        boundary_mfp=settings.get_boundary_mfp(),
        solve_collective_phonon=settings.get_solve_collective_phonon(),
        use_ave_pp=settings.get_use_ave_pp(),
        gamma_unit_conversion=settings.get_gamma_conversion_factor(),
        mesh_divisors=mesh_divs,
        coarse_mesh_shifts=settings.get_coarse_mesh_shifts(),
        is_reducible_collision_matrix=settings.get_is_reducible_collision_matrix(),
        is_kappa_star=settings.get_is_kappa_star(),
        gv_delta_q=settings.get_group_velocity_delta_q(),
        is_full_pp=settings.get_is_full_pp(),
        pinv_cutoff=settings.get_pinv_cutoff(),
        pinv_solver=settings.get_pinv_solver(),
        write_gamma=settings.get_write_gamma(),
        read_gamma=settings.get_read_gamma(),
        write_kappa=True,
        is_N_U=settings.get_is_N_U(),
        write_gamma_detail=settings.get_write_gamma_detail(),
        write_collision=settings.get_write_collision(),
        read_collision=settings.get_read_collision(),
        write_pp=settings.get_write_pp(),
        read_pp=settings.get_read_pp(),
        write_LBTE_solution=settings.get_write_LBTE_solution(),
        compression=compression,
        input_filename=input_filename,
        output_filename=output_filename)
else:
    if log_level:
        print("*" * 15 + " None of ph-ph interaction was calculated. " +
              "*" * 16)

finalize_phono3py(log_level, phono3py_conf, phono3py, interface_mode)
