Safeguard preventing from using NAC phonons wrongly for JDOS

phono3py/phonon3/interaction.py

@@ -152,7 +152,7 @@ class Interaction:
         self._g_zero = None
 
         self._phonon_done = None
-        self._is_nac_at_gamma = False  # Phonon at Gamma is calculatd with NAC.
+        self._done_nac_at_gamma = False  # Phonon at Gamma is calculatd with NAC.
         self._frequencies = None
         self._eigenvectors = None
         self._dm = None
@@ -507,7 +507,7 @@ class Interaction:
             # used.
             if (self._bz_grid.addresses[grid_point] == 0).all():
                 if self._nac_q_direction is not None:
-                    self._is_nac_at_gamma = True
+                    self._done_nac_at_gamma = True
                     self._phonon_done[grid_point] = 0
                     self.run_phonon_solver(
                         np.array(
@@ -535,30 +535,30 @@ class Interaction:
                         reciprocal_rotations=rotations,
                         is_time_reversal=False,
                     )
-            else:
-                if self._is_nac_at_gamma:
-                    if self._nac_q_direction is None:
-                        gamma_gp = get_grid_point_from_address(
-                            [0, 0, 0], self._bz_grid.D_diag
+            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.run_phonon_solver(
+                        np.array(
+                            [
+                                gamma_gp,
+                            ],
+                            dtype="int_",
                         )
-                        self._phonon_done[gamma_gp] = 0
-                        self.run_phonon_solver(
-                            np.array(
-                                [
-                                    gamma_gp,
-                                ],
-                                dtype="int_",
-                            )
-                        )
-                    else:
-                        msg = (
-                            "Phonons at Gamma has been calcualted with NAC, "
-                            "but ph-ph interaction is expected to calculate at "
-                            "non-Gamma point. Setting Interaction.nac_q_direction = "
-                            "None, can avoid raising this exception to re-run phonon "
-                            "calculation at Gamma without NAC."
-                        )
-                        raise RuntimeError(msg)
+                    )
+                else:
+                    msg = (
+                        "Phonons at Gamma has been calcualted with NAC, "
+                        "but ph-ph interaction is expected to calculate at "
+                        "non-Gamma point. Setting Interaction.nac_q_direction = "
+                        "None, can avoid raising this exception to re-run phonon "
+                        "calculation at Gamma without NAC."
+                    )
+                    raise RuntimeError(msg)
 
         reciprocal_lattice = np.linalg.inv(self._primitive.cell)
         for triplet in triplets_at_q:

phono3py/phonon3/joint_dos.py

@@ -112,6 +112,7 @@ class JointDos:
 
         self._tetrahedron_method = None
         self._phonon_done = None
+        self._done_nac_at_gamma = False  # Phonon at Gamma is calculatd with NAC.
         self._frequencies = None
         self._eigenvectors = None
 
@@ -214,8 +215,26 @@ class JointDos:
         self._frequency_points = None
         if self._phonon_done is None:
             self._allocate_phonons()
+
         gamma_gp = get_grid_point_from_address([0, 0, 0], self._bz_grid.D_diag)
-        self._phonon_done[gamma_gp] = 0
+        if (self._bz_grid.addresses[grid_point] == 0).all():
+            if self._nac_q_direction is not None:
+                self._done_nac_at_gamma = True
+                self._phonon_done[gamma_gp] = 0
+        elif self._done_nac_at_gamma:
+            if self._nac_q_direction is None:
+                self._done_nac_at_gamma = False
+                self._phonon_done[gamma_gp] = 0
+            else:
+                msg = (
+                    "Phonons at Gamma has been calcualted with NAC, "
+                    "but ph-ph interaction is expected to calculate at "
+                    "non-Gamma point. Setting Interaction.nac_q_direction = "
+                    "None, can avoid raising this exception to re-run phonon "
+                    "calculation at Gamma without NAC."
+                )
+                raise RuntimeError(msg)
+
         self.run_phonon_solver(np.array([gamma_gp, grid_point], dtype="int_"))
 
     def get_triplets_at_q(self):

test/phonon3/test_interaction.py

@@ -119,7 +119,6 @@ def test_interaction_nac_direction_phonon_NaCl(nacl_pbe: Phono3py):
     np.testing.assert_allclose(
         frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 7.41183870], rtol=0, atol=1e-6
     )
-    itr.run()
 
 
 def test_interaction_nac_direction_phonon_NaCl_second_error(nacl_pbe: Phono3py):
@@ -154,7 +153,7 @@ def test_interaction_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3p
     )
 
 
