From 8e3bc4a5b52381580e7951072044cf9f3c4dc1a5 Mon Sep 17 00:00:00 2001
From: Nakib Haider <nakib.haider.protik@gmail.com>
Date: Mon, 25 Sep 2023 15:45:57 +0200
Subject: [PATCH] Now printing out phonon full BZ energies, velocities, and
 mapping from full to irreducible BZ. These will be used later in external
 post-processing tools.

---
 src/phonon.f90 | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/src/phonon.f90 b/src/phonon.f90
index b3598d62..8f47b538 100644
--- a/src/phonon.f90
+++ b/src/phonon.f90
@@ -167,6 +167,7 @@ subroutine calculate_phonons(self, crys, sym, num, wann)
     integer(i64) :: i, iq, ii, jj, kk, l, il, s, ib, im, chunk, &
          num_active_images
     integer(i64), allocatable :: start[:], end[:]
+    integer(i64) , allocatable :: fbz2ibz_map(:)
     real(r64), allocatable :: ens_chunk(:,:)[:], vels_chunk(:,:,:)[:], &
          symmetrizers_chunk(:,:,:)[:]
     complex(r64), allocatable :: evecs_chunk(:,:,:)[:]
@@ -273,6 +274,17 @@ subroutine calculate_phonons(self, crys, sym, num, wann)
     sync all
     
     if(this_image() <= num_active_images) deallocate(symmetrizers_chunk)
+
+    !Create fbz2ibz_map
+    allocate(fbz2ibz_map(self%nwv))
+    fbz2ibz_map = -1
+    do iq = 1, self%nwv
+       do i = 1, self%nwv_irred !an irreducible point
+          do l = 1, self%nequiv(i) !number of equivalent points of i
+             if(self%ibz2fbz_map(l, i, 2) == iq) fbz2ibz_map(iq) = i
+          end do
+       end do
+    end do
     
     !Symmetrize phonon energies and velocities.
     do i = 1, self%nwv_irred !an irreducible point
@@ -312,6 +324,31 @@ subroutine calculate_phonons(self, crys, sym, num, wann)
        end do
        close(1)
     end if
+
+    !Print out full BZ phonon energies and velocities
+    if(this_image() == 1) then
+       write(numcols, "(I0)") self%numbands
+       open(1, file = "ph.ens_fbz", status = "replace")
+       do iq = 1, self%nwv
+          write(1, "(" // trim(adjustl(numcols)) // "E20.10)") &
+               self%ens(iq, :)
+       end do
+       close(1)
+
+       write(numcols, "(I0)") 3*self%numbands
+       open(1, file = "ph.vels_fbz", status = "replace")
+       do iq = 1, self%nwv
+          write(1, "(" // trim(adjustl(numcols)) // "E20.10)") &
+               self%vels(iq, :, :)
+       end do
+       close(1)
+
+       open(1, file = "ph.fbz2ibz_map", status = "replace")
+       do iq = 1, self%nwv
+          write(1, "(I10)") fbz2ibz_map(iq)
+       end do
+       close(1)
+    end if
     
     !Calculate phonon tetrahedra
     if(num%tetrahedra) then