scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Under fixing kaccum.

This commit is contained in:
Atsushi Togo 2021-06-17 23:12:05 +09:00
parent 63b17293f2
commit 77c14dcdee
1 changed files with 304 additions and 271 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

575
scripts/phono3py-kaccum
View File

@ -7,6 +7,7 @@ import h5py
from phonopy.structure.cells import get_primitive
from phonopy.structure.symmetry import Symmetry
from phonopy.phonon.tetrahedron_mesh import TetrahedronMesh
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.harmonic.force_constants import similarity_transformation
from phonopy.interface.calculator import read_crystal_structure
from phono3py.phonon.grid import (
@ -52,7 +53,7 @@ def set_T_target(temperatures,
T_target,
mean_freepath=None):
for i, t in enumerate(temperatures):
if np.abs(t - args.temperature) < epsilon:
if np.abs(t - T_target) < epsilon:
temperatures = temperatures[i:i+1]
mode_prop = mode_prop[i:i+1, :, :]
if mean_freepath is not None:
@ -127,7 +128,27 @@ def fracval(frac):
return float(x[0]) / float(x[1])
def get_grid_symmetry(primitive, symmetry, mesh, weights, qpoints):
def get_integration_weights(grid_points, bz_grid, frequencies, freq_points):
import phono3py._phono3py as phono3c
thm = TetrahedronMethod(bz_grid.microzone_lattice)
num_grid_points = len(grid_points)
num_band = frequencies.shape[1]
integration_weights = np.zeros(
(num_grid_points, num_band, num_band), dtype='double')
phono3c.integration_weights_at_grid_points(
integration_weights,
freq_points,
np.array(np.dot(thm.get_tetrahedra(), bz_grid.P.T),
dtype='int_', order='C'),
bz_grid.D_diag,
grid_points,
frequencies,
bz_grid.addresses,
bz_grid.gp_map,
bz_grid.is_dense_gp_map * 1 + 1)
def get_grid_symmetry(bz_grid, weights, qpoints):
(ir_grid_points,
weights_for_check,
grid_mapping_table) = get_ir_grid_points(bz_grid)
@ -140,12 +161,14 @@ def get_grid_symmetry(primitive, symmetry, mesh, weights, qpoints):
print("*******************************")
raise
qpoints_for_check = bz_grid.grid_address[
ir_grid_points] / mesh.astype('double')
addresses = bz_grid.grid_addresses[ir_grid_points]
D_diag = bz_grid.D_diag.astype('double')
qpoints_for_check = np.dot(bz_grid.Q, addresses / D_diag)
diff_q = qpoints - qpoints_for_check
np.testing.assert_almost_equal(diff_q, np.rint(diff_q))
return ir_grid_points, grid_address, grid_mapping_table
return ir_grid_points, grid_mapping_table
def get_gv_by_gv(gv,
@ -179,6 +202,282 @@ def get_gv_by_gv(gv,
return gv_sum2
def get_calculator(args):
"""Return calculator name."""
interface_mode = None
if args.qe_mode:
interface_mode = 'qe'
elif args.crystal_mode:
interface_mode = 'crystal'
elif args.abinit_mode:
interface_mode = 'abinit'
elif args.turbomole_mode:
interface_mode = 'turbomole'
return interface_mode
def read_files(args):
"""Read crystal structure and kappa.hdf5 files."""
interface_mode = get_calculator(args)
if len(args.filenames) > 1:
cell, _ = read_crystal_structure(args.filenames[0],
interface_mode=interface_mode)
f = h5py.File(args.filenames[1], 'r')
else:
cell, _ = read_crystal_structure(args.cell_filename,
interface_mode=interface_mode)
f = h5py.File(args.filenames[0], 'r')
return cell, f
def get_mode_property(args, f_kappa):
if args.pqj:
mode_prop = f_kappa['ave_pp'][:].reshape(
(1,) + f_kappa['ave_pp'].shape)
elif args.cv:
mode_prop = f_kappa['heat_capacity'][:]
elif args.tau:
g = f_kappa['gamma'][:]
g = np.where(g > 0, g, -1)
mode_prop = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0) # tau
elif args.gv_norm:
mode_prop = np.sqrt(
(f_kappa['group_velocity'][:, :, :] ** 2).sum(axis=2))
mode_prop = mode_prop.reshape((1,) + mode_prop.shape)
elif args.gamma:
mode_prop = f_kappa['gamma'][:]
elif args.gruneisen:
mode_prop = f_kappa['gruneisen'][:].reshape(
(1,) + f_kappa['gruneisen'].shape)
mode_prop **= 2
else:
raise RuntimeError("No property target is specified.")
return mode_prop, g
def get_init_params(args, f_kappa):
if 'mesh' in f_kappa:
mesh = np.array(f_kappa['mesh'][:], dtype='int_')
else:
mesh = np.array([int(x) for x in args.mesh.split()], dtype='int_')
if 'temperature' in f_kappa:
temperatures = f_kappa['temperature'][:]
else:
temperatures = None
weights = f_kappa['weight'][:]
return mesh, temperatures, weights
def get_parser():
"""Return args of ArgumentParser."""
parser = argparse.ArgumentParser(description="Show unit cell volume")
parser.add_argument(
"--pa", dest="primitive_matrix", default="1 0 0 0 1 0 0 0 1",
help="Primitive matrix")
parser.add_argument(
"--mesh", dest="mesh", default="1 1 1",
help="Mesh numbers")
parser.add_argument(
"-c", "--cell", dest="cell_filename",
help="Unit cell filename")
parser.add_argument(
'--gv', action='store_true',
help='Calculate for gv_x_gv (tensor)')
parser.add_argument(
'--pqj', action='store_true',
help='Calculate for Pqj (scalar)')
parser.add_argument(
'--cv', action='store_true',
help='Calculate for Cv (scalar)')
parser.add_argument(
'--tau', action='store_true',
help='Calculate for lifetimes (scalar)')
parser.add_argument(
'--gamma', action='store_true',
help='Calculate for Gamma (scalar)')
parser.add_argument(
'--gruneisen', action='store_true',
help='Calculate for mode-Gruneisen parameters squared (scalar)')
parser.add_argument(
'--gv-norm', action='store_true',
help='Calculate for |g_v| (scalar)')
parser.add_argument(
'--mfp', action='store_true',
help='Mean free path is used instead of frequency')
parser.add_argument(
'--temperature', type=float, dest='temperature',
help='Temperature to output data at')
parser.add_argument(
'--nsp', '--num-sampling-points', type=int, dest='num_sampling_points',
default=100,
help="Number of sampling points in frequency or MFP axis")
parser.add_argument(
'--average', action='store_true',
help=("Output the traces of the tensors divided by 3 "
"rather than the unique elements"))
parser.add_argument(
'--trace', action='store_true',
help=("Output the traces of the tensors "
"rather than the unique elements"))
parser.add_argument(
'--smearing', action='store_true',
help='Use smearing method (only for scalar density)')
parser.add_argument(
'--qe', '--pwscf', dest="qe_mode",
action="store_true", help="Invoke Pwscf mode")
parser.add_argument(
'--crystal', dest="crystal_mode",
action="store_true", help="Invoke CRYSTAL mode")
parser.add_argument(
'--abinit', dest="abinit_mode",
action="store_true", help="Invoke Abinit mode")
parser.add_argument(
'--turbomole', dest="turbomole_mode",
action="store_true", help="Invoke TURBOMOLE mode")
parser.add_argument(
"--noks", "--no-kappa-stars",
dest="no_kappa_stars", action="store_true",
help="Deactivate summation of partial kappa at q-stars")
parser.add_argument('filenames', nargs='*')
args = parser.parse_args()
return args
def main():
"""Calculate kappa spectrum."""
args = get_parser()
cell, f_kappa = read_files(args)
mesh, temperatures, weights = get_init_params(args, f_kappa)
primitive_matrix = np.reshape(
[fracval(x) for x in args.primitive_matrix.split()], (3, 3))
primitive = get_primitive(cell, primitive_matrix)
primitive_symmetry = Symmetry(primitive)
bz_grid = BZGrid(mesh,
lattice=primitive.cell,
symmetry_dataset=primitive_symmetry.dataset,
is_dense_gp_map=False)
grid_address = bz_grid.addresses
if args.no_kappa_stars or (weights == 1).all():
ir_grid_points = np.arange(np.prod(mesh), dtype='int_')
grid_mapping_table = np.arange(np.prod(mesh), dtype='int_')
else:
(ir_grid_points,
grid_mapping_table) = get_grid_symmetry(
bz_grid, weights, f_kappa['qpoint'][:])
################
# Set property #
################
if args.gv:
if 'gv_by_gv' in f_kappa:
gv_sum2 = f_kappa['gv_by_gv'][:]
else: # For backward compatibility. This will be removed someday.
gv = f_kappa['group_velocity'][:]
gv_sum2 = get_gv_by_gv(gv,
primitive_symmetry,
primitive,
mesh,
ir_grid_points,
grid_address)
# gv x gv is divied by primitive cell volume.
unit_conversion = primitive.volume
mode_prop = gv_sum2.reshape((1,) + gv_sum2.shape) / unit_conversion
else:
if 'mode_kappa' in f_kappa:
mode_prop = f_kappa['mode_kappa'][:]
else:
mode_prop = None
frequencies = f_kappa['frequency'][:]
conditions = frequencies > 0
if np.logical_not(conditions).sum() > 3:
sys.stderr.write("# Imaginary frequencies are found. "
"They are set to be zero.\n")
frequencies = np.where(conditions, frequencies, 0)
#######
# Run #
#######
if (args.gamma or args.gruneisen or args.pqj or
args.cv or args.tau or args.gv_norm): # noqa E129
mode_prop, g = get_mode_property(args, f_kappa)
if (args.temperature is not None and
not (args.gv_norm or args.pqj or args.gruneisen)): # noqa E129
temperatures, mode_prop = set_T_target(temperatures,
mode_prop,
args.temperature)
if args.smearing:
mode_prop_dos = GammaDOSsmearing(
mode_prop,
frequencies,
weights,
num_fpoints=args.num_sampling_points)
sampling_points, gdos = mode_prop_dos.get_gdos()
else:
mode_prop_dos = GammaDOStetrahedron(
mode_prop,
primitive,
frequencies,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
weights,
num_fpoints=args.num_sampling_points)
sampling_points, gdos = mode_prop_dos.get_gdos()
for i, gdos_t in enumerate(gdos):
total = np.dot(weights, mode_prop[i]).sum() / weights.sum()
assert np.isclose(gdos_t[-1][0], total)
show_scalar(gdos, temperatures, sampling_points, args)
else:
if args.mfp:
if 'mean_free_path' in f_kappa:
mfp = f_kappa['mean_free_path'][:]
mean_freepath = np.sqrt((mfp ** 2).sum(axis=3))
else:
mean_freepath = get_mfp(f_kappa['gamma'][:],
f_kappa['group_velocity'][:])
if args.temperature is not None:
(temperatures,
mode_prop,
mean_freepath) = set_T_target(temperatures,
mode_prop,
args.temperature,
mean_freepath=mean_freepath)
kdos, sampling_points = run_mfp_dos(mean_freepath,
mode_prop,
primitive,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
args.num_sampling_points)
show_tensor(kdos, temperatures, sampling_points, args)
else:
if args.temperature is not None and not args.gv:
temperatures, mode_prop = set_T_target(temperatures,
mode_prop,
args.temperature)
kdos, sampling_points = run_prop_dos(frequencies,
mode_prop,
primitive,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
args.num_sampling_points)
show_tensor(kdos, temperatures, sampling_points, args)
class KappaDOS(object):
def __init__(self,
mode_kappa,
@ -325,271 +624,5 @@ class GammaDOStetrahedron(GammaDOS):
self._gamma[:, i] * self._ir_grid_weights[i], iw.T)
def get_parser():
"""Return args of ArgumentParser."""
parser = argparse.ArgumentParser(description="Show unit cell volume")
parser.add_argument(
"--pa", dest="primitive_matrix", default="1 0 0 0 1 0 0 0 1",
help="Primitive matrix")
parser.add_argument(
"--mesh", dest="mesh", default="1 1 1",
help="Mesh numbers")
parser.add_argument(
"-c", "--cell", dest="cell_filename",
help="Unit cell filename")
parser.add_argument(
'--gv', action='store_true',
help='Calculate for gv_x_gv (tensor)')
parser.add_argument(
'--pqj', action='store_true',
help='Calculate for Pqj (scalar)')
parser.add_argument(
'--cv', action='store_true',
help='Calculate for Cv (scalar)')
parser.add_argument(
'--tau', action='store_true',
help='Calculate for lifetimes (scalar)')
parser.add_argument(
'--gamma', action='store_true',
help='Calculate for Gamma (scalar)')
parser.add_argument(
'--gruneisen', action='store_true',
help='Calculate for mode-Gruneisen parameters squared (scalar)')
parser.add_argument(
'--gv-norm', action='store_true',
help='Calculate for |g_v| (scalar)')
parser.add_argument(
'--mfp', action='store_true',
help='Mean free path is used instead of frequency')
parser.add_argument(
'--temperature', type=float, dest='temperature',
help='Temperature to output data at')
parser.add_argument(
'--nsp', '--num-sampling-points', type=int, dest='num_sampling_points',
default=100,
help="Number of sampling points in frequency or MFP axis")
parser.add_argument(
'--average', action='store_true',
help=("Output the traces of the tensors divided by 3 "
"rather than the unique elements"))
parser.add_argument(
'--trace', action='store_true',
help=("Output the traces of the tensors "
"rather than the unique elements"))
parser.add_argument(
'--smearing', action='store_true',
help='Use smearing method (only for scalar density)')
parser.add_argument(
'--qe', '--pwscf', dest="qe_mode",
action="store_true", help="Invoke Pwscf mode")
parser.add_argument(
'--crystal', dest="crystal_mode",
action="store_true", help="Invoke CRYSTAL mode")
parser.add_argument(
'--abinit', dest="abinit_mode",
action="store_true", help="Invoke Abinit mode")
parser.add_argument(
'--turbomole', dest="turbomole_mode",
action="store_true", help="Invoke TURBOMOLE mode")
parser.add_argument(
"--noks", "--no-kappa-stars",
dest="no_kappa_stars", action="store_true",
help="Deactivate summation of partial kappa at q-stars")
parser.add_argument('filenames', nargs='*')
args = parser.parse_args()
return args
def get_calculator(args):
"""Return calculator name."""
interface_mode = None
if args.qe_mode:
interface_mode = 'qe'
elif args.crystal_mode:
interface_mode = 'crystal'
elif args.abinit_mode:
interface_mode = 'abinit'
elif args.turbomole_mode:
interface_mode = 'turbomole'
return interface_mode
def read_files(args):
"""Read crystal structure and kappa.hdf5 files."""
interface_mode = get_calculator(args)
if len(args.filenames) > 1:
cell, _ = read_crystal_structure(args.filenames[0],
interface_mode=interface_mode)
f = h5py.File(args.filenames[1], 'r')
else:
cell, _ = read_crystal_structure(args.cell_filename,
interface_mode=interface_mode)
f = h5py.File(args.filenames[0], 'r')
return cell, f
def main():
"""Calculate kappa spectrum."""
args = get_parser()
cell, f_kappa = read_files(args)
primitive_matrix = np.reshape(
[fracval(x) for x in args.primitive_matrix.split()], (3, 3))
primitive = get_primitive(cell, primitive_matrix)
primitive_symmetry = Symmetry(primitive)
if 'mesh' in f_kappa:
mesh = np.array(f_kappa['mesh'][:], dtype='intc')
else:
mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc')
if 'temperature' in f_kappa:
temperatures = f_kappa['temperature'][:]
else:
temperatures = None
weights = f_kappa['weight'][:]
bz_grid = BZGrid(mesh,
lattice=primitive.cell,
symmetry_dataset=primitive_symmetry.dataset,
is_dense_gp_map=False)
if args.no_kappa_stars or (weights == 1).all():
grid_address = bz_grid.addresses
ir_grid_points = np.arange(np.prod(mesh), dtype='int_')
grid_mapping_table = np.arange(np.prod(mesh), dtype='int_')
else:
(ir_grid_points,
grid_address,
grid_mapping_table) = get_grid_symmetry(primitive,
primitive_symmetry,
mesh,
weights,
f_kappa['qpoint'][:])
################
# Set property #
################
if args.gv:
if 'gv_by_gv' in f_kappa:
gv_sum2 = f_kappa['gv_by_gv'][:]
else: # For backward compatibility. This will be removed someday.
gv = f_kappa['group_velocity'][:]
gv_sum2 = get_gv_by_gv(gv,
primitive_symmetry,
primitive,
mesh,
ir_grid_points,
grid_address)
# gv x gv is divied by primitive cell volume.
unit_conversion = primitive.volume
mode_prop = gv_sum2.reshape((1,) + gv_sum2.shape) / unit_conversion
else:
if 'mode_kappa' in f_kappa:
mode_prop = f_kappa['mode_kappa'][:]
else:
mode_prop = None
frequencies = f_kappa['frequency'][:]
conditions = frequencies > 0
if np.logical_not(conditions).sum() > 3:
sys.stderr.write("# Imaginary frequencies are found. "
"They are set to be zero.\n")
frequencies = np.where(conditions, frequencies, 0)
#######
# Run #
#######
if (args.gamma or args.gruneisen or args.pqj or
args.cv or args.tau or args.gv_norm):
if args.pqj:
mode_prop = f_kappa['ave_pp'][:].reshape((1,) + f_kappa['ave_pp'].shape)
elif args.cv:
mode_prop = f_kappa['heat_capacity'][:]
elif args.tau:
g = f_kappa['gamma'][:]
g = np.where(g > 0, g, -1)
mode_prop = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0) # tau
elif args.gv_norm:
mode_prop = np.sqrt(
(f_kappa['group_velocity'][:, :, :] ** 2).sum(axis=2))
mode_prop = mode_prop.reshape((1,) + mode_prop.shape)
elif args.gamma:
mode_prop = f_kappa['gamma'][:]
elif args.gruneisen:
mode_prop = f_kappa['gruneisen'][:].reshape((1,) + f_kappa['gruneisen'].shape)
mode_prop **= 2
else:
raise
if (args.temperature is not None and
not (args.gv_norm or args.pqj or args.gruneisen)):
temperatures, mode_prop = set_T_target(temperatures,
mode_prop,
args.temperature)
if args.smearing:
mode_prop_dos = GammaDOSsmearing(
mode_prop,
frequencies,
weights,
num_fpoints=args.num_sampling_points)
sampling_points, gdos = mode_prop_dos.get_gdos()
else:
mode_prop_dos = GammaDOStetrahedron(
mode_prop,
primitive,
frequencies,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
weights,
num_fpoints=args.num_sampling_points)
sampling_points, gdos = mode_prop_dos.get_gdos()
for i, gdos_t in enumerate(gdos):
total = np.dot(weights, mode_prop[i]).sum() / weights.sum()
assert np.isclose(gdos_t[-1][0], total)
show_scalar(gdos, temperatures, sampling_points, args)
else:
if args.mfp:
if 'mean_free_path' in f_kappa:
mfp = f_kappa['mean_free_path'][:]
mean_freepath = np.sqrt((mfp ** 2).sum(axis=3))
else:
mean_freepath = get_mfp(f_kappa['gamma'][:], f_kappa['group_velocity'][:])
if args.temperature is not None:
(temperatures,
mode_prop,
mean_freepath) = set_T_target(temperatures,
mode_prop,
args.temperature,
mean_freepath=mean_freepath)
kdos, sampling_points = run_mfp_dos(mean_freepath,
mode_prop,
primitive,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
args.num_sampling_points)
show_tensor(kdos, temperatures, sampling_points, args)
else:
if args.temperature is not None and not args.gv:
temperatures, mode_prop = set_T_target(temperatures,
mode_prop,
args.temperature)
kdos, sampling_points = run_prop_dos(frequencies,
mode_prop,
primitive,
mesh,
grid_address,
grid_mapping_table,
ir_grid_points,
args.num_sampling_points)
show_tensor(kdos, temperatures, sampling_points, args)
if __name__ == '__main__':
main()