-def _get_irt(ph3, mesh, nac_params=None):
+def _get_irt(ph3: Phono3py, mesh, nac_params=None):
     ph3.mesh_numbers = mesh
     itr = Interaction(
         ph3.primitive, ph3.grid, ph3.primitive_symmetry, ph3.fc3, cutoff_frequency=1e-4

test/phonon3/test_joint_dos.py

@@ -1,7 +1,10 @@
 """Tests for joint-density-of-states."""
 import numpy as np
+import pytest
 
+from phono3py import Phono3py
 from phono3py.api_jointdos import Phono3pyJointDos
+from phono3py.phonon3.joint_dos import JointDos
 
 si_freq_points = [
     0.0000000,
@@ -165,7 +168,7 @@ nacl_jdos_12_at_300K = [
 ]
 
 
-def test_jdos_si(si_pbesol):
+def test_jdos_si(si_pbesol: Phono3py):
     """Test joint-DOS by Si."""
     si_pbesol.mesh_numbers = [9, 9, 9]
     jdos = Phono3pyJointDos(
@@ -185,7 +188,7 @@ def test_jdos_si(si_pbesol):
     )
 
 
-def test_jdso_si_nomeshsym(si_pbesol):
+def test_jdso_si_nomeshsym(si_pbesol: Phono3py):
     """Test joint-DOS without considering mesh symmetry by Si."""
     si_pbesol.mesh_numbers = [9, 9, 9]
     jdos = Phono3pyJointDos(
@@ -206,7 +209,7 @@ def test_jdso_si_nomeshsym(si_pbesol):
     )
 
 
-def test_jdos_nacl(nacl_pbe):
+def test_jdos_nacl(nacl_pbe: Phono3py):
     """Test joint-DOS by NaCl."""
     nacl_pbe.mesh_numbers = [9, 9, 9]
     jdos = Phono3pyJointDos(
@@ -227,7 +230,7 @@ def test_jdos_nacl(nacl_pbe):
     )
 
 
-def test_jdos_nacl_gamma(nacl_pbe):
+def test_jdos_nacl_gamma(nacl_pbe: Phono3py):
     """Test joint-DOS at Gamma-point by NaCl."""
     nacl_pbe.mesh_numbers = [9, 9, 9]
     jdos = Phono3pyJointDos(
@@ -249,7 +252,7 @@ def test_jdos_nacl_gamma(nacl_pbe):
     )
 
 
-def test_jdos_nacl_at_300K(nacl_pbe):
+def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
     """Test joint-DOS at 300K by NaCl."""
     nacl_pbe.mesh_numbers = [9, 9, 9]
     jdos = Phono3pyJointDos(
@@ -273,3 +276,59 @@ def test_jdos_nacl_at_300K(nacl_pbe):
     np.testing.assert_allclose(
         nacl_jdos_12_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
     )
+
+
+def test_jdos_nac_direction_phonon_NaCl(nacl_pbe: Phono3py):
+    """Test JDOS of NaCl with nac_q_direction."""
+    jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
+    jdos.nac_q_direction = [1, 0, 0]
+    jdos.set_grid_point(0)
+    frequencies, _, _ = jdos.get_phonons()
+    np.testing.assert_allclose(
+        frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 7.41183870], rtol=0, atol=1e-6
+    )
+
+
+def test_jdos_nac_direction_phonon_NaCl_second_error(nacl_pbe: Phono3py):
+    """Test JDOS of NaCl with nac_q_direction.
+
+    Second setting non-gamma grid point must raise exception.
+
+    """
+    jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
+    jdos.nac_q_direction = [1, 0, 0]
+    jdos.set_grid_point(0)
+    with pytest.raises(RuntimeError):
+        jdos.set_grid_point(1)
+
+
+def test_jdos_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3py):
+    """Test JDOS of NaCl with nac_q_direction.
+
+    Second setting non-gamma grid point should not raise exception because
+    nac_q_direction = None is set, but the phonons at Gamma is updated to those without
+    NAC.
+
+    """
+    jdos = _get_jdos(nacl_pbe, [7, 7, 7], nac_params=nacl_pbe.nac_params)
+    jdos.nac_q_direction = [1, 0, 0]
+    jdos.set_grid_point(0)
+    jdos.nac_q_direction = None
+    jdos.set_grid_point(1)
+    frequencies, _, _ = jdos.get_phonons()
+    np.testing.assert_allclose(
+        frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 4.59488262], rtol=0, atol=1e-6
+    )
+
+
+def _get_jdos(ph3: Phono3py, mesh, nac_params=None):
+    ph3.mesh_numbers = mesh
+    jdos = JointDos(
+        ph3.primitive,
+        ph3.supercell,
+        ph3.grid,
+        ph3.fc2,
+        nac_params=nac_params,
+        cutoff_frequency=1e-4,
+    )
+    return jdos