1#!/usr/bin/env python
2
3import sys
4import numpy as np
5from phonopy.phonon.tetrahedron_mesh import TetrahedronMesh
6from phonopy.harmonic.force_constants import similarity_transformation
7from phonopy.interface.calculator import read_crystal_structure
8from phono3py.phonon3.triplets import (get_ir_grid_points,
9                                       get_grid_points_by_rotations,
10                                       get_grid_address)
11from phonopy.phonon.dos import NormalDistribution
12
13epsilon = 1.0e-8
14
15
16def show_tensor(kdos, temperatures, sampling_points, args):
17    for i, kdos_t in enumerate(kdos):
18        if not args.gv:
19            print("# %d K" % temperatures[i])
20
21        for f, k in zip(sampling_points[i], kdos_t):  # show kappa_xx
22            if args.average:
23                print(("%13.5f " * 3) %
24                      (f, k[0][:3].sum() / 3, k[1][:3].sum() / 3))
25            elif args.trace:
26                print(("%13.5f " * 3) % (f, k[0][:3].sum(), k[1][:3].sum()))
27            else:
28                print(("%f " * 13) % ((f,) + tuple(k[0]) + tuple(k[1])))
29
30        print('')
31        print('')
32
33
34def show_scalar(gdos, temperatures, sampling_points, args):
35    if args.pqj or args.gruneisen or args.gv_norm:
36        for f, g in zip(sampling_points, gdos[0]):
37            print("%f %e %e" % (f, g[0], g[1]))
38    else:
39        for i, gdos_t in enumerate(gdos):
40            print("# %d K" % temperatures[i])
41            for f, g in zip(sampling_points, gdos_t):
42                print("%f %f %f" % (f, g[0], g[1]))
43            print('')
44            print('')
45
46
47def set_T_target(temperatures,
48                 mode_prop,
49                 T_target,
50                 mean_freepath=None):
51    for i, t in enumerate(temperatures):
52        if np.abs(t - args.temperature) < epsilon:
53            temperatures = temperatures[i:i+1]
54            mode_prop = mode_prop[i:i+1, :, :]
55            if mean_freepath is not None:
56                mean_freepath = mean_freepath[i:i+1]
57                return temperatures, mode_prop, mean_freepath
58            else:
59                return temperatures, mode_prop
60
61
62def run_prop_dos(frequencies,
63                 mode_prop,
64                 primitive,
65                 mesh,
66                 grid_address,
67                 grid_mapping_table,
68                 ir_grid_points,
69                 num_sampling_points):
70    kappa_dos = KappaDOS(mode_prop,
71                         primitive,
72                         frequencies,
73                         mesh,
74                         grid_address,
75                         grid_mapping_table,
76                         ir_grid_points,
77                         num_sampling_points=num_sampling_points)
78    freq_points, kdos = kappa_dos.get_kdos()
79    sampling_points = np.tile(freq_points, (len(kdos), 1))
80
81    return kdos, sampling_points
82
83
84def get_mfp(g, gv):
85    g = np.where(g > 0, g, -1)
86    gv_norm = np.sqrt((gv ** 2).sum(axis=2))
87    mean_freepath = np.where(g > 0, gv_norm / (2 * 2 * np.pi * g), 0)
88    return mean_freepath
89
90
91def run_mfp_dos(mean_freepath,
92                mode_prop,
93                primitive,
94                mesh,
95                grid_address,
96                grid_mapping_table,
97                ir_grid_points,
98                num_sampling_points):
99    kdos = []
100    sampling_points = []
101    for i, mfp in enumerate(mean_freepath):
102        kappa_dos = KappaDOS(mode_prop[i:i+1, :, :],
103                             primitive,
104                             mfp,
105                             mesh,
106                             grid_address,
107                             grid_mapping_table,
108                             ir_grid_points,
109                             num_sampling_points=num_sampling_points)
110        sampling_points_at_T, kdos_at_T = kappa_dos.get_kdos()
111        kdos.append(kdos_at_T[0])
112        sampling_points.append(sampling_points_at_T)
113    kdos = np.array(kdos)
114    sampling_points = np.array(sampling_points)
115
116    return kdos, sampling_points
117
118
119def fracval(frac):
120    if frac.find('/') == -1:
121        return float(frac)
122    else:
123        x = frac.split('/')
124        return float(x[0]) / float(x[1])
125
126
127def get_grid_symmetry(primitive, mesh, weights, qpoints):
128    symmetry = Symmetry(primitive)
129    rotations = symmetry.get_pointgroup_operations()
130    (ir_grid_points,
131     weights_for_check,
132     grid_address,
133     grid_mapping_table) = get_ir_grid_points(mesh, rotations)
134
135    try:
136        np.testing.assert_array_equal(weights, weights_for_check)
137    except:
138        print("*******************************")
139        print("** Might forget --pa option? **")
140        print("*******************************")
141        raise
142
143    qpoints_for_check = grid_address[ir_grid_points] / mesh.astype('double')
144    diff_q = qpoints - qpoints_for_check
145    np.testing.assert_almost_equal(diff_q, np.rint(diff_q))
146
147    return ir_grid_points, grid_address, grid_mapping_table
148
149
150def get_gv_by_gv(gv,
151                 symmetry,
152                 primitive,
153                 mesh,
154                 grid_points,
155                 grid_address):
156    point_operations = symmetry.get_reciprocal_operations()
157    rec_lat = np.linalg.inv(primitive.get_cell())
158    rotations_cartesian = np.array(
159        [similarity_transformation(rec_lat, r)
160         for r in point_operations], dtype='double')
161
162    num_band = gv.shape[1]
163    gv_sum2 = np.zeros((gv.shape[0], num_band, 6), dtype='double')
164    for i, gp in enumerate(grid_points):
165        rotation_map = get_grid_points_by_rotations(
166            grid_address[gp],
167            point_operations,
168            mesh)
169        gv_by_gv = np.zeros((num_band, 3, 3), dtype='double')
170        for r in rotations_cartesian:
171            gvs_rot = np.dot(gv[i], r.T)
172            gv_by_gv += [np.outer(r_gv, r_gv) for r_gv in gvs_rot]
173        gv_by_gv /= len(rotation_map) // len(np.unique(rotation_map))
174        for j, vxv in enumerate(
175                ([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])):
176            gv_sum2[i, :, j] = gv_by_gv[:, vxv[0], vxv[1]]
177
178    return gv_sum2
179
180
181class KappaDOS(object):
182    def __init__(self,
183                 mode_kappa,
184                 cell,
185                 frequencies,
186                 mesh,
187                 grid_address,
188                 grid_mapping_table,
189                 ir_grid_points,
190                 grid_order=None,
191                 num_sampling_points=100):
192        self._mode_kappa = mode_kappa
193        self._tetrahedron_mesh = TetrahedronMesh(
194            cell,
195            frequencies,
196            mesh,
197            grid_address,
198            grid_mapping_table,
199            ir_grid_points)
200
201        min_freq = min(frequencies.ravel())
202        max_freq = max(frequencies.ravel()) + epsilon
203        self._frequency_points = np.linspace(min_freq,
204                                             max_freq,
205                                             num_sampling_points)
206        self._kdos = np.zeros(
207            (len(mode_kappa), len(self._frequency_points), 2, 6),
208            dtype='double')
209        self._run_tetrahedron_method()
210
211    def get_kdos(self):
212        return self._frequency_points, self._kdos
213
214    def _run_tetrahedron_method(self):
215        thm = self._tetrahedron_mesh
216        for j, value in enumerate(('J', 'I')):
217            thm.set(value=value, frequency_points=self._frequency_points)
218            for i, iw in enumerate(thm):
219                # kdos[temp, freq_points, IJ, tensor_elem]
220                # iw[freq_points, band]
221                # mode_kappa[temp, ir_gp, band, tensor_elem]
222                self._kdos[:, :, j] += np.transpose(
223                    np.dot(iw, self._mode_kappa[:, i]), axes=(1, 0, 2))
224
225
226class GammaDOS(object):
227    def __init__(self,
228                 gamma,
229                 frequencies,
230                 ir_grid_weights,
231                 num_fpoints=200):
232        self._gamma = gamma
233        self._frequencies = frequencies
234        self._ir_grid_weights = ir_grid_weights
235        self._num_fpoints = num_fpoints
236        self._set_frequency_points()
237        self._gdos = np.zeros(
238            (len(gamma), len(self._frequency_points), 2), dtype='double')
239
240    def get_gdos(self):
241        return self._frequency_points, self._gdos
242
243    def _set_frequency_points(self):
244        min_freq = np.min(self._frequencies)
245        max_freq = np.max(self._frequencies) + epsilon
246        self._frequency_points = np.linspace(min_freq,
247                                             max_freq,
248                                             self._num_fpoints)
249
250
251class GammaDOSsmearing(GammaDOS):
252    def __init__(self,
253                 gamma,
254                 frequencies,
255                 ir_grid_weights,
256                 sigma=None,
257                 num_fpoints=200):
258        GammaDOS.__init__(self,
259                          gamma,
260                          frequencies,
261                          ir_grid_weights,
262                          num_fpoints=num_fpoints)
263        if sigma is None:
264            self._sigma = (max(self._frequency_points) -
265                           min(self._frequency_points)) / 100
266        else:
267            self._sigma = 0.1
268        self._smearing_function = NormalDistribution(self._sigma)
269        self._run_smearing_method()
270
271    def _run_smearing_method(self):
272        self._dos = []
273        num_gp = np.sum(self._ir_grid_weights)
274        for i, f in enumerate(self._frequency_points):
275            dos = self._smearing_function.calc(self._frequencies - f)
276            for j, g_t in enumerate(self._gamma):
277                self._gdos[j, i, 1] = np.sum(np.dot(self._ir_grid_weights,
278                                                    dos * g_t)) / num_gp
279
280
281class GammaDOStetrahedron(GammaDOS):
282    def __init__(self,
283                 gamma,
284                 cell,
285                 frequencies,
286                 mesh,
287                 grid_address,
288                 grid_mapping_table,
289                 ir_grid_points,
290                 ir_grid_weights,
291                 num_fpoints=200):
292        GammaDOS.__init__(self,
293                          gamma,
294                          frequencies,
295                          ir_grid_weights,
296                          num_fpoints=num_fpoints)
297        self._cell = cell
298        self._mesh = mesh
299        self._grid_address = grid_address
300        self._grid_mapping_table = grid_mapping_table
301        self._ir_grid_points = ir_grid_points
302
303        self._set_tetrahedron_method()
304        self._run_tetrahedron_method()
305
306    def _set_tetrahedron_method(self):
307        self._tetrahedron_mesh = TetrahedronMesh(
308            self._cell,
309            self._frequencies,
310            self._mesh,
311            self._grid_address,
312            self._grid_mapping_table,
313            self._ir_grid_points)
314
315    def _run_tetrahedron_method(self):
316        thm = self._tetrahedron_mesh
317        for j, value in enumerate(('J', 'I')):
318            thm.set(value=value, frequency_points=self._frequency_points)
319            for i, iw in enumerate(thm):
320                # gdos[temp, freq_points, IJ]
321                # iw[freq_points, band]
322                # gamma[temp, ir_gp, band]
323                self._gdos[:, :, j] += np.dot(
324                    self._gamma[:, i] * self._ir_grid_weights[i], iw.T)
325
326
327if __name__ == '__main__':
328    """Incremental kappa with respect to frequency and the derivative"""
329
330    import h5py
331    from phonopy.structure.cells import get_primitive
332    from phonopy.structure.symmetry import Symmetry
333    import argparse
334
335    parser = argparse.ArgumentParser(description="Show unit cell volume")
336    parser.add_argument(
337        "--pa", dest="primitive_matrix", default="1 0 0 0 1 0 0 0 1",
338        help="Primitive matrix")
339    parser.add_argument(
340        "--mesh", dest="mesh", default="1 1 1",
341        help="Mesh numbers")
342    parser.add_argument(
343        "-c", "--cell", dest="cell_filename",
344        help="Unit cell filename")
345    parser.add_argument(
346        '--gv', action='store_true',
347        help='Calculate for gv_x_gv (tensor)')
348    parser.add_argument(
349        '--pqj', action='store_true',
350        help='Calculate for Pqj (scalar)')
351    parser.add_argument(
352        '--cv', action='store_true',
353        help='Calculate for Cv (scalar)')
354    parser.add_argument(
355        '--tau', action='store_true',
356        help='Calculate for lifetimes (scalar)')
357    parser.add_argument(
358        '--gamma', action='store_true',
359        help='Calculate for Gamma (scalar)')
360    parser.add_argument(
361        '--gruneisen', action='store_true',
362        help='Calculate for mode-Gruneisen parameters squared (scalar)')
363    parser.add_argument(
364        '--gv-norm', action='store_true',
365        help='Calculate for |g_v| (scalar)')
366    parser.add_argument(
367        '--mfp', action='store_true',
368        help='Mean free path is used instead of frequency')
369    parser.add_argument(
370        '--temperature', type=float, dest='temperature',
371        help='Temperature to output data at')
372    parser.add_argument(
373        '--nsp', '--num-sampling-points', type=int, dest='num_sampling_points',
374        default=100,
375        help="Number of sampling points in frequency or MFP axis")
376    parser.add_argument(
377        '--average', action='store_true',
378        help=("Output the traces of the tensors divided by 3 "
379              "rather than the unique elements"))
380    parser.add_argument(
381        '--trace', action='store_true',
382        help=("Output the traces of the tensors "
383              "rather than the unique elements"))
384    parser.add_argument(
385        '--smearing', action='store_true',
386        help='Use smearing method (only for scalar density)')
387    parser.add_argument(
388        '--qe', '--pwscf', dest="qe_mode",
389        action="store_true", help="Invoke Pwscf mode")
390    parser.add_argument(
391        '--crystal', dest="crystal_mode",
392        action="store_true", help="Invoke CRYSTAL mode")
393    parser.add_argument(
394        '--abinit', dest="abinit_mode",
395        action="store_true", help="Invoke Abinit mode")
396    parser.add_argument(
397        '--turbomole', dest="turbomole_mode",
398        action="store_true", help="Invoke TURBOMOLE mode")
399    parser.add_argument(
400        "--noks", "--no-kappa-stars",
401        dest="no_kappa_stars", action="store_true",
402        help="Deactivate summation of partial kappa at q-stars")
403    parser.add_argument('filenames', nargs='*')
404    args = parser.parse_args()
405
406    interface_mode = None
407    if args.qe_mode:
408        interface_mode = 'qe'
409    elif args.crystal_mode:
410        interface_mode = 'crystal'
411    elif args.abinit_mode:
412        interface_mode = 'abinit'
413    elif args.turbomole_mode:
414        interface_mode = 'turbomole'
415    if len(args.filenames) > 1:
416        cell, _ = read_crystal_structure(args.filenames[0],
417                                         interface_mode=interface_mode)
418        f = h5py.File(args.filenames[1], 'r')
419    else:
420        cell, _ = read_crystal_structure(args.cell_filename,
421                                         interface_mode=interface_mode)
422        f = h5py.File(args.filenames[0], 'r')
423
424    primitive_matrix = np.reshape(
425        [fracval(x) for x in args.primitive_matrix.split()], (3, 3))
426    primitive = get_primitive(cell, primitive_matrix)
427
428    if 'mesh' in f:
429        mesh = np.array(f['mesh'][:], dtype='intc')
430    else:
431        mesh = np.array([int(x) for x in args.mesh.split()], dtype='intc')
432
433    if 'temperature' in f:
434        temperatures = f['temperature'][:]
435    else:
436        temperatures = None
437    weights = f['weight'][:]
438
439    if args.no_kappa_stars or (weights == 1).all():
440        grid_address = get_grid_address(mesh)
441        ir_grid_points = np.arange(np.prod(mesh), dtype='uintp')
442        grid_mapping_table = np.arange(np.prod(mesh), dtype='uintp')
443    else:
444        (ir_grid_points,
445         grid_address,
446         grid_mapping_table) = get_grid_symmetry(primitive,
447                                                 mesh,
448                                                 weights,
449                                                 f['qpoint'][:])
450
451    ################
452    # Set property #
453    ################
454    if args.gv:
455        if 'gv_by_gv' in f:
456            gv_sum2 = f['gv_by_gv'][:]
457        else:  # For backward compatibility. This will be removed someday.
458            gv = f['group_velocity'][:]
459            symmetry = Symmetry(primitive)
460            gv_sum2 = get_gv_by_gv(gv,
461                                   symmetry,
462                                   primitive,
463                                   mesh,
464                                   ir_grid_points,
465                                   grid_address)
466
467        # gv x gv is divied by primitive cell volume.
468        unit_conversion = primitive.get_volume()
469        mode_prop = gv_sum2.reshape((1,) + gv_sum2.shape) / unit_conversion
470    else:
471        if 'mode_kappa' in f:
472            mode_prop = f['mode_kappa'][:]
473        else:
474            mode_prop = None
475
476    frequencies = f['frequency'][:]
477    conditions = frequencies > 0
478    if np.logical_not(conditions).sum() > 3:
479        sys.stderr.write("# Imaginary frequencies are found. "
480                         "They are set to be zero.\n")
481        frequencies = np.where(conditions, frequencies, 0)
482
483    #######
484    # Run #
485    #######
486    if (args.gamma or args.gruneisen or args.pqj or
487        args.cv or args.tau or args.gv_norm):
488        if args.pqj:
489            mode_prop = f['ave_pp'][:].reshape((1,) + f['ave_pp'].shape)
490        elif args.cv:
491            mode_prop = f['heat_capacity'][:]
492        elif args.tau:
493            g = f['gamma'][:]
494            g = np.where(g > 0, g, -1)
495            mode_prop = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)  # tau
496        elif args.gv_norm:
497            mode_prop = np.sqrt(
498                (f['group_velocity'][:, :, :] ** 2).sum(axis=2))
499            mode_prop = mode_prop.reshape((1,) + mode_prop.shape)
500        elif args.gamma:
501            mode_prop = f['gamma'][:]
502        elif args.gruneisen:
503            mode_prop = f['gruneisen'][:].reshape((1,) + f['gruneisen'].shape)
504            mode_prop **= 2
505        else:
506            raise
507        if (args.temperature is not None and
508            not (args.gv_norm or args.pqj or args.gruneisen)):
509            temperatures, mode_prop = set_T_target(temperatures,
510                                                   mode_prop,
511                                                   args.temperature)
512        if args.smearing:
513            mode_prop_dos = GammaDOSsmearing(
514                mode_prop,
515                frequencies,
516                weights,
517                num_fpoints=args.num_sampling_points)
518            sampling_points, gdos = mode_prop_dos.get_gdos()
519        else:
520            mode_prop_dos = GammaDOStetrahedron(
521                mode_prop,
522                primitive,
523                frequencies,
524                mesh,
525                grid_address,
526                grid_mapping_table,
527                ir_grid_points,
528                weights,
529                num_fpoints=args.num_sampling_points)
530            sampling_points, gdos = mode_prop_dos.get_gdos()
531            for i, gdos_t in enumerate(gdos):
532                total = np.dot(weights, mode_prop[i]).sum() / weights.sum()
533                assert np.isclose(gdos_t[-1][0], total)
534        show_scalar(gdos, temperatures, sampling_points, args)
535    else:
536        if args.mfp:
537            if 'mean_free_path' in f:
538                mfp = f['mean_free_path'][:]
539                mean_freepath = np.sqrt((mfp ** 2).sum(axis=3))
540            else:
541                mean_freepath = get_mfp(f['gamma'][:], f['group_velocity'][:])
542            if args.temperature is not None:
543                (temperatures,
544                 mode_prop,
545                 mean_freepath) = set_T_target(temperatures,
546                                               mode_prop,
547                                               args.temperature,
548                                               mean_freepath=mean_freepath)
549
550            kdos, sampling_points = run_mfp_dos(mean_freepath,
551                                                mode_prop,
552                                                primitive,
553                                                mesh,
554                                                grid_address,
555                                                grid_mapping_table,
556                                                ir_grid_points,
557                                                args.num_sampling_points)
558            show_tensor(kdos, temperatures, sampling_points, args)
559        else:
560            if args.temperature is not None and not args.gv:
561                temperatures, mode_prop = set_T_target(temperatures,
562                                                       mode_prop,
563                                                       args.temperature)
564            kdos, sampling_points = run_prop_dos(frequencies,
565                                                 mode_prop,
566                                                 primitive,
567                                                 mesh,
568                                                 grid_address,
569                                                 grid_mapping_table,
570                                                 ir_grid_points,
571                                                 args.num_sampling_points)
572            show_tensor(kdos, temperatures, sampling_points, args)
573