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 print_cell
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.file_IO import parse_BORN, parse_FORCE_SETS
from phonopy.units import VaspToTHz
from anharmonic.phonon3.fc3 import show_drift_fc3
from anharmonic.file_IO import parse_disp_fc2_yaml, parse_disp_fc3_yaml, \
     parse_FORCES_FC2, parse_FORCES_FC3, write_FORCES_FC3, \
     write_fc3_to_hdf5, write_fc2_to_hdf5, read_fc3_from_hdf5, \
     read_fc2_from_hdf5, write_ir_grid_points, write_grid_address, \
     write_supercells_with_displacements
from anharmonic.phonon3.triplets import get_coarse_ir_grid_points
from anharmonic.settings import Phono3pyConfParser
from anharmonic.phonon3 import Phono3py, Phono3pyJointDos, IsotopeScattering, \
     get_gruneisen_parameters

phono3py_version = "0.8.2"

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

def print_version(version):
    print " " * 42, version
    print ""

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

def print_error():
    print """  ___ _ __ _ __ ___  _ __ 
 / _ \ '__| '__/ _ \| '__|
|  __/ |  | | | (_) | |   
 \___|_|  |_|  \___/|_|
"""

def print_error_message(message):
    print message

def file_exists(filename, log_level):
    if os.path.exists(filename):
        return True
    else:
        error_text = "%s not found." % filename
        print_error_message(error_text)
        if log_level > 0:
            print_error()
        sys.exit(1)
    
# Parse options
parser = OptionParser()
parser.set_defaults(band_indices=None,
                    band_paths=None,
                    band_points=None,
                    cell_poscar=None,
                    cutoff_fc3_distance=None,
                    cutoff_frequency=1e-2,
                    cutoff_lifetime=1e-4,
                    cutoff_pair_distance=None,
                    delta_fc2=False,
                    displacement_distance=0.03,
                    delta_fc2_sets_mode=False,
                    factor=None,
                    freq_scale=None,
                    gamma_option=0,
                    grid_points=None,
                    gv_delta_q=None,
                    input_filename=None,
                    input_output_filename=None,
                    ion_clamped=False,
                    is_decay_channel=False,
                    is_nodiag=False,
                    is_displacement=False,
                    is_nosym=False,
                    is_gruneisen=False,
                    is_isotope=False,
                    is_joint_dos=False,
                    is_linewidth=False,
                    is_frequency_shift=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,                     
                    is_tetrahedron_method=False,
                    log_level=None,
                    max_freepath=None,
                    mass_variances=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,
                    frequency_pitch=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="displacement_distance", 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("--cutfc3", "--cutoff_fc3_distance",
                  dest="cutoff_fc3_distance",
                  type="float",
                  help="Cutoff distance of third-order force constants. Elements where any pair of atoms has larger distance than cut-off distance are set zero.")
parser.add_option("--cutfreq", "--cutoff_frequency", dest="cutoff_frequency",
                  type="float",
                  help="Phonon modes below this frequency are ignored.")
parser.add_option("--cutlt", "--cutoff_lifetime", dest="cutoff_lifetime",
                  type="float",
                  help="Cutoff lifetime in second used for thermal conductivity calculation")
parser.add_option("--cutpair", "--cutoff_pair_distance",
                  dest="cutoff_pair_distance",
                  type="float",
                  help="Cutoff distance between pairs of displaced atoms used for supercell creation with displacements and making third-order force constants")
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_pitch", dest="frequency_pitch", type="float",
                  help="Pitch in frequency for spectrum")
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("-i", dest="input_filename",
                  type="string",
                  help="Input filename extension")
parser.add_option("--io", dest="input_output_filename",
                  type="string",
                  help="Input and output filename extension")
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("--isotope", dest="is_isotope",
                  action="store_true",
                  help="Isotope scattering lifetime")
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("--fst", "--frequency_shift",
                  dest="is_frequency_shift",
                  action="store_true",
                  help="Calculate frequency shifts")
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("--mv", "--mass_variances", dest="mass_variances",
                  type="string",
                  help="Mass variance parameters for isotope scattering")
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="Output filename extension")
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("--thm", "--tetrahedron_method", dest="is_tetrahedron_method",
                  action="store_true",
                  help="Use tetrahedron method")
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

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

