Updated sscha code

phono3py/sscha/sscha.py

@@ -42,6 +42,7 @@
 """
 
 import numpy as np
+from phonopy.units import VaspToTHz
 from phonopy.harmonic.dynmat_to_fc import DynmatToForceConstants
 from phono3py.phonon.func import mode_length, bose_einstein
 
@@ -148,6 +149,67 @@ def _get_Fz(self, f, f_1, f_2, n_1, n_2):
 
 
 class DispCorrMatrix(object):
+    """Calculate gamma matrix"""
+
+    def __init__(self, supercell_phonon):
+        self._supercell_phonon = supercell_phonon
+        self._gamma_matrix = None
+
+    def run(self, T):
+        freqs = self._supercell_phonon.frequencies
+        eigvecs = self._supercell_phonon.eigenvectors
+        a = mode_length(freqs, T)
+        masses = np.repeat(self._supercell_phonon.supercell.masses, 3)
+        gamma = np.dot(eigvecs, np.dot(np.diag(1.0 / a ** 2), eigvecs.T))
+        self._gamma_matrix = gamma * np.sqrt(np.outer(masses, masses))
+
+    @property
+    def gamma_matrix(self):
+        return self._gamma_matrix
+
+
+class SupercellPhonon(object):
+    def __init__(self, supercell, force_constants, factor=VaspToTHz):
+        self._supercell = supercell
+        _fc2 = np.swapaxes(force_constants, 1, 2)
+        self._force_constants = np.array(
+            _fc2.reshape(-1, np.prod(_fc2.shape[-2:])),
+            dtype='double', order='C')
+        masses = np.repeat(supercell.masses, 3)
+        dynmat = self._force_constants / np.sqrt(np.outer(masses, masses))
+        eigvals, eigvecs = np.linalg.eigh(dynmat)
+        freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals) * factor
+        self._eigenvalues = np.array(eigvals, dtype='double', order='C')
+        self._eigenvectors = np.array(eigvecs, dtype='double', order='C')
+        self._frequencies = np.array(freqs, dtype='double', order='C')
+        self._dynamical_matrix = dynmat
+
+    @property
+    def eigenvalues(self):
+        return self._eigenvalues
+
+    @property
+    def eigenvectors(self):
+        return self._eigenvectors
+
+    @property
+    def frequencies(self):
+        return self._frequencies
+
+    @property
+    def dynamical_matrix(self):
+        return self._dynamical_matrix
+
+    @property
+    def supercell(self):
+        return self._supercell
+
+    @property
+    def force_constants(self):
+        return self._force_constants
+
+
+class DispCorrMatrixMesh(object):
     """Calculate gamma matrix
 
     This calculation is similar to the transformation from
@@ -167,6 +229,19 @@ def commensurate_points(self):
         return self._d2f.commensurate_points
 
     def create_gamma_matrix(self, frequencies, eigenvectors, T):
+        """
+
+        Parameters
+        ----------
+        frequencies : ndarray
+            Supercell phonon frequencies in THz (without 2pi).
+            shape=(grid_point, band), dtype='double', order='C'
+        eigenvectors : ndarray
+            Supercell phonon eigenvectors.
+            shape=(grid_point, band, band), dtype='double', order='C'
+
+        """
+
         a = mode_length(frequencies, T)
         self._d2f.create_dynamical_matrices(1.0 / a ** 2, eigenvectors)
 

test/sscha/test_sscha.py

@@ -1,19 +1,22 @@
 import numpy as np
-from phono3py.sscha.sscha import DispCorrMatrix
+from phono3py.sscha.sscha import (
+    DispCorrMatrix, DispCorrMatrixMesh, SupercellPhonon)
 from phonopy.phonon.qpoints import QpointsPhonon
 
 
-si_pbesol_gamma0_0 = [[96.22804, 0, 0], [0, 96.22804, 0], [0, 0, 96.22804]]
-si_pbesol_gamma1_34 = [[0.067368, -0.336319, -0.162847],
-                       [-0.336319, 0.067368, -0.162847],
-                       [-0.162847, -0.162847, -0.599353]]
+si_pbesol_gamma0_0 = [[3.849187e+02, 0, 0],
+                      [0, 3.849187e+02, 0],
+                      [0, 0, 3.849187e+02]]
+si_pbesol_gamma1_34 = [[1.886404, -1.549705, -1.126055],
+                       [-1.549705, 1.886404, -1.126055],
+                       [-1.126055, -1.126055, -0.006187]]
 
 
-def test_gamma_matrix(si_pbesol):
+def test_gamma_matrix_mesh(si_pbesol):
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
     dynmat = si_pbesol.dynamical_matrix
-    gmat = DispCorrMatrix(dynmat.primitive, dynmat.supercell)
+    gmat = DispCorrMatrixMesh(dynmat.primitive, dynmat.supercell)
     qpoints_phonon = QpointsPhonon(gmat.commensurate_points,
                                    dynmat,
                                    with_eigenvectors=True)
@@ -23,8 +26,26 @@ def test_gamma_matrix(si_pbesol):
     gmat.run()
     np.testing.assert_allclose(si_pbesol_gamma0_0,
                                gmat.gamma_matrix[0, 0],
-                               atol=1e-5)
+                               atol=1e-4)
     np.testing.assert_allclose(si_pbesol_gamma1_34,
                                gmat.gamma_matrix[1, 34],
-                               atol=1e-5)
+                               atol=1e-4)
     print(gmat.gamma_matrix[1, 2])
+
+
+def test_gamma_matrix(si_pbesol):
+    si_pbesol.mesh_numbers = [1, 1, 1]
+    si_pbesol.init_phph_interaction()
+    fc2 = si_pbesol.dynamical_matrix.force_constants
+    supercell = si_pbesol.phonon_supercell
+    factor = si_pbesol.unit_conversion_factor
+    supercell_phonon = SupercellPhonon(supercell, fc2, factor=factor)
+    gmat = DispCorrMatrix(supercell_phonon)
+    gmat.run(300.0)
+    np.testing.assert_allclose(si_pbesol_gamma0_0,
+                               gmat.gamma_matrix[0:3, 0:3],
+                               atol=1e-4)
+    np.testing.assert_allclose(si_pbesol_gamma1_34,
+                               gmat.gamma_matrix[1 * 3: 2 * 3, 34 * 3: 35 * 3],
+                               atol=1e-4)
+    print(gmat.gamma_matrix[1 * 3: 2 * 3, 34 * 3: 35 * 3])