phonopy/scripts/phono3py

#!/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 os
import sys
import numpy as np
from optparse import OptionParser

from phonopy.interface.vasp import read_vasp
from phonopy.structure.cells import get_supercell, Primitive, print_cell
from phonopy.structure.symmetry import Symmetry
from phonopy.harmonic.force_constants import get_force_constants, \
    set_permutation_symmetry, symmetrize_force_constants, \
    set_translational_invariance
from phonopy.hphonopy.file_IO import parse_FORCE_SETS, parse_BORN
from phonopy.units import VaspToTHz
from anharmonic.phonon3.fc3 import get_fc3, symmetrize_fc3, \
    set_translational_invariance_fc3
from anharmonic.file_IO import write_fc2_dat, write_fc3_dat,\
    parse_disp_yaml, write_FORCES_THIRD,\
    parse_DELTA_FORCES, write_supercells_with_displacements,\
    write_DELTA_FC2_SETS, parse_DELTA_FC2_SETS,\
    write_fc3_to_hdf5, write_fc2_to_hdf5,\
    read_fc3_from_hdf5, read_fc2_from_hdf5, write_ir_grid_points
from anharmonic.phonon3.displacement_fc3 import get_third_order_displacements
from anharmonic.settings import Phono3pyConfParser
from anharmonic.phonon3 import Phono3py, JointDOS, get_gruneisen_parameters    

# AA is created at http://www.network-science.de/ascii/.
def print_phono3py():
    print """        _                      _____             
  _ __ | |__   ___  _ __   ___|___ / _ __  _   _ 
 | '_ \| '_ \ / _ \| '_ \ / _ \ |_ \| '_ \| | | |
 | |_) | | | | (_) | | | | (_) |__) | |_) | |_| |
 | .__/|_| |_|\___/|_| |_|\___/____/| .__/ \__, |
 |_|                                |_|    |___/ 
"""

def print_end():
    print """                 _ 
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|
"""

def print_error(message):
    print message

# Parse options
parser = OptionParser()
parser.set_defaults(amplitude=None,
                    band_indices=None,
                    band_paths=None,
                    band_points=None,
                    cell_poscar=None,
                    delta_fc2=False,
                    factor=None,
                    delta_fc2_sets_mode=False,
                    freq_scale=None,
                    gamma_option=0,
                    grid_points=None,
                    gv_delta_q=None,
                    ion_clamped=False,
                    is_decay_channel=False,
                    is_nodiag=False,
                    is_displacement=False,
                    is_nosym=False,
                    is_gruneisen=False,
                    is_joint_dos=False,
                    is_linewidth=False,
                    is_bterta=False,
                    is_nac=False,
                    is_plusminus_displacements=False,
                    is_translational_symmetry=False,
                    is_symmetrize_fc2=False,
                    is_symmetrize_fc3_r=False,
                    is_symmetrize_fc3_q=False,                     
                    log_level=None,
                    max_freepath=None,
                    mesh_numbers=None,
                    mesh_divisors=None,
                    multiple_sigmas=None,
                    no_kappa_stars=False,
                    q_direction=None,
                    primitive_axis=None,
                    read_amplitude=False,
                    read_fc2=False,
                    read_fc2_extra=False,
                    read_fc3=False,
                    read_gamma=False,
                    freq_step=None,
                    output_filename=None,
                    qpoints=None,
                    quiet=False,
                    sigma=None,
                    supercell_dimension=None,
                    supercell_dimension_extra=None,
                    symprec=1e-5,
                    tmax=None,
                    tmin=None,
                    tstep=None,
                    temperatures=None,
                    verbose=False,
                    uplo='L',
                    write_amplitude=False,
                    write_gamma=False,
                    write_grid_points=False)
parser.add_option("--amplitude", dest="amplitude", type="float",
                  help="Distance of displacements")
