Made spectral function to run gp by gp

phono3py/phonon3/imag_self_energy.py

@@ -724,14 +724,6 @@ class ImagSelfEnergy(object):
                                      self._frequencies[self._grid_point])
 
 
-def _get_batches(tot_nelems, nelems=10):
-    nbatch = tot_nelems // nelems
-    batches = [np.arange(i * nelems, (i + 1) * nelems)
-               for i in range(nbatch)]
-    batches.append(np.arange(nelems * nbatch, tot_nelems))
-    return batches
-
-
 def _run_ise_at_frequency_points_batch(
         i, j,
         _frequency_points,
@@ -746,9 +738,13 @@ def _run_ise_at_frequency_points_batch(
         log_level=0):
     batches = _get_batches(len(_frequency_points), nelems_in_batch)
 
-    for fpts_batch in batches:
+    if log_level:
+        print("Calculations at %d frequency points are devided into "
+              "%d batches." % (len(_frequency_points), len(batches)))
+
+    for bi, fpts_batch in enumerate(batches):
         if log_level:
-            print(fpts_batch)
+            print("%d/%d: %s" % (bi + 1, len(batches), fpts_batch + 1))
 
         ise.set_frequency_points(_frequency_points[fpts_batch])
         ise.set_integration_weights(
@@ -760,3 +756,12 @@ def _run_ise_at_frequency_points_batch(
             if write_gamma_detail or return_gamma_detail:
                 detailed_gamma_at_gp[j, l, fpts_batch] = (
                     ise.get_detailed_imag_self_energy())
+
+
+def _get_batches(tot_nelems, nelems=10):
+    nbatch = tot_nelems // nelems
+    batches = [np.arange(i * nelems, (i + 1) * nelems)
+               for i in range(nbatch)]
+    if tot_nelems % nelems > 0:
+        batches.append(np.arange(nelems * nbatch, tot_nelems))
+    return batches

phono3py/phonon3/spectral_function.py

@@ -33,7 +33,8 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 import numpy as np
-from phono3py.phonon3.imag_self_energy import get_imag_self_energy
+from phono3py.phonon3.imag_self_energy import (
+    get_imag_self_energy, get_frequency_points)
 from phono3py.phonon3.real_self_energy import imag_to_real
 from phono3py.file_IO import write_spectral_function_at_grid_point
 
@@ -104,15 +105,26 @@ class SpectralFunction(object):
         self._frequency_points = None
 
     def run(self):
+        self._set_frequency_points()
+
+        self._gammas = np.zeros(
+            (len(self._grid_points), len(self._temperatures),
+             len(self._interaction.band_indices), len(self._frequency_points)),
+            dtype='double', order='C')
+        self._deltas = np.zeros_like(self._gammas)
+        self._spectral_functions = np.zeros_like(self._gammas)
+
         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()
+
+        for i, gp in enumerate(self._grid_points):
+            self._run_gamma(i, gp)
+            if self._log_level:
+                print("-" * 11 +
+                      " Real part of self energy by Kramers-Kronig relation "
+                      + "-" * 11)
+            self._run_delta(i)
+            self._run_spectral_function(i, gp)
         if self._log_level:
             print("-" * 75)
 
@@ -136,40 +148,36 @@ class SpectralFunction(object):
     def grid_points(self):
         return self._grid_points
 
-    def _run_gamma(self):
+    def _run_gamma(self, i, grid_point):
         # gammas[grid_points, sigmas, temps, band_indices, freq_points]
         fpoints, gammas = get_imag_self_energy(
             self._interaction,
-            self._grid_points,
+            [grid_point, ],
             frequency_points=self._frequency_points_in,
             frequency_step=self._frequency_step,
             num_frequency_points=self._num_frequency_points,
             temperatures=self._temperatures,
             log_level=self._log_level)
-        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)
-        for i, gp in enumerate(self._grid_points):
-            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._deltas[i, j, k] = -re_part
+        self._gammas[i] = np.array(gammas[0, 0], dtype='double', order='C')
         assert (np.abs(self._frequency_points - fpoints) < 1e-8).all()
 
-    def _run_spectral_function(self):
+    def _run_delta(self, i):
+        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._deltas[i, j, k] = -re_part
+        assert (np.abs(self._frequency_points - fpoints) < 1e-8).all()
+
+    def _run_spectral_function(self, i, grid_point):
         frequencies = self._interaction.get_phonons()[0]
-        self._spectral_functions = np.zeros_like(self._gammas)
-        for i, gp in enumerate(self._grid_points):
-            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]
-                    deltas = self._deltas[i, j, k]
-                    vals = self._get_spectral_function(gammas, deltas, freq)
-                    self._spectral_functions[i, j, k] = vals
+        for j, temp in enumerate(self._temperatures):
+            for k, bi in enumerate(self._interaction.band_indices):
+                freq = frequencies[grid_point, bi]
+                gammas = self._gammas[i, j, k]
+                deltas = self._deltas[i, j, k]
+                vals = self._get_spectral_function(gammas, deltas, freq)
+                self._spectral_functions[i, j, k] = vals
 
     def _get_spectral_function(self, gammas, deltas, freq):
         fpoints = self._frequency_points
@@ -178,3 +186,15 @@ class SpectralFunction(object):
                   + (2 * freq * gammas) ** 2)
         vals = np.where(denoms > 0, nums / denoms, 0)
         return vals
+
+    def _set_frequency_points(self):
+        if (self._interaction.get_phonons()[2] == 0).any():
+            if self._log_level:
+                print("Running harmonic phonon calculations...")
+            self._interaction.run_phonon_solver()
+        max_phonon_freq = np.amax(self._interaction.get_phonons()[0])
+        self._frequency_points = get_frequency_points(
+            max_phonon_freq=max_phonon_freq,
+            frequency_points=self._frequency_points_in,
+            frequency_step=self._frequency_step,
+            num_frequency_points=self._num_frequency_points)

test/phonon3/test_spectral_function.py

@@ -52,11 +52,13 @@ def test_SpectralFunction(si_pbesol):
     sf = SpectralFunction(si_pbesol.phph_interaction,
                           [1, 103],
                           temperatures=[300, ],
-                          num_frequency_points=10)
+                          num_frequency_points=10,
+                          log_level=1)
     sf.run()
 
     # for line in np.swapaxes(sf.spectral_functions, -2, -1).reshape(-1, 6):
     #     print(("%.7f, " * 6) % tuple(line))
+    # raise
 
     np.testing.assert_allclose(
         shifts, np.swapaxes(sf.shifts, -2, -1).ravel(), atol=1e-2)