# Title
if log_level:
    print_phono3py()
    print_version(phono3py_version)

#########################################
# Create FORCES_THIRD and FORCES_SECOND #
#########################################
if options.forces_third_mode:
    if input_filename is None:
        filename = 'disp_fc3.yaml'
    else:
        filename = 'disp_fc3.' + input_filename + '.yaml'
    file_exists(filename, log_level)
    if log_level:
        print "Displacement dataset is read from %s." % filename
    disp_dataset = parse_disp_fc3_yaml()
    write_FORCES_FC3(args, disp_dataset)
    print_end()
    exit(0)

##################
# Parse settings #
##################
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 crystal structure (POSCAR) #
###################################
if options.cell_poscar == None:
    file_exists('POSCAR', log_level)
    unitcell_filename = 'POSCAR'
else:
    file_exists(options.cell_poscar, log_level)
    unitcell_filename = options.cell_poscar
unitcell = read_vasp(unitcell_filename, settings.get_chemical_symbols())

#################################################
# Create supercells with displacements and exit #
#################################################
if options.is_displacement:
    cutoff_pair_distance = settings.get_cutoff_pair_distance()
    phono3py = Phono3py(unitcell,
                        settings.get_supercell_matrix(),
                        symprec=options.symprec)
    supercell = phono3py.get_supercell()
    phono3py.generate_displacements(
        distance=settings.get_displacement_distance(),
        cutoff_pair_distance=cutoff_pair_distance,
        is_plusminus=settings.get_is_plusminus_displacement(),
        is_diagonal=settings.get_is_diagonal_displacement())
    dds = phono3py.get_displacement_dataset()
    
    if log_level:
        print
        print "Displacement distance:", settings.get_displacement_distance()

    if output_filename is None:
        filename = 'disp_fc3.yaml'
    else:
        filename = 'disp_fc3.' + output_filename + '.yaml'
        
    num_disps, num_disp_files = write_supercells_with_displacements(
        supercell,
        dds,
        filename=filename)
        
    if log_level:
        print "Total number of displacements:", num_disps
        if cutoff_pair_distance is not None:
            print "Cutoff distance for displacements:", 
            print cutoff_pair_distance
            print "Number of displacement supercell files created:",
            print  num_disp_files
        print_end()
    sys.exit(0)

##############
# Initialize #
##############
mesh = settings.get_mesh_numbers()
mesh_divs = settings.get_mesh_divisors()
grid_points = settings.get_grid_points()
band_indices = settings.get_band_indices()
cutoff_frequency = settings.get_cutoff_frequency()
if settings.get_multiple_sigmas() is None:
    if settings.get_sigma():
        sigmas = [settings.get_sigma()]
    else:
        sigmas = []
else:
    sigmas = settings.get_multiple_sigmas()
if settings.get_is_tetrahedron_method():
    sigmas = [None] + sigmas
if settings.get_temperatures() is None:
    t_max=settings.get_max_temperature()
    t_min=settings.get_min_temperature()
    t_step=settings.get_temperature_step()
    temperature_points = [0.0, 300.0] # For spectra
    temperatures = np.arange(t_min, t_max + float(t_step) / 10, t_step)
else:
    temperature_points = settings.get_temperatures() # For spectra
    temperatures = settings.get_temperatures() # For others
if options.factor is None:
    frequency_factor_to_THz = VaspToTHz
else:
    frequency_factor_to_THz = options.factor
if settings.get_frequency_pitch() is None:
    frequency_step = 0.1
else:
    frequency_step = settings.get_frequency_pitch()
if options.freq_scale is None:
    frequency_scale_factor = 1.0
else:
    frequency_scale_factor = options.freq_scale

