Array shape of gammas returned by get_imag_self_energy was changed.

c/anharmonic_h/phonon3_h/real_self_energy.h

@@ -32,8 +32,8 @@
 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
 /* POSSIBILITY OF SUCH DAMAGE. */
 
-#ifndef __frequency_shift3_H__
+#ifndef __real_self_energy_H__
-#define __frequency_shift3_H__
+#define __real_self_energy_H__
 
 #include <phonoc_array.h>
 

phono3py/api_phono3py.py

@@ -869,18 +869,10 @@ class Phono3py(object):
     def get_phph_interaction(self):
         return self.phph_interaction
 
-    @property
-    def gammas(self):
-        return self._gammas
-
     @property
     def detailed_gammas(self):
         return self._detailed_gammas
 
-    @property
-    def frequency_points(self):
-        return self._frequency_points
-
     def init_phph_interaction(self,
                               nac_q_direction=None,
                               constant_averaged_interaction=None,
@@ -1451,6 +1443,8 @@ class Phono3py(object):
         else:
             self._frequency_points, self._gammas = vals
 
+        return vals
+
     def write_imag_self_energy(self, filename=None):
         write_imag_self_energy(
             self._gammas,
@@ -1461,8 +1455,9 @@ class Phono3py(object):
             self._temperatures,
             self._sigmas,
             scattering_event_class=self._scattering_event_class,
-            filename=filename,
+            output_filename=filename,
-            is_mesh_symmetry=self._is_mesh_symmetry)
+            is_mesh_symmetry=self._is_mesh_symmetry,
+            log_level=self._log_level)
 
     def run_real_self_energy(
             self,
@@ -1604,6 +1599,23 @@ class Phono3py(object):
             num_frequency_points=num_frequency_points,
             temperatures=temperatures,
             log_level=self._log_level)
