From 19b590c2a0bccdeb1513b768ee5d52222c2da192 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 8 Aug 2016 15:10:39 +0900 Subject: [PATCH] Update gaccum and kaccum --- anharmonic/file_IO.py | 1 + scripts/gaccum | 42 ++++++++++++++++++++++++++++++++----- scripts/kaccum | 48 +++++++++++++++++++++++++++++++++++-------- 3 files changed, 77 insertions(+), 14 deletions(-) diff --git a/anharmonic/file_IO.py b/anharmonic/file_IO.py index e6563b63..ab61f3dc 100644 --- a/anharmonic/file_IO.py +++ b/anharmonic/file_IO.py @@ -653,6 +653,7 @@ def write_kappa_to_hdf5(temperature, filename=filename) with h5py.File("kappa" + suffix + ".hdf5", 'w') as w: w.create_dataset('temperature', data=temperature) + w.create_dataset('mesh', data=mesh) if frequency is not None: w.create_dataset('frequency', data=frequency) if group_velocity is not None: diff --git a/scripts/gaccum b/scripts/gaccum index b3eba883..9573a640 100755 --- a/scripts/gaccum +++ b/scripts/gaccum @@ -1,11 +1,14 @@ #!/usr/bin/env python +import math import numpy as np from phonopy.cui.settings import fracval from phonopy.phonon.tetrahedron_mesh import TetrahedronMesh from phonopy.phonon.dos import NormalDistribution from anharmonic.phonon3.triplets import get_ir_grid_points +epsilon = 1.0e-8 + class GammaDOS: def __init__(self, gamma, @@ -25,7 +28,7 @@ class GammaDOS: def _set_frequency_points(self): min_freq = np.min(self._frequencies) - max_freq = np.max(self._frequencies) + max_freq = np.max(self._frequencies) + epsilon self._frequency_points = np.linspace(min_freq, max_freq, self._num_fpoints) @@ -125,6 +128,12 @@ if __name__ == '__main__': help='Calculate Pqj instead of Gamma') parser.add_argument('--cv', action='store_true', help='Calculate Cv instead of Gamma') + parser.add_argument('--tau', action='store_true', + help='Calculate lifetimes (tau) instead of Gamma') + parser.add_argument('--gv_norm', action='store_true', + help='Calculate |g_v| instead of Gamma') + parser.add_argument('--temperature', type=float, dest='temperature', + help='Temperature to output data at') parser.add_argument('filenames', nargs='*') args = parser.parse_args() @@ -133,7 +142,10 @@ if __name__ == '__main__': [fracval(x) for x in args.primitive_axis.split()], (3, 3))) f = h5py.File(args.filenames[1]) - mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc') + if 'mesh' in f: + mesh = np.array(f['mesh'][:], dtype='intc') + else: + mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc') frequencies = f['frequency'][:] temperatures = f['temperature'][:] weights = f['weight'][:] @@ -145,16 +157,34 @@ if __name__ == '__main__': grid_address, grid_mapping_table) = get_ir_grid_points(mesh, rotations) - assert (weights == weights_for_check).all() + if (weights != weights_for_check).any(): + print("*** gaccum exited ***") + print("Something wrong in crystal symmetry. " + "Please check --pa or --mesh") + sys.exit(1) if args.pqj: gamma = f['ave_pp'][:].reshape((1,) + f['ave_pp'].shape) elif args.cv: gamma = f['heat_capacity'][:] - print("hoge") + elif args.tau: + g = f['gamma'][:] + g = np.where(g > 0, g, -1) + gamma = np.where(g > 0, 1.0 / (2.0 * g), 0) # tau + elif args.gv_norm: + gamma = np.sqrt((f['group_velocity'][:, :, :] ** 2).sum(axis=2)) + gamma = gamma.reshape((1,) + gamma.shape) else: gamma = f['gamma'][:] + if args.temperature is not None and gamma.shape[0] > 1: + for i, t in enumerate(temperatures): + if np.abs(t - args.temperature) < epsilon: + temperatures = temperatures[i:i+1] + gamma = gamma[i:i+1,:,:] + + break + gamma_dos = GammaDOStetrahedron(gamma, primitive, frequencies, @@ -177,7 +207,9 @@ if __name__ == '__main__': print("%f %e %e" % (f, g[0], g[1])) else: for i, gdos_t in enumerate(gdos): - print("# %d K %f" % (temperatures[i], gamma[i].sum(axis=0).sum(axis=0))) + total = np.dot(weights, gamma[i]).sum() / weights.sum() + assert np.isclose(gdos_t[-1][0], total) + print("# %d K" % temperatures[i]) for f, g in zip(freq_points, gdos_t): # show kappa_xx print("%f %f %f" % (f, g[0], g[1])) print('') diff --git a/scripts/kaccum b/scripts/kaccum index 35bf1571..bc1feb82 100755 --- a/scripts/kaccum +++ b/scripts/kaccum @@ -5,6 +5,7 @@ from phonopy.phonon.tetrahedron_mesh import TetrahedronMesh from phonopy.harmonic.force_constants import similarity_transformation from anharmonic.phonon3.triplets import get_ir_grid_points, get_grid_points_by_rotations +epsilon = 1.0e-8 def fracval(frac): if frac.find('/') == -1: @@ -33,7 +34,7 @@ class KappaDOS: ir_grid_points) min_freq = min(frequencies.ravel()) - max_freq = max(frequencies.ravel()) + max_freq = max(frequencies.ravel()) + epsilon self._frequency_points = np.linspace(min_freq, max_freq, 100) self._kdos = np.zeros( (len(mode_kappa), len(self._frequency_points), 2, 6), @@ -105,6 +106,14 @@ if __name__ == '__main__': help="Mesh numbers") parser.add_argument('--gv', action='store_true', help='Calculate gv_x_gv instead of kappa') + parser.add_argument('--temperature', type=float, dest='temperature', + help='Temperature to output data at') + 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('filenames', nargs='*') args = parser.parse_args() @@ -113,7 +122,10 @@ if __name__ == '__main__': [fracval(x) for x in args.primitive_axis.split()], (3, 3))) f = h5py.File(args.filenames[1]) - mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc') + if 'mesh' in f: + mesh = np.array(f['mesh'][:], dtype='intc') + else: + mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc') frequencies = f['frequency'][:] temperatures = f['temperature'][:] weights = f['weight'][:] @@ -125,12 +137,16 @@ if __name__ == '__main__': grid_address, grid_mapping_table) = get_ir_grid_points(mesh, rotations) - assert (weights == weights_for_check).all() + if (weights != weights_for_check).any(): + print("*** kaccum exited ***") + print("Something wrong in crystal symmetry. " + "Please check --pa or --mesh") + sys.exit(1) if args.gv: if 'gv_by_gv' in f: gv_sum2 = f['gv_by_gv'][:] - else: + else: # For backward compatibility. This will be removed someday. gv = f['group_velocity'][:] gv_sum2 = get_gv_by_gv(gv, symmetry, @@ -142,7 +158,13 @@ if __name__ == '__main__': mode_kappa = gv_sum2.reshape((1,) + gv_sum2.shape) / unit_conversion else: mode_kappa = f['mode_kappa'][:] - + + if args.temperature is not None and mode_kappa.shape[0] > 1: + for i, t in enumerate(temperatures): + if np.abs(t - args.temperature) < epsilon: + temperatures = temperatures[i:i+1] + mode_kappa = mode_kappa[i:i+1,:,:] + kappa_dos = KappaDOS(mode_kappa, primitive, frequencies, @@ -155,8 +177,16 @@ if __name__ == '__main__': for i, kdos_t in enumerate(kdos): if not args.gv: - print "# %d K" % temperatures[i], mode_kappa[i].sum(axis=0).sum(axis=0) + print("# %d K" % temperatures[i]) + for f, k in zip(freq_points, kdos_t): # show kappa_xx - print ("%f " * 13) % ((f,) + tuple(k[0]) + tuple(k[1])) - print - print + if args.average: + print(("%13.5f " * 3) % + (f, k[0][:3].sum() / 3, k[1][:3].sum() / 3)) + elif args.trace: + print(("%13.5f " * 3) % (f, k[0][:3].sum(), k[1][:3].sum())) + else: + print(("%f " * 13) % ((f,) + tuple(k[0]) + tuple(k[1]))) + + print('') + print('')