Spectral function implementation updated

phono3py/api_phono3py.py

@@ -1421,7 +1421,7 @@ class Phono3py(object):
             raise RuntimeError(msg)
 
         if temperatures is None:
-            self._temperatures = [0.0, 300.0]
+            self._temperatures = [300.0, ]
         else:
             self._temperatures = temperatures
         self._grid_points = grid_points
@@ -1440,11 +1440,11 @@ class Phono3py(object):
             output_filename=output_filename,
             log_level=self._log_level)
         if keep_gamma_detail:
-            (self._gammas,
+            (self._frequency_points,
-             self._detailed_gammas,
+             self._gammas,
-             self._frequency_points) = vals
+             self._detailed_gammas) = vals
         else:
-            self._gammas, self._frequency_points = vals
+            self._frequency_points, self._gammas = vals
 
     def write_imag_self_energy(self, filename=None):
         write_imag_self_energy(
@@ -1459,6 +1459,60 @@ class Phono3py(object):
             filename=filename,
             is_mesh_symmetry=self._is_mesh_symmetry)
 
+    def run_frequency_shift(
+            self,
+            grid_points,
+            run_on_bands=False,
+            frequency_points=None,
+            frequency_step=None,
+            num_frequency_points=None,
+            temperatures=None,
+            epsilons=None,
+            write_hdf5=True,
+            output_filename=None):
+        """Frequency shift from lowest order diagram is calculated.
+
+        Parameters
+        ----------
+        epsilons : array_like
+            The value to avoid divergence. When multiple values are given
+            frequency shifts for those values are returned.
+            dtype=float, shape=(epsilons,)
+
+        """
+
+        if self._interaction is None:
+            msg = ("Phono3py.init_phph_interaction has to be called "
+                   "before running this method.")
+            raise RuntimeError(msg)
+
+        if epsilons is not None:
+            _epsilons = epsilons
+        else:
+            if len(self._sigmas) == 1 and self._sigmas[0] is None:
+                _epsilons = None
+            elif self._sigmas[0] is None:
+                _epsilons = self._sigmas[1:]
+            else:
+                _epsilons = self._sigmas
+
+        self._grid_points = grid_points
+        # (epsilon, grid_point, temperature, band)
+        self._frequency_points, self._frequency_shifts = get_frequency_shift(
+            self._interaction,
+            self._grid_points,
+            temperatures,
+            run_on_bands=run_on_bands,
+            frequency_points=frequency_points,
+            frequency_step=frequency_step,
+            num_frequency_points=num_frequency_points,
+            epsilons=_epsilons,
+            write_hdf5=write_hdf5,
+            output_filename=output_filename,
+            log_level=self._log_level)
+
+        return self._frequency_points, self._frequency_shifts
+
     def run_thermal_conductivity(
             self,
             is_LBTE=False,
@@ -1562,60 +1616,6 @@ class Phono3py(object):
                 output_filename=output_filename,
                 log_level=self._log_level)
 
-    def run_frequency_shift(
-            self,
-            grid_points,
-            run_on_bands=False,
-            frequency_points=None,
-            frequency_step=None,
-            num_frequency_points=None,
-            temperatures=None,
-            epsilons=None,
-            write_Delta_hdf5=True,
-            output_filename=None):
-        """Frequency shift from lowest order diagram is calculated.
-
-        Parameters
-        ----------
-        epsilons : array_like
-            The value to avoid divergence. When multiple values are given
-            frequency shifts for those values are returned.
-            dtype=float, shape=(epsilons,)
-
-        """
-
-        if self._interaction is None:
-            msg = ("Phono3py.init_phph_interaction has to be called "
-                   "before running this method.")
-            raise RuntimeError(msg)
-
-        if epsilons is not None:
-            _epsilons = epsilons
-        else:
-            if len(self._sigmas) == 1 and self._sigmas[0] is None:
-                _epsilons = None
-            elif self._sigmas[0] is None:
-                _epsilons = self._sigmas[1:]
-            else:
-                _epsilons = self._sigmas
-
-        self._grid_points = grid_points
-        # (epsilon, grid_point, temperature, band)
-        self._frequency_shift = get_frequency_shift(
-            self._interaction,
-            self._grid_points,
-            run_on_bands=run_on_bands,
-            frequency_points=frequency_points,
-            frequency_step=frequency_step,
-            num_frequency_points=num_frequency_points,
-            epsilons=_epsilons,
-            temperatures=temperatures,
-            write_Delta_hdf5=write_Delta_hdf5,
-            output_filename=output_filename,
-            log_level=self._log_level)
-
-        return self._frequency_shift
-
     def save(self,
              filename="phono3py_params.yaml",
              settings=None):

phono3py/phonon3/frequency_shift.py

