Update group velocity

phonopy/group_velocity/__init__.py

@@ -48,21 +48,58 @@ class GroupVelocity:
     """
     
     def __init__(self,
-                 dynmat,
-                 q_direction=None,
-                 q_length=1e-3):
-        self._dynmat = dynmat
-        self._q_direction = q_direction
+                 phonon,
+                 q_points=None,
+                 q_length=1e-4):
+        """
+        q_points is a list of sets of q-point and q-direction:
+        [[q-point, q-direction], [q-point, q-direction], ...]
+        """
+        
+        self._phonon = phonon
+        self._dynmat = phonon.get_dynamical_matrix()
+        self._reciprocal_lattice = np.linalg.inv(
+            self._phonon.get_primitive().get_cell())
+        self._q_points = q_points
         self._q_length = q_length
-
-        if self._q_direction is not None:
+        if self._q_points is not None:
             self._set_group_velocity()
 
-    def set_q_direction(self, q_direction):
-        self._q_direction = q_direction
+    def set_q_points(self, q_points):
+        self._q_points = q_points
+        self._set_group_velocity()
+
+    def set_q_length(self, q_length):
+        self._q_length = q_length
 
     def get_group_velocity(self):
         return self._group_velocity
         
     def _set_group_velocity(self):
-        pass
+        for (q, n) in self._q_points:
+            print q, n
+            self._dynmat.set_dynamical_matrix(q)
+            dm = self._dynmat.get_dynamical_matrix()
+            eigvals, eigvecs = np.linalg.eigh(dm)
+            dD = self._get_dD(np.array(q), np.array(n))
+            dD_at_q = []
+            for dD_i in dD:
+                dD_i_at_q = np.array([np.vdot(eigvec, np.dot(dD_i, eigvec)).real
+                                      for eigvec in eigvecs.T])
+                dD_at_q.append(dD_i_at_q / np.sqrt(np.abs(eigvals)) / 2)
+            print np.array(dD_at_q)
+            print 
+
+    def _get_dD(self, q, n):
+        rlat = self._reciprocal_lattice
+        nc = np.dot(n, rlat)
+        dqc = self._q_length * nc / np.linalg.norm(nc)
+        ddm = []
+        for dqc_i in np.diag(dqc):
+            dq_i = np.dot(dqc_i, np.linalg.inv(rlat))
+            self._dynmat.set_dynamical_matrix(q - dq_i)
+            dm1 = self._dynmat.get_dynamical_matrix()
+            self._dynmat.set_dynamical_matrix(q + dq_i)
+            dm2 = self._dynmat.get_dynamical_matrix()
+            ddm.append(dm2 - dm1)
+        return [ddm_i / dpc_i for (ddm_i, dpc_i) in zip(ddm, dqc)]

test/MgO-group-velocity.py

@@ -1,8 +1,25 @@
 from phonopy import Phonopy
 from phonopy.interface.vasp import read_vasp_from_strings
 from phonopy.hphonopy.file_IO import parse_FORCE_SETS_from_strings, parse_BORN_from_strings
+from phonopy.group_velocity import GroupVelocity
 import numpy as np
 
+def plot_band(phonon):
+    def append_band(bands, q_start, q_end):
+        band = []
+        for i in range(51):
+            band.append(np.array(q_start) +
+                        (np.array(q_end) - np.array(q_start)) / 50 * i)
+        bands.append(band)
+    
+    bands = []
+    append_band(bands, [0.0, 0.0, 0.0], [0.5, 0.5, 0.0])
+    # append_band(bands, [0.5, 0.0, 0.0], [0.5, 0.5, 0.0])
+    # append_band(bands, [0.5, 0.5, 0.0], [0.0, 0.0, 0.0])
+    # append_band(bands, [0.0, 0.0, 0.0], [0.5, 0.5, 0.5])
+    phonon.set_band_structure(bands)
+    phonon.plot_band_structure().show()
+
 poscar_str = """Mg O                                    
    1.00000000000000     
      4.2555564654942897    0.0000000000000000    0.0000000000000000
@@ -161,21 +178,14 @@ born_str = """14.400
 1.9715466666666668 0 0 0 1.9715466666666668 0 0 0 1.9715466666666668
 -1.9721233333333332 0 0 0 -1.9721233333333332 0 0 0 -1.9721233333333332"""
 
+#
+# initial settings
+#
 cell = read_vasp_from_strings(poscar_str)
 phonon = Phonopy(cell, np.diag([2, 2, 2]))
-symmetry = phonon.get_symmetry()
 force_sets = parse_FORCE_SETS_from_strings(force_sets_str,
                                            cell.get_number_of_atoms() * 8)
-sets_of_forces = []
-displacements = []
-for force in force_sets:
-    sets_of_forces.append(force.get_forces())
-    disp = force.get_displacement()
-    atom_number = force.get_atom_number()
-    displacements.append([atom_number,
-                          disp[0], disp[1], disp[2]])
-phonon.set_displacements(displacements)
-phonon.set_forces(sets_of_forces)
+phonon.set_force_sets(force_sets)
 phonon.set_post_process(primitive_matrix=[[0, 0.5, 0.5],
                                           [0.5, 0, 0.5],
                                           [0.5, 0.5, 0]],
@@ -184,18 +194,21 @@ born_params = parse_BORN_from_strings(born_str,
                                       phonon.get_primitive())
 phonon.set_nac_params(born_params)                                      
 
-def append_band(bands, q_start, q_end):
-    band = []
-    for i in range(51):
-        band.append(np.array(q_start) +
-                    (np.array(q_end) - np.array(q_start)) / 50 * i)
-    bands.append(band)
-
-bands = []
-append_band(bands, [0.0, 0.0, 0.0], [0.5, 0.0, 0.0])
-append_band(bands, [0.5, 0.0, 0.0], [0.5, 0.5, 0.0])
-append_band(bands, [0.5, 0.5, 0.0], [0.0, 0.0, 0.0])
-append_band(bands, [0.0, 0.0, 0.0], [0.5, 0.5, 0.5])
-phonon.set_band_structure(bands)
-phonon.plot_band_structure().show()
+print phonon.get_symmetry().get_international_table()
+#
+# Run
+#
+vg = GroupVelocity(phonon)
+vg.set_q_points([
+        # [[0,0,0], [1,1,0]],
+        [[0.1,0.1,0], [1,1,0]],
+        [[0.15,0.15,0], [1,1,0]],
+        [[0.2,0.2,0], [1,1,0]],
+        [[0.25,0.25,0], [1,1,0]],
+        [[0.3,0.3,0], [1,1,0]],
+        [[0.35,0.35,0], [1,1,0]],
+        [[0.4,0.4,0], [1,1,0]],
+        [[0.45,0.45,0], [1,1,0]],
+        [[0.5,0.5,0], [1,1,0]]])
+plot_band(phonon)