Constant ph-ph interaction strength option for Peierls approximation

anharmonic/phonon3/__init__.py

@@ -117,6 +117,7 @@ class Phono3py:
     def set_phph_interaction(self,
                              nac_params=None,
                              nac_q_direction=None,
+                             constant_strength=None,
                              frequency_scale_factor=None):
         self._interaction = Interaction(
             self._supercell,
@@ -125,6 +126,7 @@ class Phono3py:
             self._primitive_symmetry,
             fc3=self._fc3,
             band_indices=self._band_indices_flatten,
+            constant_strength=constant_strength,
             frequency_factor_to_THz=self._frequency_factor_to_THz,
             cutoff_frequency=self._cutoff_frequency,
             is_nosym=self._is_nosym,

anharmonic/phonon3/interaction.py

@@ -14,6 +14,7 @@ class Interaction:
                  symmetry,
                  fc3=None,
                  band_indices=None,
+                 constant_strength=None,
                  frequency_factor_to_THz=VaspToTHz,
                  is_nosym=False,
                  symmetrize_fc3_q=False,
@@ -29,6 +30,7 @@ class Interaction:
             self._band_indices = np.arange(num_band, dtype='intc')
         else:
             self._band_indices = np.array(band_indices, dtype='intc')
+        self._constant_strength = constant_strength
         self._frequency_factor_to_THz = frequency_factor_to_THz
 
         if cutoff_frequency is None:
@@ -60,14 +62,20 @@ class Interaction:
     def run(self, lang='C'):
         num_band = self._primitive.get_number_of_atoms() * 3
         num_triplets = len(self._triplets_at_q)
-        self._interaction_strength = np.zeros(
-            (num_triplets, len(self._band_indices), num_band, num_band),
-            dtype='double')
 
-        if lang == 'C':
-            self._run_c()
+        if self._constant_strength is None:
+            self._interaction_strength = np.zeros(
+                (num_triplets, len(self._band_indices), num_band, num_band),
+                dtype='double')
+            if lang == 'C':
+                self._run_c()
+            else:
+                self._run_py()
         else:
-            self._run_py()
+            self._interaction_strength = np.ones(
+                (num_triplets, len(self._band_indices), num_band, num_band),
+                dtype='double') * self._constant_strength
+            self.set_phonon(self._triplets_at_q.ravel())
 
     def get_interaction_strength(self):
         return self._interaction_strength

scripts/phono3py

@@ -106,6 +106,7 @@ parser.set_defaults(band_indices=None,
                     band_paths=None,
                     band_points=None,
                     cell_poscar=None,
+                    constant_pp_strength=None,
                     cutoff_fc3_distance=None,
                     cutoff_frequency=None,
                     boundary_mfp=None,
@@ -228,6 +229,9 @@ parser.add_option("--cfs", "--create_force_sets", dest="force_sets_mode",
 #                   action="store_true", help="Create DELTA_FC2_SETS")
 # parser.add_option("--dfc2", "--delta_fc2", dest="read_delta_fc2",
 #                   action="store_true", help="Read DELTA_FC2_SETS")
+parser.add_option("--constant_pp_strength", dest="constant_pp_strength",
+                  type="float",
+                  help="Constant ph-ph strength")
 parser.add_option("--factor", dest="factor", type="float",
                   help="Conversion factor to favorite frequency unit")
 parser.add_option("--fc2", dest="read_fc2", action="store_true",
@@ -1077,6 +1081,7 @@ if run_mode == "isotope":
 phono3py.set_phph_interaction(
     nac_params=nac_params,
     nac_q_direction=nac_q_direction,
+    constant_strength=options.constant_pp_strength,
     frequency_scale_factor=frequency_scale_factor)
 
 if settings.get_is_linewidth():