Add test with sigma=0.1 for RTA kappa

c/_phono3py.c

@@ -1374,6 +1374,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
   /* free(weights); */
   /* free(integration_weights); */
 
+  printf("isoiso\n");
   iso_get_isotope_scattering_strength(gamma,
                                       grid_point,
                                       mass_variances,

c/phono3py.c

@@ -42,6 +42,7 @@
 #include "imag_self_energy_with_g.h"
 #include "real_self_energy.h"
 #include "collision_matrix.h"
+#include "isotope.h"
 
 #include <stdio.h>
 
@@ -441,6 +442,65 @@ void ph3py_get_reducible_collision_matrix(double *collision_matrix,
 }
 
 
+void ph3py_get_isotope_scattering_strength(
+  double *gamma,
+  const size_t grid_point,
+  const double *mass_variances,
+  const double *frequencies,
+  const lapack_complex_double *eigenvectors,
+  const size_t num_grid_points,
+  const int *band_indices,
+  const size_t num_band,
+  const size_t num_band0,
+  const double sigma,
+  const double cutoff_frequency)
+{
+  iso_get_isotope_scattering_strength(gamma,
+                                      grid_point,
+                                      mass_variances,
+                                      frequencies,
+                                      eigenvectors,
+                                      num_grid_points,
+                                      band_indices,
+                                      num_band,
+                                      num_band0,
+                                      sigma,
+                                      cutoff_frequency);
+}
+
+/* void fc3_distribute_fc3(double *fc3,
+ *                         const int target,
+ *                         const int source,
+ *                         const int *atom_mapping,
+ *                         const size_t num_atom,
+ *                         const double *rot_cart)
+ * void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
+ *                           PHPYCONST double (*delta_fc2s)[3][3],
+ *                           const double *inv_U,
+ *                           PHPYCONST double (*site_sym_cart)[3][3],
+ *                           const int *rot_map_syms,
+ *                           const size_t num_atom,
+ *                           const size_t num_site_sym,
+ *                           const size_t num_disp)
+ * void fc3_set_permutation_symmetry_fc3(double *fc3, const size_t num_atom)
+ * void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
+ *                                               const int p2s[],
+ *                                               const int s2pp[],
+ *                                               const int nsym_list[],
+ *                                               const int perms[],
+ *                                               const size_t n_satom,
+ *                                               const size_t n_patom)
+ * void fc3_transpose_compact_fc3(double * fc3,
+ *                                const int p2s[],
+ *                                const int s2pp[],
+ *                                const int nsym_list[],
+ *                                const int perms[],
+ *                                const size_t n_satom,
+ *                                const size_t n_patom,
+ *                                const int t_type) */
+
+
+
 void ph3py_symmetrize_collision_matrix(double *collision_matrix,
                                        const long num_column,
                                        const long num_temp,

c/phono3py.h

@@ -187,6 +187,20 @@ void ph3py_get_reducible_collision_matrix(double *collision_matrix,
                                           const double temperature,
                                           const double unit_conversion_factor,
                                           const double cutoff_frequency);
+void ph3py_get_isotope_scattering_strength(
+  double *gamma,
+  const size_t grid_point,
+  const double *mass_variances,
+  const double *frequencies,
+  const lapack_complex_double *eigenvectors,
+  const size_t num_grid_points,
+  const int *band_indices,
+  const size_t num_band,
+  const size_t num_band0,
+  const double sigma,
+  const double cutoff_frequency);
+
+
 void ph3py_symmetrize_collision_matrix(double *collision_matrix,
                                        const long num_column,
                                        const long num_temp,

phono3py/api_phono3py.py

@@ -94,11 +94,8 @@ class Phono3py(object):
                  calculator=None,
                  log_level=0,
                  lapack_zheev_uplo='L'):
-        if sigmas is None:
-            self._sigmas = [None, ]
-        else:
-            self._sigmas = sigmas
-        self._sigma_cutoff = sigma_cutoff
+        self.sigmas = sigmas
+        self.sigma_cutoff = sigma_cutoff
         self._symprec = symprec
         self._frequency_factor_to_THz = frequency_factor_to_THz
         self._is_symmetry = is_symmetry
@@ -260,6 +257,27 @@ class Phono3py(object):
         """Alias to fc2"""
         return self.fc2
 
+    @property
+    def sigmas(self):
+        return self._sigmas
+
+    @sigmas.setter
+    def sigmas(self, sigmas):
+        if sigmas is None:
+            self._sigmas = [None, ]
+        elif isinstance(sigmas, float) or isinstance(sigmas, int):
+            self._sigmas = [float(sigmas), ]
+        else:
+            self._sigmas = [float(s) for s in sigmas]
+
+    @property
+    def sigma_cutoff(self):
+        return self._sigma_cutoff
+
+    @sigma_cutoff.setter
+    def sigma_cutoff(self, sigma_cutoff):
+        self._sigma_cutoff = sigma_cutoff
+
     @property
     def nac_params(self):
         """Parameters for non-analytical term correction

test/conftest.py

@@ -23,6 +23,16 @@ def si_pbesol():
                          log_level=1)
 
 
+@pytest.fixture(scope='session')
+def si_pbesol_compact_fc():
+    yaml_filename = os.path.join(current_dir, "phono3py_si_pbesol.yaml")
+    forces_fc3_filename = os.path.join(current_dir, "FORCES_FC3_si_pbesol")
+    return phono3py.load(yaml_filename,
+                         forces_fc3_filename=forces_fc3_filename,
+                         is_compact_fc=True,
+                         log_level=1)
+
+
 @pytest.fixture(scope='session')
 def si_pbesol_111():
     yaml_filename = os.path.join(current_dir, "phono3py_params_Si111.yaml")

test/other/test_isotope.py

@@ -21,3 +21,19 @@ def test_Phono3pyIsotope(si_pbesol):
         nac_params=si_pbesol.nac_params)
     iso.run([1, 103])
     np.testing.assert_allclose(si_pbesol_iso, iso.gamma[0], atol=1e-3)
+
+
+def test_Phono3pyIsotope_with_sigma(si_pbesol):
+    si_pbesol.mesh_numbers = [9, 9, 9]
+    iso = Phono3pyIsotope(
+        si_pbesol.mesh_numbers,
+        si_pbesol.phonon_primitive,
+        sigma=0.1,
+        symprec=si_pbesol.symmetry.tolerance)
+    iso.init_dynamical_matrix(
+        si_pbesol.fc2,
+        si_pbesol.phonon_supercell,
+        si_pbesol.phonon_primitive,
+        nac_params=si_pbesol.nac_params)
+    iso.run([1, 103])
+    np.testing.assert_allclose(si_pbesol_iso, iso.gamma[0], atol=1e-3)

test/phonon3/test_kappa_RTA.py

@@ -1,11 +1,28 @@
 import numpy as np
 
 si_pbesol_kappa_RTA = [107.991, 107.991, 107.991, 0, 0, 0]
+si_pbesol_kappa_RTA_with_sigmas = [109.6985, 109.6985, 109.6985, 0, 0, 0]
 
 
 def test_kappa_RTA(si_pbesol):
-    si_pbesol.mesh_numbers = [9, 9, 9]
-    si_pbesol.init_phph_interaction()
-    si_pbesol.run_thermal_conductivity(temperatures=[300, ])
-    kappa = si_pbesol.thermal_conductivity.kappa.ravel()
+    kappa = _get_kappa(si_pbesol, [9, 9, 9]).ravel()
     np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
+
+
+def test_kappa_RTA_with_sigma(si_pbesol):
+    si_pbesol.sigmas = [0.1, ]
+    kappa = _get_kappa(si_pbesol, [9, 9, 9]).ravel()
+    np.testing.assert_allclose(
+        si_pbesol_kappa_RTA_with_sigmas, kappa, atol=0.5)
+
+
+def test_kappa_RTA_compact_fc(si_pbesol_compact_fc):
+    kappa = _get_kappa(si_pbesol_compact_fc, [9, 9, 9]).ravel()
+    np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
+
+
+def _get_kappa(ph3, mesh):
+    ph3.mesh_numbers = mesh
+    ph3.init_phph_interaction()
+    ph3.run_thermal_conductivity(temperatures=[300, ])
+    return ph3.thermal_conductivity.kappa