@@ -43,14 +43,14 @@ from phono3py.phonon3.imag_self_energy import get_frequency_points
 
 def get_frequency_shift(interaction,
                         grid_points,
+                        temperatures,
                         run_on_bands=False,
                         frequency_points=None,
                         frequency_step=None,
                         num_frequency_points=None,
                         epsilons=None,
-                        temperatures=None,
                         output_filename=None,
-                        write_Delta_hdf5=True,
+                        write_hdf5=True,
                         log_level=0):
     if epsilons is None:
         _epsilons = [None, ]
@@ -134,7 +134,7 @@ def get_frequency_shift(interaction,
                               interaction.get_phonons()[0][gp][bi_set])
                         sys.stdout.flush()
 
-            if write_Delta_hdf5:
+            if write_hdf5:
                 filename = write_Delta_to_hdf5(
                     gp,
                     band_indices,
@@ -152,7 +152,7 @@ def get_frequency_shift(interaction,
                     print("\"%s\"." % filename)
                     sys.stdout.flush()
 
-    return all_deltas
+    return _frequency_points, all_deltas
 
 
 class FrequencyShift(object):

phono3py/phonon3/imag_self_energy.py

@@ -101,8 +101,8 @@ def get_imag_self_energy(interaction,
     Returns
     -------
     tuple :
-        (gammas, frequency_points) are returned. With return_gamma_detail=True,
+        (frequency_points, gammas) are returned. With return_gamma_detail=True,
-        (gammas, detailed_gammas, frequency_points) are returned.
+        (frequency_points, gammas, detailed_gammas) are returned.
 
     """
 
@@ -213,9 +213,9 @@ def get_imag_self_energy(interaction,
                 detailed_gamma.append(detailed_gamma_at_gp)
 
     if return_gamma_detail:
-        return gamma, detailed_gamma, _frequency_points
+        return _frequency_points, gamma, detailed_gamma
     else:
-        return gamma, _frequency_points
+        return _frequency_points, gamma
 
 
 def get_frequency_points(max_phonon_freq=None,

phono3py/phonon3/spectral_function.py

@@ -32,11 +32,31 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
+import sys
 import numpy as np
 from phono3py.phonon3.imag_self_energy import get_imag_self_energy
 from phono3py.phonon3.frequency_shift import imag_to_real
 
 
+def get_spectral_function(interaction,
+                          grid_points,
+                          temperatures,
+                          frequency_points=None,
+                          frequency_step=None,
+                          num_frequency_points=None,
+                          output_filename=None,
+                          write_hdf5=True,
+                          log_level=0):
+    sf = SpectralFunction(interaction,
+                          grid_points,
+                          frequency_points=frequency_points,
+                          frequency_step=frequency_step,
+                          num_frequency_points=num_frequency_points,
+                          temperatures=temperatures)
+    sf.run()
+    return sf
+
+
 class SpectralFunction(object):
     """Calculate spectral function"""
 
@@ -81,7 +101,7 @@ class SpectralFunction(object):
 
     def _run_gamma(self):
         # gammas[grid_points, sigmas, temps, freq_points, band_indices]
-        gammas, self._frequency_points = get_imag_self_energy(
+        fpoints, gammas = get_imag_self_energy(
             self._interaction,
             self._grid_points,
             frequency_points=self._frequency_points_in,
@@ -89,6 +109,7 @@ class SpectralFunction(object):
             num_frequency_points=self._num_frequency_points,
             temperatures=self._temperatures)
         self._gammas = np.array(gammas[:, 0], dtype='double', order='C')
+        self._frequency_points = fpoints
 
     def _run_delta(self):
         self._deltas = np.zeros_like(self._gammas)
@@ -114,7 +135,8 @@ class SpectralFunction(object):
 
     def _get_spectral_function(self, gammas, deltas, freq):
         fpoints = self._frequency_points
-        vals = (4 * freq ** 2 * gammas) / (
+        nums = 4 * freq ** 2 * gammas
-            (fpoints - freq ** 2 - 2 * freq * deltas) ** 2
+        denoms = ((fpoints - freq ** 2 - 2 * freq * deltas) ** 2
-            + (2 * freq * gammas) ** 2)
+                  + (2 * freq * gammas) ** 2)
+        vals = np.where(denoms > 0, nums / denoms, 0)
         return vals

test/phonon3/test_frequency_shift.py

@@ -115,10 +115,10 @@ im_part = [
 def test_frequency_shift(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
-    delta = si_pbesol.run_frequency_shift(
+    _, delta = si_pbesol.run_frequency_shift(
         [1, 103],
         temperatures=[300, ],
-        write_Delta_hdf5=False,
+        write_hdf5=False,
         run_on_bands=True)
     np.testing.assert_allclose(si_pbesol_Delta, delta[0, :, 0], atol=1e-5)
 

test/phonon3/test_spectral_function.py

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