phono3py = Phono3py(
    unitcell,
    settings.get_supercell_matrix(),
    primitive_matrix=settings.get_primitive_matrix(),
    phonon_supercell_matrix=settings.get_phonon_supercell_matrix(),
    mesh=mesh,
    band_indices=band_indices,
    sigmas=sigmas,
    cutoff_frequency=cutoff_frequency,
    frequency_factor_to_THz=frequency_factor_to_THz,
    is_symmetry=True,
    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)

supercell = phono3py.get_supercell()
primitive = phono3py.get_primitive()
phonon_supercell = phono3py.get_phonon_supercell()
phonon_primitive = phono3py.get_phonon_primitive()
symmetry = phono3py.get_symmetry()

#################
# Show settings #
#################
if log_level:
    print "Spacegroup: ", symmetry.get_international_table()
    print "---------------------------- primitive cell --------------------------------"
    print_cell(primitive)
    print "------------------------------- supercell ----------------------------------"
    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 settings.get_phonon_supercell_matrix() is not None:
        print "-------------------- primitive cell for harmonic phonon ---------------------"
        print_cell(phonon_primitive)
        print "---------------------- supercell for harmonic phonon ------------------------"
        print_cell(phonon_supercell, mapping=phonon_primitive.get_supercell_to_primitive_map())
        print "--------------- ratio (phonon supercell)/(phonon primitive) -----------------"
        for vec in np.dot(phonon_supercell.get_cell(),
                          np.linalg.inv(phonon_primitive.get_cell())):
            print "%5.2f"*3 % tuple(vec)
            
    print "----------------------------------------------------------------------------"
    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:",
        print options.is_symmetrize_fc3_r
    if options.is_symmetrize_fc3_q:
        print "FC3 symmetry of index exchange in reciprocal space:",
        print options.is_symmetrize_fc3_q
    if settings.get_cutoff_fc3_distance() is not None:
        print "FC3 cutoff distance:",
        print settings.get_cutoff_fc3_distance()
    if settings.get_phonon_supercell_matrix() is not None:
        print "Supercell for harmonic phonon is introduced."
    if settings.get_is_nac():
        print "Non-analytical term correction:", settings.get_is_nac()

#####################################################  
# Write ir-grid points and grid addresses, and exit #
#####################################################
if options.write_grid_points:
    print "---------------------------------------------------------------------------"
    if mesh is None:
        print "To write grid points, mesh numbers have to be specified."
    else:
        (ir_grid_points,
         coarse_grid_weights,
         grid_address) = get_coarse_ir_grid_points(
             primitive,
             mesh,
             mesh_divs,
             settings.get_coarse_mesh_shifts(),
             is_nosym=options.no_kappa_stars,
             symprec=options.symprec)
        write_ir_grid_points(mesh,
                             mesh_divs,
                             ir_grid_points,
                             coarse_grid_weights,
                             grid_address,
                             np.linalg.inv(primitive.get_cell()))
        gadrs_fname = write_grid_address(grid_address, mesh)

        print "Ir-grid points are written into \"ir_grid_points.yaml\"."
        print "Grid address are written into \"%s\"." % gadrs_fname

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

#######
# fc3 #
#######
if (options.is_joint_dos or
    options.is_isotope or
    settings.get_read_gamma() or
    settings.get_read_amplitude()):
    pass
else:
    if options.read_fc3: # Read fc3.hdf5
        if input_filename is None:
            filename = 'fc3.hdf5'
        else:
            filename = 'fc3.' + input_filename + '.hdf5'
        file_exists(filename, log_level)
        if log_level:
            print "Reading fc3 from %s..." % filename
        fc3 = read_fc3_from_hdf5(filename=filename)
        phono3py.set_fc3(fc3)
    else: # fc3 from FORCES_THIRD and FORCES_SECOND
        if input_filename is None:
            filename = 'disp_fc3.yaml'
        else:
            filename = 'disp_fc3.' + input_filename + '.yaml'
        file_exists(filename, log_level)
        if log_level:
            print "Displacement dataset is read from %s." % filename
            print "Solving fc3:"
        disp_dataset = parse_disp_fc3_yaml(filename=filename)
        forces_fc3 = parse_FORCES_FC3(disp_dataset)
        phono3py.produce_fc3(
            forces_fc3,
            disp_dataset,
            cutoff_distance=settings.get_cutoff_fc3_distance(),
            is_translational_symmetry=options.is_translational_symmetry,
            is_permutation_symmetry=options.is_symmetrize_fc3_r)
        if output_filename is None:
            filename = 'fc3.hdf5'
        else:
            filename = 'fc3.' + output_filename + '.hdf5'
        if log_level:
            print "Writing fc3 to %s" % filename
        write_fc3_to_hdf5(phono3py.get_fc3(), filename=filename)

    show_drift_fc3(phono3py.get_fc3())

