Refactor phonon solver call procedure

phono3py/api_phono3py.py

@@ -532,7 +532,8 @@ class Phono3py:
     @nac_params.setter
     def nac_params(self, nac_params):
         self._nac_params = nac_params
-        self._init_dynamical_matrix()
+        if self._interaction is not None:
+            self._init_dynamical_matrix()
 
     def get_nac_params(self):
         """Return NAC parameters."""
@@ -2484,11 +2485,30 @@ class Phono3py:
         )
 
     def _init_dynamical_matrix(self):
-        if self._interaction is not None:
-            self._interaction.init_dynamical_matrix(
-                self._fc2,
-                self._phonon_supercell,
-                self._phonon_primitive,
-                nac_params=self._nac_params,
-                solve_dynamical_matrices=False,
+        if self._interaction is None:
+            msg = (
+                "Phono3py.init_phph_interaction has to be called "
+                "before running this method."
             )
+            raise RuntimeError(msg)
+
+        self._interaction.init_dynamical_matrix(
+            self._fc2,
+            self._phonon_supercell,
+            self._phonon_primitive,
+            nac_params=self._nac_params,
+        )
+        gp_G = self._bz_grid.gp_Gamma
+        self.run_phonon_solver(np.array([gp_G], dtype="int_"))
+        freqs, _, _ = self.get_phonon_data()
+        if np.sum(freqs[gp_G] < self._cutoff_frequency) < 3:
+            for i, f in enumerate(freqs[gp_G, :3]):
+                if not (f < self._cutoff_frequency):
+                    freqs[gp_G, i] = 0
+                    print("=" * 26 + " Warning " + "=" * 26)
+                    print(
+                        " Phonon frequency of band index %d at Gamma "
+                        "is calculated to be %f." % (i + 1, f)
+                    )
+                    print(" But this frequency is forced to be zero.")
+                    print("=" * 61)

phono3py/phonon/grid.py

@@ -124,6 +124,7 @@ class BZGrid:
     QDinv : ndarray
     grid_matrix : ndarray
     microzone_lattice : ndarray
+    gp_Gamma : int
 
     """
 
@@ -189,6 +190,7 @@ class BZGrid:
         self._microzone_lattice = None
         self._rotations = None
         self._reciprocal_operations = None
+        self._gp_Gamma = None
 
         if reciprocal_lattice is not None:
             self._reciprocal_lattice = np.array(
@@ -290,6 +292,11 @@ class BZGrid:
         """
         return self._gp_map
 
+    @property
+    def gp_Gamma(self):
+        """Return grid point index of Gamma-point."""
+        return self._gp_Gamma
+
     @property
     def bzg2grg(self):
         """Transform grid point indices from BZG to GRG.
@@ -434,6 +441,9 @@ class BZGrid:
         self._microzone_lattice = np.dot(
             self._reciprocal_lattice, np.dot(self._QDinv, self._P)
         )
+        self._gp_Gamma = int(
+            self._grg2bzg[get_grid_point_from_address([0, 0, 0], self._D_diag)]
+        )
 
     def _set_rotations(self):
         """Rotation matrices are transformed those for non-diagonal D matrix.

phono3py/phonon3/interaction.py

@@ -50,7 +50,6 @@ from phono3py.phonon3.reciprocal_to_normal import ReciprocalToNormal
 from phono3py.phonon3.triplets import get_nosym_triplets_at_q, get_triplets_at_q
 from phono3py.phonon.grid import (
     BZGrid,
-    get_grid_point_from_address,
     get_grid_points_by_rotations,
     get_ir_grid_points,
 )
@@ -538,14 +537,11 @@ class Interaction:
             elif self._done_nac_at_gamma:
                 if self._nac_q_direction is None:
                     self._done_nac_at_gamma = False
-                    gamma_gp = get_grid_point_from_address(
-                        [0, 0, 0], self._bz_grid.D_diag
-                    )
-                    self._phonon_done[gamma_gp] = 0
+                    self._phonon_done[self._bz_grid.gp_Gamma] = 0
                     self.run_phonon_solver(
                         np.array(
                             [
-                                gamma_gp,
+                                self._bz_grid.gp_Gamma,
                             ],
                             dtype="int_",
                         )
@@ -596,7 +592,6 @@ class Interaction:
         supercell,
         primitive,
         nac_params=None,
-        solve_dynamical_matrices=True,
         decimals=None,
     ):
         """Prepare for phonon calculation on grid.
@@ -617,25 +612,6 @@ class Interaction:
             symprec=self._symprec,
         )
 
-        self._phonon_done[0] = 0
-        if solve_dynamical_matrices:
-            self.run_phonon_solver()
-        else:
-            self.run_phonon_solver(np.array([0], dtype="int_"))
-
-        if (self._bz_grid.addresses[0] == 0).all():
-            if np.sum(self._frequencies[0] < self._cutoff_frequency) < 3:
-                for i, f in enumerate(self._frequencies[0, :3]):
-                    if not (f < self._cutoff_frequency):
-                        self._frequencies[0, i] = 0
-                        print("=" * 26 + " Warning " + "=" * 26)
-                        print(
-                            " Phonon frequency of band index %d at Gamma "
-                            "is calculated to be %f." % (i + 1, f)
-                        )
-                        print(" But this frequency is forced to be zero.")
-                        print("=" * 61)
-
     def set_phonon_data(self, frequencies, eigenvectors, bz_grid_addresses):
         """Set phonons on grid."""
         if bz_grid_addresses.shape != self._bz_grid.addresses.shape:
@@ -655,20 +631,12 @@ class Interaction:
     def run_phonon_solver(self, grid_points=None):
         """Run phonon solver at BZ-grid points."""
         if grid_points is None:
-            self._get_phonons_at_all_bz_grid_points()
+            _grid_points = np.arange(len(self._bz_grid.addresses), dtype="int_")
         else:
             _grid_points = grid_points
-            self._run_phonon_solver_c(_grid_points)
+        self._run_phonon_solver_c(_grid_points)
 
-    def _get_phonons_at_all_bz_grid_points(self, expand_phonons=False):
-        """Get phonons at all BZ-grid points."""
-        if expand_phonons:
-            self._expand_phonons()
-        else:
-            grid_points = np.arange(len(self._bz_grid.addresses), dtype="int_")
-            self._run_phonon_solver_c(grid_points)
-
-    def _expand_phonons(self):
+    def run_phonon_solver_with_eigvec_rotation(self):
         """Phonons at ir-grid-points are copied by proper rotations.
 
         Some phonons that are not covered by rotations are solved.
@@ -679,9 +647,10 @@ class Interaction:
             self._phonon_done
 
         """
+        self._phonon_done[:] = 0
         ir_grid_points, _, _ = get_ir_grid_points(self._bz_grid)
         ir_bz_grid_points = self._bz_grid.grg2bzg[ir_grid_points]
-        self._run_phonon_solver_c(ir_bz_grid_points)
+        self.run_phonon_solver(grid_points=ir_bz_grid_points)
 
         d2r_map = self._get_reciprocal_rotations_in_space_group_operations()
 
@@ -784,8 +753,8 @@ class Interaction:
             )[0]
 
             if self._phonon_done[bzgp_mq] == 0:
-                self._run_phonon_solver_c(
-                    np.array(
+                self.run_phonon_solver(
+                    grid_points=np.array(
                         [
                             bzgp_mq,
                         ],

test/phonon3/test_interaction.py

@@ -153,17 +153,41 @@ def test_interaction_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3p
     )
 
 
-def _get_irt(ph3: Phono3py, mesh, nac_params=None):
+def test_phonon_solver_expand_RTA_si(si_pbesol):
+    """Test phonon solver with eigenvector rotation of Si.
+
+    Eigenvectors can be different but frequencies must be almost the same.
+
+    """
+    itr = _get_irt(si_pbesol, [4, 4, 4])
+    freqs, _, phonon_done = itr.get_phonons()
+    assert (phonon_done == 1).all()
+    itr = _get_irt(si_pbesol, [4, 4, 4], solve_dynamical_matrices=False)
+    itr.run_phonon_solver_with_eigvec_rotation()
+    freqs_expanded, _, _ = itr.get_phonons()
+    np.testing.assert_allclose(freqs, freqs_expanded, rtol=0, atol=1e-6)
+
+
+def _get_irt(ph3: Phono3py, mesh, nac_params=None, solve_dynamical_matrices=True):
     ph3.mesh_numbers = mesh
     itr = Interaction(
         ph3.primitive, ph3.grid, ph3.primitive_symmetry, ph3.fc3, cutoff_frequency=1e-4
     )
     if nac_params is None:
-        itr.init_dynamical_matrix(ph3.fc2, ph3.phonon_supercell, ph3.phonon_primitive)
+        itr.init_dynamical_matrix(
+            ph3.fc2,
+            ph3.phonon_supercell,
+            ph3.phonon_primitive,
+        )
     else:
         itr.init_dynamical_matrix(
-            ph3.fc2, ph3.phonon_supercell, ph3.phonon_primitive, nac_params=nac_params
+            ph3.fc2,
+            ph3.phonon_supercell,
+            ph3.phonon_primitive,
+            nac_params=nac_params,
         )
+    if solve_dynamical_matrices:
+        itr.run_phonon_solver()
     return itr