parser.add_option("--bi", "--band_indices", dest="band_indices",
                  type="string",
                  help="Band indices where life time is calculated")
parser.add_option("--band", dest="band_paths",
                  action="store", type="string",
                  help="Band structure paths calculated for Gruneisen parameter")
parser.add_option("--band_points", dest="band_points",
                  type="int",
                  help="Number of points calculated on a band segment in the band structure Gruneisen parameter calculation")
parser.add_option("--br", "--bterta", dest="is_bterta",
                  action="store_true",
                  help="Calculate thermal conductivity in BTE-RTA")
parser.add_option("-c", "--cell", dest="cell_poscar",
                  action="store", type="string",
                  help="Read unit cell", metavar="FILE")
parser.add_option("-d", "--disp", dest="is_displacement",
                  action="store_true",
                  help="As first stage, get least displacements")
parser.add_option("--decay", dest="is_decay_channel",
                  action="store_true", help="Calculate decay channels")
parser.add_option("--dim",
                  dest="supercell_dimension",
                  type="string",
                  help="Supercell dimension")
parser.add_option("--dim2",
                  dest="supercell_dimension_extra",
                  type="string",
                  help="Supercell dimension for extra fc2")
parser.add_option("--cf3", "--create_f3",
                  dest="forces_third_mode",
                  action="store_true",
                  help="Create FORCES_THIRD")
parser.add_option("--cdfc2", "--create_delta_fc2",
                  dest="delta_fc2_sets_mode",
                  action="store_true",
                  help="Create DELTA_FC2_SETS")
parser.add_option("--dfc2", "--delta_fc2",
                  dest="read_delta_fc2",
                  action="store_true",
                  help="Read DELTA_FC2_SETS")
parser.add_option("--factor", dest="factor", type="float",
                  help="Conversion factor to favorite frequency unit")
parser.add_option("--fc2",
                  dest="read_fc2",
                  action="store_true",
                  help="Read second order force constants")
parser.add_option("--fc2_extra",
                  dest="read_fc2_extra",
                  action="store_true",
                  help="Read extra second order force constants")
parser.add_option("--fc3",
                  dest="read_fc3",
                  action="store_true",
                  help="Read third order force constants")
parser.add_option("--freepath", dest="max_freepath", type="float",
                  help="Maximum mean free path of phonon in meter")
parser.add_option("--freq_scale", dest="freq_scale", type="float",
                  help="Scale factor is multiplied to frequencies only, i.e., changes frequencies but assumed not to change the physical unit")
parser.add_option("--freq_step", dest="freq_step", type="float",
                  help="Pitch of frequency calculated")
parser.add_option("--gamma_option", dest="gamma_option", type="int",
                  help="Option for the test of imaginary part of self energy")
parser.add_option("--gp", "--grid_points",
                  dest="grid_points",
                  type="string",
                  help="Fixed grid points where damping functions are calculated ")
parser.add_option("--gruneisen", dest="is_gruneisen",
                  action="store_true",
                  help="Calculate phonon Gruneisen parameter")
parser.add_option("--gv_delta_q", dest="gv_delta_q", type="float",
                  help="Delta-q distance used for group velocity calculation")
parser.add_option("--ion_clamped", dest="ion_clamped",
                  action="store_true",
                  help="Atoms are clamped under applied strain in Gruneisen parameter calculation")
parser.add_option("--jdos",
                  dest="is_joint_dos",
                  action="store_true",
                  help="Calculate joint density of states")
parser.add_option("--lw", "--linewidth",
                  dest="is_linewidth",
                  action="store_true",
                  help="Calculate linewidths")
parser.add_option("--md", "--mesh_divisors",
                  dest="mesh_divisors",
                  type="string",
                  help="Divisors for mesh numbers")
parser.add_option("--mesh",
                  dest="mesh_numbers",
                  type="string",
                  help="Mesh numbers")