##############
# phonon fc2 #
##############
if options.read_fc2:
    if input_filename is None:
        filename = 'fc2.hdf5'
    else:
        filename = 'fc2.' + input_filename + '.hdf5'
    file_exists(filename, log_level)
    if log_level:
        print "Reading fc2 from %s..." % filename
    phonon_fc2 = read_fc2_from_hdf5(filename=filename)
    phono3py.set_fc2(phonon_fc2)
else:
    if settings.get_phonon_supercell_matrix() is None:
        if phono3py.get_fc2() is None:
            if input_filename is None:
                filename = 'disp_fc3.yaml'
            else:
                filename = 'disp_fc3.' + input_filename + '.yaml'
            if log_level:
                print "Displacement dataset is read from %s." % filename
                print "Solving fc2:"
            file_exists(filename, log_level)
            disp_dataset = parse_disp_fc3_yaml(filename=filename)
            forces_fc3 = parse_FORCES_FC3(disp_dataset)
            phono3py.produce_fc2(forces_fc3, disp_dataset)
    else:
        if input_filename is None:
            filename = 'disp_fc2.yaml'
        else:
            filename = 'disp_fc2.' + input_filename + '.yaml'
        if log_level:
            print "Displacement dataset is read from %s." % filename
            print "Solving fc2:"
        file_exists(filename, log_level)
        disp_dataset = parse_disp_fc2_yaml(filename=filename)
        forces_fc2 = parse_FORCES_FC2(disp_dataset)
        phono3py.produce_fc2(forces_fc2, disp_dataset)

    if output_filename is None:
        filename = 'fc2.hdf5'
    else:
        filename = 'fc2.' + output_filename + '.hdf5'
    if log_level:
        print "Writing fc2 to %s..." % filename
    write_fc2_to_hdf5(phono3py.get_fc2(), filename=filename)
    
show_drift_force_constants(phono3py.get_fc2(), name='fc2')

if settings.get_is_nac():
    file_exists('BORN', log_level)
    nac_params = parse_BORN(phonon_primitive)
    nac_q_direction = settings.get_nac_q_direction()
else:
    nac_params = None
    nac_q_direction = None

##############################
# Phonon Gruneisen parameter #
##############################
if options.is_gruneisen:
    fc2 = phono3py.get_fc2()
    fc3 = phono3py.get_fc3()
    if len(fc2) != len(fc3):
        print_error_message("Supercells used for fc2 and fc3 have to be same.")
        if log_level:
            print_error()
        sys.exit(1)
    
    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):

        print_error_message("An option of --mesh, --band, or --qpoints "
                            "has to be specified.")
        if log_level:
            print_error()
        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()

    if output_filename is None:
        filename = 'gruneisen3.yaml'
    else:
        filename = 'gruneisen3.' + output_filename + '.yaml'
    gr.write_yaml(filename=filename)

    if log_level:
        print_end()

