phonopy/scripts/phono3py

721 lines
28 KiB
Plaintext
Raw Normal View History

2012-12-11 16:19:23 +08:00
#!/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.forces import Forces
from phonopy.harmonic.force_constants import get_force_constants, set_permutation_symmetry
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix
from phonopy.hphonopy.file_IO import parse_FORCE_SETS
from phonopy.units import THzToCm, VaspToTHz
from phonopy.hphonopy.file_IO import parse_BORN
from anharmonic.fc3 import get_fc3, symmetrize_fc3
from anharmonic.file_IO import write_fc2_dat, write_fc3_dat,\
parse_disp_yaml, parse_FORCES_THIRD, write_FORCES_THIRD,\
parse_DELTA_FORCES, write_supercells_with_displacements,\
parse_QPOINTS3, parse_fc3, parse_fc2,\
write_DELTA_FC2_SETS, parse_DELTA_FC2_SETS,\
write_FC2_FOURTH_SETS, parse_FC2_FOURTH_SETS
from anharmonic.displacement_fc3 import get_third_order_displacements
from anharmonic.settings import Phono3pyConfParser
from anharmonic import Phono3py, JointDOS, get_gruneisen_parameters
2012-12-11 16:19:23 +08:00
# 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,
cell_poscar=None,
delta_fc2=False,
factor=None,
delta_fc2_sets_mode=False,
fc2_fourth_sets_mode=False,
freq_scale=None,
gamma_option=0,
grid_points=None,
is_cm=False,
is_decay_channel=False,
is_nodiag=False,
is_displacement=False,
is_nosym=False,
is_gruneisen=False,
is_joint_dos=False,
2013-02-08 17:34:46 +08:00
is_linewidth=False,
2012-12-11 16:19:23 +08:00
is_lifetime=False,
is_nac=False,
is_plusminus_displacements=False,
is_read_triplets=False,
is_translational_symmetry=False,
is_symmetrize_fc2=False,
is_symmetrize_fc3_r=False,
is_symmetrize_fc3_q=False,
is_Peierls=False,
loglevel=None,
mesh_numbers=None,
q_direction=None,
primitive_axis=None,
read_fc2=False,
read_fc2_extra=False,
read_fc3=False,
r2q_TI_index=None,
omega_step=None,
output_filename=None,
qpoints=None,
2012-12-11 16:19:23 +08:00
sigma=None,
supercell_dimension=None,
supercell_dimension_extra=None,
symprec=1e-5,
tmax=None,
tmin=None,
tstep=None,
temperature=None,
verbose=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("-c", "--cell", dest="cell_poscar",
action="store", type="string",
help="Read unit cell", metavar="FILE")
parser.add_option("--cm", dest="is_cm",
action="store_true", help="Convert THz to cm^-1")
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("-o", dest="output_filename",
type="string",
help="Filename of output of damping function")
parser.add_option("--nodiag", dest="is_nodiag",
action="store_true",
help="Set displacements parallel to axes")
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("--create_fc2_fourth",
dest="fc2_fourth_sets_mode",
action="store_true",
help="Create FC2_FOURTH_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("--delta_fc2",
dest="read_delta_fc2",
action="store_true",
help="Read DELTA_FC2_SETS")
parser.add_option("--fc2_fourth",
dest="read_fc2_fourth_sets",
action="store_true",
help="Read FC2_FOURTH_SETS")
parser.add_option("--freq_scale", dest="freq_scale", type="float",
2013-02-06 23:29:22 +08:00
help="Scale factor is multiplied to frequencies only, i.e., changes frequencies but assumed not to change the physical unit")
2012-12-11 16:19:23 +08:00
parser.add_option("-g", "--grid_points",
dest="grid_points",
type="string",
help="Fixed grid points where damping functions are 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("--gruneisen", dest="is_gruneisen",
action="store_true",
help="Calculate phonon Gruneisen parameter")
parser.add_option("--jdos",
dest="is_joint_dos",
action="store_true",
help="Calculate joint density of states")
parser.add_option("--lt", "--lifetime", dest="is_lifetime",
action="store_true",
2013-02-08 17:34:46 +08:00
help="Calculate lifetimes")
parser.add_option("--lw", "--linewidth", dest="is_linewidth",
action="store_true",
help="Calculate linewidths")
2012-12-11 16:19:23 +08:00
parser.add_option("--mesh",
dest="mesh_numbers",
type="string",
help="Mesh numbers")
parser.add_option("--nac", dest="is_nac",
action="store_true",
help="Non-analytical term correction")
parser.add_option("--nosym", dest="is_nosym",
action="store_true",
help="No symmetrization of triplets")
parser.add_option("--omega_step", dest="omega_step", type="float",
help="Pitch of frequency calculated")
parser.add_option("--pa", "--primitive_axis", dest="primitive_axis",
action="store", type="string",
help="Same as PRIMITIVE_AXIS tags")
parser.add_option("--peierls", dest="is_Peierls",
action="store_true",
help="Peierls approximation")
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")
2012-12-11 16:19:23 +08:00
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("--r2q", dest="r2q_TI_index", type="int",
help="Index to set translational invariance in transforming fc3 from real to reciprocal space")
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("--tsym", dest="is_translational_symmetry",
action="store_true",
help="Impose translational invariance condition")
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("--tolerance", dest="symprec", type="float",
help="Symmetry tolerance to search")
parser.add_option("--tp",
dest="is_read_triplets",
action="store_true",
help="Read triplets")
parser.add_option("-v", "--verbose", dest="verbose",
action="store_true",
help="Detailed run-time information is displayed")
parser.add_option("--loglevel", dest="loglevel", type="int",
help="Log level")
parser.add_option("--temp", dest="temperature", type="float",
help="Temperature")
(options, args) = parser.parse_args()
option_list = parser.option_list
# Log level
if options.verbose:
log_level = 1
elif not options.loglevel==None:
log_level = options.loglevel
else:
log_level = 0
# 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)
# Create FC2_FOURTH_SETS
if options.fc2_fourth_sets_mode:
displacements = parse_disp_yaml('disp.yaml')
write_FC2_FOURTH_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 > 0:
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 > 0:
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 not settings.get_supercell_matrix_extra()==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 options.is_nac:
print "Non-analytical term correction:", options.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)
# fc2
if options.read_fc2 or options.read_delta_fc2:
if log_level:
print "----- Read fc2 -----"
sys.stdout.flush()
fc2_with_dim = parse_fc2(supercell.get_number_of_atoms())
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)
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 = parse_fc2(supercell_extra.get_number_of_atoms(),
filename='fc2_extra.dat')
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')
if options.is_symmetrize_fc2:
set_permutation_symmetry(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 > 0:
print_end()
sys.exit(1)
else:
nac_params = None
nac_q_direction = None
# fc3
if not options.is_joint_dos:
if options.read_fc3: # Read FORCES_THIRD
if log_level:
print "----- Read fc3 -----"
sys.stdout.flush()
fc3 = parse_fc3(supercell.get_number_of_atoms())
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 in range(num_atom):
for j in range(num_atom):
for k in range(3):
for l in range(3):
for m in range(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:", 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)
if options.read_fc2_fourth_sets:
displacements = parse_disp_yaml('disp.yaml')
parse_FC2_FOURTH_SETS(displacements)
count = 0
for first_disp in displacements['first_atoms']:
count += 1
print count
print first_disp['fc2'][0, 0]
for first_disp in displacements['first_atoms']:
for second_disp in first_disp['second_atoms']:
for disp in second_disp['delta_fc2']:
count += 1
print count
print disp[0, 0]
sys.exit(0)
#============================
# Phonon Gruneisen parameter
#============================
if options.is_gruneisen:
mesh = settings.get_mesh_numbers()
if log_level:
print "------ Phonon Gruneisen parameter ------"
if mesh is not None:
print "Mesh:", mesh
sys.stdout.flush()
gr = get_gruneisen_parameters(fc2,
fc3,
supercell,
primitive,
2013-02-06 00:17:30 +08:00
factor=VaspToTHz,
is_ion_clamped=False,
symprec=options.symprec)
qpoints = settings.get_qpoints()
if mesh is not None:
gr.set_sampling_mesh(mesh,
is_gamma_center=True)
else:
gr.set_qpoints(qpoints)
gr.run()
2013-02-06 00:17:30 +08:00
gr.write_yaml()
2012-12-11 16:19:23 +08:00
#===========================
# phonon-phonon interaction
#===========================
2013-02-10 14:59:14 +08:00
else:
2012-12-11 16:19:23 +08:00
mesh = settings.get_mesh_numbers()
2013-02-08 22:44:52 +08:00
if options.grid_points is None:
grid_points = None
2012-12-11 16:19:23 +08:00
else:
2013-02-06 00:17:30 +08:00
grid_points = np.array([int(x)
for x in options.grid_points.split()])
2012-12-11 16:19:23 +08:00
if log_level:
print "------ Phonon-phonon ------"
print "Mesh:", mesh
2013-02-08 22:44:52 +08:00
if grid_points is not None:
print "Grid points to be calculated:", grid_points
2012-12-11 16:19:23 +08:00
sys.stdout.flush()
if options.factor==None:
factor = VaspToTHz
else:
factor = options.factor
if options.is_cm:
freq_factor = THzToCm
else:
freq_factor = 1.0
if settings.get_sigma()==None:
sigma = 0.2 * freq_factor
else:
sigma = settings.get_sigma()
if settings.get_omega_step()==None:
omega_step = 0.1 * freq_factor
else:
omega_step = settings.get_omega_step()
if options.freq_scale == None:
freq_scale = 1.0
else:
freq_scale = options.freq_scale
2013-02-08 17:34:46 +08:00
if settings.get_supercell_matrix_extra() is None:
supercell_dm = supercell
primitive_dm = primitive
else:
supercell_dm = supercell_extra
primitive_dm = primitive_extra
2012-12-11 16:19:23 +08:00
2013-02-08 17:34:46 +08:00
if options.is_joint_dos:
jdos = JointDOS(supercell_dm,
primitive_dm,
mesh,
fc2,
nac_params,
sigma,
omega_step,
factor,
freq_factor,
freq_scale,
options.is_nosym,
options.symprec,
log_level)
2012-12-11 16:19:23 +08:00
jdos.get_jointDOS(grid_points,
2013-02-08 17:34:46 +08:00
filename=options.output_filename)
2012-12-11 16:19:23 +08:00
else:
phono3py = Phono3py(supercell,
primitive,
mesh,
fc3,
factor,
freq_factor,
freq_scale,
options.is_nosym,
options.is_symmetrize_fc3_q,
options.is_read_triplets,
options.r2q_TI_index,
options.is_Peierls,
options.symprec,
log_level)
2013-02-08 17:34:46 +08:00
phono3py.set_dynamical_matrix(fc2,
supercell_dm,
primitive_dm,
nac_params,
nac_q_direction)
2012-12-11 16:19:23 +08:00
if options.is_decay_channel:
phono3py.get_decay_channels(grid_points,
settings.get_band_indices(),
options.temperature)
2013-02-08 17:34:46 +08:00
elif settings.get_is_linewidth():
phono3py.get_linewidth(
2013-02-06 23:29:22 +08:00
grid_points,
settings.get_band_indices(),
2013-02-08 17:34:46 +08:00
sigma=sigma,
t_max=settings.get_max_temperature(),
t_min=settings.get_min_temperature(),
t_step=settings.get_temperature_step(),
2013-02-06 23:29:22 +08:00
gamma_option=options.gamma_option,
filename=options.output_filename)
2012-12-11 16:19:23 +08:00
2013-02-08 17:34:46 +08:00
elif settings.get_is_lifetime():
2013-02-10 00:30:08 +08:00
phono3py.get_thermal_conductivity(
2013-02-08 22:44:52 +08:00
sigma=sigma,
t_max=settings.get_max_temperature(),
t_min=settings.get_min_temperature(),
t_step=settings.get_temperature_step(),
gamma_option=options.gamma_option,
filename=options.output_filename)
2012-12-11 16:19:23 +08:00
else:
2013-02-07 22:09:36 +08:00
phono3py.get_damping_function(
2013-02-08 17:34:46 +08:00
grid_points,
settings.get_band_indices(),
2013-02-07 22:09:36 +08:00
sigma=sigma,
omega_step=omega_step,
temperature=options.temperature,
filename=options.output_filename,
gamma_option=options.gamma_option)
2012-12-11 16:19:23 +08:00
if log_level:
print_end()