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

Update gaccum and kaccum

This commit is contained in:
Atsushi Togo 2016-08-08 15:10:39 +09:00
parent a394674b94
commit 19b590c2a0
3 changed files with 77 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
anharmonic/file_IO.py
View File

@ -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:

42
scripts/gaccum
View File

@ -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('')

48
scripts/kaccum
View File

@ -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('')