+        self._spectral_function.run()
+
+        # if write_hdf5:
+        #     filename = write_spectral_function_to_hdf5(
+        #         gp,
+        #         band_indices,
+        #         _temperatures,
+        #         all_deltas[i, j],
+        #         mesh,
+        #         frequency_points=_frequency_points,
+        #         frequencies=frequencies,
+        #         filename=output_filename)
+
+        #     if self._log_level:
+        #         print("Spectral functions were stored in \"%s\"." % filename)
+        #         sys.stdout.flush()
+
 
     def run_thermal_conductivity(
             self,

phono3py/cui/phono3py_argparse.py

@@ -172,9 +172,6 @@ def get_parser(fc_symmetry=False,
         "--fs2f2", "--force-sets-to-forces-fc2",
         dest="force_sets_to_forces_fc2_mode", default=False,
         action="store_true", help="Create FORCES_FC2 from FORCE_SETS")
-    parser.add_argument(
-        "--fst", "--frequency-shift", dest="is_frequency_shift", default=False,
-        action="store_true", help="Calculate frequency shifts")
     parser.add_argument(
         "--full-pp", dest="is_full_pp", action="store_true", default=False,
         help=("Calculate full ph-ph interaction for RTA conductivity."
@@ -327,6 +324,10 @@ def get_parser(fc_symmetry=False,
         "--reducible-colmat", dest="is_reducible_collision_matrix",
         action="store_true", default=False,
         help="Solve reducible collision matrix")
+    parser.add_argument(
+        "--rse", dest="is_real_self_energy", action="store_true",
+        default=False,
+        help="Calculate real part of self energy")
     parser.add_argument(
         "--scattering-event-class", dest="scattering_event_class", type=int,
         default=None,
@@ -339,6 +340,10 @@ def get_parser(fc_symmetry=False,
     parser.add_argument(
         "--sigma-cutoff", dest="sigma_cutoff_width", type=float, default=None,
         help="Cutoff width of smearing function (ratio to sigma value)")
+    parser.add_argument(
+        "--spf", dest="is_spectral_function", action="store_true",
+        default=False,
+        help="Calculate spectral function")
     parser.add_argument(
         "--stp", "--show-num-triplets", dest="show_num_triplets",
         action="store_true", default=False,

phono3py/cui/phono3py_script.py

@@ -1108,7 +1108,7 @@ def main(**argparse_control):
     #####################################################
     # Run frequency shift calculation of bubble diagram #
     #####################################################
-    if run_mode == "real_self_energy":
+    elif run_mode == "real_self_energy":
         phono3py.run_real_self_energy(
             updated_settings['grid_points'],
             updated_settings['temperature_points'],
@@ -1119,7 +1119,7 @@ def main(**argparse_control):
     #######################################################
     # Run spectral function calculation of bubble diagram #
     #######################################################
-    if run_mode == "spectral_function":
+    elif run_mode == "spectral_function":
         phono3py.run_spectral_function(
             updated_settings['grid_points'],
             updated_settings['temperature_points'],
@@ -1130,7 +1130,7 @@ def main(**argparse_control):
     ####################################
     # Run lattice thermal conductivity #
     ####################################
-    if run_mode == "conductivity-RTA" or run_mode == "conductivity-LBTE":
+    elif run_mode == "conductivity-RTA" or run_mode == "conductivity-LBTE":
         phono3py.run_thermal_conductivity(
             is_LBTE=settings.is_lbte,
             temperatures=updated_settings['temperatures'],

phono3py/cui/settings.py

@@ -53,7 +53,6 @@ class Phono3pySettings(Settings):
         'ion_clamped': False,
         'is_bterta': False,
         'is_compact_fc': False,
-        'is_frequency_shift': False,
         'is_full_pp': False,
         'is_gruneisen': False,
         'is_imag_self_energy': False,
@@ -62,7 +61,9 @@ class Phono3pySettings(Settings):
         'is_kappa_star': True,
         'is_lbte': False,
         'is_N_U': False,
+        'is_real_self_energy': False,
         'is_reducible_collision_matrix': False,
+        'is_spectral_function': False,
         'is_symmetrize_fc2': False,
         'is_symmetrize_fc3_q': False,
         'is_symmetrize_fc3_r': False,
@@ -141,9 +142,6 @@ class Phono3pySettings(Settings):
     def set_is_compact_fc(self, val):
         self._v['is_compact_fc'] = val
 
-    def set_is_frequency_shift(self, val):
-        self._v['is_frequency_shift'] = val
-
     def set_is_full_pp(self, val):
         self._v['is_full_pp'] = val
 
@@ -168,9 +166,15 @@ class Phono3pySettings(Settings):
     def set_is_N_U(self, val):
         self._v['is_N_U'] = val
 
+    def set_is_real_self_energy(self, val):
+        self._v['is_real_self_energy'] = val
+
     def set_is_reducible_collision_matrix(self, val):
         self._v['is_reducible_collision_matrix'] = val
 
+    def set_is_spectral_function(self, val):
+        self._v['is_spectral_function'] = val
+
     def set_is_symmetrize_fc2(self, val):
         self._v['is_symmetrize_fc2'] = val
 
@@ -343,10 +347,6 @@ class Phono3pyConfParser(ConfParser):
             if self._args.is_gruneisen:
                 self._confs['gruneisen'] = '.true.'
 
-        if 'is_frequency_shift' in self._args:
-            if self._args.is_frequency_shift:
-                self._confs['frequency_shift'] = '.true.'
-
         if 'is_full_pp' in self._args:
             if self._args.is_full_pp:
                 self._confs['full_pp'] = '.true.'
@@ -375,10 +375,18 @@ class Phono3pyConfParser(ConfParser):
             if self._args.is_N_U:
                 self._confs['N_U'] = '.true.'
 
+        if 'is_real_self_energy' in self._args:
+            if self._args.is_real_self_energy:
+                self._confs['real_self_energy'] = '.true.'
+
         if 'is_reducible_collision_matrix' in self._args:
             if self._args.is_reducible_collision_matrix:
                 self._confs['reducible_collision_matrix'] = '.true.'
 
+        if 'is_spectral_function' in self._args:
+            if self._args.is_spectral_function:
+                self._confs['spectral_function'] = '.true.'
+
         if 'is_symmetrize_fc2' in self._args:
             if self._args.is_symmetrize_fc2:
                 self._confs['symmetrize_fc2'] = '.true.'
@@ -504,9 +512,9 @@ class Phono3pyConfParser(ConfParser):
                     'write_gamma_detail', 'write_gamma',
                     'write_collision', 'write_phonon', 'write_pp',
                     'write_LBTE_solution', 'full_pp', 'ion_clamped',
-                    'bterta', 'compact_fc', 'frequency_shift',
+                    'bterta', 'compact_fc', 'real_self_energy',
                     'gruneisen', 'imag_self_energy', 'isotope',
-                    'joint_dos', 'lbte', 'N_U',
+                    'joint_dos', 'lbte', 'N_U', 'spectral_function',
                     'reducible_collision_matrix', 'symmetrize_fc2',
                     'symmetrize_fc3_q', 'symmetrize_fc3_r', 'kappa_star'):
                 if confs[conf_key].lower() == '.true.':
@@ -679,10 +687,6 @@ class Phono3pyConfParser(ConfParser):
         if 'compact_fc' in params:
             self._settings.set_is_compact_fc(params['compact_fc'])
 
-        # Calculate frequency_shifts
-        if 'frequency_shift' in params:
-            self._settings.set_is_frequency_shift(params['frequency_shift'])
-
         # Calculate full ph-ph interaction strength for RTA conductivity
         if 'full_pp' in params:
             self._settings.set_is_full_pp(params['full_pp'])
@@ -763,6 +767,10 @@ class Phono3pyConfParser(ConfParser):
             self._settings.set_pp_conversion_factor(
                 params['pp_conversion_factor'])
 
+        # Calculate real_self_energys
+        if 'real_self_energy' in params:
+            self._settings.set_is_real_self_energy(params['real_self_energy'])
+
         # Read phonon-phonon interaction amplitudes from hdf5
         if 'read_amplitude' in params:
             self._settings.set_read_amplitude(params['read_amplitude'])
@@ -800,6 +808,10 @@ class Phono3pyConfParser(ConfParser):
         if 'sigma_cutoff_width' in params:
             self._settings.set_sigma_cutoff_width(params['sigma_cutoff_width'])
 
+        # Calculate spectral_functions
+        if 'spectral_function' in params:
+            self._settings.set_is_spectral_function(params['spectral_function'])
+
         # Subtract residual forces to create FORCES_FC2 and FORCES_FC3
         if 'subtract_forces' in params:
             self._settings.set_subtract_forces(params['subtract_forces'])

phono3py/cui/show_log.py

@@ -178,7 +178,10 @@ def show_phono3py_settings(phono3py, settings, updated_settings, log_level):
 
     if settings.is_lbte and settings.read_collision is not None:
         pass
-    elif settings.is_frequency_shift or settings.is_bterta:
+    elif (settings.is_real_self_energy or
+          settings.is_imag_self_energy or
+          settings.is_spectral_function or
+          settings.is_bterta):
         if len(temperatures) > 5:
             text = (" %.1f " * 5 + "...") % tuple(temperatures[:5])
             text += " %.1f" % temperatures[-1]
@@ -222,7 +225,10 @@ def show_phono3py_settings(phono3py, settings, updated_settings, log_level):
     if frequency_scale_factor is not None:
         print("Frequency scale factor: %8.5f" % frequency_scale_factor)
 
-    if settings.is_joint_dos or settings.is_imag_self_energy:
+    if (settings.is_joint_dos or
+        settings.is_imag_self_energy or
+        settings.is_real_self_energy or
+        settings.is_spectral_function):
         if frequency_step is not None:
             print("Frequency step for spectrum: %s" % frequency_step)
         if num_frequency_points is not None:

phono3py/file_IO.py

@@ -419,6 +419,8 @@ def write_imag_self_energy_at_grid_point(gp,
         w.write("%15.7f %20.15e\n" % (freq, g))
     w.close()
 
+    return gammas_filename
+
 
 def write_joint_dos(gp,
                     mesh,
@@ -478,34 +480,6 @@ def _write_joint_dos_at_t(grid_point,
         return jdos_filename
 
 
-def write_linewidth_at_grid_point(gp,
-                                  band_indices,
-                                  temperatures,
-                                  gamma,
-                                  mesh,
-                                  sigma=None,
-                                  filename=None,
-                                  is_mesh_symmetry=True):
-
-    lw_filename = "linewidth"
-    lw_filename += "-m%d%d%d-g%d-" % (mesh[0], mesh[1], mesh[2], gp)
-    if sigma is not None:
-        lw_filename += ("s%f" % sigma).rstrip('0') + "-"
-
-    for i in band_indices:
-        lw_filename += "b%d" % (i + 1)
-
-    if filename is not None:
-        lw_filename += ".%s" % filename
-    elif not is_mesh_symmetry:
-        lw_filename += ".nosym"
-    lw_filename += ".dat"
-
-    with open(lw_filename, 'w') as w:
-        for v, t in zip(gamma.sum(axis=1) * 2 / gamma.shape[1], temperatures):
-            w.write("%15.7f %20.15e\n" % (t, v))
-
-
 def write_real_self_energy(gp,
                            band_indices,
                            frequency_points,
@@ -516,7 +490,7 @@ def write_real_self_energy(gp,
                            filename=None,
                            is_mesh_symmetry=True):
 
-    fst_filename = "real_self_energy"
+    fst_filename = "deltas"
     suffix = _get_filename_suffix(mesh,
                                   grid_point=gp,
                                   band_indices=band_indices,
@@ -539,19 +513,20 @@ def write_real_self_energy(gp,
     return fst_filename
 
 
-def write_Delta_to_hdf5(grid_point,
+def write_real_self_energy_to_hdf5(grid_point,
-                        band_indices,
+                                   band_indices,
-                        temperatures,
+                                   temperatures,
-                        deltas,
+                                   deltas,
-                        mesh,
+                                   mesh,
-                        epsilon,
+                                   epsilon,
-                        frequency_points=None,
+                                   frequency_points=None,
-                        frequencies=None,
+                                   frequencies=None,
-                        filename=None):
+                                   filename=None):
-    """Wirte frequency shifts (currently only bubble) in hdf5
+    """Wirte real part of self energy (currently only bubble) in hdf5
 
     deltas : ndarray
-        Frequency shifts.
+        Real part of self energy.
+
         With frequency_points:
             shape=(temperature, frequency_points, band_index),
             dtype='double', order='C'
@@ -559,7 +534,7 @@ def write_Delta_to_hdf5(grid_point,
             shape=(temperature, band_index), dtype='double', order='C'
 
     """
-    full_filename = "Delta"
+    full_filename = "deltas"
     suffix = _get_filename_suffix(mesh, grid_point=grid_point)
     _band_indices = np.array(band_indices, dtype='intc')
 
@@ -585,6 +560,43 @@ def write_Delta_to_hdf5(grid_point,
     return full_filename
 
 
+def write_spectral_function_to_hdf5(grid_point,
+                                    band_indices,
+                                    temperatures,
+                                    spectral_functions,
+                                    mesh,
+                                    frequency_points=None,
+                                    frequencies=None,
+                                    filename=None):
+    """Wirte spectral functions (currently only bubble) in hdf5
+
+    spectral_functions : ndarray
+        Spectral functions.
+        shape=(temperature, band_index, frequency_points),
+        dtype='double', order='C'
+
+    """
+    full_filename = "spectral"
+    suffix = _get_filename_suffix(mesh, grid_point=grid_point)
+    _band_indices = np.array(band_indices, dtype='intc')
+
+    full_filename += suffix
+    full_filename += ".hdf5"
+
+    with h5py.File(full_filename, 'w') as w:
+        w.create_dataset('grid_point', data=grid_point)
+        w.create_dataset('mesh', data=mesh)
+        w.create_dataset('band_index', data=_band_indices)
+        w.create_dataset('delta', data=spectral_functions)
+        w.create_dataset('temperature', data=temperatures)
+        if frequency_points is not None:
+            w.create_dataset('frequency_points', data=frequency_points)
+        if frequencies is not None:
+            w.create_dataset('frequency', data=frequencies)
+
+    return full_filename
+
+
 def write_collision_to_hdf5(temperature,
                             mesh,
                             gamma=None,

phono3py/phonon3/imag_self_energy.py

@@ -131,7 +131,7 @@ def get_imag_self_energy(interaction,
         interaction, with_detail=(write_gamma_detail or return_gamma_detail))
     gamma = np.zeros(
         (len(grid_points), len(_sigmas), len(temperatures),
-         len(_frequency_points), len(interaction.band_indices)),
+         len(interaction.band_indices), len(_frequency_points)),
         dtype='double', order='C')
     detailed_gamma = []
 
@@ -179,6 +179,7 @@ def get_imag_self_energy(interaction,
 
             ise.set_sigma(sigma)
 
+            # Run one by one at frequency points
             for k, freq_point in enumerate(_frequency_points):
                 ise.set_frequency_points([freq_point])
                 ise.set_integration_weights(
@@ -187,7 +188,7 @@ def get_imag_self_energy(interaction,
                 for l, t in enumerate(temperatures):
                     ise.set_temperature(t)
                     ise.run()
-                    gamma[i, j, l, k] = ise.get_imag_self_energy()[0]
+                    gamma[i, j, l, :, k] = ise.get_imag_self_energy()[0]
                     if write_gamma_detail or return_gamma_detail:
                         detailed_gamma_at_gp[j, l, k] = (
                             ise.get_detailed_imag_self_energy()[0])
@@ -271,8 +272,9 @@ def write_imag_self_energy(imag_self_energy,
                            temperatures,
                            sigmas,
                            scattering_event_class=None,
-                           filename=None,
+                           output_filename=None,
-                           is_mesh_symmetry=True):
+                           is_mesh_symmetry=True,
+                           log_level=0):
     for gp, ise_sigmas in zip(grid_points, imag_self_energy):
         for sigma, ise_temps in zip(sigmas, ise_sigmas):
             for t, ise in zip(temperatures, ise_temps):
@@ -280,17 +282,20 @@ def write_imag_self_energy(imag_self_energy,
                     pos = 0
                     for j in range(i):
                         pos += len(band_indices[j])
-                    write_imag_self_energy_at_grid_point(
+                    filename = write_imag_self_energy_at_grid_point(
                         gp,
                         bi,
                         mesh,
                         frequency_points,
-                        ise[:, pos:(pos + len(bi))].sum(axis=1) / len(bi),
+                        ise[pos:(pos + len(bi))].sum(axis=0) / len(bi),
                         sigma=sigma,
                         temperature=t,
                         scattering_event_class=scattering_event_class,
-                        filename=filename,
+                        filename=output_filename,
                         is_mesh_symmetry=is_mesh_symmetry)
+                    if log_level:
+                        print("Imaginary parts of self-energies were "
+                              "written to \"%s\"." % filename)
 
 
 def average_by_degeneracy(imag_self_energy, band_indices, freqs_at_gp):

phono3py/phonon3/real_self_energy.py

@@ -37,7 +37,8 @@ import numpy as np
 from phonopy.units import Hbar, EV, THz
 from phonopy.phonon.degeneracy import degenerate_sets
 from phono3py.phonon.func import bose_einstein
-from phono3py.file_IO import write_real_self_energy, write_Delta_to_hdf5
+from phono3py.file_IO import (
+    write_real_self_energy, write_real_self_energy_to_hdf5)
 from phono3py.phonon3.imag_self_energy import get_frequency_points
 
 
@@ -135,7 +136,7 @@ def get_real_self_energy(interaction,
                         sys.stdout.flush()
 
             if write_hdf5:
-                filename = write_Delta_to_hdf5(
+                filename = write_real_self_energy_to_hdf5(
                     gp,
                     band_indices,
                     _temperatures,

phono3py/phonon3/spectral_function.py

@@ -62,9 +62,17 @@ class SpectralFunction(object):
         self._frequency_points = None
 
     def run(self):
+        if self._log_level:
+            print("-" * 28 + " Spectral function " + "-" * 28)
         self._run_gamma()
+        if self._log_level:
+            print("-" * 11 +
+                  " Real part of self energy by Kramers-Kronig relation "
+                  + "-" * 11)
         self._run_delta()
         self._run_spectral_function()
+        if self._log_level:
+            print("-" * 75)
 
     @property
     def spectral_functions(self):
@@ -87,10 +95,7 @@ class SpectralFunction(object):
         return self._grid_points
 
     def _run_gamma(self):
-        if self._log_level > 0:
+        # gammas[grid_points, sigmas, temps, band_indices, freq_points]
-            print("Start calculation of imaginary part of self energies")
-
-        # gammas[grid_points, sigmas, temps, freq_points, band_indices]
         fpoints, gammas = get_imag_self_energy(
             self._interaction,
             self._grid_points,
@@ -108,8 +113,8 @@ class SpectralFunction(object):
             for j, temp in enumerate(self._temperatures):
                 for k, bi in enumerate(self._interaction.band_indices):
                     re_part, fpoints = imag_to_real(
-                        self._gammas[i, j, :, k], self._frequency_points)
+                        self._gammas[i, j, k], self._frequency_points)
-                    self._deltas[i, j, :, k] = -re_part
+                    self._deltas[i, j, k] = -re_part
         assert (np.abs(self._frequency_points - fpoints) < 1e-8).all()
 
     def _run_spectral_function(self):
@@ -119,10 +124,10 @@ class SpectralFunction(object):
             for j, temp in enumerate(self._temperatures):
                 for k, bi in enumerate(self._interaction.band_indices):
                     freq = frequencies[gp, bi]
-                    gammas = self._gammas[i, j, :, k]
+                    gammas = self._gammas[i, j, k]
-                    deltas = self._deltas[i, j, :, k]
+                    deltas = self._deltas[i, j, k]
                     vals = self._get_spectral_function(gammas, deltas, freq)
-                    self._spectral_functions[i, j, :, k] = vals
+                    self._spectral_functions[i, j, k] = vals
 
     def _get_spectral_function(self, gammas, deltas, freq):
         fpoints = self._frequency_points

test/phonon3/test_imag_self_energy.py

@@ -72,48 +72,48 @@ detailed_gamma = [0.00000000, 0.00653193, 0.02492913, 0.01682092, 0.01001680,
 def test_imag_self_energy_npoints(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    si_pbesol.run_imag_self_energy(
+    _fpoints, _gammas = si_pbesol.run_imag_self_energy(
         [1, 103],
         [300, ],
         num_frequency_points=10)
     np.testing.assert_allclose(
-        gammas, si_pbesol.gammas.ravel(), atol=1e-2)
+        gammas, np.swapaxes(_gammas, -1, -2).ravel(), atol=1e-2)
     np.testing.assert_allclose(
-        freq_points, si_pbesol.frequency_points.ravel(),
+        freq_points, _fpoints.ravel(),
         atol=1e-5)
 
 
 def test_imag_self_energy_freq_points(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    si_pbesol.run_imag_self_energy(
+    _fpoints, _gammas = si_pbesol.run_imag_self_energy(
         [1, 103],
         [300, ],
         frequency_points=freq_points)
     np.testing.assert_allclose(
-        gammas, si_pbesol.gammas.ravel(), atol=1e-2)
+        gammas, np.swapaxes(_gammas, -1, -2).ravel(), atol=1e-2)
     np.testing.assert_allclose(
-        freq_points, si_pbesol.frequency_points.ravel(), atol=1e-5)
+        freq_points, _fpoints.ravel(), atol=1e-5)
 
 
 def test_imag_self_energy_detailed(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    si_pbesol.run_imag_self_energy(
+    _fpoints, _gammas, _detailed_gammas = si_pbesol.run_imag_self_energy(
         [1, ],
         [300, ],
         frequency_points=freq_points,
         keep_gamma_detail=True)
     np.testing.assert_allclose(
         detailed_gamma,
-        si_pbesol.detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4)),
+        _detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4)),
         atol=1e-2)
 
 
 def test_imag_self_energy_scat_class1(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    si_pbesol.run_imag_self_energy(
+    _fpoints, _gammas = si_pbesol.run_imag_self_energy(
         [1, 103],
         [300, ],
         frequency_points=freq_points,
@@ -121,13 +121,13 @@ def test_imag_self_energy_scat_class1(si_pbesol):
     # for line in si_pbesol.gammas.reshape(-1, 6):
     #     print(("%10.8f, " * 6) % tuple(line))
     np.testing.assert_allclose(
-        gammas_class1, np.array(si_pbesol.gammas).ravel(), atol=1e-2)
+        gammas_class1, np.swapaxes(_gammas, -1, -2).ravel(), atol=1e-2)
 
 
 def test_imag_self_energy_scat_class2(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    si_pbesol.run_imag_self_energy(
+    _fpoints, _gammas = si_pbesol.run_imag_self_energy(
         [1, 103],
         [300, ],
         frequency_points=freq_points,
@@ -135,4 +135,4 @@ def test_imag_self_energy_scat_class2(si_pbesol):
     # for line in si_pbesol.gammas.reshape(-1, 6):
     #     print(("%10.8f, " * 6) % tuple(line))
     np.testing.assert_allclose(
-        gammas_class2, np.array(si_pbesol.gammas).ravel(), atol=1e-2)
+        gammas_class2, np.swapaxes(_gammas, -1, -2).ravel(), atol=1e-2)

test/phonon3/test_spectral_function.py

@@ -57,6 +57,8 @@ def test_SpectralFunction(si_pbesol):
     # for line in sf.spectral_functions.reshape(-1, 6):
     #     print(("%.7f, " * 6) % tuple(line))
     # raise
-    np.testing.assert_allclose(shifts, sf.shifts.ravel(), atol=1e-2)
+    np.testing.assert_allclose(
-    np.testing.assert_allclose(spec_funcs, sf.spectral_functions.ravel(),
+        shifts, np.swapaxes(sf.shifts, -2, -1).ravel(), atol=1e-2)
-                               atol=1e-4, rtol=1e-2)
+    np.testing.assert_allclose(
+        spec_funcs, np.swapaxes(sf.spectral_functions, -2, -1).ravel(),
+        atol=1e-2, rtol=1e-2)