if log_level:
    print "--------------------------- Settings ----------------------------"
    print "Mesh sampling: [ %d %d %d ]" % tuple(mesh)
    if mesh_divs is not None and settings.get_is_bterta():
        print "Mesh divisors: [ %d %d %d ]" % tuple(mesh_divs)
    if band_indices is not None and not settings.get_is_bterta():
        print "Band indices:", band_indices
    if sigmas:
        print "Sigma:",
        for sigma in sigmas:
            if sigma:
                print sigma,
            else:
                print "Tetrahedron-method",
        print
    if (settings.get_is_linewidth() or
        settings.get_is_frequency_shift() or
        settings.get_is_bterta()):
        print "Temperature:",
        if len(temperatures) > 5:
            print ((" %.1f " * 5) % tuple(temperatures[:5])), "...",
            print " %.1f" % temperatures[-1]
        else:
            print ("%.1f  " * len(temperatures)) % tuple(temperatures)
    elif options.is_isotope or options.is_joint_dos:
        pass
    else:
        print "Temperatures:",
        print ("%.1f  " * len(temperature_points)) % tuple(temperature_points)
    if grid_points is not None:
        print "Grid point to be calculated:",
        if len(grid_points) > 8:
            for i, gp in enumerate(grid_points):
                if i % 10 == 0:
                    print
                    print " ",
                print gp,
            print
        else:
            for gp in grid_points:
                print gp,
            print
    if cutoff_frequency:
        print "Cutoff frequency:", cutoff_frequency
    if log_level > 1:
        print "Frequency factor to THz:", frequency_factor_to_THz
        print "Frequency step for spectrum:", frequency_step
        print "Frequency scale factor:", frequency_scale_factor
    sys.stdout.flush()
    
#############
# Joint DOS #
#############
if options.is_joint_dos:
    joint_dos = Phono3pyJointDos(
        supercell_dm,
        primitive_dm,
        mesh,
        fc2,
        nac_params=nac_params,
        sigmas=sigmas,
        frequency_step=frequency_step,
        frequency_factor_to_THz=frequency_factor_to_THz,
        frequency_scale_factor=frequency_scale_factor,
        is_nosym=options.is_nosym,
        symprec=options.symprec,
        output_filename=output_filename,
        log_level=log_level)
    joint_dos.run(grid_points)
    if log_level:
        print_end()
    
######################
# Isotope scattering #
######################
if options.is_isotope:
    if log_level:
        print "Cutoff frequency:", settings.get_cutoff_frequency()
    mass_variances = settings.get_mass_variances()
    iso = IsotopeScattering(
        mesh,
        mass_variances,
        frequency_factor_to_THz=frequency_factor_to_THz,
        symprec=options.symprec,
        cutoff_frequency=settings.get_cutoff_frequency(),
        lapack_zheev_uplo=options.uplo)
    iso.set_dynamical_matrix(fc2,
                             supercell,
                             primitive,
                             nac_params=nac_params,
                             frequency_scale_factor=frequency_scale_factor)
    if settings.get_band_indices() is None:
        band_indices = None
    else:
        band_indices = np.hstack(self._band_indices).astype('intc')
        
        for sigma in sigmas:
            if log_level:
                print "Sigma:", sigma
            iso.set_sigma(sigma)
            print iso.run(gp, band_indices)
    if log_level:
        print_end()
            
#########
# Ph-ph #
#########
phono3py.set_phph_interaction(nac_params=nac_params,
                              nac_q_direction=nac_q_direction,
                              frequency_scale_factor=frequency_scale_factor)

if settings.get_is_linewidth():
    phono3py.run_linewidth(
        grid_points,
        temperatures=temperatures)
    phono3py.write_linewidth(filename=output_filename)
elif settings.get_is_frequency_shift():
    phono3py.get_frequency_shift(
        grid_points,
        epsilon=sigma,
        temperatures=temperatures,
        output_filename=output_filename)
elif settings.get_is_bterta():
    phono3py.run_thermal_conductivity(
        temperatures=temperatures,
        sigmas=sigmas,
        mass_variances=settings.get_mass_variances(),
        grid_points=grid_points,
        mesh_divisors=mesh_divs,
        coarse_mesh_shifts=settings.get_coarse_mesh_shifts(),
        cutoff_lifetime=settings.get_cutoff_lifetime(),
        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(),
        input_filename=input_filename,
        output_filename=output_filename)
else:
    phono3py.run_imag_self_energy(
        grid_points,
        frequency_step=frequency_step,
        temperatures=temperature_points)
    phono3py.write_imag_self_energy(filename=output_filename)

if log_level:
    print_end()