Sampling frequency points for jdos

Handling of nac_q_direction in JointDOS was written similarly to Interaction.

phono3py/api_jointdos.py

@@ -98,6 +98,7 @@ class Phono3pyJointDos:
         self._filename = output_filename
         self._log_level = log_level
 
+        self._bz_grid = None
         self._joint_dos = None
         self._num_frequency_points_in_batch = num_points_in_batch
         self._frequency_step = frequency_step
@@ -108,6 +109,7 @@ class Phono3pyJointDos:
         )
 
         if mesh is not None:
+            self.mesh_numbers = mesh
             self.initialize(mesh)
 
     @property
@@ -161,6 +163,7 @@ class Phono3pyJointDos:
         self._jdos = JointDos(
             self._primitive,
             self._supercell,
+            self._bz_grid,
             self._fc2,
             nac_params=self._nac_params,
             cutoff_frequency=self._cutoff_frequency,
@@ -176,7 +179,6 @@ class Phono3pyJointDos:
             print("Generating grid system ... ", end="", flush=True)
 
         self.mesh_numbers = mesh_numbers
-        self._jdos.bz_grid = self._bz_grid
 
         if self._log_level:
             if self._bz_grid.grid_matrix is None:
@@ -203,7 +205,10 @@ class Phono3pyJointDos:
 
         self._jdos.run_phonon_solver()
         frequencies, _, _ = self._jdos.get_phonons()
+        self._jdos.run_phonon_solver_at_gamma()
         max_phonon_freq = np.max(frequencies)
+        self._jdos.run_phonon_solver_at_gamma(is_nac=True)
+
         self._frequency_points = get_frequency_points(
             max_phonon_freq=max_phonon_freq,
             sigmas=self._sigmas,

phono3py/cui/phono3py_script.py

@@ -829,12 +829,13 @@ def run_gruneisen_then_exit(phono3py, settings, output_filename, log_level):
 
 
 def run_jdos_then_exit(
-    phono3py, settings, updated_settings, output_filename, log_level
+    phono3py: Phono3py, settings, updated_settings, output_filename, log_level
 ):
     """Run joint-DOS calculation."""
     joint_dos = Phono3pyJointDos(
         phono3py.phonon_supercell,
         phono3py.phonon_primitive,
+        phono3py.grid,
         phono3py.fc2,
         mesh=settings.mesh_numbers,
         nac_params=phono3py.nac_params,

phono3py/phonon3/interaction.py

@@ -929,9 +929,9 @@ class Interaction:
         num_grid = len(self._bz_grid.addresses)
         self._phonon_done = np.zeros(num_grid, dtype="byte")
         self._frequencies = np.zeros((num_grid, num_band), dtype="double", order="C")
-        itemsize = self._frequencies.itemsize
+        complex_dtype = "c%d" % (np.dtype("double").itemsize * 2)
         self._eigenvectors = np.zeros(
-            (num_grid, num_band, num_band), dtype=("c%d" % (itemsize * 2)), order="C"
+            (num_grid, num_band, num_band), dtype=complex_dtype, order="C"
         )
         gp_Gamma = self._bz_grid.gp_Gamma
         self.run_phonon_solver_at_gamma()

phono3py/phonon3/joint_dos.py

@@ -56,8 +56,9 @@ class JointDos:
 
     def __init__(
         self,
-        primitive,
-        supercell,
+        primitive: Primitive,
+        supercell: Supercell,
+        bz_grid: BZGrid,
         fc2,
         nac_params=None,
         nac_q_direction=None,
@@ -77,6 +78,7 @@ class JointDos:
         self._grid_point = None
         self._primitive = primitive
         self._supercell = supercell
+        self._bz_grid = bz_grid
         self._fc2 = fc2
         self._nac_params = nac_params
         self.nac_q_direction = nac_q_direction
@@ -99,7 +101,6 @@ class JointDos:
 
         self._num_band = len(self._primitive) * 3
         self._reciprocal_lattice = np.linalg.inv(self._primitive.cell)
-        self._init_dynamical_matrix()
 
         self._tetrahedron_method = None
         self._phonon_done = None
@@ -114,7 +115,8 @@ class JointDos:
         self._g_zero = None
         self._ones_pp_strength = None
         self._temperature = None
-        self._bz_grid = None
+
+        self._init_dynamical_matrix()
 
     @property
     def dynamical_matrix(self) -> DynamicalMatrix:
@@ -216,10 +218,6 @@ class JointDos:
         """Setter and getter of BZGrid."""
         return self._bz_grid
 
-    @bz_grid.setter
-    def bz_grid(self, bz_grid: BZGrid):
-        self._bz_grid = bz_grid
-
     @property
     def temperature(self):
         """Setter and getter of temperature."""
@@ -241,8 +239,6 @@ class JointDos:
         self._grid_point = grid_point
         self._set_triplets()
         self._joint_dos = 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)
         if (self._bz_grid.addresses[grid_point] == 0).all():
@@ -276,8 +272,7 @@ class JointDos:
             _grid_points = np.arange(len(self._bz_grid.addresses), dtype="int_")
         else:
             _grid_points = grid_points
-        if self._phonon_done is None:
-            self._allocate_phonons()
+
         run_phonon_solver_c(
             self._dm,
             self._frequencies,
@@ -291,6 +286,29 @@ class JointDos:
             self._lapack_zheev_uplo,
         )
 
+    def run_phonon_solver_at_gamma(self, is_nac=False):
+        """Run phonon solver at Gamma point.
+
+        Run phonon solver at Gamma point with/without NAC. When `self._nac_q_direction`
+        is None, always without NAC. `self._nac_q_direction` will be unchanged in any
+        case.
+
+        Parameters
+        ----------
+        is_nac : bool, optional
+            With NAC when is_nac is True and `self._nac_q_direction` is not None,
+            otherwise without NAC. Default is False.
+
+        """
+        self._phonon_done[self._bz_grid.gp_Gamma] = 0
+        if is_nac:
+            self.run_phonon_solver(np.array([self._bz_grid.gp_Gamma], dtype="int_"))
+        else:
+            _nac_q_direction = self._nac_q_direction
+            self._nac_q_direction = None
+            self.run_phonon_solver(np.array([self._bz_grid.gp_Gamma], dtype="int_"))
+            self._nac_q_direction = _nac_q_direction
+
     def run(self):
         """Calculate joint-density-of-states."""
         self.run_phonon_solver()
@@ -402,6 +420,7 @@ class JointDos:
             frequency_scale_factor=self._frequency_scale_factor,
             symprec=self._symprec,
         )
+        self._allocate_phonons()
 
     def _set_triplets(self):
         if not self._is_mesh_symmetry:
@@ -421,8 +440,8 @@ class JointDos:
         num_grid = len(self._bz_grid.addresses)
         num_band = self._num_band
         self._phonon_done = np.zeros(num_grid, dtype="byte")
-        self._frequencies = np.zeros((num_grid, num_band), dtype="double")
-        itemsize = self._frequencies.itemsize
+        self._frequencies = np.zeros((num_grid, num_band), dtype="double", order="C")
+        complex_dtype = "c%d" % (np.dtype("double").itemsize * 2)
         self._eigenvectors = np.zeros(
-            (num_grid, num_band, num_band), dtype=("c%d" % (itemsize * 2))
+            (num_grid, num_band, num_band), dtype=complex_dtype, order="C"
         )

test/phonon3/test_joint_dos.py

@@ -169,6 +169,40 @@ nacl_jdos_12_at_300K = [
     0.0000000,
     0.0000000,
 ]
+nacl_freq_points_gamma_at_300K = [
+    0.0000000,
+    1.6322306,
+    3.2644613,
+    4.8966919,
+    6.5289225,
+    8.1611531,
+    9.7933838,
+    11.4256144,
+    13.0578450,
+    14.6900756,
+]
+nacl_jdos_gamma_at_300K = [
+    0.0000000,
+    0.0000000,
+    6.3607672,
+    0.4210009,
+    2.9647113,
+    2.3994749,
+    0.9360874,
+    5.2286115,
+    0.1977176,
+    22.0282005,
+    0.0000000,
+    32.0059314,
+    0.0000000,
+    13.9738865,
+    0.0000000,
+    2.1095895,
+    0.0000000,
+    0.3079461,
+    0.0000000,
+    0.0213677,
+]
 
 
 @pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
@@ -297,6 +331,38 @@ def test_jdos_nacl_at_300K(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
     )
 
 
+def test_jdos_nacl_nac_gamma_at_300K_npoints(nacl_pbe: Phono3py):
+    """Real part of self energy spectrum of NaCl.
+
+    * at 10 frequency points sampled uniformly.
+    * at q->0
+
+    """
+    nacl_pbe.mesh_numbers = [9, 9, 9]
+    jdos = Phono3pyJointDos(
+        nacl_pbe.phonon_supercell,
+        nacl_pbe.phonon_primitive,
+        nacl_pbe.fc2,
+        mesh=nacl_pbe.mesh_numbers,
+        nac_params=nacl_pbe.nac_params,
+        nac_q_direction=[1, 0, 0],
+        num_frequency_points=10,
+        temperatures=[
+            300,
+        ],
+        log_level=1,
+    )
+    jdos.run([nacl_pbe.grid.gp_Gamma])
+    # print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
+    np.testing.assert_allclose(
+        nacl_freq_points_gamma_at_300K, jdos.frequency_points, atol=1e-5
+    )
+    # print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
+    np.testing.assert_allclose(
+        nacl_jdos_gamma_at_300K[2:], jdos.joint_dos.ravel()[2:], rtol=1e-2, atol=1e-5
+    )
+
+
 @pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
 def test_jdos_nac_direction_phonon_NaCl(
     nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
@@ -434,10 +500,10 @@ def _get_jdos(ph3: Phono3py, mesh, nac_params=None, store_dense_gp_map=False):
     jdos = JointDos(
         ph3.primitive,
         ph3.supercell,
+        bz_grid,
         ph3.fc2,
         nac_params=nac_params,
         store_dense_gp_map=store_dense_gp_map,
         cutoff_frequency=1e-4,
     )
-    jdos.bz_grid = bz_grid
     return jdos