parser.add_option("--multiple_sigmas", dest="multiple_sigmas",
                  type="string",
                  help="Multiple sigmas for smearing width used for limited functions")
parser.add_option("--nac", dest="is_nac",
                  action="store_true",
                  help="Non-analytical term correction")
parser.add_option("--nodiag", dest="is_nodiag",
                  action="store_true",
                  help="Set displacements parallel to axes")
parser.add_option("--noks", "--no_kappa_stars", dest="no_kappa_stars",
                  action="store_true",
                  help="Deactivate summation of partial kappa at q-stars"),
parser.add_option("--nosym", dest="is_nosym",
                  action="store_true",
                  help="No symmetrization of triplets")
parser.add_option("-o", dest="output_filename",
                  type="string",
                  help="Filename of output of damping function")
parser.add_option("--pa", "--primitive_axis", dest="primitive_axis",
                  action="store", type="string",
                  help="Same as PRIMITIVE_AXIS tags")
parser.add_option("--pm", dest="is_plusminus_displacements",
                  action="store_true",
                  help="Set plus minus displacements")
parser.add_option("--qpoints", dest="qpoints", type="string",
                  help="Calculate at specified q-points")
parser.add_option("--q_direction",
                  dest="q_direction",
                  type="string",
                  help="q-vector direction at q->0 for non-analytical term correction")
parser.add_option("-q", "--quiet", dest="quiet",
                  action="store_true",
                  help="Print out smallest information")
parser.add_option("--read_amplitude", dest="read_amplitude",
                  action="store_true",
                  help="Read phonon-phonon interaction amplitudes")
parser.add_option("--read_gamma", dest="read_gamma",
                  action="store_true",
                  help="Read Gammas from files")
parser.add_option("--sigma", dest="sigma", type="float",
                  help="Smearing width for DOS")
parser.add_option("--sym_fc2", dest="is_symmetrize_fc2",
                  action="store_true",
                  help="Symmetrize fc2 by index exchange")
parser.add_option("--sym_fc3r", dest="is_symmetrize_fc3_r",
                  action="store_true",
                  help="Symmetrize fc3 in real space by index exchange")
parser.add_option("--sym_fc3q", dest="is_symmetrize_fc3_q",
                  action="store_true",
                  help="Symmetrize fc3 in reciprocal space by index exchange")
parser.add_option("--tmax", dest="tmax", type="string",
                  help="Maximum calculated temperature")
parser.add_option("--tmin", dest="tmin", type="string",
                  help="Minimum calculated temperature")
parser.add_option("--tstep", dest="tstep", type="string",
                  help="Calculated temperature step")
parser.add_option("--tsym", dest="is_translational_symmetry",
                  action="store_true",
                  help="Impose translational invariance condition")
parser.add_option("--tolerance", dest="symprec", type="float",
                  help="Symmetry tolerance to search")
parser.add_option("-v", "--verbose", dest="verbose",
                  action="store_true",
                  help="Detailed run-time information is displayed")
parser.add_option("--loglevel", dest="log_level", type="int",
                  help="Log level")
parser.add_option("--ts", dest="temperatures", 
                  type="string", help="Temperatures for damping functions")
parser.add_option("--uplo",
                  dest="uplo",
                  type="string",
                  help="Lapack zheev UPLO")
parser.add_option("--wgp", "--write_grid_points", dest="write_grid_points",
                  action="store_true",
                  help="Write grid address of irreducible grid points for specified mesh numbers to ir_grid_address.yaml")
parser.add_option("--write_amplitude", dest="write_amplitude",
                  action="store_true",
                  help="Write phonon-phonon interaction amplitudes")
parser.add_option("--write_gamma", dest="write_gamma",
                  action="store_true",
                  help="Write gamma")
(options, args) = parser.parse_args()
option_list = parser.option_list

# Log level
log_level = 1
if options.verbose:
    log_level = 2
if options.quiet:
    log_level = 0
if not options.log_level==None:
    log_level=options.log_level
    
# Title
if log_level:
    print_phono3py()

# Create FORCES_THIRD
if options.forces_third_mode:
    displacements = parse_disp_yaml('disp.yaml')
    write_FORCES_THIRD(args, displacements)
    print_end()
    exit(0)

# Create DELTA_FC2_SETS
if options.delta_fc2_sets_mode:
    displacements = parse_disp_yaml('disp.yaml')
    write_DELTA_FC2_SETS(args, displacements)
    print_end()
    exit(0)

# Import input files
if len(args) > 0:
    phono3py_conf = Phono3pyConfParser(filename=args[0],
                                       options=options,
                                       option_list=option_list)
    settings = phono3py_conf.get_settings()

else:
    phono3py_conf = Phono3pyConfParser(options=options,
                                       option_list=option_list)
    settings = phono3py_conf.get_settings()

# Read POSCAR
if options.cell_poscar == None:
    if os.path.exists('POSCAR'):
        unitcell_filename = 'POSCAR'
    else:
        print_error("POSCAR could not be found.")
        if log_level:
            print_end()
        sys.exit(1)
else:
    if os.path.exists(options.cell_poscar):
        unitcell_filename = options.cell_poscar
    else:
        print_error("The file \'%s\' could not be found." %
                    options.cell_poscar)
        if log_level:
            print_end()
        sys.exit(1)
    
unitcell = read_vasp(unitcell_filename,
                     settings.get_chemical_symbols())

# Supercell and Symmetry 
supercell = get_supercell(unitcell, settings.get_supercell_matrix())
symmetry = Symmetry(supercell, options.symprec)

if settings.get_supercell_matrix_extra() is not None:
    supercell_extra = get_supercell(unitcell, settings.get_supercell_matrix_extra())
    symmetry_extra = Symmetry(supercell_extra, options.symprec)

# Log
if log_level:
    if options.is_translational_symmetry:
        print "Translational symmetry:", options.is_translational_symmetry
    if options.is_symmetrize_fc2:
        print "FC2 symmetry of index exchange:", options.is_symmetrize_fc2
    if options.is_symmetrize_fc3_r:
        print "FC3 symmetry of index exchange in real space:", options.is_symmetrize_fc3_r
    if options.is_symmetrize_fc3_q:
        print "FC3 symmetry of index exchange in reciprocal space:", options.is_symmetrize_fc3_q
    if not options.supercell_dimension_extra == None:
        print "Extra supercell for fc2 is supplied."
    if settings.get_is_nac():
        print "Non-analytical term correction:", settings.get_is_nac()
    print "Spacegroup: ", symmetry.get_international_table()

###############################################################
# Create supercells with displacements and exit (pre-process) #
###############################################################
if options.is_displacement:
    dds = get_third_order_displacements(
        supercell,
        symmetry,
        is_plusminus=settings.get_is_plusminus_displacement(),
        is_diagonal=settings.get_is_diagonal_displacement())
    write_supercells_with_displacements(supercell,
                                        dds,
                                        options.amplitude)

#########################################
# Calculate third-order force constants #
#########################################
else:
    primitive = Primitive(
        supercell,
        np.dot(np.linalg.inv(settings.get_supercell_matrix()),
               settings.get_primitive_matrix()),
        options.symprec)

    if not settings.get_supercell_matrix_extra()==None:
        primitive_extra = Primitive(
            supercell_extra,
            np.dot(np.linalg.inv(settings.get_supercell_matrix_extra()),
                   settings.get_primitive_matrix()),
            options.symprec)



    if log_level:
        print "------------------------ primitive cell for fc ---------------------------"
        print_cell(primitive)
        print "-------------------------- supercell for fc ------------------------------"
        print_cell(supercell, mapping=primitive.get_supercell_to_primitive_map())
        print "----------------- ratio (supercell for fc)/(primitive) -------------------"
        for vec in np.dot(supercell.get_cell(), np.linalg.inv(primitive.get_cell())):
            print "%5.2f"*3 % tuple(vec)
    if log_level and (not settings.get_supercell_matrix_extra()==None):
        print "------------------------- primitive cell extra ----------------------------"
        print_cell(primitive_extra)
        print "--------------------------- supercell extra -------------------------------"
        print_cell(supercell_extra, mapping=primitive_extra.get_supercell_to_primitive_map())
        print "--------------- ratio (supercell extra)/(primitive  extra) ----------------"
        for vec in np.dot(supercell_extra.get_cell(),
                           np.linalg.inv(primitive_extra.get_cell())):
            print "%5.2f"*3 % tuple(vec)

    # Write ir-grid points
    if options.write_grid_points:
        print "---------------------------------------------------------------------------"
        mesh = settings.get_mesh_numbers()
        if mesh is None:
            print "To write grid points, mesh numbers have to be specified."
        else:
            print "Ir-grid points are written into ir_grid_points.yaml."
            mesh_divs = settings.get_mesh_divisors()
            coarse_mesh_shifts = settings.get_coarse_mesh_shifts()
            write_ir_grid_points(primitive, mesh, mesh_divs, coarse_mesh_shifts)
        sys.exit(0)
            
    # fc2
    if options.read_fc2 or options.read_delta_fc2:
        if log_level:
            print  "----- Read fc2 -----"
            sys.stdout.flush()
        if os.path.exists('fc2.hdf5'):
            fc2_with_dim = read_fc2_from_hdf5()
        elif os.path.exists('force_constants.hdf5'):
            fc2_with_dim = read_fc2_from_hdf5(filename='force_constants.hdf5')
        else:
            print "fc2.hdf5 not found"
            if log_level:
                print_end()
            sys.exit(0)
            
    else:
        if log_level:
            print "----- Solve fc2 -----"
            sys.stdout.flush()

        forces_second = parse_FORCE_SETS(supercell.get_number_of_atoms())
        fc2_with_dim = get_force_constants(forces_second,
                                           symmetry,
                                           supercell)

        write_fc2_dat(fc2_with_dim)
        write_fc2_to_hdf5(fc2_with_dim)

    if settings.get_supercell_matrix_extra()==None:
        fc2 = fc2_with_dim
    else:
        # fc2 extra (FORCE_SETS_EXTRA)
        if options.read_fc2_extra:
            if log_level:
                print "----- Read fc2 extra -----"
                sys.stdout.flush()

            fc2 = read_fc2_from_hdf5(filename='fc2_extra.hdf5')
        else:                
            if log_level:
                print "----- Solve fc2 extra -----"
                sys.stdout.flush()
    
            forces_second_extra = parse_FORCE_SETS(
                supercell_extra.get_number_of_atoms(),
                filename="FORCE_SETS_EXTRA")
            fc2 = get_force_constants(forces_second_extra,
                                      symmetry_extra,
                                      supercell_extra)
            write_fc2_dat(fc2, filename='fc2_extra.dat')
            write_fc2_to_hdf5(fc2, filename='fc2_extra.hdf5')

    if options.is_symmetrize_fc2:
        set_permutation_symmetry(fc2)
        # symmetrize_force_constants(fc2)

    if options.is_translational_symmetry:
        set_translational_invariance(fc2)

    if settings.get_is_nac():
        if os.path.exists('BORN'):
            if settings.get_supercell_matrix_extra()==None:
                nac_params = parse_BORN(primitive)
            else:
                nac_params = parse_BORN(primitive_extra)

            nac_q_direction = settings.get_q_direction()
        else:
            print_error("BORN not found")
            if log_level:
                print_end()
            sys.exit(1)
    else:
        nac_params = None
        nac_q_direction = None

    # fc3
    if (options.is_joint_dos or
        settings.get_read_gamma() or
        settings.get_read_amplitude()):
        fc3 = None
    else:
        if options.read_fc3: # Read FORCES_THIRD
            if log_level:
                print  "----- Read fc3 -----"
                sys.stdout.flush()
            fc3 = read_fc3_from_hdf5()
            if options.is_translational_symmetry:
                set_translational_invariance_fc3(fc3)
        else: # fc3 from FORCES_THIRD and FORCES_SECOND
            if log_level:
                print  "----- Solve fc3 -----"
                sys.stdout.flush()

            displacements = parse_disp_yaml('disp.yaml')
            if options.read_delta_fc2:
                parse_DELTA_FC2_SETS(displacements)
            else:
                parse_DELTA_FORCES(displacements)
            fc3 = get_fc3(supercell,
                          displacements,
                          fc2_with_dim,
                          symmetry,
                          options.is_translational_symmetry,
                          verbose=log_level)

        num_atom = fc3.shape[0]
        maxval1 = 0
        maxval2 = 0
        maxval3 = 0
        for i, j, k, l, m in list(np.ndindex((num_atom, num_atom, 3, 3, 3))):
            val1 = fc3[:, i, j, k, l, m].sum()
            val2 = fc3[i, :, j, k, l, m].sum()
            val3 = fc3[i, j, :, k, l, m].sum()
            if abs(val1) > abs(maxval1):
                maxval1 = val1
            if abs(val2) > abs(maxval2):
                maxval2 = val2
            if abs(val3) > abs(maxval3):
                maxval3 = val3
        print "max drift:", maxval1, maxval2, maxval3 
        

        # Symmetrize fc3_r
        if options.is_symmetrize_fc3_r:
            if log_level:
                print "----- Symmetrize fc3 real space -----"
            symmetrize_fc3(fc3)
            
        # Write fc3
        if not options.read_fc3:
            if log_level:
                print "----- Write fc3.dat -----"
            write_fc3_dat(fc3)
            write_fc3_to_hdf5(fc3)
    
    #============================
    # Phonon Gruneisen parameter
    #============================
    if options.is_gruneisen:
        mesh = settings.get_mesh_numbers()
        band_paths = settings.get_bands()
        qpoints = settings.get_qpoints()
        ion_clamped = settings.get_ion_clamped()

        if (mesh is None and
            band_paths is None and
            qpoints is None):
            if log_level:
                print "An option of --mesh, --band, or --qpoints has to be specified."
                sys.exit(1)

        if log_level:
            print "------ Phonon Gruneisen parameter ------"
            if mesh is not None:
                print "Mesh:", mesh
            elif band_paths is not None:
                print "Paths in reciprocal reduced coordinates:"
                for path in band_paths:
                    print ("[%5.2f %5.2f %5.2f] --> [%5.2f %5.2f %5.2f]" % 
                           (tuple(path[0]) + tuple(path[-1])))
            if ion_clamped:
                print "To be calculated with ion clamped."
                
            sys.stdout.flush()

        gr = get_gruneisen_parameters(fc2,
                                      fc3,
                                      supercell,
                                      primitive,
                                      nac_params=nac_params,
                                      nac_q_direction=nac_q_direction,
                                      ion_clamped=ion_clamped,
                                      factor=VaspToTHz,
                                      symprec=options.symprec)
        if mesh is not None:
            gr.set_sampling_mesh(mesh,
                                 is_gamma_center=True)
        elif band_paths is not None:
            gr.set_band_structure(band_paths)
        elif qpoints is not None:
            gr.set_qpoints(qpoints)
        gr.run()
        gr.write_yaml()

    #===========================
    # phonon-phonon interaction
    #===========================
    elif settings.get_mesh_numbers() is not None:
        mesh = settings.get_mesh_numbers()
        if options.grid_points is None:
            grid_points = None
        else:
            grid_points = np.array(
                [int(x) for x in options.grid_points.replace(',', ' ').split()])
            
        if options.factor is None:
            factor = VaspToTHz
        else:
            factor = options.factor

        freq_factor = 1.0

        if settings.get_sigma() is None:
            sigma = 0.2 * freq_factor
        else:
            sigma = settings.get_sigma()

        if settings.get_multiple_sigmas() is None:
            multiple_sigmas = [sigma]
        else:
            multiple_sigmas = settings.get_multiple_sigmas()

        if settings.get_omega_step() is None:
            freq_step = 0.1 * freq_factor
        else:
            freq_step = settings.get_omega_step()

        if options.freq_scale is None:
            freq_scale = 1.0
        else:
            freq_scale = options.freq_scale

        if settings.get_supercell_matrix_extra() is None:
            supercell_dm = supercell
            primitive_dm = primitive
        else:
            supercell_dm = supercell_extra
            primitive_dm = primitive_extra

        if settings.get_temperatures() is None:
            temperatures = [None]
        else:
            temperatures = settings.get_temperatures()

        if log_level:
            print "------ Settings ------"
            print "Mesh sampling: [ %d %d %d ]" % tuple(mesh)
            if grid_points is not None:
                print "Grid points to be calculated:", grid_points
            sys.stdout.flush()

        if options.is_joint_dos:
            jdos = JointDOS(supercell_dm,
                            primitive_dm,
                            mesh,
                            fc2,
                            nac_params=nac_params,
                            sigma=sigma,
                            frequency_step=freq_step,
                            factor=factor,
                            frequency_factor=freq_factor,
                            frequency_scale_factor=freq_scale,
                            is_nosym=options.is_nosym,
                            symprec=options.symprec,
                            log_level=log_level)
            jdos.get_jointDOS(grid_points,
                              filename=options.output_filename)
        else:
            phono3py = Phono3py(supercell,
                                primitive,
                                mesh,
                                fc3,
                                band_indices=settings.get_band_indices(),
                                frequency_factor_to_THz=factor,
                                is_nosym=options.is_nosym,
                                symmetrize_fc3_q=options.is_symmetrize_fc3_q,
                                symprec=options.symprec,
                                log_level=log_level,
                                lapack_zheev_uplo=options.uplo)

            phono3py.set_dynamical_matrix(fc2,
                                          supercell_dm,
                                          primitive_dm,
                                          nac_params=nac_params,
                                          nac_q_direction=nac_q_direction,
                                          frequency_scale_factor=freq_scale)

            if settings.get_is_linewidth():
                phono3py.get_linewidth(
                    grid_points,
                    sigmas=multiple_sigmas,
                    t_max=settings.get_max_temperature(),
                    t_min=settings.get_min_temperature(),
                    t_step=settings.get_temperature_step(),
                    filename=options.output_filename)

            elif settings.get_is_bterta():
                phono3py.get_thermal_conductivity(
                    sigmas=multiple_sigmas,
                    t_max=settings.get_max_temperature(),
                    t_min=settings.get_min_temperature(),
                    t_step=settings.get_temperature_step(),
                    grid_points=grid_points,
                    mesh_divisors=settings.get_mesh_divisors(),
                    coarse_mesh_shifts=settings.get_coarse_mesh_shifts(),
                    no_kappa_stars=settings.get_no_kappa_stars(),
                    gv_delta_q=settings.get_group_velocity_delta_q(),
                    write_gamma=settings.get_write_gamma(),
                    read_gamma=settings.get_read_gamma(),
                    write_amplitude=settings.get_write_amplitude(),
                    read_amplitude=settings.get_read_amplitude(),
                    filename=options.output_filename)
            else:
                phono3py.get_imag_self_energy(
                    grid_points,
                    frequency_step=freq_step,
                    sigmas=multiple_sigmas,
                    temperatures=temperatures,
                    filename=options.output_filename)
    else:
        pass

if log_level:
    print_end()